Chapter 33 Support Vector Machine
Definition: Support vector machine is a representation of the training data as points in space separated into categories by a clear gap that is as wide as possible. New examples are then mapped into that same space and predicted to belong to a category based on which side of the gap they fall.
SVMs are supervised learning models with associated learning algorithms that analyze data for classification and regression analysis.
Advantages: Effective in high dimensional spaces and uses a subset of training points in the decision function so it is also memory efficient.
Disadvantages: The algorithm does not directly provide probability estimates, these are calculated using an expensive five-fold cross-validation.
Parameters of SVM:
* Type of kernel
* Gamma value
* C value
### Support Vector Machine
# Classification
library(e1071)
# Generate data
set.seed(10111)
x = matrix(rnorm(40), 20, 2)
y = rep(c(-1, 1), c(10, 10))
x[y == 1,] = x[y == 1,] + 1
# Plot generated data
plot(x, col = y + 3, pch = 19)
dat = data.frame(x, y = as.factor(y))
svmfit = svm(y ~ ., data = dat, kernel = "linear", cost = 10, scale = FALSE)
print(svmfit)
plot(svmfit, dat)
make.grid = function(x, n = 75) {
grange = apply(x, 2, range)
x1 = seq(from = grange[1,1], to = grange[2,1], length = n)
x2 = seq(from = grange[1,2], to = grange[2,2], length = n)
expand.grid(X1 = x1, X2 = x2)
}
xgrid = make.grid(x)
xgrid[1:10,]
ygrid = predict(svmfit, xgrid)
plot(xgrid, col = c("red","blue")[as.numeric(ygrid)], pch = 20, cex = .2)
points(x, col = y + 3, pch = 19)
points(x[svmfit$index,], pch = 5, cex = 2)
beta = drop(t(svmfit$coefs)%*%x[svmfit$index,])
beta0 = svmfit$rho
plot(xgrid, col = c("red", "blue")[as.numeric(ygrid)], pch = 20, cex = .2)
points(x, col = y + 3, pch = 19)
points(x[svmfit$index,], pch = 5, cex = 2)
abline(beta0 / beta[2], -beta[1] / beta[2])
abline((beta0 - 1) / beta[2], -beta[1] / beta[2], lty = 2)
abline((beta0 + 1) / beta[2], -beta[1] / beta[2], lty = 2)
### 2. Non-Linear SVM Classifier
load(file = "ESL.mixture.rda")
names(ESL.mixture)
rm(x, y)
attach(ESL.mixture)
plot(x, col = y + 1)
dat = data.frame(y = factor(y), x)
fit = svm(factor(y) ~ ., data = dat, scale = FALSE, kernel = "radial", cost = 5)
xgrid = expand.grid(X1 = px1, X2 = px2)
ygrid = predict(fit, xgrid)
plot(xgrid, col = as.numeric(ygrid), pch = 20, cex = .2)
points(x, col = y + 1, pch = 19)
func = predict(fit, xgrid, decision.values = TRUE)
func = attributes(func)$decision
xgrid = expand.grid(X1 = px1, X2 = px2)
ygrid = predict(fit, xgrid)
plot(xgrid, col = as.numeric(ygrid), pch = 20, cex = .2)
points(x, col = y + 1, pch = 19)
contour(px1, px2, matrix(func, 69, 99), level = 0, add = TRUE)
contour(px1, px2, matrix(func, 69, 99), level = 0.5, add = TRUE, col = "blue", lwd = 2)Support vector machine
В практических примерах ниже показано как:
обучить модель классификатора на опорных векторах
обучить модель SVM с различными формами ядерной функции
строить ROC-кривые встроенными функциями R
Модели: SVM
Данные: Khan {ISLR}
Подробные комментарии к коду лабораторных см. в [1], глава 9.
library('e1071') # SVM
library('ROCR') # ROC-кривые
library('ISLR') # данные по экспрессии генов
my.seed <- 1
Классификатор на опорных векторах
Cгенерированные данные: два линейно неразделимых класса
# создаём наблюдения
set.seed(my.seed)
x <- matrix(rnorm(20*2), ncol = 2)
y <- c(rep(-1, 10), rep(1, 10))
x[y == 1, ] <- x[y == 1, ] + 1
# данные не разделяются линейно
plot(x, pch = 19, col = (3 - y))
# таблица с данными, отклик -- фактор
dat <- data.frame(x = x, y = as.factor(y))
# классификатор на опорных векторах с линейной границей
svmfit <- svm(y ~ ., data = dat, kernel = "linear", cost = 10, scale = FALSE)
# на графике опорные наблюдения показаны крестиками
plot(svmfit, dat)
# список опорных векторов
svmfit$index
## [1] 1 2 5 7 14 16 17
# сводка по модели
summary(svmfit)
##
## Call:
## svm(formula = y ~ ., data = dat, kernel = "linear", cost = 10,
## scale = FALSE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 10
## gamma: 0.5
##
## Number of Support Vectors: 7
##
## ( 4 3 )
##
##
## Number of Classes: 2
##
## Levels:
## -1 1
# уменьшаем штрафной параметр
svmfit <- svm(y ~ ., data = dat, kernel = "linear", cost = 0.1, scale = FALSE)
plot(svmfit, dat)
svmfit$index
## [1] 1 2 3 4 5 7 9 10 12 13 14 15 16 17 18 20
# делаем перекрёстную проверку, изменяя штраф (аргумент cost)
set.seed(1)
tune.out <- tune(svm, y ~ ., data = dat, kernel = "linear",
ranges = list(cost = c(0.001, 0.01, 0.1, 1, 5, 10, 100)))
summary(tune.out)
##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost
## 0.1
##
## - best performance: 0.1
##
## - Detailed performance results:
## cost error dispersion
## 1 1e-03 0.70 0.4216370
## 2 1e-02 0.70 0.4216370
## 3 1e-01 0.10 0.2108185
## 4 1e+00 0.15 0.2415229
## 5 5e+00 0.15 0.2415229
## 6 1e+01 0.15 0.2415229
## 7 1e+02 0.15 0.2415229
# лучшая модель -- с минимальной ошибкой
bestmod <- tune.out$best.model
summary(bestmod)
##
## Call:
## best.tune(method = svm, train.x = y ~ ., data = dat, ranges = list(cost = c(0.001,
## 0.01, 0.1, 1, 5, 10, 100)), kernel = "linear")
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 0.1
## gamma: 0.5
##
## Number of Support Vectors: 16
##
## ( 8 8 )
##
##
## Number of Classes: 2
##
## Levels:
## -1 1
# генерируем контрольные данные
xtest <- matrix(rnorm(20*2), ncol = 2)
ytest <- sample(c(-1,1), 20, rep = TRUE)
xtest[ytest == 1, ] <- xtest[ytest == 1, ] + 1
testdat <- data.frame(x = xtest, y = as.factor(ytest))
# делаем прогноз по лучшей модели
ypred <- predict(bestmod, testdat)
# матрица неточностей
table(predict = ypred, truth = testdat$y)
## truth
## predict -1 1
## -1 11 1
## 1 0 8
# прогноз по модели с cost = 0.01
svmfit <- svm(y ~ ., data = dat, kernel = "linear", cost = .01, scale = FALSE)
ypred <- predict(svmfit, testdat)
# матрица неточностей
table(predict = ypred, truth = testdat$y)
## truth
## predict -1 1
## -1 11 2
## 1 0 7
Сгенерированные данные: два линейно разделимых класса
# создаём наблюдения
x[y == 1, ] <- x[y == 1, ] + 0.5
plot(x, col = (y+5)/2, pch = 19)
# таблица с данными, отклик -- фактор
dat <- data.frame(x = x, y = as.factor(y))
# очень большой cost (маленький зазор, высокая точность классификации)
svmfit <- svm(y ~ ., data = dat, kernel = "linear", cost = 1e5)
summary(svmfit)
##
## Call:
## svm(formula = y ~ ., data = dat, kernel = "linear", cost = 1e+05)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 1e+05
## gamma: 0.5
##
## Number of Support Vectors: 3
##
## ( 1 2 )
##
##
## Number of Classes: 2
##
## Levels:
## -1 1
plot(svmfit, dat)
svmfit <- svm(y ~ ., data = dat, kernel = "linear", cost = 1)
summary(svmfit)
##
## Call:
## svm(formula = y ~ ., data = dat, kernel = "linear", cost = 1)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 1
## gamma: 0.5
##
## Number of Support Vectors: 7
##
## ( 4 3 )
##
##
## Number of Classes: 2
##
## Levels:
## -1 1
plot(svmfit,dat)
Машина опорных векторов
Сгенерированные данные: нелинейная граница между классами
# создаём наблюдения
set.seed(1)
x <- matrix(rnorm(200*2), ncol = 2)
x[1:100, ] <- x[1:100, ] + 2
x[101:150, ] <- x[101:150, ] - 2
y <- c(rep(1, 150), rep(2, 50))
# таблица с данными, отклик -- фактор
dat <- data.frame(x = x, y = as.factor(y))
plot(x, col = y, pch = 19)
# обучающая выборка
train <- sample(200, 100)
# SVM с радиальным ядром и маленьким cost
svmfit <- svm(y ~ ., data = dat[train, ], kernel = "radial",
gamma = 1, cost = 1)
plot(svmfit, dat[train, ])
summary(svmfit)
##
## Call:
## svm(formula = y ~ ., data = dat[train, ], kernel = "radial",
## gamma = 1, cost = 1)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: radial
## cost: 1
## gamma: 1
##
## Number of Support Vectors: 37
##
## ( 17 20 )
##
##
## Number of Classes: 2
##
## Levels:
## 1 2
# SVM с радиальным ядром и большим cost
svmfit <- svm(y ~ ., data = dat[train, ], kernel = "radial",
gamma = 1, cost = 1e5)
plot(svmfit, dat[train, ])
# перекрёстная проверка
set.seed(1)
tune.out <- tune(svm, y ~ ., data = dat[train, ], kernel = "radial",
ranges = list(cost = c(0.1, 1, 10, 100, 1000),
gamma = c(0.5, 1, 2, 3, 4)))
summary(tune.out)
##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost gamma
## 1 2
##
## - best performance: 0.12
##
## - Detailed performance results:
## cost gamma error dispersion
## 1 1e-01 0.5 0.27 0.11595018
## 2 1e+00 0.5 0.13 0.08232726
## 3 1e+01 0.5 0.15 0.07071068
## 4 1e+02 0.5 0.17 0.08232726
## 5 1e+03 0.5 0.21 0.09944289
## 6 1e-01 1.0 0.25 0.13540064
## 7 1e+00 1.0 0.13 0.08232726
## 8 1e+01 1.0 0.16 0.06992059
## 9 1e+02 1.0 0.20 0.09428090
## 10 1e+03 1.0 0.20 0.08164966
## 11 1e-01 2.0 0.25 0.12692955
## 12 1e+00 2.0 0.12 0.09189366
## 13 1e+01 2.0 0.17 0.09486833
## 14 1e+02 2.0 0.19 0.09944289
## 15 1e+03 2.0 0.20 0.09428090
## 16 1e-01 3.0 0.27 0.11595018
## 17 1e+00 3.0 0.13 0.09486833
## 18 1e+01 3.0 0.18 0.10327956
## 19 1e+02 3.0 0.21 0.08755950
## 20 1e+03 3.0 0.22 0.10327956
## 21 1e-01 4.0 0.27 0.11595018
## 22 1e+00 4.0 0.15 0.10801234
## 23 1e+01 4.0 0.18 0.11352924
## 24 1e+02 4.0 0.21 0.08755950
## 25 1e+03 4.0 0.24 0.10749677
# матрица неточностей для прогноза по лучшей модели
table(true = dat[-train, "y"],
pred = predict(tune.out$best.model, newdata = dat[-train, ]))
## pred
## true 1 2
## 1 74 3
## 2 7 16
ROC-кривые
# функция построения ROC-кривой: pred -- прогноз, truth -- факт
rocplot <- function(pred, truth, ...){
predob = prediction(pred, truth)
perf = performance(predob, "tpr", "fpr")
plot(perf,...)}
# последняя оптимальная модель
svmfit.opt <- svm(y ~ ., data = dat[train, ],
kernel = "radial", gamma = 2, cost = 1, decision.values = T)
# количественные модельные значения, на основе которых присваивается класс
fitted <- attributes(predict(svmfit.opt, dat[train, ],
decision.values = TRUE))$decision.values
# график для обучающей выборки
par(mfrow = c(1, 2))
rocplot(fitted, dat[train, "y"], main = "Training Data")
# более гибкая модель (gamma выше)
svmfit.flex = svm(y ~ ., data = dat[train, ], kernel = "radial",
gamma = 50, cost = 1, decision.values = T)
fitted <- attributes(predict(svmfit.flex, dat[train, ],
decision.values = T))$decision.values
rocplot(fitted, dat[train,"y"], add = T, col = "red")
# график для тестовой выборки
fitted <- attributes(predict(svmfit.opt, dat[-train, ],
decision.values = T))$decision.values
rocplot(fitted, dat[-train, "y"], main = "Test Data")
fitted <- attributes(predict(svmfit.flex, dat[-train, ],
decision.values = T))$decision.values
rocplot(fitted, dat[-train, "y"], add = T, col = "red")
SVM с несколькими классами
# генерируем данные
set.seed(1)
x <- rbind(x, matrix(rnorm(50*2), ncol = 2))
y <- c(y, rep(0, 50))
x[y == 0, 2] <- x[y == 0, 2] + 2
dat <- data.frame(x = x, y = as.factor(y))
# график и модель по методу "один против одного"
par(mfrow = c(1, 1))
plot(x, col = (y + 1))
svmfit = svm(y ~ ., data = dat, kernel = "radial", cost = 10, gamma = 1)
plot(svmfit, dat)
Анализ данных по уровню экспрессии генов
# данные по образцам тканей четырёх типов саркомы
names(Khan)
## [1] "xtrain" "xtest" "ytrain" "ytest"
dim(Khan$xtrain) # обучающая выборка, предикторы
## [1] 63 2308
dim(Khan$xtest) # тестовая выборка, предикторы
## [1] 20 2308
length(Khan$ytrain) # обучающая выборка, отклик
## [1] 63
length(Khan$ytest) # тестовая выборка, отклик
## [1] 20
table(Khan$ytrain)
##
## 1 2 3 4
## 8 23 12 20
table(Khan$ytest)
##
## 1 2 3 4
## 3 6 6 5
dat <- data.frame(x = Khan$xtrain, y = as.factor(Khan$ytrain))
# SVM с линейным ядром
out <- svm(y ~ ., data = dat, kernel = "linear", cost = 10)
summary(out)
##
## Call:
## svm(formula = y ~ ., data = dat, kernel = "linear", cost = 10)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 10
## gamma: 0.0004332756
##
## Number of Support Vectors: 58
##
## ( 20 20 11 7 )
##
##
## Number of Classes: 4
##
## Levels:
## 1 2 3 4
# матрица неточностей
table(out$fitted, dat$y)
##
## 1 2 3 4
## 1 8 0 0 0
## 2 0 23 0 0
## 3 0 0 12 0
## 4 0 0 0 20
# тестовые данные
dat.te <- data.frame(x = Khan$xtest, y = as.factor(Khan$ytest))
# прогноз на тестовой выборке
pred.te <- predict(out, newdata = dat.te)
# матрица неточностей
table(pred.te, dat.te$y)Sources
https://www.datacamp.com/community/tutorials/support-vector-machines-r