Rで機械学習するならチューニングもグリッドサーチ関数orオプションでお手軽に

ちょっと調べてみたらタイトルの件について言及してる記事があまり多くなかったので、ざっくり書いてみます。なお、この記事はid:shakezoさんの

• 実務でRandomForestを使ったときに聞かれたこと

へのオマージュです。というか、実は僕もこの記事を読んでから「多分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 のチューニングのしかた（２）」を参考にさせていただきました。というか、そちらの記事の方が

チューニングの手順としては、

1. グリッドサーチで大雑把に検索する。
2. 最適なパラメータがありそうなところを絞って再びグリッドサーチを行う。

という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で話そうかな。

*1:実はMatlabスクリプトをGitHubの某所に転がしてあります。。。

*2:これが10万レコードで50次元以上の実データとかだと見たくもない数字になることも。。。

*3:言わずと知れた多言語対応SVMパッケージです

*4:でも肝心の決定木を実行する関数が見当たらない。。。

*5:そもそも他のパッケージを山ほど依存関係で引っ張ってきている

*6:{kernlab}のksvm()関数のラッパーらしい