時系列性を持つデータに対するLightGBM予測モデルをhyperoptでパラメータチューニングした：結果の再現とn_estimatorsを整数型にするのにハマった件

はじめに

python機械学習の6章には、GridSearchによるハイパーパラメータチューニングの方法は載っていました。
やはり、ベイズ最適化によるハイパーパラメータチューニングもおさえておきたいです。Kaggleなどでよく使われているhyperoptを実際に使ってみました。この際、以下のことにハマってしまいました。

• 乱数シードを固定し、結果を再現する
• RandomForest、LightGBMのn_estimatorsの様に整数型として扱う

自分の備忘録を兼ねて、hyperoptの使い方を上記の点に重きを置きながら記事を書きたいと思います。

この記事を書くにあたっては、以下の記事を参考にさせてもらいました。筆者の皆さんにここでお礼を申し上げます。
hyperoptって何してんの？
hyperoptで再現性を担保するために必要なパラメータ設定について
Hyperparameters tunning with Hyperopt

hyperoptについて

ベイズ最適化は、シーケンシャルモデルベース最適化（SMBO）とも呼ばれています。SMBOのうち、Tree-structured Parzen Estimator Approach（TPE）というロジックを実装して最適なパラメータを探索してくれるライブラリにhyperoptがあります。
hyperoptドキュメント
TPEのアルゴリズムなどについて、図など可視化した説明と共に解説してあるものはHyperparameters tunning with Hyperoptを参照ください（英語です。悪しからず）。

私はanacondaを使っているので、以下のコマンドでhyperoptライブラリをインストールします。

$ conda install -c conda-forge hyperopt

hyperoptを使ったモデルの比較対象：GridSearchCVを使ったモデルも作ります。

hyperoptによるパラメータチューニングの結果を評価するため、比較対象としてGrid Searchでパラメータを探索したモデルも作成します。Grid Seachではまず広い範囲を粗く探索してあたりをつけて、その後であたりをつけた近傍領域を細かく探索することにします。

この記事で扱うデータ

この記事では以下の2つのデータを扱います。

• Python機械学習6章で良性腫瘍と悪性腫瘍に関するthe Breast Cancer Wisconsin datase　：　こちらのデータのダウンロードなどについては拙著良性腫瘍と悪性腫瘍に関するthe Breast Cancer Wisconsin dataseの準備を参照ください。
• KaggleのRecruit Restraunt Visotors Prediction

the Breast Cancer Wisconsin dataseを使ってサポートベクターによる分類モデルを作成

the Breast Cancer Wisconsin dataseを使い、サポートベクターによる分類モデル（SVC）を作ることにします。
基底関数としては動径基底関数（RBF）を使うことにします。正則化パラメータCと動径基底関数のパラメータγの２つを調整するパラメータとします。

データの準備

#データの読み込み
df = pd.read_csv('./rawdata/wdbc.data', header=None)

X = df.loc[:,2:].values
y = df.loc[:,1]

# ラベルを0,1で置き換える
from sklearn.preprocessing import LabelEncoder

le = LabelEncoder()
y=le.fit_transform(y)

# 訓練データとテストデータに分割
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.20, stratify=y, random_state=1)

Grid SearchによるSVCモデル

はじめに、パラメータを0.0001～1000の範囲で粗く探索します。

param_range = [0.0001, 0.001, 0.01, 0.1, 1, 10, 100, 1000]

param_grid = [{'svc__C':param_range, 'svc__kernel':['rbf'], 'svc__gamma':param_range}]

gs = GridSearchCV(estimator=pipe_svc,
                 param_grid=param_grid,
                 scoring='accuracy',
                 cv=st_kfold)

gs.fit(X_train, y_train)
print(gs.best_score_)
print(gs.best_params_)
print(gs.score(X_test, y_test))

0.9781159420289856
{'svc_C': 10, 'svcgamma': 0.01, 'svc_kernel': 'rbf'}
0.9736842105263158

次に、C=10、γ=0.01の近傍で探索します。

param_C = [5,10,15,20,30,40,50,60,70,80]
param_gamma = [0.006, 0.008, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06]

param_grid2 = [{'svc__C':param_C, 'svc__kernel':['rbf'], 'svc__gamma':param_gamma}]

gs2 = GridSearchCV(estimator=pipe_svc,
                 param_grid=param_grid2,
                 scoring='accuracy',
                 cv=st_kfold)

gs2.fit(X_train, y_train)
print(gs2.best_score_)
print(gs2.best_params_)
print(gs2.score(X_test, y_test))

0.9846376811594203
{'svc_C': 30, 'svcgamma': 0.008, 'svc_kernel': 'rbf'}
0.9736842105263158

テストデータに対する精度は0.974でした。上の粗く探索した結果から、訓練データに対しては精度が上がりましたが、テストデータに対しては精度は変わりませんでした。。。

hyperoptによるSVCモデル

つぎに、hyperoptを使って、正則化パラメータCと動径基底関数のパラメータγを最適化します。
hyperoptをインポートします。また、最適化の対象とする目的関数として、クロスバリデーションによるスコアの平均値を使うことにします。

from hyperopt import hp, tpe, Trials, fmin
from sklearn.model_selection import cross_val_score

探索範囲の定義

探索範囲を指定します。0.0001～1000のように、桁を大きく動かすものは、logをとった値が一様とあるように探索したいですね。この場合は、hp.loguniformを使います。
詳しくはhyperoptのドキュメントDefining a Search Spaceに説明があります。

hyperopt_parameters = {
    'C': hp.loguniform('C', np.log(0.0001), np.log(1000)),
    'gamma': hp.loguniform('gamma', np.log(0.0001), np.log(1000)),
    'kernel': hp.choice('kernel', ['rbf'])
}

目的関数の定義

目的関数を定義します。hyperoptでは、目的関数の値を最小化するようにパラメータを探索します。精度の様に大きいほどよい値を目的とする場合は、-1をかけて目的関数を定義します。

def objective(args):

    svcH = SVC(**args, random_state=1)
    st_kfold = StratifiedKFold(n_splits=10, random_state=1, shuffle=True)

    pipe_svcH = make_pipeline(StandardScaler(), svcH)

    score_array = cross_val_score(estimator=pipe_svcH, X=X_train, y=y_train, cv=st_kfold)

    score_mean = np.mean(score_array)

    return -1*score_mean

実行

繰り返す回数を指定して、fminを呼び出して最適なパラメータ探索を実行します。
結果に再現性を持たせるには、rstate=np.random.RandomState(0)を指定して乱数シードを指定します。　これをやらないと、結果が毎回かわってしまいます。。。

# iterationする回数
max_evals = 500
# 試行の過程を記録するインスタンス
trials = Trials()

best = fmin(
    # 最小化する値を定義した関数
    objective,
    # 探索するパラメータのdictもしくはlist
    space=hyperopt_parameters,
    # どのロジックを利用するか、基本的にはtpe.suggestでok
    algo=tpe.suggest,
    max_evals=max_evals,
    trials=trials,
    # 試行の過程を出力
    verbose=1,
    # 結果再現のため乱数を固定
    rstate=np.random.RandomState(0)
)

print(best)
print(trials.best_trial['result'])

{'C': 26.271132480263706, 'gamma': 0.009415740834817117, 'kernel': 0}
{'loss': -0.9846376811594203, 'status': 'ok'}

GridSearchCVを使って、2回目の細かく探索した結果に近い値が選ばれています。500回繰り返しているので、計算に若干（私のPCで約43秒）かかっていますが、一度で同じ精度の結果が得られました。

探索されたパラメータを使って、訓練データ全体を使ってパラメータを推定し、テストデータに対する精度を求めます。

svc_best = SVC(C=best['C'], gamma=best['gamma'], kernel='rbf', random_state=1)
pipe_best = make_pipeline(StandardScaler(), svc_best)
pipe_best.fit(X_train, y_train)
print(pipe_best.score(X_test, y_test))

0.9736842105263158
精度はGridSearchCVの結果と全く同じでした。この検証では、精度上では優位性を確認できませんでした。ですが、広めの探索範囲をとっておいて、一度の探索でこの結果が得られているので、効率的に探せていると考えられます。

Kaggle：Recruit Restaurant Visitor Forecastingについて LightGBMで”簡単な”予測モデルを作る

次に、KaggleのRecruit Restaurant Visitor Forecastingデータを使って、レストランを訪れる客数を予測するモデルを作ります。
このコンペ自体に関しては、例えばkaggle リクルート レストラン客数予測チャレンジを参照ください。

ここでは、予測するモデル手法としてはLightGBMを使います。過去データから未来を予測する時系列性があるデータです。クロスバリデーションを行うときは、Time Series Cross Validationを使うことにします。

データの準備

ここではリクルートのAirレジユーザの来客数実績から、曜日、休みの前日かなどの情報だけを使って来客数を予測する”簡単な”モデルを作ります。よって

• air_visit_data.csv : Airレジを使っているユーザID、来客日、来客数のデータ
• date_info.csv　：　日付、曜日、休日フラグのデータ

df_air_visit_data = pd.read_csv('./raw_data/air_visit_data.csv/air_visit_data.csv')
df_air_visit_data['visit_date'] = df_air_visit_data['visit_date'].apply(lambda x: datetime.datetime.strptime(x, '%Y-%m-%d'))

df_date_info = pd.read_csv('./raw_data/date_info.csv/date_info.csv')
df_date_info['calendar_date'] = df_date_info['calendar_date'].apply(lambda x: datetime.datetime.strptime(x, '%Y-%m-%d'))

#ユーザー数の増加に影響されずに人数を予測するため、平均値を求める
df_air_visit_mean = df_air_visit_data.groupby('visit_date', as_index=False).mean()

df_air_visit_mean = pd.merge(df_air_visit_mean, df_date_info.rename(columns={'calendar_date':'visit_date'}),\
                                  how='inner', on='visit_date')

#前日が休日
df_air_visit_mean['prev_holiday_flg'] = df_air_visit_mean['holiday_flg'].shift(-1).fillna(0)
#曜日をダミー変数化
df_air_visit_mean_dummy = pd.get_dummies(df_air_visit_mean, columns=['day_of_week'])

X = df_air_visit_mean_dummy.drop(['visit_date', 'visitors'], axis=1)
y = df_air_visit_mean_dummy['visitors']

num_train = int(X.shape[0]*0.8)
print('number of train', num_train)

X_train = X[:num_train]
X_test = X[num_train:]
print('shape of X_train', X_train.shape)
print('shape of X_test ', X_test.shape)

y_train = y[:num_train]
y_test = y[num_train:]

number of train 382
shape of X_train (382, 9)
shape of X_test (96, 9)

このデータを使って来客数を予測するモデルとして、LightGBMを用います。予測ですのでLGBMRegressorを使います。
lightgbm.LGBMRegressor 公式ドキュメント
また、LightGBMはハイパーパラメータの種類が多いですが、チューニングの仕方に関してはLightGBM公式のParameters Tuningを参考にしています。

from sklearn.model_selection import TimeSeriesSplit
from sklearn.model_selection import GridSearchCV

from sklearn.metrics import mean_squared_error 
from sklearn.model_selection import cross_val_score

#LightGBMのインポート
import lightgbm as lgb

Grid SearchによるLightGBMモデル

ここでも、GridSearchCVを使って、LightGBMのハイパーパラメータをチューニングします。
時系列性を考慮したクロスバリデーションを定義しておき、まずは粗くハイパーパラメータを探索します。
予測ですので、評価しひょとしてはroot_mean_squared_errorを使います。

#時系列を考慮したクロスバリデーション
tscv = TimeSeriesSplit()
tscv.split(X_train, y_train)

param_grid1 = {
    'num_leaves':[5,10,20,30,40],
    'max_depth':[5,10,15,20],
    'n_estimators':[10,25,50,100,200],
    'reg_lambda':[0, 0.001, 0.01, 0.1, 1]

}

lgb_estimator = lgb.LGBMRegressor(random_state=0)
gsearch1 = GridSearchCV(estimator=lgb_estimator, param_grid=param_grid1, cv=tscv, scoring = 'neg_root_mean_squared_error')
lgb_model1 = gsearch1.fit(X=X_train, y=y_train)

lgb_model1.best_params_

{'max_depth': 5, 'n_estimators': 25, 'num_leaves': 10, 'reg_lambda': 0.1}

次に、このハイパーパラメータの値の近傍で細かく探索します。

param_grid2 = {
    'num_leaves':[7,8,9,10,11,12,13,14],
    'max_depth':[7,8,9,10,11,12,13],
    'n_estimators':[30,40,50,60,70],
    'reg_lambda':[0,0.005,0.0075,0.01,0.02,0.003,0.005]

}
gsearch2 = GridSearchCV(estimator=lgb_estimator, param_grid=param_grid2, cv=tscv, scoring = 'neg_root_mean_squared_error')
lgb_model2 = gsearch2.fit(X=X_train, y=y_train)
lgb_model2.best_params_

{'max_depth': 7, 'n_estimators': 30, 'num_leaves': 8, 'reg_lambda': 0.01}

テストデータに対して適用します。

np.sqrt( mean_squared_error(y_test,lgb_model2.predict(X_test) ) )

1.5863437047685356

hyperoptによるLightGBMモデル

ここでもまず探索範囲を定義します。

params_hyperopt1 = {
    'num_leaves':hp.quniform('num_leaves',5,40,2),
    'max_depth':hp.quniform('max_depth',5,20,2),
    'n_estimators':hp.quniform('n_estimators',10,200,10),
    'reg_lambda':hp.loguniform('reg_lambda', np.log(0.001), np.log(1))
}

num_leaves（葉の数）など整数値だけを離散的に取るものについては、hp.quniformを使います。reg_lambdaについては0.001～1の間で探索します。

目的関数を以下の様に定義します。ここで、n_estimatorsのように整数型で渡さないといけないものについては、args = {'n_estimators':int(args['n_estimators']),…}のようにします。

def objective(args):
    # ポイント！ int(args['ラベル'])を使って、整数型にする。
    args = {
        'num_leaves':int(args['num_leaves']),
        'max_depth':int(args['max_depth']),
        'n_estimators':int(args['n_estimators']),
        'reg_lambda':args['reg_lambda']

    }

    lgb_hyp = lgb.LGBMRegressor(**args, random_state=0)

    # cvとしてtscvを指定する
    score_array = cross_val_score(estimator=lgb_hyp, X=X_train, y=y_train, cv=tscv, scoring = 'neg_root_mean_squared_error')
    score_mean = np.mean(score_array)

    return -1*score_mean

実行します。

# iterationする回数
max_evals = 2000
# 試行の過程を記録するインスタンス
trials1 = Trials()

best1 = fmin(
    # 最小化する値を定義した関数
    objective,
    # 探索するパラメータのdictもしくはlist
    space=params_hyperopt1,
    # どのロジックを利用するか、基本的にはtpe.suggestでok
    algo=tpe.suggest,
    max_evals=max_evals,
    trials=trials1,
    # 試行の過程を出力
    verbose=1,
    # 結果再現のため乱数を固定
    rstate=np.random.RandomState(0)
)

best1 = {
    'num_leaves':int(best1['num_leaves']),
    'max_depth':int(best1['max_depth']),
    'reg_lambda':best1['reg_lambda'],
    'n_estimators':int(best1['n_estimators'])

}
best1

{'num_leaves': 36,
'max_depth': 14,
'reg_lambda': 0.9067818858763805,
'n_estimators': 30}

テストデータに対して予測し、評価します。

svc_hyp_best1 = lgb.LGBMRegressor(**best1, random_state=0)
svc_hyp_best1
svc_hyp_best1.fit(X_train, y_train)
print( np.sqrt(mean_squared_error(y_test, svc_hyp_best1.predict(X_test))) )

1.5923797093822105

テストデータに対しては、GridSearchCVで作ったモデルより、誤差が大きくなっていますね。Cross Validationで訓練データ内では検証していますが、hyperoptを使うと過学習がやはりおこりやすいのでしょうか。

n_estimatorsが整数型であることを考慮せず目的関数を定義した場合

上の様に目的関数objectiveを定義すればエラーが起こらずに最適化を実行できます。しかし、私は、はじめ次のように目的関数を定義していました。

# エラーが起きた時の目的関数objective
def objective(args):    
    lgb_hyp = lgb.LGBMRegressor(**args, random_state=0)

    score_array = cross_val_score(estimator=lgb_hyp, X=X_train, y=y_train, cv=tscv, scoring = 'neg_root_mean_squared_error')
    score_mean = np.mean(score_array)

    return -1*score_mean

これで
実行すると以下のようなエラーメッセージが出ました。

lightgbm.basic.LightGBMError: Parameter num_iterations should be of type int, got "95.0"

  FitFailedWarning)

---------------------------------------------------------------------------
AllTrialsFailed                           Traceback (most recent call last)
<timed exec> in <module>

# 中略

~\anaconda3\envs\ana0\lib\site-packages\hyperopt\base.py in best_trial(self)
    620         ]
    621         if not candidates:
--> 622             raise AllTrialsFailed
    623         losses = [float(t["result"]["loss"]) for t in candidates]
    624         if len(losses) == 0:

AllTrialsFailed: 

n_estimatorsに95.0という浮動小数点が渡されているため、怒られています。上に書いたように``args = args{ n_estimators':int(args['n_estimators'], ... } '''を冒頭に書くようにしましょう。
ただ、この方法はちょっと面倒くさいですよね。もっと効率的なやり方を知っている方がいらっしゃたら、是非教えてほしいです。

提出用のテストデータに対して予測する

上のテストデータは、コンペでは訓練データとして与えらているデータを分割して用意したものでした。訓練データを時系列に並べて、最後の20%の期間を評価用として使っていました。
コンペでは、訓練データ期間の後の2017/4/23～2017/5/31を提出用のテストデータとして用意してあり、店ごとの訪問客数を予測して、予測誤差を競います。
上の２つのモデルを使って、提出用のテストデータに対して予測を行い、Kaggleへ提出（Submit）して実際に評価してもらいました。
提出するデータの作成はこの記事の本題からはそれるので、ここにはコードは示しません。もし興味があればこのGitHubの002_Recruit_Restaurant_Visitors/002_measure_rmse_for_test_data.ipynb を参照してください。

評価指標は root means squared logarithmic error で与えられます。誤差ですので、小さい値ほど良いということになります。結果は以下の通りです。
GridSearchCV　0.60582
hyperopt     0.60593

どちらも大差なく、順位は約1700位でした。airレジデータだけ使って、日付と休日だけを考慮し、店の位置や季節性は考慮していないので、こんなものでしょうか。

まとめ

hyperoptを使ってサポートベクターの分類モデルとLightGBMによる回帰モデルのハイパーパラメータをチューニングしてみました。GridSearchCVなら、初めに粗く探索して、それから細かく探索して見つけるようなハイパーパラメータを、hyperoptなら１回で探索してくれて、精度も遜色ないモデルが作れました。
これからも使っていきたいと改めて思いました。

n_estimatorsを整数型としてきちんと渡すときなど、ハマったポイントも記載しました。パラメータ探索範囲の指定方法や、目的関数の与え方などの作法を覚える必要はあります。しかし、それほど難しいわけでもないので、一度やり方を覚えるといろんな場面で使えそうです。

コードはあまりきれいではありませんが、GitHubに公開しています。