Загружаем нужные пакеты:

library("knitr")
# opts_chunk$set(cache=FALSE)

library("ggplot2") # графики
library("kernlab") # Support Vector Machines
library("dplyr") # манипуляции с данными
library("caret") # построение моделей

Другой популярный пакет для SVM в R — это e1071.

Ядрёная функция!

Ядерная функция, или просто ядро, — это скалярное произведение в преобразованном пространстве.

\[ K(x,x')=(\phi(x),\phi(x')) \]

Гауссовское ядро: \[ K(x, x') = \exp(-\sigma |x-x'|^2) \]

x <- 1:6
y <- 1:6
k1 <- vanilladot()
k2 <- rbfdot(sigma = 1)
k1(x, y)
##      [,1]
## [1,]   91
k2(x, y)
##      [,1]
## [1,]    1

Скалярное произведение трудно интерпретировать “в лоб”. Лучшая известная мне интерпретация: скалярное произведение — это произведение длин векторов на косинус угла между ними, \((y,y')=|y|\cdot |y'| \cos(y,y')\).

Скалярное произведение содержит в себе всю информацию о геометрии пространства. Если знать чему равно скалярное произведение между любыми векторами, то можно посчитать длину любого вектора и угол между любыми векторами.

Теория SVM

\[ \min_{w, w_0} \frac{1}{2} w'w + C \sum_{i=1}^n \xi_i \]

Уравнение разделяющей поверхности: \(w'x = w_0\), края полосы: \(w'x=w_0+1\) и \(w'x=w_0-1\). Нарушителями считаются наблюдения, которые попали на нейтральную полосу или на чужую территорию. Здесь \(\xi_i = |w| \cdot d_i\), где \(d_i\) — “заступ” наблюдения за черту “своих”.

SVM пытается разделить гиперплоскостью наши данные в более многомерном пространстве, чем исходное, видео. Это новое пространство называется спрямляющим.

Для гауссовского ядра это спрямляющее пространство будет бесконечномерным, а все векторы превращаются в векторы единичной длины.

Оцениваем SVM с заданными параметрами

Загружаем данные по стоимости квартир в Москве:

filename <- "../datasets/flats_moscow.txt"
flats <- read.table(filename, header = TRUE)

str(flats)
## 'data.frame':    2040 obs. of  11 variables:
##  $ n       : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ price   : int  81 75 128 95 330 137 98 88 225 140 ...
##  $ totsp   : int  58 44 70 61 104 76 59 55 80 86 ...
##  $ livesp  : int  40 28 42 37 60 50 39 36 56 51 ...
##  $ kitsp   : num  6 6 6 6 11 9 6 6 9 10 ...
##  $ dist    : num  12.5 13.5 14.5 13.5 10.5 11 7.5 9 9 12.7 ...
##  $ metrdist: int  7 7 3 7 7 7 10 5 5 10 ...
##  $ walk    : int  1 1 1 1 0 1 0 1 1 1 ...
##  $ brick   : int  1 0 1 0 1 1 0 1 1 0 ...
##  $ floor   : int  1 1 1 1 1 1 1 0 1 1 ...
##  $ code    : int  3 6 3 1 3 8 8 4 3 5 ...

Для ksvm нужно указать, что зависимая переменная — факторная, а не количественная. Более того, названия категорий должны быть валидными именами для переменных.

flats$brick <- as.factor(flats$brick)
levels(flats$brick) <- list(no = "0", yes = "1") 
flats$walk <- as.factor(flats$walk)

Разделим выборку на две части — 75% для обучения и 25% для оценки качества обучения. Указание зависимой переменной важно, т.к. R случайно отберет 75% единичек и 75% нулей в обучающую выборку.

set.seed(777)

train_index <- createDataPartition(y = flats$brick, p = 0.75, list = FALSE)
train_flats <- flats[train_index, ]
test_flats <- flats[-train_index, ]

Команда set.seed обеспечивает воспроизводимость эксперимента на другом компьютере.

Находим оптимальную разделяющую гиперплоскость:

m1 <- ksvm(brick ~ price + totsp + livesp + kitsp + dist + metrdist + walk + floor + code, 
           data = train_flats, 
           kernel = "rbfdot", 
           kpar = list(sigma = 0.05), 
           C = 5)

Настраивать можно:

Есть встроенные методы подбора параметра \(\sigma\) для Гауссовского ядра:

m1 <- ksvm(brick ~ price + totsp + livesp + kitsp + dist + metrdist + walk + floor + code, 
           data = train_flats, 
           kernel = "rbfdot", 
           C = 5)
m1
## Support Vector Machine object of class "ksvm" 
## 
## SV type: C-svc  (classification) 
##  parameter : cost C = 5 
## 
## Gaussian Radial Basis kernel function. 
##  Hyperparameter : sigma =  0.110127445667694 
## 
## Number of Support Vectors : 786 
## 
## Objective Function Value : -2946.599 
## Training error : 0.157413

Слабо добраться до оптимального автоматически найденного \(\sigma\)?

kernelf(m1) # описание ядра: тип и параметры
## Gaussian Radial Basis kernel function. 
##  Hyperparameter : sigma =  0.110127445667694
kpar(kernelf(m1))$sigma # из описания извлекаем параметры, из них --- сигму
## [1] 0.1101274

Вектора, лежащие на границе разделяющей полосы, называются опорными.

Найденная оптимальная гиперплоскость не даёт нам простой интепретации зависимости. Из объекта m1 можно извлечь опорные вектора, но в реальной практической задаче их как минимум сотни, и что с ними делать в такой ситуации мне не известно.

Прогнозы и прогнозные вероятности

Строим прогнозы для тестовой части выборки:

test_flats$brick.pred <- predict(m1, test_flats)
table(test_flats$brick.pred, test_flats$brick)
##      
##        no yes
##   no  308  75
##   yes  37  89

Также можно спрогнозировать расстояние до разделяющей гиперплоскости…

…

Скомбинировав SVM и logit, можно в каком-то смысл оценить \(P(y_i=1)\) с помощью SVM. Сначала нужно построить прогноз по SVM для всех объектов в выборке. (тут точнее…) Используя расстояние от каждого прогноза до разделяющей гиперплоскости мы оценим модель похожую на logit…

Этот подход уже автоматизирован:

m1 <- ksvm(brick ~ price + totsp + livesp + kitsp + dist + metrdist + walk + floor + code, 
           data = train_flats,
           kernel = "rbfdot",
           kpar = list(sigma = 0.05),
           C = 5,
           prob.model = TRUE)

probs <- predict(m1, test_flats, type = "probabilities")
head(probs)
##             no       yes
## [1,] 0.4157523 0.5842477
## [2,] 0.7558480 0.2441520
## [3,] 0.5577324 0.4422676
## [4,] 0.3340514 0.6659486
## [5,] 0.2358058 0.7641942
## [6,] 0.3530213 0.6469787

Используя этот подход можно построить и кривые предельных эффектов.

В случае, если SVM идеально разделяет игреки на обучающей выборке оценить вероятности не получится. Дело здесь в том, что при идеальной классификации не существует оптимума правдоподобия в задаче оценки логит-модели.

Няшные графики

Красивые графики, как мне кажется, можно получить только для случая двух объясняющих переменных

Пример из документации:

x <- rbind(matrix(rnorm(120), ncol = 2), matrix(rnorm(120, mean = 3), ncol = 2))
y <- matrix(c(rep(1, 60), rep(-1, 60)))
 
svp <- ksvm(x, y, type = "C-svc")
plot(svp, data = x)

Закрашенные треугольники и кружочки — это опорные вектора, т.е. те точки, которые оказались на границе разделяющей полосы в спрямляющем пространстве.

Выбор между точностью подгонки и простотой модели

Во многих моделях есть параметр, отвечающий за простоту модели. Чем проще модель, тем хуже модель описывает выборку, по которой она оценивалась.

В линейной регрессии параметр сложности модели — это количество регрессоров, \(k\). Чем больше регрессоров, тем ниже сумма квадратов остатков, \(RSS\).

На гистограмме параметр сложности — число столбцов.

В ridge regression и LASSO параметр простоты — это \(\lambda\). Действительно, LASSO минимизирует \[ \sum_{i=1}^n (y_i-\hat{y}_i)^2+\lambda \sum_{j=1}^k |\hat{\beta}_j| \] Значит, чем больше параметр \(\lambda\), тем больше стремление алгоритма LASSO занулить некоторые \(\hat{\beta}_j\).

В SVM за “простоту” отвечают \(\sigma\) и \(C\).

Этот параметр сложности не так легко оценить, как коэффициенты модели. Если наивно попытаться выбрать параметр сложности так, чтобы модель как можно лучше описывала бы выборку, по которой она оценивалась… То ничего хорошего не выйдет. Окажется, что оптимальное количество регрессоров равно плюс бесконечности, а оптимальное \(\lambda\) в LASSO равно нулю. В регрессии одно из решений этой проблемы — это проверка гипотез о значимости коэффициентов.

Есть и универсальный способ выбора сложности модели — кросс валидация (перекрёстная проверка, cross validation). Её идея состоит в том, что надо оценивать качество прогнозов не по той же выборке, на которой оценивалась модель, а на новых наблюдениях.

k-кратная кросс-валидация

  • Разбиваем случайным образом всю выборку на k частей. Сразу мораль: используйте set.seed() для воспроизводимости эксперимента.
  • Прогнозы для первой части выборки строим, оценивая модель по наблюдениям всех остальных частей. Прогнозы для второй части выборки строим, оценивая модель по наблюдениям всех частей кроме второй. И так далее.
  • Получаем по одному прогнозу для каждого наблюдения.
  • Считаем сумму квадратов ошибок прогнозов или другой показатель их качества.

Популярные значения k:

  • 10-кратная кроссвалидация.
  • k равно числу наблюдений. Т.е. модель оценивается по всем наблюдениям кроме одного. Для этого невключенного наблюдения считается ошибка прогноза. Исключая по очереди то одно, то другое наблюдение, получаем ошибку прогноза для каждого наблюдения. По этим ошибкам считаем сумму квадратов. При этом подходе алгоритм не является случайным и зачастую есть готовые формулы для суммы квадратов ошибок прогнозов.

Кросс-валидация на примере SVM

Создаём табличку перебираемых \(C\) и \(\sigma\):

C <- c(1, 10, 100)
sigma <- c(0.1, 1, 10)
d <- expand.grid(C, sigma) # Случайно не декартово ли произведение это? :)
colnames(d) <- c("C", "sigma")
head(d)
##     C sigma
## 1   1   0.1
## 2  10   0.1
## 3 100   0.1
## 4   1   1.0
## 5  10   1.0
## 6 100   1.0

Пробуем 10-кратную кросс-валидацию для конкретных \(C\) и \(\sigma\):

set.seed(33222233) # любимый seed Фрекен Бок

m1 <- ksvm(brick ~ price + totsp + livesp + kitsp + dist + metrdist + walk + floor + code,
           data = flats,
           kernel = rbfdot,
           kpar = list(sigma = 1),
           C = 1,
           cross = 10)

cross(m1) # cross validation error
## [1] 0.2210784

Для удобства оформляем это в функцию двух переменных:

f_cross <- function(sig, C) {
  model <- ksvm(brick ~ price + totsp + livesp + kitsp + dist + metrdist + walk + floor + code,
                data = flats,
                kernel = rbfdot,
                kpar = list(sigma = sig),
                C = C,
                cross = 10)
  
  return(cross(model))
}

f_cross(1, 1) # тестируем сделанную функцию
## [1] 0.2279412

Применяем её ко всем возможным \(C\) и \(\sigma\):

d <- mutate(d, cr = f_cross(sigma, C))

head(d)
##     C sigma    cr
## 1   1   0.1 0.475
## 2  10   0.1 0.475
## 3 100   0.1 0.475
## 4   1   1.0 0.475
## 5  10   1.0 0.475
## 6 100   1.0 0.475

Вкусная Морковка

Есть замечательный пакет caret созвучный с carrot :) В нем реализовано и деление выборки на тестовую и обучающую части и подбор параметров с помощью кросс-валидации для кучи методов, в том числе для SVM. Морковка очень вкусная :)

library(caret)
ctrl <- trainControl(classProbs = TRUE)
fit <- train(brick ~ price + totsp + livesp + kitsp + dist + metrdist + walk + floor + code,
             data = train_flats,
             method = "svmRadial",
             trControl = ctrl)
fit
## Support Vector Machines with Radial Basis Function Kernel 
## 
## 1531 samples
##   10 predictor
##    2 classes: 'no', 'yes' 
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 1531, 1531, 1531, 1531, 1531, 1531, ... 
## Resampling results across tuning parameters:
## 
##   C     Accuracy   Kappa    
##   0.25  0.7680351  0.4340640
##   0.50  0.7733471  0.4488185
##   1.00  0.7748407  0.4539361
## 
## Tuning parameter 'sigma' was held constant at a value of 0.09918708
## Accuracy was used to select the optimal model using  the largest value.
## The final values used for the model were sigma = 0.09918708 and C = 1.

Можно прогнозировать с помощью наилучшей модели

head(predict(fit, test_flats , type = "prob"))
##          no       yes
## 1 0.3985826 0.6014174
## 2 0.7497911 0.2502089
## 3 0.4992112 0.5007888
## 4 0.2976260 0.7023740
## 5 0.2369349 0.7630651
## 6 0.3789946 0.6210054

Недостатки SVM в задаче классификации

  • Оптимальная разделяющая гиперплоскость не даёт простой формулы, описывающей зависимость \(y\) от регрессоров. Нет оценок, которые бы легко интерпретировались.
  • Напрямую SVM применима только для двух классов. Есть способы расширить её на много классов.

Как применить SVM, если \(y\) принимает больше двух значений?

Если \(y\) принимает больше двух значений, то можно использовать метод “все против всех”:

  • Для каждой пары возможных значений \(y\), \(y=a\) и \(y=b\), составить выборку, содержащую только эти значения \(y\).
  • На каждой из этих выборок запустить SVM, которая будет выбирать из двух значений \(y\) одно.
  • В результате мы получили кворум SVM-ок, голосующих каждая за свой прогноз \(y\). Выбираем прогноз простым большинством голосов.

Естественно, это уже автоматизировано. Есть и другие методы.

SVM на практике

Чайники используют следующую процедуру оценивания SVM:

  • Оценить несколько SVM с параметрами “от фонаря”.
  • Выбрать наилучшую.

Суровые челябинские жители советают делать так:

  • Нормировать все переменные. Например, можно вычитать среднее и делить на корень из дисперсии. Можно все переменные привести к диапазону \([0;1]\). Гуру постигшие абсолют рекомендуют диапазон \([0;10]\), как более понятный простым смертным.
  • Поделить выборку на две части, обучающую и тестовую.
  • По обучающей части оценить SVM с гауссовским ядром, rbfdot в пакете kernlab.
  • Параметры \(C\) и \(\sigma\) подобрать с помощью кросс-валидации.
  • Оценить прогнозную силу полученной модели по тестовой выборке.
  • Если набор данных слишком большой и кросс-валидация занимает не меряно времени, то её можно делать на небольшой случайной подвыборке обучающей части.

Ссылки: