Python 勾配ブースティングにおけるパラメータと調整方法について

また機械学習ネタです。機械学習の醍醐味である予測モデル作製において勾配ブースティング(Gradient Boosting)について今回は勉強したいと思います。以前のメルカリ価格予測（http://rautaku.hatenablog.com/entry/2017/12/22/195649/ ）でもLightGBMという勾配ブースティングを使ったんですが使ってみたものの結局中身がどういう動きをしているのかだとかパラメーターなどもどういうものがあるのかよく理解できていなかったのと、あまり検索してもパラメータとその調整についてわかりやすく書いているものがあまり見つからなかったので。このページ(https://www.analyticsvidhya.com/blog/2016/02/complete-guide-parameter-tuning-gradient-boosting-gbm-python/)　がわかりやすかったので参考にして実際のデータも使いながら書いていきます。

勾配ブースティングとは

まず勾配ブースティングとは複数の弱学習器を組み合わせたアンサンブル学習の一種で、その中でも1つずつ順番に弱学習器を構築していく手法です。一つずつ順番作っていくので、新しい弱学習器を作るときには、バギングと呼ばれるすべての弱学習器を独立に学習する方法と比べて並列計算ができず時間がかかります。時間はかかるのですが、過学習（高バリアンス）と学習不足（高バイアス）のバランスをうまく取ることができます。イメージとしては、ブースティングは最適化、バギングは平均化。　　

各ステップごとに弱学習器を構築して、前のステップで間違って識別されたものへのウェイトを重くして，次のステップで間違ったものをうまく識別できるようにしていきます。視覚的に見るとこのような感じです。

左から見ていくと、最初の弱学習結果で２個の青丸と６個の赤丸をうまく分けることができていますが、一方で３個の青丸が間違えてしまっています。そこで間違えた３個のプラスを重み付けして次の弱学習器ではこの３個を間違わないように分類して見ます。そうすると次は、赤丸２個と青丸１個が間違えてしまいました。今度はこれらを重み付けして分類してみます。３回目では赤丸２個と青丸１個が間違えてしまっています。

このような作業を繰り返していき、最終的な境界線を描いていきます。今回の絵では３回試行した後での境界線でも青丸１個が間違えてしまっています。実際ぽくていいのではないかなと思います。この辺の最終モデルのちょうど良いバランス（適切なバリアンスとバイアス）を探すために分類回数や数値など適切なパラメータ調整が必要になります。

パラメータについて

GBMの調整できるパラメータはハイパーパラメータと呼ばれるそうですが、このハイパーパラメータは大まかにいうと３種類に分けることができます。 1. Tree-Specific Parameters: ツリー固有パラメータ
ツリーの形状や動作などを規定するパラメータ 2. Boosting Parameters: ブースティングパラメータ
勾配動作などを規定するパラメータ 3. Miscellaneous Parameters: 他パラメータ
上記２枠以外のパラメータ

ツリー固有パラメータ

ツリー固有パラメータは以下のようなものがあります（イメージ図を書いてみましたが間違いがあればご指摘下さい><）

ブースティングパラメータ

ブースティングパラメータは以下のようなものがあります

他パラメータ

他パラメータは以下のようなものがあります。今回はこれらのパラメータの調整は行いませんが参考に。

実データを用いたパラメータ調整

各パラメータを見たところで、実際のデータでのGBMのパラメータ調整の影響を確認していきたいと思います。
今回は参考サイトのと同じ、『インドでの貸付者情報に基づくローン支払い』に関するデータを使います。生データはこちらから

インドのクレジットカード会社かどこかのデータぽいです。ローンを組んだ人の収入情報や住所、年齢、その他個人情報からその人がちゃんとローンを返したかという情報が入っています。クレジットカード会社側はこのような個人情報から顧客にどのくらい貸付できるかどのくらいで返せるだろうかなど試算しているんですね。怖い怖い。。一つデータの中に聞き慣れない単語が入っていました。"EMI"という単語ですがこれは、Equated Monthly Installmentの略で、インドでのクレジットカードを利用した定額月分割払いの総称らしいです。言わば月払いのローンのことです。インドはクレジットカードの需要がかなり高いらしく、EMIという制度が大流行したそうです（ちょっと調べていると、借金が膨らみすぎる人が多発したらしく今は禁止されてるような内容もありました。）「EMI Liya Hai To Chukana Padega(EMIで借りた金はちゃんと返せよ)」という映画もあったそうです。ポスターがインパクトありすぎたので貼っておきますね。

データ前処理

今回、前処理については参考サイトですでに行われたものを用います。<a　href="https://github.com/aarshayj/Analytics_Vidhya/tree/master/Articles/Parameter_Tuning_GBM_with_Example">こちら(github)からダウンロードします。
前処理前のデータで一つだけ遊んでみました。地域毎の支払率の可視化です。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import folium
from folium import plugins
%matplotlib inline

train = pd.read_csv("train.csv", encoding="iso-8859-1")
num = 50
city_list = train["City"].value_counts().head(num).index
disbursed_df = pd.DataFrame(index=range(num), columns=["City", "ratio", "lat", "lon"])
disbursed_df["City"] = city_list
disbursed_ratio_list = []
for c in city_list:
    city_count = train.loc[train["City"]==c, "City"].count()
    disbursed_count = train.loc[train["City"]==c, "Disbursed"].sum()
    disbursed_ratio = disbursed_count/city_count
    disbursed_ratio_list.append(disbursed_ratio)
disbursed_df["ratio"] = disbursed_ratio_list

# http://openweathermap.org/help/city_list.txt: get each city's lattitude and longtitude
world_city_df = pd.read_table("http://openweathermap.org/help/city_list.txt", encoding="iso-8859-1")
india_city_df = world_city_df.loc[world_city_df["countryCode"]=="IN"]

def addLatLon(disbursed_df, india_city_df):
    disbursed_df_city = list(disbursed_df["City"])
    disbursed_df_lat = []
    disbursed_df_lon = []
    for d in disbursed_df_city:
        lat_v = india_city_df[india_city_df["nm"]==d]["lat"].values
        lon_v = india_city_df[india_city_df["nm"]==d]["lon"].values
        try:
            lat_v = lat_v[0]
            lon_v = lon_v[0]
        except:
            lat_v = np.nan
            lon_v = np.nan
        disbursed_df_lat.append(lat_v)
        disbursed_df_lon.append(lon_v)
    disbursed_df["lat"] = disbursed_df_lat
    disbursed_df["lon"] = disbursed_df_lon
    return disbursed_df

disbursed_df = addLatLon(disbursed_df, india_city_df)
latlon_df = disbursed_df.dropna().loc[:, ["lat","lon","ratio"]]

# For Folium Heatmap, make tuple list
latlon_list = []
for l in latlon_df.values:
    latlon_list.append(tuple(l))

# Make Heatmap by using Folium
map = folium.Map(location=[disbursed_df["lat"].median(), disbursed_df["lon"].median()], zoom_start=5)
map.add_child(plugins.HeatMap(latlon_list))
map

Foliumというパッケージを使って、地域別のローン支払い率のヒートマップを地図上に可視化してみました。インドの知識はほとんどないのでこのマッピングを見て何も感じるものがないのですが、都市部が高いとかあるのではないでしょうか。また今度日本のデータ系で何か可視化してみたいです。

まぁ確かにインドのデータだとわかったところで、本題のパラメータ調整とその影響の調査に戻りたいと思います。

まずはデータのインポートやモデルフィット関数の準備です。 modelfit関数の中身は、分類アルゴリズム（今回はGradientBoostingClassifier）を使ってデータを分類して、精度評価として正確度評価(Accuracy)、AUC評価(Area under the ROC curve)、交差検定でスコア評価し、特徴量の影響度グラフと共に出力します。

import pandas as pd
import numpy as np
from sklearn.ensemble import GradientBoostingClassifier  #GBM algorithm
from sklearn import cross_validation, metrics   #Additional scklearn functions
from sklearn.grid_search import GridSearchCV   #Perforing grid search
import matplotlib.pylab as plt
%matplotlib inline
from matplotlib.pylab import rcParams
import seaborn as sns
rcParams['figure.figsize'] = 10, 4

train = pd.read_csv('train_modified.csv')
target = 'Disbursed'
IDcol = 'ID'

def modelfit(alg, dtrain, predictors, performCV=True, printFeatureImportance=True, cv_folds=5):
    #Fit the algorithm on the data
    alg.fit(dtrain[predictors], dtrain['Disbursed'])
        
    #Predict training set:
    dtrain_predictions = alg.predict(dtrain[predictors])
    dtrain_predprob = alg.predict_proba(dtrain[predictors])[:,1]
    
    #Perform cross-validation:
    if performCV:
        cv_score = cross_validation.cross_val_score(alg, dtrain[predictors], dtrain['Disbursed'], cv=cv_folds, scoring='roc_auc')
    
    #Print model report:
    print ("Model Report")
    print ("Accuracy : {:.4f}".format(metrics.accuracy_score(dtrain['Disbursed'].values, dtrain_predictions)))
    print ("AUC Score (Train): {:.4f}".format(metrics.roc_auc_score(dtrain['Disbursed'], dtrain_predprob)))
    
    if performCV:
        print ("CV Score : Mean - {:.6f} | Std - {:.6f} | Min - {:.6f} | Max - {:.6f}".format(np.mean(cv_score),np.std(cv_score),np.min(cv_score),np.max(cv_score)))
        
    #Print Feature Importance:
    if printFeatureImportance:
        feat_imp = pd.Series(alg.feature_importances_, predictors).sort_values(ascending=False)
        sns.set_palette("husl")
        sns.barplot(feat_imp.head(10).index, feat_imp.head(10).values)
        plt.title('Top10 Feature Importances')
        plt.ylabel('Feature Importance Score')
        plt.xticks(rotation=60)
        plt.show()

準備ができたところで、まずは標準値のまま勾配ブースティングをかけて見ます。random_stateは"10"で固定します。

#Choose all predictors except target & IDcols
predictors = [x for x in train.columns if x not in [target, IDcol]]
gbm0 = GradientBoostingClassifier(random_state=10)
modelfit(gbm0, train, predictors)

Model Report
Accuracy : 0.9856
AUC Score (Train): 0.8584
CV Score : Mean - 0.830063 | Std - 0.010523 | Min - 0.815374 | Max - 0.843712

交差検定スコアとして、0.830063を得ました。これをパラメータ調整により上回ることを目標とします。

パラメータ調整アプローチ

前述のように調整すべきパラメータは2種類あります。ツリー固有パラメータとブースティングパラメータです。 learning_rateには最適値というものはなくツリー数(n_estimatiors)が十分にあるのであれば小さければ小さいほど良いですが、あまりに小さすぎると計算負荷が大きく処理に時間がかかってしまいます。

以下の点に気をつけるべきです。
・比較的高めのlearnig_rateを設定する。標準値は0.1だが0.05~0.2あたりを試すべき
・設定したlearnig_rateに対してツリー数を設定する。40〜70くらい。いろんなパターンを試すためにも処理スピードにも気をつける・learning_rateとツリー数を決めたらツリー固有パラメータを調整する。・learning_rateを低くしたなら、より堅牢なモデルを得る為、それに反比例してツリー数を増やす。

ツリー固有パラメータ調整のためのlearning_rateとツリー数(n_estimators)の検証

learning_rateを0.1に固定し、n_estimatorsを20~80まで10個間隔で検証して見ます。

他パラメータは以下で固定します。

min_samples_split = 500 : 総サンプル数の0.5~1%が一般的で、今回のようなクラスに属するサンプル数に偏りがある不均衡データの場合低い方（0.5%）を使います。
min_samples_leaf = 50 : 直感で選びます。過学習を防ぐために使うのと不均衡データなのであまりに小さい数は採用しません。
max_depth = 8 : サンプル数に応じて5~8ぐらいで選択。今回は87000×49と結構大きいので8。
max_features = ‘sqrt’: sqrtで始めるのが一般的である。
subsample = 0.8 : これも一般的に使われる値

GridSearchCVを使ってn_estimatorsを振った際の結果変化を検証します。

predictors = [x for x in train.columns if x not in [target, IDcol]]
param_test1 = {'n_estimators':list(range(20,81,10))}
gsearch1 = GridSearchCV(estimator = GradientBoostingClassifier(learning_rate=0.1, min_samples_split=500,min_samples_leaf=50,max_depth=8,max_features='sqrt',subsample=0.8,random_state=10), 
param_grid = param_test1, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch1.fit(train[predictors],train[target])
mean_list = []
for s in gsearch1.grid_scores_:
    mean_list.append(s[1])
plt.figure(figsize=(4,2))
plt.plot(list(range(20,81,10)), mean_list)
plt.title("GridSearchCV Score")
plt.xlabel("n_estimators")
plt.show()
gsearch1.grid_scores_, gsearch1.best_params_

([mean: 0.83234, std: 0.01247, params: {'n_estimators': 20},
  mean: 0.83432, std: 0.01280, params: {'n_estimators': 30},
  mean: 0.83564, std: 0.01126, params: {'n_estimators': 40},
  mean: 0.83666, std: 0.01144, params: {'n_estimators': 50},
  mean: 0.83682, std: 0.01193, params: {'n_estimators': 60},
  mean: 0.83612, std: 0.01242, params: {'n_estimators': 70},
  mean: 0.83564, std: 0.01238, params: {'n_estimators': 80}],
 {'n_estimators': 60})

ツリー数60がベストスコアを出しました。この先はn_estimators=60として進めていきます。

from mpl_toolkits.mplot3d import Axes3D
import time
start_time = time.time()
param_test2 = {'max_depth':list(range(5,16,2)), 'min_samples_split':list(range(200,1001,200))}
gsearch2 = GridSearchCV(estimator = GradientBoostingClassifier(learning_rate=0.1, n_estimators=60, max_features='sqrt', subsample=0.8, random_state=10), 
param_grid = param_test2, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch2.fit(train[predictors],train[target])
x=[]
y=[]
z=[]
for s in gsearch2.grid_scores_:
    x.append(s[0]["max_depth"])
    y.append(s[0]["min_samples_split"])
    z.append(s[1])
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(x, y, z)
plt.show()
print('Showing results took {} secs.'.format(time.time() - start_time))
gsearch2.grid_scores_, gsearch2.best_params_, gsearch2.best_score_

Showing results took 593.9939701557159 secs.





([mean: 0.83094, std: 0.01075, params: {'max_depth': 5, 'min_samples_split': 200},
  mean: 0.83225, std: 0.01140, params: {'max_depth': 5, 'min_samples_split': 400},
  mean: 0.83265, std: 0.01255, params: {'max_depth': 5, 'min_samples_split': 600},
  mean: 0.83239, std: 0.01121, params: {'max_depth': 5, 'min_samples_split': 800},
  mean: 0.83203, std: 0.01215, params: {'max_depth': 5, 'min_samples_split': 1000},
  mean: 0.83260, std: 0.00662, params: {'max_depth': 7, 'min_samples_split': 200},
  mean: 0.83504, std: 0.01166, params: {'max_depth': 7, 'min_samples_split': 400},
  mean: 0.83517, std: 0.00963, params: {'max_depth': 7, 'min_samples_split': 600},
  mean: 0.83528, std: 0.01013, params: {'max_depth': 7, 'min_samples_split': 800},
  mean: 0.83471, std: 0.01230, params: {'max_depth': 7, 'min_samples_split': 1000},
  mean: 0.83094, std: 0.01107, params: {'max_depth': 9, 'min_samples_split': 200},
  mean: 0.83384, std: 0.01082, params: {'max_depth': 9, 'min_samples_split': 400},
  mean: 0.83522, std: 0.01034, params: {'max_depth': 9, 'min_samples_split': 600},
  mean: 0.83252, std: 0.00916, params: {'max_depth': 9, 'min_samples_split': 800},
  mean: 0.83356, std: 0.01169, params: {'max_depth': 9, 'min_samples_split': 1000},
  mean: 0.82438, std: 0.01229, params: {'max_depth': 11, 'min_samples_split': 200},
  mean: 0.82908, std: 0.01067, params: {'max_depth': 11, 'min_samples_split': 400},
  mean: 0.82865, std: 0.01038, params: {'max_depth': 11, 'min_samples_split': 600},
  mean: 0.83004, std: 0.01081, params: {'max_depth': 11, 'min_samples_split': 800},
  mean: 0.83232, std: 0.01079, params: {'max_depth': 11, 'min_samples_split': 1000},
  mean: 0.81861, std: 0.01130, params: {'max_depth': 13, 'min_samples_split': 200},
  mean: 0.82383, std: 0.00630, params: {'max_depth': 13, 'min_samples_split': 400},
  mean: 0.82776, std: 0.00844, params: {'max_depth': 13, 'min_samples_split': 600},
  mean: 0.83018, std: 0.01263, params: {'max_depth': 13, 'min_samples_split': 800},
  mean: 0.82855, std: 0.01021, params: {'max_depth': 13, 'min_samples_split': 1000},
  mean: 0.81289, std: 0.01212, params: {'max_depth': 15, 'min_samples_split': 200},
  mean: 0.82010, std: 0.01238, params: {'max_depth': 15, 'min_samples_split': 400},
  mean: 0.82385, std: 0.01026, params: {'max_depth': 15, 'min_samples_split': 600},
  mean: 0.82805, std: 0.01025, params: {'max_depth': 15, 'min_samples_split': 800},
  mean: 0.82918, std: 0.01099, params: {'max_depth': 15, 'min_samples_split': 1000}],
 {'max_depth': 7, 'min_samples_split': 800},
 0.8352763290658064)

'max_depth': 7, 'min_samples_split': 800がベストスコアでした。どんどん行きましょう、min_samples_leaf、max_featuresを順に確認して行きます。

param_test3 = {'min_samples_leaf':list(range(30,71,10))}
gsearch3 = GridSearchCV(estimator = GradientBoostingClassifier(learning_rate=0.1, n_estimators=60, min_samples_split=800, max_depth=7,max_features='sqrt', subsample=0.8, random_state=10), 
param_grid = param_test3, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch3.fit(train[predictors],train[target])
mean_list = []
for s in gsearch3.grid_scores_:
    mean_list.append(s[1])
plt.figure(figsize=(4,2))
plt.plot(list(range(30,71,10)), mean_list)
plt.title("GridSearchCV Score")
plt.xlabel("min_samples_leaf")
plt.show()
gsearch3.grid_scores_, gsearch3.best_params_, gsearch3.best_score_

([mean: 0.83565, std: 0.01055, params: {'min_samples_leaf': 30},
  mean: 0.83599, std: 0.01158, params: {'min_samples_leaf': 40},
  mean: 0.83488, std: 0.00898, params: {'min_samples_leaf': 50},
  mean: 0.83634, std: 0.01103, params: {'min_samples_leaf': 60},
  mean: 0.83609, std: 0.00992, params: {'min_samples_leaf': 70}],
 {'min_samples_leaf': 60},
 0.8363437200261575)

param_test4 = {'max_features':list(range(7,20,2))}
gsearch4 = GridSearchCV(estimator = GradientBoostingClassifier(learning_rate=0.1, n_estimators=60,max_depth=7, min_samples_split=800, min_samples_leaf=60, subsample=0.8, random_state=10),
param_grid = param_test4, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch4.fit(train[predictors],train[target])
mean_list = []
for s in gsearch4.grid_scores_:
    mean_list.append(s[1])
plt.figure(figsize=(4,2))
plt.plot(list(range(7,20,2)), mean_list)
plt.title("GridSearchCV Score")
plt.xlabel("max_features")
plt.show()
gsearch4.grid_scores_, gsearch4.best_params_, gsearch4.best_score_

([mean: 0.83442, std: 0.01179, params: {'max_features': 7},
  mean: 0.83606, std: 0.01068, params: {'max_features': 9},
  mean: 0.83704, std: 0.00725, params: {'max_features': 11},
  mean: 0.83873, std: 0.00966, params: {'max_features': 13},
  mean: 0.83671, std: 0.01033, params: {'max_features': 15},
  mean: 0.83629, std: 0.00860, params: {'max_features': 17},
  mean: 0.83582, std: 0.00988, params: {'max_features': 19}],
 {'max_features': 13},
 0.8387261096220344)

各パラメータのGridSearchCVによる調整した値が算出できました。
n_estimators = 60
max_depth = 7
min_samples_split = 800
min_samples_leaf = 60
max_features = 13
これらをGradientBoostingClassifierのパラメータに設定して再度評価してみましょう。標準値でのCVスコアは0.830063でしたのでこれを上回れば良いです。

#Choose all predictors except target & IDcols
predictors = [x for x in train.columns if x not in [target, IDcol]]
gbm1 = GradientBoostingClassifier(learning_rate=0.1, n_estimators=60,max_depth=7,min_samples_split=800, min_samples_leaf=60, subsample=0.8, random_state=10, max_features=13)
modelfit(gbm1, train, predictors)

Model Report
Accuracy : 0.9854
AUC Score (Train): 0.8871
CV Score : Mean - 0.838726 | Std - 0.009660 | Min - 0.829550 | Max - 0.851014

CVスコア=0.838726でした。無事上回りましたので良かったです！
という形で勾配ブースティングのパラメータ調整の影響を実感することができました。より良いスコアを出そうとしたらsubsampleやlearning_rateとツリー数のバランス、warm_startなどを調整していくことで、どんどんスコアを上げていくことができると思います。今回は以上にしておきます！

データとかITとか統計とか好きなこと

データや統計に関すること＆ITなどが好きで仕事にも生かせてます。音楽とかyoutubeも好き。