Математическое моделирование, лабораторная 9
Мелихова И.С.
Математическое моделирование
Практика 9
Модели: SVM
Данные: Auto {ISLR}
Машины опорных векторов
Пример на данных по автомобилям: Auto {ISLR}.
Построим график разброса наблюдений в пространств предикторов, класс показан цветом (‘Yes’, ‘No’).
Зависимая переменная high.mpg, объсняющие: displacement, horsepower.
Переменная high.mpg – высокое значение mpg (сколько автомобиль проходит на галлоне топлива).
Перменная high.mpg принимает значение Yes, если mpg > 23, иначе - Nо.
library('e1071') # SVM
library('ROCR') # ROC-кривые
library('ISLR') # данные по экспрессии генов
library('GGally')
my.seed <- 1
attach(Auto)
head(Auto)str(Auto)## 'data.frame': 392 obs. of 9 variables:
## $ mpg : num 18 15 18 16 17 15 14 14 14 15 ...
## $ cylinders : num 8 8 8 8 8 8 8 8 8 8 ...
## $ displacement: num 307 350 318 304 302 429 454 440 455 390 ...
## $ horsepower : num 130 165 150 150 140 198 220 215 225 190 ...
## $ weight : num 3504 3693 3436 3433 3449 ...
## $ acceleration: num 12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
## $ year : num 70 70 70 70 70 70 70 70 70 70 ...
## $ origin : num 1 1 1 1 1 1 1 1 1 1 ...
## $ name : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...# новая переменная
high.mpg <- ifelse(mpg < 23, 'No', 'Yes')
# присоединяем к таблице данных
Auto <- data.frame(Auto, high.mpg)
high.mpg <- as.factor(high.mpg)
ggp <- ggpairs(Auto[, c('high.mpg', 'displacement', 'horsepower')],
mapping = ggplot2::aes(color = high.mpg))
print(ggp, progress = F)
# таблица с данными, отклик — фактор
dat <- data.frame(displacement, horsepower , high.mpg)
# обучающая выборка
train <- sample(1:nrow(dat), nrow(dat)/2)
# SVM с радиальным ядром и маленьким cost
svmfit <- svm(high.mpg ~ ., data = dat[train, ], kernel = "radial",
gamma = 1, cost = 1)
plot(svmfit, dat[train, ])
summary(svmfit)##
## Call:
## svm(formula = high.mpg ~ ., 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: 56
##
## ( 27 29 )
##
##
## Number of Classes: 2
##
## Levels:
## No Yes# SVM с радиальным ядром и большим cost
svmfit <- svm(high.mpg ~ ., data = dat[train, ], kernel = "radial",
gamma = 1, cost = 1e5)
plot(svmfit, dat[train, ])
summary(svmfit)##
## Call:
## svm(formula = high.mpg ~ ., data = dat[train, ], kernel = "radial",
## gamma = 1, cost = 1e+05)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: radial
## cost: 1e+05
## gamma: 1
##
## Number of Support Vectors: 52
##
## ( 22 30 )
##
##
## Number of Classes: 2
##
## Levels:
## No Yes# перекрёстная проверка
set.seed(my.seed)
tune.out <- tune(svm, high.mpg ~ ., 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 1
##
## - best performance: 0.09631579
##
## - Detailed performance results:
## cost gamma error dispersion
## 1 1e-01 0.5 0.10157895 0.05830794
## 2 1e+00 0.5 0.10157895 0.05333160
## 3 1e+01 0.5 0.10131579 0.04688266
## 4 1e+02 0.5 0.10131579 0.04688266
## 5 1e+03 0.5 0.11157895 0.05704048
## 6 1e-01 1.0 0.10157895 0.05830794
## 7 1e+00 1.0 0.09631579 0.04962463
## 8 1e+01 1.0 0.10131579 0.04688266
## 9 1e+02 1.0 0.11157895 0.05704048
## 10 1e+03 1.0 0.11657895 0.06697815
## 11 1e-01 2.0 0.10157895 0.05830794
## 12 1e+00 2.0 0.10131579 0.04688266
## 13 1e+01 2.0 0.10657895 0.05021116
## 14 1e+02 2.0 0.11157895 0.05704048
## 15 1e+03 2.0 0.11631579 0.07466762
## 16 1e-01 3.0 0.10157895 0.05830794
## 17 1e+00 3.0 0.10131579 0.04688266
## 18 1e+01 3.0 0.11157895 0.05704048
## 19 1e+02 3.0 0.11631579 0.07466762
## 20 1e+03 3.0 0.11684211 0.07933821
## 21 1e-01 4.0 0.11710526 0.07338337
## 22 1e+00 4.0 0.10131579 0.04688266
## 23 1e+01 4.0 0.11157895 0.05704048
## 24 1e+02 4.0 0.11631579 0.07466762
## 25 1e+03 4.0 0.09631579 0.07280385Построим матрицу неточностей для прогноза по лучшей модели и рассчитаем MSE.
# матрица неточностей для прогноза по лучшей модели
matrix <- table(true = dat[-train, "high.mpg"],
pred = predict(tune.out$best.model, newdata = dat[-train, ]))
bestmod <- tune.out$best.model
summary(bestmod)##
## Call:
## best.tune(method = svm, train.x = high.mpg ~ ., data = dat[train,
## ], ranges = list(cost = c(0.1, 1, 10, 100, 1000), gamma = c(0.5,
## 1, 2, 3, 4)), kernel = "radial")
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: radial
## cost: 1
## gamma: 1
##
## Number of Support Vectors: 56
##
## ( 27 29 )
##
##
## Number of Classes: 2
##
## Levels:
## No Yes#MSE
sum(diag(matrix))/sum(matrix) ## [1] 0.8979592ROC-кривые
# функция построения ROC-кривой: pred -- прогноз, truth -- факт
rocplot <- function(pred, truth, ...){
predob = prediction(pred, truth)
perf = performance(predob, "tpr", "fpr")
plot(perf,...)}
# последняя оптимальная модель
svmfit.opt <- svm(high.mpg ~ ., data = dat[train, ],
kernel = "radial", gamma = 2, cost = 0.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, "high.mpg"], main = "Training Data")
# более гибкая модель (gamma выше)
svmfit.flex = svm(high.mpg ~ ., 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,"high.mpg"], add = T, col = "red")
# график для тестовой выборки
fitted <- attributes(predict(svmfit.opt, dat[-train, ],
decision.values = T))$decision.values
rocplot(fitted, dat[-train, "high.mpg"], main = "Test Data")
fitted <- attributes(predict(svmfit.flex, dat[-train, ],
decision.values = T))$decision.values
rocplot(fitted, dat[-train, "high.mpg"], add = T, col = "red")
Построенные ROC-кривые показывают большое количество точных предсказаний, что говорит о правильности выбора объясняющих перменных.