とりあえずニューラルネットを使ってみる
今後ニューラルネット(以下、NN)について調べていくが、とりあえず、nnet()関数でモデルを作成してみる。
NNであるが、複雑な非線形モデルと言ってよい。説明変数と目的変数が線形な関係でないとき、線形モデル(重回帰、主成分回帰、PLSなど)ではうまく近似できないが、NNを用いると近似することができる。
データはオリジナルなものを使用。説明変数2つ、目的変数1つの簡単な例。
head( NNexample )
X1 X2 Y 1 0.0 0 39.87536 2 0.5 0 39.99583 3 1.0 0 40.11665 4 1.5 0 40.23488 5 2.0 0 40.34901 6 2.5 0 40.45954
str( NNexample )
'data.frame': 121 obs. of 3 variables: $ X1: num 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 ... $ X2: num 0 0 0 0 0 0 0 0 0 0 ... $ Y : num 39.9 40 40.1 40.2 40.3 ...
plot( NNexample )
scatterplot3dパッケージのscatterplot3d()関数を用い、3Dプロットしてみる。
library( scatterplot3d )
scatterplot3d( NNexample, angle=280 )
引数angleで眺める角度の調整が行える。
説明変数(X1、X2)と目的変数(Y)の関係は明らかに非線形である。
後の重ね描きで必要なのでオブジェクトを作成しておく。
plotNN <- scatterplot3d( NNexample, angle=280 )
データに線形モデルである重回帰モデルをあてはめてみる。
目的変数がYで、説明変数がX1とX2。
res_lr <- lm( Y~X1+X2, data=NNexample )
あてはめ結果とプロットを重ね描きしてみる。
scatterplot3d( NNexample,angle=280, main="重回帰であてはめ" )
ドットを重ねる。
plotNN$points3d(NNexample$X1, NNexample$X2, predict(res_lr), col="red", pch=20)
回帰モデルに関しては平面も重ねられるようなのでついでに。
plotNN$plane3d(res_lr, col="red")
このように、線形モデルではこのデータの現象を捉えることはできない。今回は1次式だったが、2次式をあてはめても同様である。
NNをあてはめてみる。nnetパッケージのnnet()関数を使用。
library(nnet)
関数は、nnet(式, データ, size=ユニット数, linout=出力層への関数(TRUE:logistic、FALSE:linear))と指定する。
引数lineoutは、出力層への最後の関数を線形とするかロジスティックとするかの選択のよう。判別の問題を扱うときはロジスティックとするのだと思われる。今回は、目的変数(Y)は連続なので線形とする。
引数sizeはユニット数であるが、これの最適数の決定方法に関して私は知らない。
今回は説明変数が2つなので、ユニット数が幾つでも大して変わらないような気がするし、3Dプロットを使ってあてはまりが正しいか簡単に確認できると思う。
まあ試しで、データを3:1で学習とテスト用に分け、学習データでユニットが1から15までのモデルを作成し、テストデータで評価してみる。
評価の指標は、予測残差の平方和をテストデータ数で割り一行あたりの平均としたもの。小さいほうが良い。
乱数(二項分布)を使い、0と1を約3:1の割合でフラグを立てる。
flag <- rbinom( nrow(NNexample), 1, 0.25)
NNexample2 <- data.frame( NNexample, flag)
フラグが0を学習データとする。
NNexampleTrain <- subset( NNexample2, flag=="0", c(X1,X2,Y) )
head( NNexampleTrain )
X1 X2 Y 1 0.0 0 39.87536 2 0.5 0 39.99583 4 1.5 0 40.23488 5 2.0 0 40.34901 8 3.5 0 40.68306 9 4.0 0 40.80690
nrow( NNexampleTrain )
[1] 92
フラグが1をテストデータとする。
NNexampleTest <- subset( NNexample2, flag=="1", c(X1,X2,Y) )
head( NNexampleTest )
X1 X2 Y 3 1.0 0.0 40.11665 6 2.5 0.0 40.45954 7 3.0 0.0 40.56929 14 1.0 0.5 42.04056 21 4.5 0.5 42.60309 27 2.0 1.0 42.64893
nrow( NNexampleTest )
[1] 29
各ユニット数での評価結果を格納するオブジェクト。
SSEMean <- 1:15
各ユニット数での結果を繰り返し計算。
predict()関数は学習データ自身の予測値だけでなく、テストデータへの予測値も返せるようである。
for(i in 1:15) {
res_temp <- nnet(Y~X1+X2, data=NNexampleTrain, size=i, linout=TRUE)
pred_temp <- predict( res_temp, NNexampleTest )
SSEMean[i] <- sum( (NNexampleTest$Y - pred_temp)^2 ) / nrow(NNexampleTest)
}
結果。
data.frame(SSEMean)
SSEMean 1 1.529610362 2 1.529610923 3 1.529614121 4 0.333588701 5 0.006297643 6 0.002305100 7 0.006254313 8 0.008703751 9 0.605880589 10 1.530035349 11 0.003378762 12 0.009971504 13 1.529620984 14 0.049907408 15 0.016740579
数回、以上のプログラム(データの分割から各ユニット数でのモデル作成)を試してみた。
ユニット数が5以上もあれば、安定して小さな値になるようである。まあ、今回のデータは外れ値もないきれいなデータなので、もちろん一般的な結論でもなんでもない。
では、ユニット数は5としてモデルを作成する。
res_nn <- nnet(Y~X1+X2, data=NNexample, size=5, linout=TRUE)
あてはめ結果とプロットを重ね描きしてみる。
scatterplot3d(NNexample,angle=280, main="ニューラルネットであてはめ")
ドットを重ねる。
plotNN$points3d(NNexample$X1, NNexample$X2, predict(res_nn), col="blue", pch=20)
良くあてはまっている。
残差を確認してみる。
重回帰とニューラルで、横軸を予測値、縦軸を実測値でプロットする。
par(mfrow=c(1,2))
plot(predict(res_lr), NNexample$Y, main="重回帰", xlab="予測", ylab="実測")
abline(a=0, b=1, col="red", lwd=2)
plot(predict(res_nn), NNexample$Y, main="ニューラル", xlab="予測", ylab="実測")
abline(a=0, b=1, col="red", lwd=2)
重回帰は使えないであろうことが分かり、ニューラルネットの有効性が分かる。
回帰モデルのパラメータ推定に用いられる通常の最小二乗法は「SSE を最小化」する。
NNのパラメータ推定ではペナルティ付き最小二乗法が用いられる。「SSE + λ *(パラメータの二乗和)を最小化」する。
SSEは残差平方和。λ(λ> 0)はペナルティ。
λは正の値なので、ペナルティ付き最小二乗法は、パラメータが大きく推定されるのを避ける手法だということが分かる。
λが大きいと、推定されるパラメータの絶対値は全体的に小さくなり、λが小さいと、推定されるパラメータの絶対値は全体的に大きくなる。
推定するパラメータの絶対値が大きすぎるとオーバーフィットの危険性が増す。NNモデルはあらゆる応答を近似することができ、オーバーフィットの危険性を心配しなくてはいけないので、ペナルティ付き最小二乗法はこれの回避のために導入されていることになる。
パラメータの二乗和は、中間層の数やユニットの数を増やすと増えるので、それに掛け合わされているλは、それ自身の値に意味はないのではかと思う。
λの決定は、λの値を変化させながらクロスバリデーションを行い、予測誤差が最小になる値に決定されるのだと思われる。
λに関してもそうだが、nnet()の出力などいろいろと不明なので、後々調べて行く。