ちょっと調べてみたらタイトルの件について言及してる記事があまり多くなかったので、ざっくり書いてみます。なお、この記事はid:shakezoさんの
へのオマージュです。というか、実は僕もこの記事を読んでから「多分Rなら専用の関数なんかもあるだろうし簡単にできるはず」と思って以前よりも積極的にやるようになったのでした（笑）。
総論：何で機械学習するのにチューニングが必要なの？
どんな機械学習でも、何かしらのチューニングパラメータを持っています。例えばソフトマージンSVMならマージンパラメータCがありますし、非線形ガウシアンカーネルSVMならさらにカーネルパラメータのσとかが入ります。SMO（逐次最大最適化）アルゴリズムを利用するのであれば、さらにさらにtoleranceとかも入ってきます。
しかも、ちょっといじってみればすぐ分かると思うんですがこの辺のチューニングパラメータを変えるだけで機械学習の結果はガンガン変わります。試しに2次元SVMとかで分離超平面を書いてみると、Cとσを変えるだけで全く変わってしまいます。
下の図2つは僕が以前独習用にMatlabで書いたSMOアルゴリズムによる非線形ガウシアンカーネルSVM*1の実行例なんですが。。。
ほぼ同じ学習データに対して異なるCとσの組み合わせで分離超平面をコンタープロットで描いてみると、こんなに違います。上は綺麗に汎化されていますが、下は明らかにオーバーフィッティングです。
ということで、機械学習するならパラメータチューニングは必須。そのRでのやり方を、代表的なパッケージを取り上げながらちょろっと見てみようと思います。僕の勉強も兼ねて。。。
{randomForest}パッケージの場合
ランダムフォレストの原理とか利点とかは全てid:shakezoさんの記事を読んで下さいということで、まずこのブログでも頻出の{randomForest}パッケージでやる場合のパラメータチューニングのところを取り上げます。
{randomForest}パッケージおさらい
面倒なので、以前の記事で使ったサンプルデータをそのまま転用します。その時同様とりあえずsample_dとかいう名前にしておきます。
既に何度か色々な記事で触れてますが、ランダムフォレストはrandomForest(){randomForest}関数で簡単にやれます。とりあえず何もいじらずにデフォルトのままやってみましょう。
> sample_d.rf<-randomForest(cv~.,sample_d) > print(sample_d.rf) Call: randomForest(formula = cv ~ ., data = sample_d) Type of random forest: classification Number of trees: 500 No. of variables tried at each split: 2 OOB estimate of error rate: 6.37% Confusion matrix: No Yes class.error No 1399 101 0.06733333 Yes 90 1410 0.06000000 > importance(sample_d.rf) MeanDecreaseGini a1 21.554057 a2 11.912213 a3 2.550909 a4 219.898301 a5 82.449264 a6 735.583208 a7 2.543989
まぁ見たまんまです。OOB estimate of error rateが6.37%と、この時点でもそれほど悪くない気はします*2。
tuneRF()でグリッドサーチをかけて最適化する
で、今回のお題なんですが。推定パラメータをグリッドサーチでチューニングして最適化した方が良いですよ！ということなので、やってみましょう。id:shakezoさんの記事にもあるように、
RandomForestの主要なパラメータは次の2つです。
- 作成する決定木の数
- １つ１つの決定木を作成する際に使用する特徴量の数
これらのうち、randomForest(){randomForest}関数では作成する決定木の数はntree、1つ1つの決定木を作成する際に使用する特徴量の数はmtryで指定できます。このうち、ntreeは学習結果のrandomForest.formulaクラスデータをplot()関数で図示することでntreeを増やすごとに収束度合いがどう変わるかを見ることができるので、事後の判定になります。
よって、ここでの焦点はmtry一択。これをグリッドサーチする方法は簡単、tuneRF(){randomForest}関数を使うだけです。
> sample_d.tune<-tuneRF(sample_d[,-8],sample_d[,8],doBest=T) mtry = 2 OOB error = 6.17% Searching left ... mtry = 1 OOB error = 9% -0.4594595 0.05 Searching right ... mtry = 4 OOB error = 6.63% -0.07567568 0.05
こんな感じで、doBestをTrueにしておくとmtryの候補を図示して選んでくれます。なおtuneRF()関数でもntreeを指定できますが、デフォルトで問題ないでしょう。
ということで、実際にmtry=2を指定してやり直してみました。
> sample_d.rf2<-randomForest(cv~.,sample_d,mtry=2) > print(sample_d.rf2) Call: randomForest(formula = cv ~ ., data = sample_d, mtry = 2) Type of random forest: classification Number of trees: 500 No. of variables tried at each split: 2 OOB estimate of error rate: 6.33% Confusion matrix: No Yes class.error No 1400 100 0.06666667 Yes 90 1410 0.06000000 > importance(sample_d.rf2) MeanDecreaseGini a1 21.687817 a2 12.408733 a3 2.508693 a4 196.222452 a5 76.332301 a6 776.395852 a7 2.674710
OOB estimate of error rateが6.33%とほんの少しですが、改善されました。ついでなので、ntreeも変えてもう一度やってみましょう。
> sample_d.rf3<-randomForest(cv~.,sample_d,mtry=2,ntree=2000) > print(sample_d.rf3) Call: randomForest(formula = cv ~ ., data = sample_d, mtry = 2, ntree = 2000) Type of random forest: classification Number of trees: 2000 No. of variables tried at each split: 2 OOB estimate of error rate: 6.3% Confusion matrix: No Yes class.error No 1403 97 0.06466667 Yes 92 1408 0.06133333 > importance(sample_d.rf3) MeanDecreaseGini a1 21.806920 a2 12.199403 a3 2.467463 a4 201.304439 a5 79.126512 a6 767.086724 a7 2.765361 > plot(sample_d.rf3)
OOB estimate of error rateが6.3%ぴったりまで落ちました（まぁほんのちょっとですが）。学習の収束状況をプロットしてみると、実はntree = 500ぐらいでもう十分に収束していて、2000まで増やす必要はなかったらしいということも分かります。なので、計算負荷を考えてntree = 500が妥当という結論に至りました。
・・・ということで、tuneRF()でグリッドサーチしてチューニングをかけられるというのを見てみました。実データだとチューニングなしだと全然error rateが落ちないよ～みたいなこともバンバン出てくると思うので、重要ですね。
{e1071}パッケージの場合
{e1071}は汎用パッケージなので色々ありますが、代表例としてLIBSVM*3のR実装であるsvm()関数でやってみます。モデルの推定自体はこんな感じです。
> sample_d.libsvm<-svm(cv~.,sample_d) > summary(sample_d.libsvm) Call: svm(formula = cv ~ ., data = sample_d) Parameters: SVM-Type: C-classification SVM-Kernel: radial cost: 1 gamma: 0.1428571 Number of Support Vectors: 468 ( 216 252 ) Number of Classes: 2 Levels: No Yes
svm(){e1071}に対してチューニングしたければ、tune.svm()関数を使います。
> t<-tune.svm(cv~.,data=sample_d) > summary(t) Error estimation of ‘svm’ using 10-fold cross validation: 0.06468537 > t$best.parameters dummyparameter 1 0 > t$best.performance [1] 0.06468537 > t$best.model Call: best.svm(x = cv ~ ., data = sample_d) Parameters: SVM-Type: C-classification SVM-Kernel: radial cost: 1 gamma: 0.1428571 Number of Support Vectors: 468
最初にデフォルトで勝手にチューニングさせてしまっているので、あまり面白くないですね（汗）。普通にirisデータとかでもっと遊んでみても良いだろうとは思います。
なお、こちらの説明はid:hoxo_mさんの「SVM のチューニングのしかた（２）」を参考にさせていただきました。というか、そちらの記事の方が
チューニングの手順としては、
- グリッドサーチで大雑把に検索する。
- 最適なパラメータがありそうなところを絞って再びグリッドサーチを行う。
という2段階のグリッドサーチを行います。
というように王道のチューニング手順をきちんと踏んで実践している（最適パラメータ探索にプロットも使っている）例を紹介しているので、お薦めです。
ちなみに{e1071}は他にもk最近傍・ニューラルネットワーク・決定木*4・ランダムフォレストといった手法を実装していて、その全てにtune.XXX()関数を備えているので、どの手法でもチューニングが手軽にできます。
{caret}パッケージの場合
{caret}は膨大な機械学習手法*5をカバーしていますが、基本的にはメインとなるtrain()関数の中のtuneLength, tuneGridという2つの引数でチューニング周りはコントロールします。例えばmethod="svmRadial"*6の場合はこんな感じになります。
> sample_d.c_svm<-train(cv~.,data=sample_d,method="svmRadial",trace=T,tuneLength=10) > print(sample_d.c_svm) 3000 samples 7 predictors 2 classes: 'No', 'Yes' No pre-processing Resampling: Bootstrap (25 reps) Summary of sample sizes: 3000, 3000, 3000, 3000, 3000, 3000, ... Resampling results across tuning parameters: C Accuracy Kappa Accuracy SD Kappa SD 0.25 0.935 0.87 0.00551 0.011 0.5 0.936 0.871 0.00532 0.0107 1 0.935 0.87 0.00561 0.0112 2 0.934 0.867 0.00624 0.0125 4 0.933 0.865 0.00596 0.0119 8 0.932 0.863 0.00587 0.0118 16 0.931 0.861 0.00602 0.0121 32 0.93 0.861 0.00596 0.0119 64 0.93 0.861 0.00596 0.0119 128 0.93 0.861 0.00596 0.0119 Tuning parameter 'sigma' was held constant at a value of 0.1179734 Accuracy was used to select the optimal model using the largest value. The final values used for the model were C = 0.5 and sigma = 0.118.
tuneLength=10としたので、10通りのパラメータに対してチューニングテストが行われています。2の倍数に沿ってマージンパラメータCを動かしてる感じでしょうか。
で、この結果を見れば分かるようにもうちょっとグリッドサーチの範囲を絞ればパフォーマンスが良くなる気がするので、creatGrid()関数でグリッドサーチの範囲を明示的に決めて、train()関数の引数tuneGridとして与えてみることにしました。
> t.grid<-createGrid("svmRadial",data=sample_d,len=4) > print(t.grid) .sigma .C 1 0.01 0.25 2 0.10 0.25 3 1.00 0.25 4 10.00 0.25 5 0.01 0.50 6 0.10 0.50 7 1.00 0.50 8 10.00 0.50 9 0.01 1.00 10 0.10 1.00 11 1.00 1.00 12 10.00 1.00 13 0.01 2.00 14 0.10 2.00 15 1.00 2.00 16 10.00 2.00 > sample_d.c_svm<-train(cv~.,data=sample_d,method="svmRadial",trace=F,tuneGrid=t.grid) > print(sample_d.c_svm) 3000 samples 7 predictors 2 classes: 'No', 'Yes' No pre-processing Resampling: Bootstrap (25 reps) Summary of sample sizes: 3000, 3000, 3000, 3000, 3000, 3000, ... Resampling results across tuning parameters: C sigma Accuracy Kappa Accuracy SD Kappa SD 0.25 0.01 0.935 0.87 0.00546 0.0109 0.25 0.1 0.936 0.872 0.00574 0.0115 0.25 1 0.932 0.864 0.00544 0.0109 0.25 10 0.931 0.862 0.00512 0.0103 0.5 0.01 0.935 0.87 0.00546 0.0109 0.5 0.1 0.936 0.872 0.00545 0.0109 0.5 1 0.932 0.865 0.00431 0.0087 0.5 10 0.932 0.865 0.00431 0.0087 1 0.01 0.935 0.87 0.00546 0.0109 1 0.1 0.937 0.873 0.00516 0.0103 1 1 0.932 0.865 0.00431 0.0087 1 10 0.932 0.865 0.00431 0.0087 2 0.01 0.935 0.87 0.00546 0.0109 2 0.1 0.935 0.871 0.00454 0.00912 2 1 0.932 0.865 0.00431 0.0087 2 10 0.932 0.865 0.00431 0.0087 Accuracy was used to select the optimal model using the largest value. The final values used for the model were C = 1 and sigma = 0.1.
何となくところどころローカルミニマムみたいな組み合わせがあるようですが、わずかながらさらにaccuracyが向上しました。あんまりグリッドサーチ範囲を広げると計算時間がかかるので、その辺は適当に。
他にもmethod="nnet"、すなわちニューラルネットワークの場合だとこんな感じになります。
> sample_d.c_nnet<-train(cv~.,data=sample_d,method="nnet",tuneLength=4,maxit=100,trace=F) > print(sample_d.c_nnet) 3000 samples 7 predictors 2 classes: 'No', 'Yes' No pre-processing Resampling: Bootstrap (25 reps) Summary of sample sizes: 3000, 3000, 3000, 3000, 3000, 3000, ... Resampling results across tuning parameters: size decay Accuracy Kappa Accuracy SD Kappa SD 1 0 0.937 0.874 0.00796 0.0159 1 1e-04 0.937 0.874 0.00717 0.0143 1 0.00316 0.937 0.874 0.0069 0.0138 1 0.1 0.936 0.872 0.00703 0.0141 3 0 0.934 0.869 0.00704 0.0141 3 1e-04 0.935 0.87 0.00648 0.013 3 0.00316 0.933 0.866 0.0144 0.0288 3 0.1 0.936 0.872 0.0066 0.0132 5 0 0.931 0.863 0.00653 0.0131 5 1e-04 0.933 0.866 0.00691 0.0138 5 0.00316 0.932 0.864 0.0067 0.0134 5 0.1 0.935 0.87 0.00702 0.0141 7 0 0.928 0.856 0.0217 0.0429 7 1e-04 0.932 0.863 0.00673 0.0135 7 0.00316 0.932 0.864 0.00673 0.0135 7 0.1 0.934 0.868 0.00722 0.0145 Accuracy was used to select the optimal model using the largest value. The final values used for the model were size = 1 and decay = 0.
tuneLength=4とすることで、4通りのパラメータ（nnetの場合は個々の括りごとに4通り）をグリッドサーチで試した結果が示されています。nnetだと収束に時間がかかるケースがあるので、これまたグリッドサーチ範囲の決め方は色々気を付けた方が良いかもです。
おまけ
特に{e1071}, {caret}は手法ごとにチューニングのやり方がまちまちらしいので、この辺はまたナレッジがたまったら紹介します。。。そのうちどこかのTokyoRで話そうかな。
