Как и обещал, продолжаю публикацию разборов задач, которые я прорешал за время работы с ребятами из MLClass.ru. В этот раз мы разберем метод главных компонент на примере известной задачи распознавания цифр Digit Recognizer с платформы Kaggle. Статья будет полезна новичкам, которые еще только начинают изучать анализ данных. Кстати, еще не поздно записаться на курс Прикладной анализ данных, получив возможность максимально быстро прокачаться в данной области.
Вступление
Данная работа является естественным продолжением исследования, изучающего зависимость качества модели от размера выборки. В ней для была показана возможность уменьшения количества используемых объектов в обучающей выборке с целью получения приемлемых результатов в условиях ограниченных вычислительных ресурсов. Но, кроме количества объектов, на размер данных влияет и количество используемых признаков. Рассмотрим эту возможность на тех же данных. Используемые данные были подробно изучены в предыдущей работе, поэтому просто загрузим тренировочную выборку в R.
library(readr)
library(caret)
library(ggbiplot)
library(ggplot2)
library(dplyr)
library(rgl)
data_train <- read_csv("train.csv")
##
|================================================================================| 100% 73 MB
Как мы уже знаем данные имеют 42000 объектов и 784 признака, представляющие собой значение яркости каждого из пикселей составляющего изображение цифры. Разобъём выборку на тренировочную и тестовую в соотношении 60/40.
set.seed(111)
split <- createDataPartition(data_train$label, p = 0.6, list = FALSE)
train <- slice(data_train, split)
test <- slice(data_train, -split)
Теперь удалим признаки, имеющие константное значение.
zero_var_col <- nearZeroVar(train, saveMetrics = T)
train <- train[, !zero_var_col$nzv]
test <- test[, !zero_var_col$nzv]
dim(train)
## [1] 25201 253
В итоге осталось 253 признака.
Теория
Метод главных компонент (PCA) преобразует базовые признаки в новые, каждый из которых является линейной комбинацией изначальных таким образом, что разброс данных (то есть среднеквадратичное отклонение от среднего значения) вдоль них максимален. Метод применяется для визуализации данных и для уменьшения размерности данных (сжатия).
PCA
Для большей наглядности случайным образом отберём из тренировочной выборки 1000 объектов и изобразим их в пространстве первых двух признаков.
train_1000 <- train[sample(nrow(train), size = 1000),]
ggplot(data = train_1000, aes(x = pixel152, y = pixel153, color = factor(label))) + geom_point()
Очевидно, что объекты перемешаны и выделить среди них группы объектов принадлежащих одному классу проблематично. Проведём преобразование данных по методу главных компонент и изобразим в пространстве первых двух компонент. Замечу, что компоненты расположены в убывающем порядке в зависимости от разброса данных, который приходится вдоль них.
pc <- princomp(train_1000[, -1], cor=TRUE, scores=TRUE)
ggbiplot(pc, obs.scale = 1, var.scale = 1, groups = factor(train_1000$label),
ellipse = TRUE, circle = F, var.axes = F) +
scale_color_discrete(name = '') +
theme(legend.direction = 'horizontal', legend.position = 'top')
Очевидно, что даже в пространстве всего лишь двух признаков уже можно выделить явные группы объектов. Теперь рассмотрим те же данные, но уже в пространстве первых трёх компонент.
plot3d(pc$scores[,1:3], col= train_1000$label + 1, size = 0.7, type = "s")
Выделение различных классов ещё больше упростилось. Теперь выберем количество компонент, которое будем использовать для дальнейшей работы. Для этого посмотрим на соотношение дисперсии и количество компонент объясняющие её, но уже используя всю тренировочную выборку.
pc <- princomp(train[, -1], cor=TRUE, scores=TRUE)
variance <- pc$sdev^2/sum(pc$sdev^2)
cumvar <- cumsum(variance)
cumvar <- data.frame(PC = 1:252, CumVar = cumvar)
ggplot(data = cumvar, aes(x = PC, y = CumVar)) + geom_point()
variance <- data.frame(PC = 1:252, Var = variance*100)
ggplot(data = variance[1:10,], aes(x = factor(PC), y = Var)) + geom_bar(stat = "identity")
sum(variance$Var[1:70])
## [1] 92.69142
Для того, чтобы сохранить более 90 процентов информации, содержащейся в данных достаточно всего лишь 70 компонент, т.е. мы от 784 признаков пришли к 70 и, при этом, потеряли менее 10 процентов вариации данных!
Преобразуем тренировочную и тестовую выборки в пространство главных компонент.
train <- predict(pc) %>% cbind(train$label, .) %>% as.data.frame(.) %>% select(1:71)
colnames(train)[1]<- "label"
train$label <- as.factor(train$label)
test %<>% predict(pc, .) %>% cbind(test$label, .) %>% as.data.frame(.) %>% select(1:71)
colnames(test)[1]<- "label"
Для выбора параметров моделей я использую пакет caret, предоставляющий возможность выполнять параллельные вычисления, используя многоядерность современных процессоров.
library("doParallel")
cl <- makePSOCKcluster(2)
registerDoParallel(cl)
KNN
Теперь приступим к созданию предсказывающих моделей используя преобразованные данные. Создадим первую модель используя метод k ближайших соседей (KNN). В этой модели есть только один параметр — количество ближайших объектов, используемых для классификации объекта. Подбирать этот параметр будем с помощью десятикратной кросс-проверки (10-fold cross-validation (CV)) с разбиением выборки на 10 частей. Оценка производится на случайно отобранной части изначальных объектов. Для оценки качества моделей будем использовать метрику Accuracy, представляющий собой процент точно предсказанных классов объектов.
set.seed(111)
train_1000 <- train[sample(nrow(train), size = 1000),]
ctrl <- trainControl(method="repeatedcv",repeats = 3)
Для начала определим область поиска значений параметра.
knnFit <- train(label ~ ., data = train_1000, method = "knn", trControl = ctrl,tuneLength = 20)
knnFit
## k-Nearest Neighbors
##
## 1000 samples
## 70 predictor
## 10 classes: '0', '1', '2', '3', '4', '5', '6', '7', '8', '9'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 3 times)
## Summary of sample sizes: 899, 901, 900, 901, 899, 899, ...
## Resampling results across tuning parameters:
##
## k Accuracy Kappa Accuracy SD Kappa SD
## 5 0.8749889 0.8608767 0.03637257 0.04047629
## 7 0.8679743 0.8530101 0.03458659 0.03853048
## 9 0.8652707 0.8500155 0.03336461 0.03713965
## 11 0.8529954 0.8363199 0.03692823 0.04114777
## 13 0.8433141 0.8255274 0.03184725 0.03548771
## 15 0.8426833 0.8248052 0.04097424 0.04568565
## 17 0.8423694 0.8244683 0.04070299 0.04540152
## 19 0.8340150 0.8151256 0.04291349 0.04788273
## 21 0.8263450 0.8065723 0.03914363 0.04369889
## 23 0.8200042 0.7995067 0.03872017 0.04320466
## 25 0.8156764 0.7946582 0.03825163 0.04269085
## 27 0.8093227 0.7875839 0.04299301 0.04799252
## 29 0.8010018 0.7783100 0.04252630 0.04747852
## 31 0.8019849 0.7794036 0.04327120 0.04827493
## 33 0.7963572 0.7731147 0.04418378 0.04930341
## 35 0.7936906 0.7701616 0.04012802 0.04478789
## 37 0.7889930 0.7649252 0.04163075 0.04644193
## 39 0.7863463 0.7619669 0.03947693 0.04404655
## 41 0.7829758 0.7582087 0.03482612 0.03889550
## 43 0.7796388 0.7544879 0.03745359 0.04179976
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was k = 5.
Теперь сократим её и получим точное значение.
grid <- expand.grid(k=2:5)
knnFit <- train(label ~ ., data = train_1000, method = "knn", trControl = ctrl, tuneGrid=grid)
knnFit
## k-Nearest Neighbors
##
## 1000 samples
## 70 predictor
## 10 classes: '0', '1', '2', '3', '4', '5', '6', '7', '8', '9'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 3 times)
## Summary of sample sizes: 900, 901, 901, 899, 901, 899, ...
## Resampling results across tuning parameters:
##
## k Accuracy Kappa Accuracy SD Kappa SD
## 2 0.8699952 0.8553199 0.03055544 0.03402108
## 3 0.8799832 0.8664399 0.02768544 0.03082014
## 4 0.8736591 0.8593777 0.02591618 0.02888557
## 5 0.8726753 0.8582703 0.02414173 0.02689738
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was k = 3.
Наилучший показатель модель имеет при значении параметра k равному 3. Используя это значение получим предсказание на тестовых данных. Построим Confusion Table и вычислим Accuracy.
library(class)
prediction_knn <- knn(train, test, train$label, k=3)
table(test$label, prediction_knn)
## prediction_knn
## 0 1 2 3 4 5 6 7 8 9
## 0 1643 0 6 1 0 1 2 0 0 0
## 1 0 1861 4 1 2 0 0 0 0 0
## 2 7 7 1647 3 0 0 1 11 0 0
## 3 1 0 9 1708 2 19 4 6 1 3
## 4 0 4 0 0 1589 0 10 7 0 6
## 5 3 2 1 20 1 1474 13 0 6 2
## 6 0 0 0 1 2 3 1660 0 0 0
## 7 0 6 3 0 2 0 0 1721 0 13
## 8 0 1 0 11 1 16 12 4 1522 20
## 9 0 0 1 3 3 5 1 23 5 1672
sum(diag(table(test$label, prediction_knn)))/nrow(test)
## [1] 0.9820227
Random Forest
Вторая модель — это Random Forest. У этой модели будем выбирать параметр mtry — количество используемых признаков при получении каждого из используемых в ансамбле деревьев. Для выбора наилучшего значения данного параметра пойдём тем же путём, что и ранее.
rfFit <- train(label ~ ., data = train_1000, method = "rf", trControl = ctrl,tuneLength = 3)
rfFit
## Random Forest
##
## 1000 samples
## 70 predictor
## 10 classes: '0', '1', '2', '3', '4', '5', '6', '7', '8', '9'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 3 times)
## Summary of sample sizes: 901, 900, 900, 899, 902, 899, ...
## Resampling results across tuning parameters:
##
## mtry Accuracy Kappa Accuracy SD Kappa SD
## 2 0.8526986 0.8358081 0.02889351 0.03226317
## 36 0.8324051 0.8133909 0.03442843 0.03836844
## 70 0.8026823 0.7802912 0.03696172 0.04118363
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 2.
grid <- expand.grid(mtry=2:6)
rfFit <- train(label ~ ., data = train_1000, method = "rf", trControl = ctrl,tuneGrid=grid)
rfFit
## Random Forest
##
## 1000 samples
## 70 predictor
## 10 classes: '0', '1', '2', '3', '4', '5', '6', '7', '8', '9'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 3 times)
## Summary of sample sizes: 898, 900, 900, 901, 900, 898, ...
## Resampling results across tuning parameters:
##
## mtry Accuracy Kappa Accuracy SD Kappa SD
## 2 0.8553016 0.8387134 0.03556811 0.03967709
## 3 0.8615798 0.8457973 0.03102887 0.03458732
## 4 0.8669329 0.8517297 0.03306870 0.03690844
## 5 0.8739532 0.8595897 0.02957395 0.03296439
## 6 0.8696883 0.8548470 0.03203166 0.03568138
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 5.
Выбираем mtry равным 4 и получаем Accuracy на тестовых данных. Замечу, что пришлось обучать модель на части от доступных тренировочных данных, т.к. для использования всех данных требуется больше оперативной памяти. Но, как показано в предыдущей работе, это не сильно повлияет на конечный результат.
library(randomForest)
rfFit <- randomForest(label ~ ., data = train[sample(nrow(train), size = 15000),], mtry = 4)
prediction_rf<-predict(rfFit,test)
table(test$label, prediction_rf)
## prediction_rf
## 0 1 2 3 4 5 6 7 8 9
## 0 1608 0 6 3 4 1 20 2 9 0
## 1 0 1828 9 9 3 5 3 2 9 0
## 2 12 9 1562 16 15 5 6 19 31 1
## 3 12 1 26 1625 2 33 12 14 22 6
## 4 0 6 11 1 1524 0 22 7 6 39
## 5 12 3 3 48 12 1415 10 1 15 3
## 6 13 4 8 0 4 11 1623 0 3 0
## 7 3 14 25 2 13 3 0 1653 4 28
## 8 4 10 12 64 8 21 12 5 1428 23
## 9 4 4 10 39 38 10 0 39 7 1562
sum(diag(table(test$label, prediction_rf)))/nrow(test)
## [1] 0.9421989
SVM
И, наконец, Support Vector Machine. В этой модели будет использоваться Radial Kernel и подбираются уже два параметра: sigma (регуляризационный параметр) и C (параметр, определяющий форму ядра).
svmFit <- train(label ~ ., data = train_1000, method = "svmRadial", trControl = ctrl,tuneLength = 5)
svmFit
## Support Vector Machines with Radial Basis Function Kernel
##
## 1000 samples
## 70 predictor
## 10 classes: '0', '1', '2', '3', '4', '5', '6', '7', '8', '9'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 3 times)
## Summary of sample sizes: 901, 900, 898, 900, 900, 901, ...
## Resampling results across tuning parameters:
##
## C Accuracy Kappa Accuracy SD Kappa SD
## 0.25 0.7862419 0.7612933 0.02209354 0.02469667
## 0.50 0.8545924 0.8381166 0.02931921 0.03262332
## 1.00 0.8826064 0.8694079 0.02903226 0.03225475
## 2.00 0.8929180 0.8808766 0.02781461 0.03090255
## 4.00 0.8986322 0.8872208 0.02607149 0.02898200
##
## Tuning parameter 'sigma' was held constant at a value of 0.007650572
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were sigma = 0.007650572 and C = 4.
grid <- expand.grid(C = 4:6, sigma = seq(0.006, 0.009, 0.001))
svmFit <- train(label ~ ., data = train_1000, method = "svmRadial", trControl = ctrl,tuneGrid=grid)
svmFit
## Support Vector Machines with Radial Basis Function Kernel
##
## 1000 samples
## 70 predictor
## 10 classes: '0', '1', '2', '3', '4', '5', '6', '7', '8', '9'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 3 times)
## Summary of sample sizes: 901, 900, 900, 899, 901, 901, ...
## Resampling results across tuning parameters:
##
## C sigma Accuracy Kappa Accuracy SD Kappa SD
## 4 0.006 0.8943835 0.8824894 0.02999785 0.03335171
## 4 0.007 0.8970537 0.8854531 0.02873482 0.03194698
## 4 0.008 0.8984139 0.8869749 0.03068411 0.03410783
## 4 0.009 0.8990838 0.8877269 0.03122154 0.03469947
## 5 0.006 0.8960834 0.8843721 0.03061547 0.03404636
## 5 0.007 0.8960703 0.8843617 0.03069610 0.03412880
## 5 0.008 0.8990774 0.8877134 0.03083329 0.03427321
## 5 0.009 0.8990838 0.8877271 0.03122154 0.03469983
## 6 0.006 0.8957534 0.8840045 0.03094360 0.03441242
## 6 0.007 0.8963971 0.8847267 0.03081294 0.03425451
## 6 0.008 0.8990774 0.8877134 0.03083329 0.03427321
## 6 0.009 0.8990838 0.8877271 0.03122154 0.03469983
##
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were sigma = 0.009 and C = 4.
library(kernlab)
svmFit <- ksvm(label ~ ., data = train,type="C-svc",kernel="rbfdot",kpar=list(sigma=0.008),C=4)
prediction_svm <- predict(svmFit, newdata = test)
table(test$label, prediction_svm)
## prediction_svm
## 0 1 2 3 4 5 6 7 8 9
## 0 1625 0 5 1 0 3 13 0 6 0
## 1 1 1841 6 6 4 1 0 3 5 1
## 2 8 4 1624 5 7 1 1 13 11 2
## 3 2 0 18 1684 0 23 2 6 12 6
## 4 1 3 3 0 1567 0 9 7 5 21
## 5 8 3 2 24 6 1465 7 0 6 1
## 6 2 1 2 1 5 5 1649 0 1 0
## 7 3 8 15 3 3 0 0 1695 3 15
## 8 1 6 10 10 5 9 3 4 1530 9
## 9 3 1 5 13 14 9 0 21 3 1644
sum(diag(table(test$label, prediction_svm)))/nrow(test)
## [1] 0.9717245
Ансамбль моделей
Создадим четвёртую модель, которая представляет собой ансамбль из трёх моделей, созданных ранее. Эта модель предсказывает то значение, за которое «голосует» большинство из использованных моделей.
all_prediction <- cbind(as.numeric(levels(prediction_knn))[prediction_knn],
as.numeric(levels(prediction_rf))[prediction_rf],
as.numeric(levels(prediction_svm))[prediction_svm])
predictions_ensemble <- apply(all_prediction, 1, function(row) {
row %>% table(.) %>% which.max(.) %>% names(.) %>% as.numeric(.)
})
table(test$label, predictions_ensemble)
## predictions_ensemble
## 0 1 2 3 4 5 6 7 8 9
## 0 1636 0 5 1 0 1 8 0 2 0
## 1 1 1851 3 5 3 0 0 1 4 0
## 2 7 6 1636 3 6 0 0 11 7 0
## 3 6 0 14 1690 1 18 4 8 7 5
## 4 0 5 4 0 1573 0 12 6 2 14
## 5 5 1 2 21 5 1478 7 0 3 0
## 6 3 1 2 0 5 3 1651 0 1 0
## 7 1 11 12 2 1 0 0 1704 0 14
## 8 1 5 11 17 4 13 4 3 1514 15
## 9 4 2 4 21 11 5 0 20 1 1645
sum(diag(table(test$label, predictions_ensemble)))/nrow(test)
## [1] 0.974939
Итоги
На тестовой выборке получены следующие результаты:
| Model | Test Accuracy |
| KNN | 0.981 |
| Random Forest | 0.948 |
| SVM | 0.971 |
| Ensemble | 0.974 |
Лучший показатель Accuracy имеет модель использующая метод k ближайших соседей (KNN).
Оценка моделей на сайте Kaggle приведена в следующей таблице.
| Model | Kaggle Accuracy |
| KNN | 0.97171 |
| Random Forest | 0.93286 |
| SVM | 0.97786 |
| Ensemble | 0.97471 |
И лучшие результаты здесь у SVM.
Eigenfaces
Ну и напоследок, уже из чистого любопытства, посмотрим наглядно на произведённые методом главных компонент преобразования. Для этого, во-первых получим изображение цифр в первоначальном виде.
set.seed(100)
train_1000 <- data_train[sample(nrow(data_train), size = 1000),]
colors<-c('white','black')
cus_col<-colorRampPalette(colors=colors)
default_par <- par()
number_row <- 28
number_col <- 28
par(mfrow=c(5,5),pty='s',mar=c(1,1,1,1),xaxt='n',yaxt='n')
for(i in 1:25)
{
z<-array(as.matrix(train_1000)[i,-1],dim=c(number_row,number_col))
z<-z[,number_col:1]
image(1:number_row,1:number_col,z,main=train_1000[i,1],col=cus_col(256))
}
par(default_par)
И изображение этих же цифр, но уже после того, как мы использовали метод PCA и оставили первые 70 компонент. Получившиеся объекты принято называть eigenfaces
zero_var_col <- nearZeroVar(train_1000, saveMetrics = T)
train_1000_cut <- train_1000[, !zero_var_col$nzv]
pca <- prcomp(train_1000_cut[, -1], center = TRUE, scale = TRUE)
restr <- pca$x[,1:70] %*% t(pca$rotation[,1:70])
restr <- scale(restr, center = FALSE , scale=1/pca$scale)
restr <- scale(restr, center = -1 * pca$center, scale=FALSE)
restr <- as.data.frame(cbind(train_1000_cut$label, restr))
test <- data.frame(matrix(NA, nrow = 1000, ncol = ncol(train_1000)))
zero_col_number <- 1
for (i in 1:ncol(train_1000)) {
if (zero_var_col$nzv[i] == F) {
test[, i] <- restr[, zero_col_number]
zero_col_number <- zero_col_number + 1
}
else test[, i] <- train_1000[, i]
}
par(mfrow=c(5,5),pty='s',mar=c(1,1,1,1),xaxt='n',yaxt='n')
for(i in 1:25)
{
z<-array(as.matrix(test)[i,-1],dim=c(number_row,number_col))
z<-z[,number_col:1]
image(1:number_row,1:number_col,z,main=test[i,1],col=cus_col(256))
}
par(default_par)
В следующий раз мы рассмотрим одну из задач Text Mining'а, ну а пока можете присоединиться к курсу по Практическому анализу данных — рекомендую!
This entry passed through the Full-Text RSS service - if this is your content and you're reading it on someone else's site, please read the FAQ at http://ift.tt/jcXqJW.
Комментариев нет:
Отправить комментарий