Математическое моделирование. Практика 9
Трапезникова Анастасия
28 апреля, 2021
Упражнение 9
Необходимо построить модель на основе SVM для указанной в варианте зависимой переменной.
Данные взять из упражнения №3
Для модели:
- Отложить 25% наблюдений в тестовую выборку (ядро генератора случайных чисел указано в варианте к упражнению №3).
- На обучающей выборке (оставшихся 75% наблюдений) сравнить несколько видов ядер SVM по точности модели (AUC) методом сеточного поиска.
- Для оптимальной формы ядерной функции на обучающей выборке подобрать оптимальное значение настроечных параметров по минимальной ошибке с перекрёстной проверкой (функция tune).
- Подогнать лучшую модель на всей обучающей выборке. Построить ROC-кривую и рассчитать матрицу неточностей, чувствительность и специфичность.
- Сделать прогноз по лучшей модели на тестовую выборку, оценить его качество точность по матрице неточностей, чувствительность и специфичность, построить ROC-кривую.
- Сравнить результаты, которые дал SVM, с результатами, полученными в упражнении 3. Какой из методов оказался лучше?
Варианты
| Вариант | Ядро для set.seed() |
Данные | Зависимая переменная | Объясняющие переменные | Методы |
|---|---|---|---|---|---|
| 9 | 345 | Default{ISLR} – долги по кредитным картам | default (Yes – наличие признака, No – отсутствие) |
все остальные | LDA, QDA |
Как сдавать: прислать на почту преподавателя ссылки:
* на html-отчёт с видимыми блоками кода (блоки кода с параметром echo = T), размещённый на rpubs.com.
* на код, генерирующий отчёт, в репозитории на github.com. В текст отчёта включить постановку задачи и ответы на вопросы задания.
Решение
Подключаем необходимые пакеты и набор данных Default
library('ISLR')
library('MASS')
library('e1071') # SVM
data(Default)Зададим ядро генератора случайных чисел и объём обучающей выборки.
my.seed <- 345
train.percent <- 0.75Исходные данные: набор Default
# ?Default # справка по набору данных
head(Default)## default student balance income
## 1 No No 729.5265 44361.625
## 2 No Yes 817.1804 12106.135
## 3 No No 1073.5492 31767.139
## 4 No No 529.2506 35704.494
## 5 No No 785.6559 38463.496
## 6 No Yes 919.5885 7491.559str(Default)## 'data.frame': 10000 obs. of 4 variables:
## $ default: Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ student: Factor w/ 2 levels "No","Yes": 1 2 1 1 1 2 1 2 1 1 ...
## $ balance: num 730 817 1074 529 786 ...
## $ income : num 44362 12106 31767 35704 38463 ...Отбираем наблюдения в обучающую выборку
set.seed(my.seed)
inTrain <- sample(seq_along(Default$default),nrow(Default)*train.percent)
df_train <- Default[inTrain, ]
df_test <- Default[-inTrain, ]На обучающей выборке (оставшихся 75% наблюдений) сравнить несколько видов ядер SVM по точности модели (AUC) методом сеточного поиска.
Сеточный поиск:
# обучающие данные
dat.tr <- data.frame(x = df_train[,-1], y = df_train[,1])
# тестовые данные
dat.te <- data.frame(x = df_test[,-1], y = df_test[,1])
# параметры алгоритма
kernel.grid <- c('linear', 'polynomial')
cost.grid <- seq(1, 20, by = 0.5)
AUC <- matrix(0, length(kernel.grid), length(cost.grid))
colnames(AUC) <- paste0('cost = ', cost.grid)
rownames(AUC) <- paste0('kernel = ', kernel.grid)
# SVM
for (i in 1:length(kernel.grid)) {
print(paste0('Starting ', kernel.grid[i], ' kernel'))
for (j in 1:length(cost.grid)) {
out <- svm(y ~ ., data = dat.tr, kernel = kernel.grid[i], cost = cost.grid[j])
# прогноз на тестовой выборке
pred.te <- predict(out, newdata = dat.te)
# матрица неточностей
tbl <- table(pred.te, dat.te$y)
AUC[i, j] <- sum(diag(tbl)) / sum(tbl)
}
}## [1] "Starting linear kernel"
## [1] "Starting polynomial kernel"AUC## cost = 1 cost = 1.5 cost = 2 cost = 2.5 cost = 3 cost = 3.5
## kernel = linear 0.9648 0.9648 0.9648 0.9648 0.9648 0.9648
## kernel = polynomial 0.9700 0.9704 0.9700 0.9700 0.9700 0.9700
## cost = 4 cost = 4.5 cost = 5 cost = 5.5 cost = 6 cost = 6.5
## kernel = linear 0.9648 0.9648 0.9648 0.9648 0.9648 0.9648
## kernel = polynomial 0.9704 0.9704 0.9704 0.9704 0.9708 0.9704
## cost = 7 cost = 7.5 cost = 8 cost = 8.5 cost = 9 cost = 9.5
## kernel = linear 0.9648 0.9648 0.9648 0.9648 0.9648 0.9648
## kernel = polynomial 0.9704 0.9704 0.9704 0.9704 0.9704 0.9704
## cost = 10 cost = 10.5 cost = 11 cost = 11.5 cost = 12
## kernel = linear 0.9648 0.9648 0.9648 0.9648 0.9648
## kernel = polynomial 0.9704 0.9704 0.9704 0.9704 0.9704
## cost = 12.5 cost = 13 cost = 13.5 cost = 14 cost = 14.5
## kernel = linear 0.9648 0.9648 0.9648 0.9648 0.9648
## kernel = polynomial 0.9704 0.9704 0.9704 0.9704 0.9704
## cost = 15 cost = 15.5 cost = 16 cost = 16.5 cost = 17
## kernel = linear 0.9648 0.9648 0.9648 0.9648 0.9648
## kernel = polynomial 0.9704 0.9704 0.9704 0.9708 0.9704
## cost = 17.5 cost = 18 cost = 18.5 cost = 19 cost = 19.5
## kernel = linear 0.9648 0.9648 0.9648 0.9648 0.9648
## kernel = polynomial 0.9704 0.9704 0.9704 0.9704 0.9704
## cost = 20
## kernel = linear 0.9648
## kernel = polynomial 0.9704# View(t(AUC))
# Максимальное значение AUC (area under ROC curve) равно 0.9708 при значении параметров:
# kernel = 'polynomial'
# cost = 6Для оптимальной формы ядерной функции на обучающей выборке подобрать оптимальное значение настроечных параметров по минимальной ошибке с перекрёстной проверкой (функция tune).
# делаем перекрёстную проверку, изменяя штраф (аргумент cost)
set.seed(my.seed)
tune.out <- tune(svm, y ~ ., data = dat.tr,
kernel = "polynomial",
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
## 10
##
## - best performance: 0.0268
##
## - Detailed performance results:
## cost error dispersion
## 1 1e-03 0.03266667 0.007200823
## 2 1e-02 0.03160000 0.007491229
## 3 1e-01 0.02826667 0.007484634
## 4 1e+00 0.02693333 0.006645894
## 5 5e+00 0.02693333 0.006993471
## 6 1e+01 0.02680000 0.006926778
## 7 1e+02 0.02706667 0.006741799# лучшая модель -- с минимальной ошибкой
bestmod <- tune.out$best.model
summary(bestmod)##
## Call:
## best.tune(method = svm, train.x = y ~ ., data = dat.tr, ranges = list(cost = c(0.001,
## 0.01, 0.1, 1, 5, 10, 100)), kernel = "polynomial")
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: polynomial
## cost: 10
## degree: 3
## coef.0: 0
##
## Number of Support Vectors: 467
##
## ( 246 221 )
##
##
## Number of Classes: 2
##
## Levels:
## No Yes# делаем прогноз по лучшей модели
ypred <- predict(bestmod, dat.te)
# матрица неточностей
table(Прогноз = ypred, Факт = dat.te$y)## Факт
## Прогноз No Yes
## No 2407 69
## Yes 5 19# прогноз по модели с cost = 0.01
svmfit <- svm(y ~ ., data = dat.tr, kernel = "polynomial", cost = .01, scale = FALSE)
ypred <- predict(svmfit, dat.te)
# матрица неточностей
table(Прогноз = ypred, Факт = dat.te$y)## Факт
## Прогноз No Yes
## No 1292 24
## Yes 1120 64Подогнать лучшую модель на всей обучающей выборке. Построить ROC-кривую и рассчитать матрицу неточностей, чувствительность и специфичность. Сделать прогноз по лучшей модели на тестовую выборку, оценить его качество точность по матрице неточностей, чувствительность и специфичность, построить ROC-кривую. 6 Сравнить результаты, которые дал SVM, с результатами, полученными в упражнении 3. Какой из методов оказался лучше?
ROC-кривые
# функция построения ROC-кривой: pred -- прогноз, truth -- факт
rocplot <- function(pred, truth, ...){
require(ROCR)
# вызов функции непосредственно из пакета ROCR
predob = ROCR::prediction(pred, truth)
perf = performance(predob, "tpr", "fpr")
plot(perf,...)}
# последняя оптимальная модель
svmfit.opt <- svm(y ~ ., data = dat.te, kernel = "polynomial", gamma = 2, cost = 1, probability = T)
# матрица неточностей на обучающей (p = 0.5)
table(Прогноз = predict(svmfit.opt, dat.tr), Факт = dat.tr$y)## Факт
## Прогноз No Yes
## No 7237 182
## Yes 18 63# прогноз вероятностей, на основе которых присваивается класс
fitted.prob <- predict(svmfit.opt, dat.tr, type = "prob",probability = TRUE)
fitted.prob <- attr(fitted.prob, "probabilities")[, 2]
# график для обучающей выборки
par(mfrow = c(1, 2))
# ROC-кривая для первой модели
rocplot(fitted.prob, dat.tr[, "y"], main = "Training Data")## Loading required package: ROCRpredob = prediction(fitted.prob, dat.tr[, "y"])
perf = performance(predob, "tpr", "fpr")
# более гибкая модель (gamma выше)
svmfit.flex = svm(y ~ ., data = dat.tr, kernel = "radial", gamma = 10, cost = 1, probability = T)
# матрица неточностей на обучающей (p = 0.5)
table(Прогноз = predict(svmfit.flex, dat.tr), Факт = dat.tr$y)## Факт
## Прогноз No Yes
## No 7231 164
## Yes 24 81# прогноз вероятностей, на основе которых присваивается класс
fitted.prob <- predict(svmfit.flex, dat.tr, type = "prob",probability = TRUE)
fitted.prob <- attr(fitted.prob, "probabilities")[, 2]
# ROC-кривая для второй модели
rocplot(fitted.prob, dat.tr[, "y"], add = T, col = "red")
# график для тестовой выборки
fitted.prob <- predict(svmfit.opt, dat.te[, ], type = "prob",probability = TRUE)
fitted.prob <- attr(fitted.prob, "probabilities")[, 2]
# матрица неточностей на тестовой (p = 0.5)
table(Прогноз = predict(svmfit.opt, dat.tr), Факт = dat.tr$y)## Факт
## Прогноз No Yes
## No 7237 182
## Yes 18 63rocplot(fitted.prob, dat.te[, "y"], main = "Test Data")
fitted.prob <- predict(svmfit.flex, dat.te, type = "prob",probability = TRUE)
fitted.prob <- attr(fitted.prob, "probabilities")[, 2]
# матрица неточностей на тестовой (p = 0.5)
table(Прогноз = predict(svmfit.flex, dat.te), Факт = dat.te$y)## Факт
## Прогноз No Yes
## No 2406 70
## Yes 6 18rocplot(fitted.prob, dat.te[, "y"], add = T, col = "red")
Источники
- Джеймс Г., Уиттон Д., Хасти Т., Тибширани Р. Введение в статистическое обучение с примерами на языке R / пер. с англ. С.Э. Мастицкого. – М.: ДМК Пресс, 2016 – 450 с. Репозиторий с примерами к книге на русском языке: https://github.com/ranalytics/islr-ru