Rによるニューラルネットワーク―neuralnet関数を利用した例―
個人的なメモとしてneuralnet関数の使い方についてまとめてみます。
関数の基本的な使い方については以下の本が詳しいのではないかと思います。
本記事ではニューラルネットワークを構築して回帰問題に取り組みます。
目次
データセット
UCI Machine Learning Repository: Concrete Compressive Strength Data Set
- データセットの中身は以下のとおり。
- Number of instances 1030
- Number of Attributes 9
- Attribute breakdown 8 quantitative input variables, and 1 quantitative output variable
- input variables: Cement, Slag, Ash, Water, Superplasticizer, Coarse Aggregate, Fine Aggregate, Age
- Output Variable: Concrete compressive strength
- コンクリートの圧縮強度に影響を及ぼすと考えられている8つの特徴量を使用して、圧縮強度をニューラルネットワークを構築して予測します。
使用するパッケージ
> library(readr) > library(neuralnet) > library(caret)
データの読込み
> data <- read.csv("C:/Users/admin/data/Concrete_data.csv") > str(data) 'data.frame': 1030 obs. of 9 variables: $ cement : num 540 540 332 332 199 ... $ slag : num 0 0 142 142 132 ... $ ash : num 0 0 0 0 0 0 0 0 0 0 ... $ water : num 162 162 228 228 192 228 228 228 228 228 ... $ superplasticizer: num 2.5 2.5 0 0 0 0 0 0 0 0 ... $ coarseagg : num 1040 1055 932 932 978 ... $ fineagg : num 676 676 594 594 826 ... $ age : num 28 28 270 365 360 90 365 28 28 28 ... $ strength : num 80 61.9 40.3 41 44.3 ...
特徴量の正規化
機械学習では特徴量を正規化して構築したモデルに投入することがお作法らしいです。特徴量を正規化することで変数の値（桁数）の大きさが予測精度に及ぼす影響を緩和しているということなんでしょうか？とりあえず、ここではお作法どおり各特徴量が区間[0, 1]の値を有するように正規化（置換）します。
> norm <- function(x){ > return((x-min(x)) / (max(x)-min(x))) > } > data.norm <- as.data.frame(lapply(data, norm)) > summary(data.norm$strength) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.0000 0.2664 0.4001 0.4172 0.5457 1.0000 # 正規化前の要約統計量 > summary(data$strength) Min. 1st Qu. Median Mean 3rd Qu. Max. 2.33 23.71 34.44 35.82 46.14 82.60
訓練データとテストデータの作成
Yeh(1998)と同様にデータの75%を訓練データに利用して、25%をテストデータに利用する。
smp.size <- floor(0.75*nrow(data.norm)) set.seed(1030) train.ind <- sample(seq_len(nrow(data.norm)), smp.size) train <- data.norm[train.ind,] test <- data.norm[-train.ind,]
neuralnet関数によるモデルの構築
neuralnet関数の引数のうちhiddenは隠れ層におけるユニット数を指定します。デフォルトでは1となっており、ここではとりあえずデフォルトどおり1を設定します。act.fctは活性化関数を指定する引数です。デフォルトではロジスティックシグモイド関数となっており、これもとりあえずデフォルトにしておきます。learningrateは学習率でありデフォルトはNULLですが、ここでは0.01に設定します。thresholdはイテレーション（学習）を止める閾値であり、デフォルトの0.01を設定します。
n <- names(train) formula.bpn <- as.formula(paste("strength~",paste(n[!n%in%c("strength")],collapse="+"))) bpn <- neuralnet(formula = formula.bpn, data = train, hidden = c(1), # 隠れ層におけるユニット数 act.fct="logistic", # 活性化関数にロジスティックシグモイド関数 learningrate = 0.01, # 学習率 threshold = 0.01, # 学習を止める閾値 linear.output=TRUE )
推定結果
neuralnet関数ではplot関数で推定結果を図示することができます。引数のrep="best"によって誤差（2乗誤差）が最も小さなものをプロットします。
plot(bpn, rep="best")
上図の左側の〇はコンクリートの圧縮強度を予測するために使用した8つの特徴量を示しており、隠れ層へ向かう黒線に書かれている数値は各特徴量のユニットが隠れ層のユニットに及ぼす影響度を表すウエイトを示しています。次に青色の線は各ステップにおけるバイアス項の値を示しており、これは回帰モデルにおける定数項みたいなものらしいです。
上図の下部にあるErrorとは損失関数（2乗誤差）の値を意味しており、ここでは正解ラベルの値と予測値の差の2乗和を2で割り算した値が5.053034であることを意味しています。この損失関数はneuralnet関数の引数（err.fct）で指定することができ、デフォルトでは"sse"（2乗誤差）が指定されています。sseの具体的な計算式は以下のように計算されます。
> bpn$err.fct structure(function (x, y) { 1/2 * (y - x)^2 }, type = "sse")
err.fctには他にも交差エントロピーを指定することも可能です。Stepsは学習が収束するまでのステップ数であり、ここでは4063回となっています。
上図で示されているErrorやStepsなどの値は以下のコマンドから取得できます。
> bpn$result.matrix 1 error 5.053033984221 reached.threshold 0.008933137079 steps 4063.000000000000 Intercept.to.1layhid1 -3.010942534867 cement.to.1layhid1 4.695003963273 slag.to.1layhid1 2.944619408533 ash.to.1layhid1 1.093491879204 water.to.1layhid1 -1.388814333100 superplasticizer.to.1layhid1 1.251426794179 coarseagg.to.1layhid1 0.274430837516 fineagg.to.1layhid1 0.473816307600 age.to.1layhid1 10.011792413222 Intercept.to.strength -0.013639063787 1layhid.1.to.strength 0.691225238214
次に、目的変数の予測に利用した個々の特徴量が予測に役立つものかを確認する手軽な手法として一般化ウエイトと特徴量の散布図を確認することがあげられます*1。この散布図は横軸に特徴量、縦軸に一般化ウエイトをプロットしており、散布図の縦軸（一般化ウエイト）の値がゼロ近辺に特徴量がプロットされているとき、その特徴量は目的変数の予測に役立っていないと判断できます。ちなみに、今回のデータによる散布図は以下のとおりです。
# 一般化ウエイトと特徴量の散布図を作成 par(mfrow=c(2,4)) gwplot(bpn,selected.covariate="cement") gwplot(bpn,selected.covariate="slag") gwplot(bpn,selected.covariate="ash") gwplot(bpn,selected.covariate="water") gwplot(bpn,selected.covariate="superplasticizer") gwplot(bpn,selected.covariate="coarseagg") gwplot(bpn,selected.covariate="fineagg") gwplot(bpn,selected.covariate="age")
上図をみると縦軸の値がゼロのところで横軸と水平に値が分布している図は見当たらないので、予測に利用した特徴量はこのまま含めてもよいと判断できそうです。
この一般化ウエイトの値を特徴量ごとに確認するには以下のようにします。
> print(head(bpn$generalized.weights[[1]])) [,1] [,2] [,3] [,4] [,5] 895 3.5893379669 2.2511662022 0.8359762738 -1.0617507576 0.9567177664 838 3.1450228972 1.9725000310 0.7324928850 -0.9303193164 0.8382880936 482 0.5892127611 0.3695433158 0.1372308468 -0.1742931709 0.1570513343 698 5.2561940936 3.2965874500 1.2241961033 -1.5548182178 1.4010088544 955 3.9371617592 2.4693148337 0.9169863209 -1.1646394179 1.0494282341 253 3.0937397406 1.9403361864 0.7205487598 -0.9151494074 0.8246188578 [,6] [,7] [,8] 895 0.20980280999 0.36223331770 7.654031159 838 0.18383185071 0.31739337136 6.706557998 482 0.03444047178 0.05946291357 1.256458120 698 0.30723333964 0.53045119811 11.208494084 955 0.23013369264 0.39733543609 8.395742909 253 0.18083426441 0.31221791335 6.597199982
訓練データによる予測とRMSE
ここではパラメタを未調整の状態で構築したニューラルネットワークによる予測値を計算し、その予測精度の指標であるRMSEを計算してみます。
> pred <- compute(bpn, test[,1:8]) > pred.strength <- pred$net.result > rmse.bpn <- sqrt(sum((test$strength - pred.strength)^2)/nrow(test)) > print(rmse.bpn) [1] 0.1215388084
RMSEの値は0.1215388084だと分かりました。これをチューニング後のモデルのRMSEと比べて、チューニングによって予測精度が高まっているかを以降で確認します。
neuralnet関数におけるパラメタのチューニング
caretパッケージのtrain関数を使用します。neuralnet関数を使用する場合は、隠れ層の数とそのユニット数がチューニングするパラメタとなります。
隠れ層の数やユニット数の組合せごとにRMSEの値を計算し、最も小さい値を有する組合せで再度ニューラルネットワークを構築します。
ここでは隠れ層が1つの場合と2つの場合を考え、1つ目の隠れ層のユニットが1から3つの場合と2つ目の隠れ層のユニットが0から3つの場合の合計12通りの組合せからRMSEを計算してみます*2。
model <- train(form=formula.bpn, data=train, method="neuralnet", tuneGrid = expand.grid(.layer1=c(1:3), .layer2=c(0:3), .layer3=c(0)), learningrate = 0.01, threshold = 0.01 ) > (model) Neural Network 772 samples 8 predictor No pre-processing Resampling: Bootstrapped (25 reps) Summary of sample sizes: 772, 772, 772, 772, 772, 772, ... Resampling results across tuning parameters: layer1 layer2 RMSE Rsquared 1 0 0.11673050316 0.6861291084 1 1 0.12088359121 0.6591349899 1 2 0.11698356020 0.6850522478 1 3 0.11656958967 0.6869186251 2 0 0.08534486470 0.8324254050 2 1 0.08344890188 0.8401588492 2 2 0.08380815321 0.8387608912 2 3 0.08371587352 0.8390426413 3 0 0.08594628891 0.8324012544 3 1 0.08061781486 0.8499824496 3 2 0.08074193620 0.8498935825 3 3 0.07930272754 0.8556659702 Tuning parameter 'layer3' was held constant at a value of 0 RMSE was used to select the optimal model using the smallest value. The final values used for the model were layer1 = 3, layer2 = 3 and layer3 = 0. > plot(model)
上図をみると、隠れ層1のユニットが3つで、隠れ層2のユニットが3つの場合のRMSEの値が約0.079と12通りの組合せの中で最も小さいことが分かります。これは、パラメタのチューニングによってモデルの精度が高くなったことを意味します。したがって、チューニング後のパラメタを用いて再度モデルを構築してみたいと思います。
チューニング後のパラメタを利用してモデルを再構築
bpn2 <- neuralnet(formula = formula.bpn, data = train, hidden = c(3,3), # 隠れ層1：ニューロン数2 act.fct="logistic", learningrate = 0.01, # learning rate threshold = 0.01, linear.output=TRUE ) plot(bpn2, rep="best")
上図をみると、最初のモデルと比べて、2乗誤差で表現されるErrorの値が減少していることが分かります。
テストデータによる予測とRMSE
結果は以下のとおり。
pred2 <- compute(bpn2, test[,1:8]) pred2.strength <- pred2$net.result rmse.bpn2 <- sqrt(sum((test$strength - pred2.strength)^2)/nrow(test)) print(rmse.bpn2) [1] 0.09096336435
RMSEの値は0.09096336435という結果になりました。
*1:一般化ウエイトについてはGünther and Fritsch(2010)を参照してください。
*2:本当はより多くの組合せを試して見たかったのですが、PCのスペックが低いので今回の12通りが限界でした
