XGBoostでグリッドサーチとクロスバリデーション1

気付きがあったので書いておきます。

多分みんな知ってるんだと思う。

Pythonでsklearn.datasetsにあるload.boston()の回帰をXGBoostを用いて行います。
XGBoostでは、DMatrixという目的変数と目標値が格納されたデータを作成します。testデータ2割、まずはdefaultで。

import xgboost as xgb
import matplotlib.pyplot as plt
from sklearn.datasets import load_boston
from sklearn import metrics as met
import sklearn as skl
import numpy as np

dataset = load_boston()
x = dataset.data
t = dataset.target

x_train, x_test, t_train, t_test = skl.model_selection.train_test_split(x, t, test_size=0.2, random_state=0)

dtrain = xgb.DMatrix(x_train, label=t_train, feature_names=dataset.feature_names)
dtest = xgb.DMatrix(x_test, label=t_test, feature_names=dataset.feature_names)
# 変数名マッチしてないと怒られる

param = {}
param['eta'] = 0.3
param['gamma'] = 0
param['max_depth'] = 6
# max_depthだけ分岐を許す

param['min_child_weight'] = 1
param['subsample'] = 1
param['colsample_bytree'] = 1

param['objective'] = 'reg:squarederror'
param['eval_metric'] = 'rmse'
# default

num_round = 10000
evals = [(dtest, 'eval'),(dtrain, 'train')]
# early_stoppingのためのevaluationはtrainデータでやるみたい。evalデータは値を見せるだけ。
evals_result = {}

bst = xgb.train(param, dtrain, num_boost_round=num_round, evals=evals, evals_result=evals_result, early_stopping_rounds=100)
train_pred = bst.predict(dtrain, ntree_limit=bst.best_ntree_limit)
test_pred = bst.predict(dtest, ntree_limit=bst.best_ntree_limit)

print('train_RMSE_score_is_{:.8f}, test_RMSE_score_is_{:.8f}'.format(np.sqrt(met.mean_squared_error(t_train, train_pred)), np.sqrt(met.mean_squared_error(t_test, test_pred))))
# 0.00116996, 4.62175620,  当然だがeval_metricが出してくる値と同じ
print('train_R2_score_is_{:.4f}, test_R2_score_is_{:.4f}'.format(met.r2_score(t_train, train_pred), met.r2_score(t_test, test_pred)))
# 1.0000, 0.7377 かなりoverfitting?
print('Best Score:{:.4f}, Iteration:{:d}, Ntree_Limit:{:d}'.format(bst.best_score, bst.best_iteration, bst.best_ntree_limit))
# Best Score:0.0012, Iteration:165, Ntree_Limit:166


plt.rcParams['figure.facecolor'] = 'white'

fig = plt.figure(figsize=(13, 5))
ax1 = plt.subplot(121)
ax1 = plt.plot(evals_result['eval']['rmse'], label = 'eval')
ax1 = plt.plot(evals_result['train']['rmse'], label = 'train')
plt.legend()
ax2 = plt.subplot(122)
ax2 = xgb.plot_importance(bst, ax=ax2)

eval_resultはeval_matricのTree数発展で並んだlist. dictになってるからkeyはドキュメント参照。
Feature importanceは目的変数が何回分岐に使われているかというF scoreを返す。

xgb.plot_tree(bst, num_trees=X)

でX番目のTreeの構造見れる。今回はmax_depth=6を入れてるから分岐が最大6つ発生している。
ちなみにnum_trees=0において、max_depth=6の木の構造は3回目の分岐までmax_depthが3の木と同じ。再現性あり。(それ以降はひとつ前の木に依存するので異なる。)

Tree0の構造は同じ

で、上の結果はtrain setのrmseが0.0012に対してtest setが4.6とかなりのoverfittingになっているのでハイパーパラメーターを検討する。

ハイパーパラメーターの調整

ハイパーパラメーターを調整する際にグリッドサーチとクロスバリデーション(CV)という手法がよく行われています。

ここ気付きポイントなのですが、

クロスバリデーションはパラメーターを決めるための指標ではなく、指定したパラメーターがどれだけ汎化性能があるかを示すための数字なので、グリッドサーチと組み合わせないと全く意味がない。

ということ。

xgb.cvを使ってCVだけをやってみます。

params = {}
num_round = 10000

bst_cv = xgb.cv(param, dtrain, num_boost_round=num_round, early_stopping_rounds=100, nfold=5, metrics='rmse')


print(bst_cv)
# bst_cvに対するメゾットがないので使いづらいがとりあえず行数=最適なTree数だと思う。

bst = xgb.train(param, dtrain, num_boost_round=num_round)
train_pred = bst.predict(dtrain)
test_pred = bst.predict(dtest)

print('train_RMSE_score_is_{:.4f}, test_RMSE_score_is_{:.4f}'.format(np.sqrt(met.mean_squared_error(t_train, train_pred)), np.sqrt(met.mean_squared_error(t_test, test_pred))))
print('train_R2_score_is_{:.4f}, test_R2_score_is_{:.4f}'.format(met.r2_score(t_train, train_pred), met.r2_score(t_test, test_pred)))

CVの結果、nfold回の学習が行われ、その過程における各Tree数でのmetricsの評価の平均値がbst_cvに表れています。xgb.trainと同じことをやっていますが、modelは作られないのでxgb.trainを実行する必要があります。

CVにおける各学習で、trainデータとtestデータの分け方によってmetric最小を取るTree数が違いますが、metricの平均値は82本のTreeで最小になっていますね。

パラメーターはすべてdefaultで指定しており、CVによってパラメーターが変わるということはありません。なのでこのままtrainを回してもCVをしない時と全く同じ結果になります。

nfold-CV時のTree数を使って学習をしたいときは、それらしいメゾットがない気がするのでlen(bst_cv)でDataFrameの行数を代入してやればいいと思います。

ちなみにその結果は以下のようになり、Tree数を約半分に減らしてtrainデータに対するscoreが下がっているにも関わらず、testデータに対するscoreはほとんど変わっていません。

# train_RMSE_score_is_0.0324, test_RMSE_score_is_4.6219
# train_R2_score_is_1.0000, test_R2_score_is_0.7377

Scikit-learn AIPによるXGBoostとクロスバリデーション

Scikit-learnっぽくXGBoostやるよ、みたいな関数があります。web検索してもこっちのコードばっかりヒットするが、GridSearchCVが使えるからみんなこっち使ってるんだと思う。

import xgboost as xgb
import matplotlib.pyplot as plt
from sklearn.datasets import load_boston
from sklearn import metrics as met
import numpy as np
from sklearn.model_selection import train_test_split, GridSearchCV

dataset = load_boston()
x = dataset.data
t = dataset.target

x_train, x_test, t_train, t_test = train_test_split(x, t, test_size=0.2, random_state=0)

# sklearn AIPなのでDMatrix作る必要ない


mod = xgb.XGBRegressor(max_depth = 6, learning_rate = 0.3, n_estimators=10000, objective='reg:squarederror', gamma=0, min_child_weight=1, subsample=1, colsample_bytree=1)
# モデルを作成。xgb.trainにおけるparamの部分

params = {'max_depth': [6], 'learning_rate': [0.3]}
gscv = GridSearchCV(mod, params, cv = 5, refit=True,  scoring='r2', verbose = 1, n_jobs=1)
# GridSearchCVをこれからやるよっていう宣言だけ。グリッドサーチは.fitで実行される。
# 今回はparamsは固定

evallist = [(x_train, t_train)]
gscv.fit(x_train, t_train, eval_metric='rmse', eval_set=evallist, early_stopping_rounds=100)
# 学習を行う。evallistの値に対してRMSEで評価を行い、100round後も変化がなければ終了。

print(gscv.best_estimator_)

train_pred = gscv.predict(x_train)
test_pred = gscv.predict(x_test)

print('train_RMSE_score_is_{:.4f}, test_RMSE_score_is_{:.4f}'.format(np.sqrt(met.mean_squared_error(t_train, train_pred)), np.sqrt(met.mean_squared_error(t_test, test_pred))))
# 0.0012, 4.6218
print('train_R2_score_is_{:.4f}, test_R2_score_is_{:.4f}'.format(met.r2_score(t_train, train_pred), met.r2_score(t_test, test_pred)))
# 1.0000, 0.7377
print('Best Score:{:.4f}'.format(gscv.best_score_))
# Best Score:0.8928 (R2)

XGBRegressorでモデル作ったりGridSrearchCVでグリッドサーチとCVのパラメーター指定したりするが実際に実行するのは.fitの部分。

つまりXGBRegressorはXGBにおけるparam = {}の部分、
GridSearchCV＋.fitはセットで、xgb.cvからxgb.trainまでを行っていることになります。

GridSearchCVのrefitをTrueにすると, CVの結果による最適なパラメーターを用いて、データ全体に対するmodelの作成を行います。
従ってrefitはxgb.trainの部分を行うかどうかを指定します。

CVだけやっても意味がないの本質はここにあって、
CVによって汎化性能が高い最適なパラメーターを探索することがCVの意義であるので、パラメーターを振らない状態でやっても出てくるデータは同じということです。

結果はxgb.cvの時と全く一緒。

まとめ

CVによって汎化性能が高い最適なパラメーターを探索することがCVの意義であるので、パラメーターを振らない状態でやっても出てくるデータは同じ

とにかくグリッドサーチとクロスバリデーションは一緒にやれ

文章にすると当たり障りのない感じになってしまう。。

そうなるとxgb.cvで求めた汎化性能が高いTree数とxgb.trainで求めたtrainデータに対してscoreがいいTree数のどちらを採用すればよいのかとか、またいろいろ考えることが出てきますが、長くなったので一旦切ります。

次はグリッドサーチとクロスバリデーションを一緒にやった結果を書く予定です。