【XGB】交差検証法を使った勾配ブースティング決定木の実装|python

今回は、kaggleなどのデータ分析コンペでもよく使われる「勾配ブースティング決定木アルゴリズム」の解説を行います。

このコンテンツでわかること

・実際中身でどのような計算をしているのか

・コード例

勾配ブースティング決定木(gradient boosting decision tree)

勾配ブースティング決定木(以下GBDT)とは、決定木の集合です。

であれば、ランダムフォレスト(random forest)との違いが気になります。

ランダムフォレストは、並列に決定木を作成する(アンサンブルのバギングをベースに少しずつ異なる決定木を集める)のに対して、GBDTは「逐次的に決定木を追加していく」ことに特徴があります。

・まず、一本目の決定木が作成され、予測値が出ます。

・予測値と教師データの目的変数との誤差から、「どこの特徴量のどの重みをどのように変えるか」を計算し、誤差を小さくするように次の木を作ります。

・これを繰り返します。(ハイパーパラメータで定めた分だけ繰り返します)

・予測値は、予測対象のデータがそれぞれの決定木で属するノードの重みの和で計算されます。

特徴としては、大きく2点あります。

1つ目：変数間の相互作用が反映される

決定木アルゴリズムの一種なので、分岐を繰り返すことによる説明変数の相互作用が反映されます。

2つ目：正規化が不要な場合が多い

決定木アルゴリズムは、大きいか小さいかだけを問題にしており、変数間のスケールの違いが予測値に致命的な影響を与えることは少ないです。

＊ただし弱学習器の数を増やし過ぎてしまうと、過学習が起きてしまいます。

アルゴリズムについて(xgboostを例に)

GBDTのライブラリは、xgboostやLightGBM,catboostなどありますが、今回はxgboostのアルゴリズムについて解説いたします。

基本的な決定木の作り方は、【機械学習】決定木の仕組みと実装方法について|pythonに記載してあります。

w：重み

m：決定木の添字。1~Mの値をとります。

i：レコードの添字。1~Nの値をとります。

$$l(y_{i},\hat{y})$$

上は、目的変数\(y_i\)と予測値を与えられた時の目的関数です。

そして、決定木をfmとし、決定木に対する正則化項を\(Ω(fm)\)とした時、目的関数の和は以下のように表すことができます。

$$L = \sum_{i=1}^{N}l(y_{I},\hat{y})+\sum_{m=1}^{M}Ω(f_{m})$$

正則化項とは、関数に対して罰則を与えるもので、パラメーターが大きくなるほど大きな値を返すことになります。

目的関数を大きく減少させる重みが採用されるので、罰則が大きいと最適ではありません。

GBDTの正則化項は以下のように表すことができます。

$$Ω(f_{m})=\sqrt{T}+\sqrt{2}λ\sum_{j}w{j}^2+α\sum_{j}|w{i}|$$

見た目ほど難しいものではありません。

特に第2項と第3項ですが、絶対値と平方が使われており、「重みが大きくなる」と正則化項(罰則)が大きくなることがわかります。

Tはノードの数です。

とにかく大切なのは、どの特徴量のどの値で分岐させるかは、目的関数の減少が最も大きい時のノードへの重みの付け方で決定するということです。

なぜそのような重みの決定をするかといえば、「過学習を防ぐため」です。

勾配とは

では、目的関数をどのように最小にするのでしょうか。

一旦簡略化のため先ほどの正則化項を除いた目的関数を考えます。

重み\(w\)も考えると、レコード\(i\)の集合の目的関数の和は以下のようになります。

$$L_{j}=\sum_{I\in{I_{i}}}l(y_{I},\hat{y}_{I}+w_{j})$$

これを最小化したいわけですが、少しテクニカルな話になります。

勾配と二階微分値を使って、2次関数に近似させます。

$$勾配g_{i}=\frac{\partial l}{\partial \hat{y}_{i}},二回微分値h_{i}=\frac{\partial^2l}{∂\hat{y}^2}$$

$$\tilde{L}_{j}=\sum_{i \in {I_{i}}}l(y_{i},\hat{y}_{i})+g_{i}w_{i}+\frac{1}{2}h_{i}w_{i}$$

このようなwについての二次関数に近似できました。

また、第1項はwを変数にもたないので取り除けます。定数ですね。

ここからは高校数学でもわかりますね。

「微分してイコールゼロ」というやつです。

$$w_{j}^{*}=- \frac{\sum_{i \in{I}}g_{i}}{\sum_{i \in{I}}h_{i}},\tilde{L}_{j}=-\frac{1}{2} \frac{(\sum_{i\in{I}}g_{i})^2}{\sum_{i\in{I}}h_{i}}$$

CODE|乳がんデータ

scilitlearnの乳がんデータを使って、targetを予測してみましょう。

乳房塊の微細針吸引物（FNA）のデジタル化画像から計算されており、画像中に存在する細胞核の特徴を捉えたものです。

データセットの中では悪性(malignant)は0、良性(benign)は1で表されており、targetカラムで表されております。

from sklearn.model_selection import train_test_split
from sklearn.datasets import load_breast_cancer
import numpy as np
import pandas as pd
data_breast_cancer = load_breast_cancer()

# Pandasによるデータの表示
df_target = pd.DataFrame(data_breast_cancer["target"], columns=["target"])
df_data = pd.DataFrame(data_breast_cancer["data"], columns=data_breast_cancer["feature_names"])
df = pd.concat([df_target, df_data], axis=1)
df.head()

まずはデータのロードとデータの可視化です。

df.head()で先頭5行が表示されたと思います。

説明変数をXに格納し、目的変数をyに格納していきます。

import xgboost as xgb
from sklearn.model_selection import KFold
from sklearn.model_selection import cross_val_score

y = df.loc[:,["target"]]
X = df.loc[:, "mean radius":]

xgb_reg_grid = xgb.XGBRegressor(random_state = 0)

# 訓練データとテストデータに分ける
k_fold = KFold(n_splits = 5, shuffle = True, random_state = 0)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5, random_state=0)
xgb_scores = cross_val_score(xgb_reg_grid,X,y,cv=k_fold,scoring="r2")

xgb_reg_grid にXGBのインスタンスを格納します。 xgb.XGBRegressor(random_state = 0) のrandom_state=0は再現性を確保するために必要な引数です。

train_test_splitで教師データとテストデータに分割します。

from sklearn.model_selection import GridSearchCV
params = {"booster":["gbtree"],
          "n_estimators":[10,20,30,40,50,100],
          "max_depth":[2,3,4,5,6,7,8],
          "learning_rate":[0.1,0.25,0.5,0.75,1.0],
          "colsample_bytree":[0.1,0.25,0.5,0.75,1.0],
          "random_state":[0]
         }
grid = GridSearchCV(estimator=xgb_reg_grid, param_grid = params,cv = k_fold,scoring="r2" )

grid.fit(X_train,y_train)
print(grid.best_params_)
print(grid.best_score_)

今回は、パラメータのチューニングにグリッドサーチ(grid search)を使用します。

与えられたパラメーターを総当たりで確認し、最も評価が高いパラメータをbest_params_として保持します。

ベイズ最適化やランダムサーチに対して、計算量が多く時間がかかるが、選択肢の中から抜け漏れなく最適なパラメータを出せるという特徴があります。

初手からグリッドサーチをするのはあまりお勧めできません。データ量が多いと探索に時間がかかってしまいます。

ランダムサーチなどで当たりをつけてから、的を絞ったグリッドサーチを行うのがセオリーです。

また、ガウス過程回帰を使って探索回数のコスパをよくするベイズ最適化という手法もあります。

詳しくは、【kaggle】ベイズ最適化とXGBでtitanicの予測問題を解く|pythonをご覧ください。

for文で簡単に実装できますが、便利なメソッドがあるので使います。

まずparamsに辞書型でパラメータの名前と値の候補を格納します。

GridSearchCVというメソッドの引数である、param_gridに代入しています。

→{'booster': 'gbtree', 'colsample_bytree': 0.1, 'learning_rate': 0.25, 'max_depth': 2, 'n_estimators': 40, 'random_state': 0}
0.8323395038069155

他の機械学習のライブラリと同じように、インスタンスの生成→fit→predictで予測値が出せます。

print(grid.best_params_)でパラメーターチューニングの結果、最適だった組み合わせが出ました。

y_test_pred = grid.predict(X_test)
y_test_pred = np.expand_dims(y_test_pred,1)

#確率が格納されているので、0.5を閾値に確率に直す
pred = np.where(y_test_pred > 0.5, 1, 0)
pred[:5]

最後にテストデータを入れて、予測値をみてみましょう。

np.expand_dims() は、第2引数の axis で指定した場所の直前に dim=1 を挿入しています。

新しい軸が追加されて、配列が拡大されています。

→array([[0],
       [1],
       [1],
       [1],
       [1],

出力してみればわかりますが、y_test_predには、1(benign)になる確率が格納されています。

0.5を閾値に0と1に振り分けて、きちんとした予測値を出力させましょう。

症例の少ない病気にかかるかどうかなど、2値の分布が均等でない場合は閾値を適当に0.5と決めてしまうのは早計です。

しっかりedaなどを行い、分布を確認しておきましょう。

ROU曲線やAUCなどで閾値の分析ができます。いづれ扱います。

【多変量解析】ROC曲線とAUCによる判別分析|pythonをご覧ください。

他にも機械学習コンテンツを用意しております。

【分類タスク】ロジスティック回帰の使い方|python

【機械学習】決定木の仕組みと実装方法について|python

【Sequential】Kerasを使ったニューラルネットワーク|python

【判別問題】サポートベクトルマシン(SVM)の仕組み|python

FOLLOW ME ！

@blue_statistics