xgboostでの予測
はじめに
こんにちは。はんぺんです。
最近、xgboostを使う機会があったので、それについてまとめようと思います。
使うデータはKaggleにあるBike Sharing Demandです。
コードはGitにあげています。
環境は以下の通りです
- macOS High Sierra 10.13.6
- Python 3.6
- Anaconda 5.6
コードの解説
Pythonで書きました。
少しずつ解説していきます。
インポート
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
%matplotlib inline import numpy as np import pandas as pd import matplotlib.pyplot as plt import xgboost as xgb from sklearn.model_selection import GridSearchCV from sklearn.model_selection import train_test_split import sklearn.metrics as metrics import pickle from matplotlib import cm import seaborn as sns from mpl_toolkits.mplot3d import Axes3D |
データ取得
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 |
df = pd.read_csv('../CSV/train.csv', parse_dates=['datetime']) df_test = pd.read_csv('../CSV/test.csv', parse_dates=['datetime']) # 学習のために型変換 def add_features(df): df['year'] = df['datetime'].dt.year df['month'] = df['datetime'].dt.month df['day'] = df['datetime'].dt.day df['dayofweek'] = df['datetime'].dt.dayofweek df['hour'] = df['datetime'].dt.hour add_features(df) add_features(df_test) # log1pの変換によってデータの分散が抑えられる df["count"] = df["count"].map(np.log1p) col = ['count', 'season', 'holiday', 'workingday', 'weather', 'temp', 'atemp', 'humidity', 'windspeed', 'year', 'month', 'day', 'dayofweek', 'hour'] col_test = ['season', 'holiday', 'workingday', 'weather', 'temp', 'atemp', 'humidity', 'windspeed', 'year', 'month', 'day', 'dayofweek', 'hour'] df = df.loc[:,col] df_test = df_test.loc[:,col_test] |
ここで、データおログ変換することでデータの分散が抑えられるので変換します。
データ分割
|
1 2 3 4 |
#目的変数の指定 objective = 'count' X_train_validate, X_test, y_train_validate, y_test = train_test_split(df.drop(objective, axis=1), df[objective].ravel()) |
ここでは、train-validationとtestにデータをわけます。
パラメータチューニング
max_depth & learning_rate
|
1 2 3 4 5 6 7 8 9 10 11 12 |
dtrain = xgb.DMatrix(X_train_validate.as_matrix(), label=y_train_validate.tolist()) fit_params = { 'eval_metric': 'rmse', 'eval_set': [[X_train_validate,y_train_validate]] } #グリッドサーチの範囲 params = { 'learning_rate': list(np.arange(0.05, 0.41, 0.05)), 'max_depth': list(np.arange(3, 11, 1)) } |
今回はまずmax_depthとlearning_rateをグリッドサーチで探索します。探索の範囲は,max_depthは3~10、learning_rateは0.5~0.4くらいで私はいつもやっています。
刻み幅は時間との兼ね合いもありますが、時間があるのなら0.2とかもっと細かくしてもよいです。
この二つを決めればとりあえずOKみたいなところある思ってます。(違ったらごめんなさい)
|
1 2 3 4 5 |
def GSfit(params): regressor = xgb.XGBRegressor(n_estimators=100) grid = GridSearchCV(regressor, params, cv=3, fit_params=fit_params, scoring='neg_mean_squared_error',verbose=2) grid.fit(X_train_validate,y_train_validate) return grid |
グリッドサーチのための関数を定義します。
クロスバリデーションはcv=3(データを３つに分ける)で行います。
また、sklearnのGridSearchの仕様でrmseではなくmseしかつかえないのでそれを使います。
ルートついてるかついてないかの違いなので、影響ないと考えています。
|
1 2 3 4 5 |
#Grid search grid = GSfit(params) grid_best_params = grid.best_params_ grid_scores_df = pd.DataFrame(grid.cv_results_) grid_scores_df.to_csv('../CSV/grid_scores_' + name_param + '_' + str(num_hue) + '.csv', index=False) |
グリッドサーチを行って結果を保存します。
データやPCの性能によって早さは異なりますが今回は私のMacbookproで1分ほどで終わりました。
この結果から、max_depth=7, learning_rate=0.1が最適なパラメータであることがわかります。
なんとなく好きなので3Dでもプロットしてみましたが、ヒートマップの方が見やすそうです。
n_estimators
|
1 2 3 4 5 |
#best n by CV cv=xgb.cv(grid_best_params, dtrain, num_boost_round=200, nfold=3) n_best = cv[cv['test-rmse-mean'] == cv['test-rmse-mean'].min()]['test-rmse-mean'].index[0] grid_best_params['n_estimators'] = n_best + 1 pd.io.json.json_normalize(grid_best_params).to_csv('../CSV/param_best_part_' +'.csv', index=False) |
これは何個の木を生成して回帰をするかというパラメータです。基本的に多くするほど学習エラーは下がりますが、ある値でバリデーションエラーは悪化して行く傾向にあります。
経験からすると70 〜130くらいになることが多いです。(今回は175と多めでした)
train-errorとvalidation-errorの可視化結果が以下のようになります。
これではよくわかりませんが、拡大すると、
175が最適な値であることがわかります。
以上のクロスバリデーションの結果より、n_estimators=175を採用します。
学習
|
1 2 3 4 5 |
#fit by best params regressor = xgb.XGBRegressor(learning_rate=grid_best_params['learning_rate'], max_depth=grid_best_params['max_depth'], n_estimators=grid_best_params['n_estimators']) regressor.fit(X_train_validate, y_train_validate, verbose=False) |
|
1 2 |
#save model pickle.dump(regressor, open('model.pkl', 'wb')) |
予測
|
1 2 3 4 |
# prediction result = regressor.predict(X_test) X_test[objective] = y_test X_test[objective + '_pred'] = result |
そしてtestデータに実際の値と予測値のカラムを追加して、評価を行なっていきます。
誤差率計算
|
1 2 3 4 5 6 7 8 9 |
# Metric Use By Kaggle def compute_rmsle(y_true, y_pred): if type(y_true) != np.ndarray: y_true = np.array(y_true) if type(y_pred) != np.ndarray: y_pred = np.array(y_pred) return(np.average((np.log1p(y_pred) - np.log1p(y_true))**2)**.5) |
|
1 |
print("RMSLE: {0}".format(compute_rmsle(X_test['count'].map(np.expm1),X_test['count_pred'].map(np.expm1)))) |
可視化
Feature importance
|
1 |
xgb.plot_importance(regressor) |
Xgboostには説明変数の中でどれが大事だったかを示す、Feature importanceというものがあります。
これを見ながらデータの特徴を考えて説明変数を追加して行くというのが一般的です。
散布図
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
plt.rcParams["font.size"] = 12 #フォントサイズの指定 fig = plt.figure(figsize=(6,6)) #図の生成 ax = fig.add_subplot(1,1,1) #軸の指定 plt.grid(True) ax.scatter(X_test[objective], X_test[objective + '_pred']) #データのプロット x = np.arange(0, 10, 0.1) ax.plot(x,x, color='red') ax.set_xlim([0, 8]) #x軸の範囲 ax.set_ylim([0, 8]) #y軸の範囲 ax.set_title('name') #グラフのタイトル ax.set_xlabel(objective) #x軸の名前 ax.set_ylabel(objective + '_pred') #y軸の名前 ax.legend(loc='upper left') #凡例の表示 plt.savefig('../Graph/' + 'name' + '_' + '_scatter.png', bbox_inches="tight") plt.show() |
実際の値と予測値がどれくらいあっているかを直感的にわかるようにグラフ化ました。
X軸に実績値、Y軸に予測値を取っています。
y=xの赤い線に点が乗っているほど性能の高いモデルということが言えます。
さいごに
いかがでしたでしょうか。
今回はxgboostでKaggleの自電車需要のデータを題材にして説明を行いました。
xgboostは非常に強力な機械学習のパッケージで、Kaggleのみならず業務データの分析においても非常に高い精度が出ることが多く、使えて損はないと思います。
私ももっと使いこなせるように頑張ろうと思います。
ありがとうございました。
ソースコードは以下においてあります。
コメントを残す