Загружаем нужные пакеты:
library(knitr)
opts_chunk$set(cache = FALSE)
library(ggplot2) # графики
library(kernlab) # Support Vector Machines
# install.packages('kernlab') # может быть нужно установить пакеты...
require(plyr)
require(caret)
Другой популярный пакет для SVM в R — это e1071.
Ядерная функция, или просто ядро, — это скалярное произведение в преобразованном пространстве.
\[ K(x,x')=(\phi(x),\phi(x')) \]
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 пытается разделить гиперплоскостью наши данные в более многомерном пространстве, чем исходное, видео. Это новое пространство называется спрямляющим.
Для гауссовского ядра это спрямляющее пространство будет бесконечномерным.
Загружаем данные по стоимости квартир в Москве:
filename <- "../datasets/flats_moscow.txt"
h <- read.table(filename, header = TRUE)
str(h)
## '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 нужно указать, что зависимая переменная — факторная, а не количественная. Более того, названия категорий должны быть валидными именами для переменных.
h$brick <- as.factor(h$brick)
levels(h$brick) <- list(no = "0", yes = "1")
h$walk <- as.factor(h$walk)
Разделим выборку на две части — 75% для обучения и 25% для оценки качества обучения. Указание зависимой переменной важно, т.к. R случайно отберет 75% единичек и 75% нулей в обучающую выборку.
set.seed(777)
train.index <- createDataPartition(y = h$brick, p = 0.75, list = FALSE)
train.h <- h[train.index, ]
test.h <- h[-train.index, ]
Команда set.seed обеспечивает воспроизводимость эксперимента на другом компьютере.
Находим оптимальную разделяющую гиперплоскость
m1 <- ksvm(brick ~ price + totsp + livesp + kitsp + dist + metrdist + walk +
floor + code, data = train.h, kernel = "rbfdot", kpar = list(sigma = 0.05),
C = 5)
Настраивать можно:
kernel="rbfdot"kpar=list(sigma=0.05)C=5Есть встроенные методы подбора параметра \( \sigma \) для Гауссовского ядра:
m1 <- ksvm(brick ~ price + totsp + livesp + kitsp + dist + metrdist + walk +
floor + code, data = train.h, kernel = "rbfdot", C = 5)
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
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.119322729770998
##
## Number of Support Vectors : 787
##
## Objective Function Value : -2887
## Training error : 0.152188
Слабо добраться до оптимального автоматически найденного \( \sigma \)?
kernelf(m1) # описание ядра: тип и параметры
## Gaussian Radial Basis kernel function.
## Hyperparameter : sigma = 0.119322729770998
kpar(kernelf(m1))$sigma # из описания извлекаем параметры, из них --- сигму
## [1] 0.1193
Вектора, лежащие на границе разделяющей полосы, называются опорными.
Найденная оптимальная гиперплоскость не даёт нам простой интепретации зависимости. Из объекта m1 можно извлечь опорные вектора, но в реальной практической задаче их как минимум сотни, и что с ними делать в такой ситуации мне не известно.
Строим прогнозы для тестовой части выборки:
test.h$brick.pred <- predict(m1, test.h)
table(test.h$brick.pred, test.h$brick)
##
## no yes
## no 305 69
## yes 40 95
Также можно спрогнозировать расстояние до разделяющей гиперплоскости…
…
Скомбинировав SVM и logit, можно в каком-то смысл оценить \( P(y_i=1) \) с помощью SVM. Сначала нужно построить прогноз по SVM для всех объектов в выборке. (тут точнее…) Используя расстояние от каждого прогноза до разделяющей гиперплоскости мы оценим модель похожую на logit…
Этот подход уже автоматизирован:
m1 <- ksvm(brick ~ price + totsp + livesp + kitsp + dist + metrdist + walk +
floor + code, data = train.h, kernel = "rbfdot", kpar = list(sigma = 0.05),
C = 5, prob.model = TRUE)
probs <- predict(m1, test.h, type = "probabilities")
head(probs)
## no yes
## [1,] 0.4158 0.5842
## [2,] 0.7558 0.2442
## [3,] 0.5577 0.4423
## [4,] 0.3341 0.6659
## [5,] 0.2358 0.7642
## [6,] 0.3530 0.6470
Используя этот подход можно построить и кривые предельных эффектов.
Красивые графики, как мне кажется, можно получить только для случая двух объясняющих переменных
Пример из документации:
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")
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
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). Её идея состоит в том, что надо оценивать качество прогнозов не по той же выборке, на которой оценивалась модель, а на новых наблюдениях.
set.seed() для воспроизводимости эксперимента.Популярные значения k:
Создаём табличку перебираемых \( 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(data = h, brick ~ price + totsp + livesp + kitsp + dist + metrdist +
walk + floor + code, kernel = rbfdot, kpar = list(sigma = 1), C = 1, cross = 10)
cross(m1) # cross validation error
## [1] 0.2211
Для удобства оформляем это в функцию двух переменных
f_cross <- function(sig, C) {
model <- ksvm(data = h, brick ~ price + totsp + livesp + kitsp + dist +
metrdist + walk + floor + code, kernel = rbfdot, kpar = list(sigma = sig),
C = C, cross = 10)
return(cross(model))
}
f_cross(1, 1) # тестируем сделанную функцию
## [1] 0.2279
Применяем её ко всем возможным \( C \) и \( \sigma \)
d <- ddply(d, ~C + sigma, transform, cr = f_cross(sigma, C))
head(d)
## C sigma cr
## 1 1 0.1 0.2186
## 2 1 1.0 0.2250
## 3 1 10.0 0.2941
## 4 10 0.1 0.2211
## 5 10 1.0 0.2363
## 6 10 10.0 0.2877
Есть замечательный пакет caret созвучный с carrot :) В нем реализовано и деление выборки на тестовую и обучающую части и подбор параметров с помощью кросс-валидации для кучи методов, в том числе для SVM. Морковка очень вкусная :)
require(caret)
ctrl <- trainControl(classProbs = TRUE)
fit <- train(brick ~ price + totsp + livesp + kitsp + dist + metrdist + walk +
floor + code, data = train.h, method = "svmRadial", trControl = ctrl)
fit
## Support Vector Machines with Radial Basis Function Kernel
##
## 1531 samples
## 10 predictors
## 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 Accuracy SD Kappa SD
## 0.2 0.8 0.4 0.02 0.04
## 0.5 0.8 0.4 0.02 0.05
## 1 0.8 0.5 0.02 0.04
##
## Tuning parameter 'sigma' was held constant at a value of 0.1014
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were sigma = 0.1 and C = 1.
Можно прогнозировать с помощью наилучшей модели
head(predict(fit, test.h, type = "prob"))
## no yes
## 1 0.3847 0.6153
## 2 0.7522 0.2478
## 3 0.4840 0.5160
## 4 0.2794 0.7206
## 5 0.2209 0.7791
## 6 0.3718 0.6282
Если \( y \) принимает больше двух значений, то можно использовать метод “все против всех”:
Естественно, это уже автоматизировано. Есть и другие методы.
Чайники используют следующую процедуру оценивания SVM:
Суровые челябинские жители советают делать так:
rbfdot в пакете kernlabA Practical Guide to Support Vector Classiffication, отсюда взята инструкция по практике SVM
A User's Guide to Support Vector Machines — картинки, помогающие понять разницу гауссовского ядра для разных \( C \) и \( \sigma \) с примерами на питоне. Питон — это еще один подходящий язык програмирования для анализа данных.
Пакет caret: официальный сайт, статья в JSS, и введение побольше. Пакет достаточно активно обновляется и это радует :)