Grid Searchによる勾配ブースティングのハイパーパラメータ調整

## はじめに
本日は以前に行った[Jリーグの観客動員数予測](https://signate.jp/competitions/27)のデータコンペの修正の続きに取り組んでいこうと思います。自分が初めて解いたときは[こちら](https://qiita.com/flystaslingan40/private/95115459b0384f03f90e)（修正後は[こちら](https://qiita.com/flystaslingan40/items/7fe67fb47d88e4811301)）
さて、前回は目的変数の正規分布化により、スコアを100程度縮めたわけですが、本日はコンペにて用いたGradient Boosting（勾配ブースティング）のハイパーパラメータ調整をGrid Search(グリッドサーチ)を用いて行うことにします。分類に比べて回帰にアンサンブル学習を用いる時の記事が少ないと感じたので回帰を題材とした記事を書いてみようと思います。

##Grid Searchとは？
さて本日利用するGrid Searchについて説明します。これはハイパーパラメータの値を適切なものに近づかせるために用いる手法です。その名の通りGrid(格子点)をSearch(探索)するものです。つまり複数のハイパーパラメータについて組み合わせを全て試し、最も評価精度の良いものを発見する手法です。[SlideShare](https://www.slideshare.net/hoxo_m/ss-77421091)から取ってきた下図が非常にわかりやすいです。
![-22-638.jpg](https://qiita-image-store.s3.amazonaws.com/0/287300/1a33d24b-0472-7bf0-d629-6723119aed30.jpeg)
scikit learnにはGrid Searchを行う関数GridSearchCVがある。この引数は下図にて表される。

|パラメータ名|説明|
|:---:|:---|
|estimator|チューニングを行うモデル|
|param_grid|パラメタ候補値を「パラメタ名, 候補値リスト」の辞書で与える|
|n_jobs|同時実行数(-1にするとコア数で同時実行)|
|refit|Trueだと最良だったパラメタを使い学習データ全体で再学習する|
|cv|Cross validationの分割数(デフォルト値は3)|
|verbose|ログ出力レベル|

##勾配ブースティングのハイパーパラメータの概要
勾配ブースティングはアンサンブル学習の一つで、ある決定木の間違いを修正するような決定木を逐次的につくることで､汎化性能の向上を目指すアルゴリズムです。最適化のために勾配降下法を用いることがその名の由来です。
さて、勾配ブースティングのハイパーパラメータは以下の通りです。 [scikit learnの公式ドキュメント](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingRegressor.html)を参考に作りました。

|パラメータ名|デフォルト値|説明|
|:---------:|:----:|:----|
|loss|ls|損失関数を最適化するもの<br>{'ls', 'lad', 'huber', 'quantile'}から選択|
|learning_rate|0.1|合計でいくつの決定木を学習させるか。<br>n_estimatorsとトレードオフの関係|
|n_estimators|100|後に学習させる決定木が､前の決定木の間違いをどれだけ強く修正するか|
|subsample|1.0|各ステップの決定木の構築に関係するデータの割合。<br>1.0より小さいと[確率的勾配ブースティング](https://www.rco.recruit.co.jp/career/engineer/blog/32/)になる。|
|criterion|friedman_mse|分割の品質を測定する。基本的にデフォルト値が良い。|
|min_samples_split|2|内部ノードを分割するために必要なサンプルの最小数|
|min_samples_leaf|1|葉として必要なサンプルの最小数。<br>左右の枝分かれの葉がこの数より多いようにする|
|min_weight_fraction_leaf|0|min_samples_leafの割合表示|
|max_depth|3|各木構造の最大深度。<br>勾配ブースティングの特徴を反映させるには小さく設定するべき値。|
|init|None|初期予測を計算するもの|
|random_state|None|乱数のシード値の指定|
|max_features|None|最適な分割を探す際に考慮する機能の数<br>大きな値ほど過学習が起きやすい。<br>autoの時はmax_features=n_features<br>sqrtの時はmax_features=sqrt(n_features)<br>log2の時はmax_features=log2(n_features)<br>Noneの時はmax_features=n_features|
|alpha|0.9|[分位点回帰](http://kosugitti.net/archives/171)、[Huber損失](https://ja.wikipedia.org/wiki/Huber%E6%90%8D%E5%A4%B1)のα|
|verbose|0|1の場合は冗長出力を有効化する。|
|max_leaf_nodes|None|葉の個数の最大値、Noneを指定すると無制限になる。|
|warm_start|False|Trueで前のフィッティングに対して、追加ツリーをフィットさせる。|
|presort|auto|フィッティングを早めるために事前にソートを行うかどうか(boolかauto)|
|validation_fraction|0.1|早期停止のための妥当性確認として設定される訓練データの割合|
|n_iter_no_change|None|検証スコアが改善されない時に早期終了するかの決定|
|tol|0.1|早期終了のための許容差を決定する|

##パラメータ調整
[前回](https://qiita.com/flystaslingan40/items/7fe67fb47d88e4811301)の点数は3,268.64733点で、これは勾配ブースティングのハイパーパラメータをデフォルトの値にしたものであり、このスコアを越すことを目標にしようと思いました。
import文は以下の通りです。

```python
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.grid_search import GridSearchCV
from mpl_toolkits.mplot3d import Axes3D
% matplotlib inline
```

#### 1. n_estimators、learning_rateの調整
n_estimatorsは20~100で調整しlearning_rateは0.05~0.2で調整します。他パラメータはまず以下のように固定します。<br><b>min_samples_split = 5</b>: 総サンプル数の0.5~1%が一般的、今回はデータに偏りがあまりないので5<br><b>min_samples_leaf = 50</b>: 直感で選びます。過学習を防ぐ目的のためあまり小さな値を採用しません<br><b>max_depth = 5</b>: サンプル数に応じて5~8ぐらいで選択。今回はデータ数が少ないので5<br><b>max_features = ‘sqrt’</b>: 最初はsqrtが一般的<br><b>subsample = 0.8</b>: 一般的な値

```python
from sklearn.grid_search import GridSearchCV
from mpl_toolkits.mplot3d import Axes3D
#　ホールドアウト法
Xg_train_J1, Xg_test_J1, yg_train_J1, yg_test_J1 = train_test_split(X_J1, y_J1, test_size=0.3, random_state=0)
mod1J1 =GradientBoostingRegressor(
    min_samples_split = 5, 
    min_samples_leaf = 50, 
    max_depth = 5, 
    max_features = 'sqrt', 
    subsample = 0.8)
paramJ1_1 = {'n_estimators': list(range(20, 101, 10)),
             'learning_rate': list(np.arange(0.05, 0.20, 0.01))}
# scoringを neg_mean_squared_error(負のMSE)とする
gsearch1 = GridSearchCV(estimator = mod1J1, 
                        param_grid = paramJ1_1, 
                        cv = 5, 
                        n_jobs=4, 
                        scoring = 'neg_mean_squared_error')
gsearch1.fit(Xg_train_J1, yg_train_J1)
x = []
y = []
z = []
for i in gsearch1.grid_scores_:
    x.append(i[0]['n_estimators'])
    y.append(i[0]['learning_rate'])
    z.append(i[1])
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(x, y, z)
plt.show()
# test精度の平均が最も高かった組み合わせを出力
print(gsearch1.best_params_)
```
この結果、

```
{'learning_rate': 0.19, 'n_estimators': 100}
```
と出力されました。Axes3Dによる図も出力したがよくわからないグラフになっていたので割愛します。
また、RMSEを確認するため

```python
# そのときのtest精度の平均を出力(RMSE)
print((- gsearch1.best_score_) ** (1/2))
# ホールド・アウト法によるtest精度を出力(RMSE)
print((- gsearch1.score(Xg_test_J1, yg_test_J1)) ** (1/2))
```
と入力すると結果として

```
0.298537282383042
0.27801005403635193
```
と出力されました。

#### 2. min_samples_splitとmax_depthの調整
min_samples_splitとmax_depthをともに3~7の自然数で調整します。上記とほぼ同様のコードを書くことで調整を行いました。すると

```
{'max_depth': 6, 'min_samples_split': 5}
```
と出力され、test精度の平均(RMSE)、ホールド・アウト法によるtest精度(RMSE)はそれぞれ

```
0.2942748176491349
0.2784377684697284
```
でした。　

####3. min_samples_leafの調整
続いて、min_samples_leafを1~71の10飛びの値で調整した。今回は変数が少ないためAxes3Dを用いずmatplotlib.pyplotを利用し、図を表示した。図は下記のようになる。
<img width="280" alt="スクリーンショット 2018-12-26 18.30.01.png" src="https://qiita-image-store.s3.amazonaws.com/0/287300/449fdb9b-a013-659e-7a43-53c1e506d687.png">

```
{'min_samples_leaf': 11}
```
と出力されました。ただ、グラフを見る限り1~11の間の値がmin_samples_leafに最適な値であるように考えられたのでmin_samples_leafを1~20の間で推移させ調整を行いました。図は下図のようになりました。
<img width="303" alt="スクリーンショット 2018-12-26 18.33.23.png" src="https://qiita-image-store.s3.amazonaws.com/0/287300/5556adfd-4d1f-6376-ae94-fb933be02ae9.png">

```
{'min_samples_leaf': 4}
```
またtest精度の平均(RMSE)、ホールド・アウト法によるtest精度(RMSE)はそれぞれ

```
0.2125592561762786
0.1994161961345939
```
でした。　

####4. max_featuresの調整
続いて、max_featuresを7~19の2飛びの値で調整した。
<img width="309" alt="スクリーンショット 2018-12-26 18.36.07.png" src="https://qiita-image-store.s3.amazonaws.com/0/287300/047695cb-baee-fc7d-2b59-62bb6d6aeb4c.png">

```
{'max_features': 13}
```
と出力されました。
またtest精度の平均(RMSE)、ホールド・アウト法によるtest精度(RMSE)はそれぞれ

```
0.20940111752331694
0.19811875953543406
```
でした。　
この手順でJ２リーグのデータも調整することで完了です。お疲れ様でした。
全コードは[こちら](https://github.com/kylasakerskog/SignateJLeague/blob/master/JLeagueGridSearch.ipynb)

##まとめ
このGridSearchによるパラメータ調整の結果、Scoreが
<b>3,268.64733点</b>→<b>3,245.83137点</b>
と20点程度、スコアを伸ばすことができました。ちなみに順位は
<b>62位/614位</b>
でした。結構時間がかかった割にスコアが伸びにくいなと思いました。