Neural Network
Небожатко Екатерина
library(neuralnet)
library(MASS)
library(randomForest)
library(nnet)
data <- Boston
set.seed(500)
index <- sample(1:nrow(data),round(0.75*nrow(data)))
train <- data[index,]
test <- data[-index,]Регрессия
Линейная регрессия
lm.fit <- glm(medv~., data=train)
summary(lm.fit)##
## Call:
## glm(formula = medv ~ ., data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -14.9143 -2.8607 -0.5244 1.5242 25.0004
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 43.469681 6.099347 7.127 5.50e-12 ***
## crim -0.105439 0.057095 -1.847 0.065596 .
## zn 0.044347 0.015974 2.776 0.005782 **
## indus 0.024034 0.071107 0.338 0.735556
## chas 2.596028 1.089369 2.383 0.017679 *
## nox -22.336623 4.572254 -4.885 1.55e-06 ***
## rm 3.538957 0.472374 7.492 5.15e-13 ***
## age 0.016976 0.015088 1.125 0.261291
## dis -1.570970 0.235280 -6.677 9.07e-11 ***
## rad 0.400502 0.085475 4.686 3.94e-06 ***
## tax -0.015165 0.004599 -3.297 0.001072 **
## ptratio -1.147046 0.155702 -7.367 1.17e-12 ***
## black 0.010338 0.003077 3.360 0.000862 ***
## lstat -0.524957 0.056899 -9.226 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 23.26491)
##
## Null deviance: 33642 on 379 degrees of freedom
## Residual deviance: 8515 on 366 degrees of freedom
## AIC: 2290
##
## Number of Fisher Scoring iterations: 2pr.lm <- predict(lm.fit,test)
MSE.lm <- sum((pr.lm - test$medv)^2)/nrow(test)
MSE.lm## [1] 21.62976Нейронная сеть
Для нейронной сети необходимо стандартизировать данные.
scaled <- as.data.frame(scale(data))
train_ <- scaled[index,]
test_ <- scaled[-index,]
test.error <- NULL
train.error <- NULLcrossvalidate <- function(data,hidden_l=c(5))
{
scaled <- as.data.frame(scale(data))
cv.error <- NULL
k <- 5
for(j in 1:k)
{
print(k)
index <- sample(1:nrow(data),round(0.90*nrow(data)))
train.cv <- scaled[index,]
test.cv <- scaled[-index,]
nn <- neuralnet(f,data=train.cv,hidden=hidden_l,linear.output=T, stepmax = 10^6)
pr.nn <- compute(nn,test.cv[,1:13])
pr.nn <- pr.nn$net.result*sd(data$medv) + mean(data$medv)
test.cv.r <- (test.cv$medv)*sd(data$medv) + mean(data$medv)
cv.error[j] <- sum((test.cv.r - pr.nn)^2)/nrow(test.cv)
}
return(mean(cv.error))
}Используем нейронную сеть с одним слоем. Наиболее популярное правило выбирать количество нейронов между 1 и количеством признаков. Для выбора количества нейронов будем использовать 5-fold cross-validation. Выберем то количество, которое минимизирует ошибку.
n <- names(train)
f <- as.formula(paste("medv ~", paste(n[!n %in% "medv"], collapse = " + ")))
# set.seed(100)
# for(i in 1:13)
# {
# nn <- neuralnet(f,data=scaled,hidden=c(i),linear.output=T, stepmax = 10^6)
# print(i)
# train.error[i] <- sum(((as.data.frame(nn$net.result)*sd(data$medv) + mean(data$medv)) - (scaled$medv*sd(data$medv) + mean(data$medv)))^2)/nrow(scaled)
# test.error[i] <- crossvalidate(data,hidden_l=c(i))
# }
test.error <- c(14.197351092, 19.118568862, 15.323856797, 9.900850077, 13.013030110, 10.294646035, 18.273162987, 14.455016860, 17.336128854,16.800375573, 18.927995691, 12.791386809, 19.295322392)
train.error <- c(15.579924258, 11.583278409, 6.542277976, 7.211608826, 5.299934423 , 5.344445163, 3.685781822, 4.148210662, 3.280870741,2.679354351, 2.202861799, 2.153608806 , 1.886667069)Функция neuralnet используется для обучения нейронной сети. По непонятным причиным, она не принимает формулу в виде “y~ …”.
Важные аргуенты
hidden - вектор, определяющий количество нейроннов в скрытых слоях;
stepmax - максимальное количество шагов, в обучении нейронной сети;
algorithm - какой алгоритм будет использоваться для построения сети;
linear.output - логический. Когда TRUE - регрессия, FALSE - классификация;
act.fct - выбор функции активации. ’logistic’ или ’tanh.
Функция compute вычисляет результат по построенной нейронной сети для тестовой выборки.
plot(train.error,main='MSE vs hidden neurons',xlab="Hidden neurons",ylab='Train error MSE',type='l',col='red',lwd=2)
plot(test.error,main='MSE vs hidden neurons',xlab="Hidden neurons",ylab='Test error MSE',type='l',col='blue',lwd=2)
opt_n <- which(min(test.error) == test.error)
which(min(train.error) == train.error)## [1] 13opt_n## [1] 4Оптимальное количество нейронов - 4.
Сравним результаты линейной регрессии и нейронной сети.
nn <- neuralnet(f,data=train_,hidden=opt_n,linear.output=T, stepmax = 10^6)
pr.nn <- compute(nn,test_[,1:13])
pr.nn_ <- pr.nn$net.result*sd(data$medv) + mean(data$medv)
test.r <- (test_[,14])*sd(data$medv) + mean(data$medv)
par(mfrow=c(1,2))
plot(test$medv,pr.nn_,col='red',main='Real vs predicted NN',pch=18,cex=0.7)
abline(0,1,lwd=2)
legend('bottomright',legend='NN',pch=18,col='red', bty='n')
plot(test$medv,pr.lm,col='blue',main='Real vs predicted lm',pch=18, cex=0.7)
abline(0,1,lwd=2)
legend('bottomright',legend='LM',pch=18,col='blue', bty='n', cex=.95)
plot(test$medv,pr.nn_,col='red',main='Real vs predicted NN',pch=18,cex=0.7)
points(test$medv,pr.lm,col='blue',pch=18,cex=0.7)
abline(0,1,lwd=2)
legend('bottomright',legend=c('NN','LM'),pch=18,col=c('red','blue'))
Теперь посчитаем 5-fold CV error для регрессии.
library(boot)
set.seed(200)
k = 5
lm.fit <- glm(medv~.,data=data)
cv.lm <- cv.glm(data,lm.fit,K=k)$delta[1]
cv.lm## [1] 23.72094885-fold CV error для нейронной сети.
cv.nn <- test.error[opt_n]
cv.nn## [1] 9.900850077Классификация
Нейронная сеть
wines <- read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data", sep = ",")
colnames(wines) <- c("type", "alcohol", "malic", "ash", "alcalinity", "magnesium", "phenols", "flavanoids", "nonflavanoids","proanthocyanins", "color", "hue", "dilution", "proline")
head(wines)## type alcohol malic ash alcalinity magnesium phenols flavanoids
## 1 1 14.23 1.71 2.43 15.6 127 2.80 3.06
## 2 1 13.20 1.78 2.14 11.2 100 2.65 2.76
## 3 1 13.16 2.36 2.67 18.6 101 2.80 3.24
## 4 1 14.37 1.95 2.50 16.8 113 3.85 3.49
## 5 1 13.24 2.59 2.87 21.0 118 2.80 2.69
## 6 1 14.20 1.76 2.45 15.2 112 3.27 3.39
## nonflavanoids proanthocyanins color hue dilution proline
## 1 0.28 2.29 5.64 1.04 3.92 1065
## 2 0.26 1.28 4.38 1.05 3.40 1050
## 3 0.30 2.81 5.68 1.03 3.17 1185
## 4 0.24 2.18 7.80 0.86 3.45 1480
## 5 0.39 1.82 4.32 1.04 2.93 735
## 6 0.34 1.97 6.75 1.05 2.85 1450Создадим dummy переменные для переменной type, которую будем предсказывать.
data <- cbind(wines[, 2:14], class.ind(as.factor(wines$type)))
names(data) <- c(names(wines)[2:14],"l1","l2","l3")
data[, 1:13] <- data.frame(lapply(data[, 1:13], scale))Построим 2х слойную нейронную сеть с 9 и 5 нейронами в скрытых слоях соответственно.
set.seed(123)
index <- sample(1:nrow(data),round(0.75*nrow(data)))
train <- data[index,]
test <- data[-index,]
n <- names(data)
f <- as.formula(paste("l1 + l2 + l3 ~", paste(n[!n %in% c("l1","l2","l3")], collapse = " + ")))
nn <- neuralnet(f,data = train,hidden = c(9,5), linear.output = FALSE)
pr.nn <- compute(nn, test[, 1:13])
pr.nn_ <- pr.nn$net.result
original_values <- max.col(test[, 14:16])
pr.nn_2 <- max.col(pr.nn_)
test.nn <- mean(pr.nn_2 == original_values)
test.nn## [1] 0.9772727273
Neural Network
10-fold CV error.
set.seed(10)
k <- 10
outs <- NULL
proportion <- 0.95
for(i in 1:k)
{
index <- sample(1:nrow(train), round(proportion*nrow(train)))
train_cv <- train[index, ]
test_cv <- train[-index, ]
nn_cv <- neuralnet(f,
data = train_cv,
hidden = c(9,5),
linear.output = FALSE)
pr.nn <- compute(nn_cv, test_cv[, 1:13])
pr.nn_ <- pr.nn$net.result
original_values <- max.col(test_cv[, 14:16])
pr.nn_2 <- max.col(pr.nn_)
outs[i] <- mean(pr.nn_2 == original_values)
}
cv.nn <- mean(outs)
cv.nn## [1] 0.9857142857Random forest
set.seed(11)
wines$type <- as.factor(wines$type)
rf <- randomForest(type~., data = wines, mtry = 4, importance = TRUE)
test.rf <- mean(rf$predicted == wines$type)
test.rf## [1] 0.9775280899cv.rf <- 1-mean(rf$err.rate)
cv.rf## [1] 0.9771086488