# '20/02/16更新：グラフのフォントサイズ調整やタイトル等の見栄え関連のコード更新
本記事では、決定木アルゴリズムの勾配ブースティング法のひとつLight Gradient Boosting Machine（LightGBM, lgbm）の回帰分析について、実務向けコードの雛形を記載しています。このライブラリはMicrosoft社製ですが、下記リンク先にあるように、ライセンスはMIT licenseでありフリーです。
▼インストール方法
conda install -c conda-forge lightgbm
▼以下、本記事で使用した分析データと、出力データ等のプログラムの仕様について説明します。
1. 分析するデータはcsv形式
csvファイルであれば分析できるように汎用性を持たせている。parameter以下で設定できる。本記事では、scikit-learn(サイキットラーン)に含まれているボストン住宅価格に関する回帰分析用のデータセットを使用した。Python scikit-learn付属のボストン市の住宅価格データ（Boston house prices dataset）をcsvファイル化する - HK29’s blog
2. 学習データとテストデータに分割する
下図は学習データの行列散布図である。scikit-learnの「import train_test_split」を使用していて、非常に簡便に書ける。例：X_train, X_test, y_train, y_test = train_test_split(X_df, Y_df, test_size=0.3)
下図はテストデータの行列散布図である。上記でtest_size=0.3と指定したように、全データの30％がテストデータとして割り当てられているため、視覚的にも量が少ないことがわかる。しかし、上図と凡そ同様な分布をしている様子がわかる。
3. 評価指標(メトリック)と繰り返し回数(イタレーション)の関係をグラフ化する
下図は平均絶対誤差 (Mean Absolute Error, MAE)です。各実測値と予測値の差を絶対値として、平均化した指標です。値が小さい程、精度が高いことを意味します。
下図は平均二乗誤差 (Mean Squared Error, MSE)です。各実測値と予測値の差を2乗して、平均化した指標です。値が小さい程、精度が高いことを意味します。
下図は二乗平均平方根誤差 (Root Mean Squared Error, RMSE)です。上記のMSEをRoot（ルート,平方根）とった指標です。値が小さい程、精度が高いことを意味します。
下図は、決定係数 (R2)です。yxを実測値、yaを実測値の平均値、fxを予測値とした時に次のように定義される。R2=1－Σ(yx－fx)^2／Σ(yx－ya)^2
式からわかるように、簡単に言うと差分に対して評価している指標である。値が1に近い程、精度が高いことを意味します。
以上が一般的な指標ですが、次のような散布図で視覚的にも適合度を確認する。
4. 設計変数（特徴量）について、目的変数Yに対する影響度の高い順にグラフ化する
上位の要因を設計因子として設計最適化と制御管理項目とすれば良い。
5. 作成した回帰モデルは、シリアル化により「.joblig」で保存する
作成した回帰モデルの再利用方法は、Python joblibで保存した回帰モデルを読み込んで実行する - HK29’s blogを参照。
■本プログラム
特記事項：LightGBMには枝刈り機能があり「early_stopping_rounds」を使用すれば良いです。そのため、最大試行回数「num_boost_round」は大きくとっていて構わない。
#!/usr/bin/env python3 # coding: utf-8 # LightGBM regressor import os, glob import pandas as pd import numpy as np import scipy as sp from scipy import optimize import math from sklearn.preprocessing import StandardScaler from sklearn.preprocessing import MinMaxScaler from sklearn.model_selection import train_test_split from sklearn.metrics import r2_score import lightgbm as lgbm import matplotlib.pyplot as plt import seaborn as sns from joblib import dump, load import datetime now = datetime.datetime.now() now = now.strftime("%y%m%d") ######################################## csvファイルからデータをロードし前処理する def my_load_csv(target_Yname): #check_flag, check_name, check_duplicates_columns = check_list df = pd.read_csv(file_path, skiprows=start_row - 1, encoding="utf-8") #"shift-jis") print("df -> " + str(df)) print(df.dtypes) # pandasの列毎の型確認。その列に文字列があるかで決まる。但し、object型は要素毎に異なる可能性がある。 df.replace('none', '0', inplace=True) # 置換。 inplace=Trueでdfが置換される(df=と新規に代入しなくてよい） df.replace('FALSE', '0', inplace=True) df_buf1 = df.iloc[:, :out_start_column] df_buf2 = df.iloc[:, out_start_column:] # ilocとreplaceは同時に使用できなかったため、一旦分割して処理する df_buf2.replace('0', np.nan, inplace=True) #文字列0をnumpy形式のNaNにする df_buf2.replace(0, np.nan, inplace=True) #数値0をnumpy形式のNaNにする #columns_list = df_buf2.columns #df_buf2[columns_list] = pd.to_numeric(df_buf2[columns_list], errors='coerce') #文字列をNaNへ置換する df = pd.concat([df_buf1, df_buf2], axis=1) #データフレームの列を連結して、ひとつのdfへ戻す df.to_csv("data_nan.csv", sep=",", header=True, index=False) print('check NaN -> ' + str(df.isnull().any())) # ひとつでもNaNがある列はTrue、そうでなければFalseと表示される。 print(df.dtypes) if check_duplicates_columns: #指定した列名があった場合、重複行を削除する df2 = df.drop_duplicates(check_duplicates_columns, keep='first') #最初の行が重複とみなされずに残す else: df2 = df print("df2 -> " + str(df2)) if check_flag: df3 = df2.dropna(subset=[target_Yname]) # NaNやNoneなどの欠損値を除外する場合 else: df3 = df2[df2[target_Yname].str.contains(check_name, na=False)] # 指定文字列の行を残す場合 print("df3 -> " + str(df3)) if target_value2ln: # 目的変数Yを対数に変換する場合 df4 = df3.assign(tmp = np.log(df3[target_Yname])) target_Yname = target_Yname + '_ln' tmp_columns = df4.columns.values print("tmp_columns -> " + str(tmp_columns)) tmp_columns[-1] = target_Yname df4.columns = tmp_columns print("df4 -> " + str(df4)) df4.to_csv("data_df4.csv", sep=",", header=True, index=False) DF = df4 else: DF = df3 x_DF = DF.iloc[:, input_columns[0]:input_columns[1]] print("x_DF -> " + str(x_DF)) y = DF.loc[:,target_Yname] x_DF_columns = x_DF.columns print("x_DF_columns -> " + str(x_DF_columns)) DF.to_csv("data_DF.csv", sep=",", header=True, index=False) print('check2 NaN -> ' + str(DF.isnull().any())) return x_DF, y, x_DF_columns ######################################## 説明変数Xを標準化する場合 def my_StandardScaler(X_df): scaler = StandardScaler() x_sc = scaler.fit_transform(X_df) # sklearnの上記のようなfit関連の実行後、データの書式がnumpy形式となりカラム(列)名が無くなる。 # そのため、再度、pandas形式に変換する(次の2行)。これはグラフの軸名で使用するため。 x_sc_df = pd.DataFrame(x_sc) # pndas形式のデータフレーム型へ変換 x_DF_columns = X_df.columns # カラム名(列名)を抽出 x_sc_df.columns = x_DF_columns # カラム名(列名)を挿入 return x_sc_df ######################################## 説明変数Xを正規化する場合 def my_MinMaxScaler(X_df): mscaler = MinMaxScaler() x_msc = mscaler.fit_transform(X_df) x_msc_df = pd.DataFrame(x_msc) x_DF_columns = X_df.columns x_msc_df.columns = x_DF_columns return x_msc_df ######################################## lgbm用のユーザー定義メトリックの関数 def my_r2_score(preds, data): y_true = data.get_label() y_pred = np.where(preds < 0, 0, preds) R2 = r2_score(y_true, y_pred) # 決定係数R2 return 'R2', R2, True ######################################## LightGBM def my_lgbm(X_df, Y_df, my_feature_names, target_Yname, title): # 訓練データとテストデータに分割 X_train, X_test, y_train, y_test = train_test_split(X_df, Y_df, test_size=0.3) # LightGBM用の変数にデータセット lgbm_train = lgbm.Dataset(X_train, y_train) lgbm_test = lgbm.Dataset(X_test, y_test, reference=lgbm_train) my_params = { 'objective': 'regression', 'learning_rate': 0.1, # 学習率 'max_depth': -1, # 木の数 (負の値で無制限) 'num_leaves': 9, # 枝葉の数 'metric': ('mean_absolute_error', 'mean_squared_error', 'rmse'), #メトリック https://lightgbm.readthedocs.io/en/latest/Parameters.html#metric-parameters 'drop_rate': 0.15, 'verbose': 0 } evaluation_results = {} # to record evaluation results for plotting # lgbm.trainが、scikit-learnのmodel.fitに相当する regr = lgbm.train(my_params, lgbm_train, num_boost_round = 500, # 最大試行数 early_stopping_rounds=15, # この数分、連続でメトリックに変化なければ終了する valid_sets = [lgbm_train, lgbm_test], valid_names=['train', 'test'], feature_name = list(my_feature_names), evals_result = evaluation_results, feval = my_r2_score, verbose_eval = 4) ##### 以下、データ分析の可視化など my_fontsize = 15 # メトリックの変化の過程をグラフ化 # 平均絶対誤差 ax = lgbm.plot_metric(evaluation_results, metric='l1', ylabel='mean_absolute_error') my_x_label = ax.get_xlabel() # xラベルをゲット my_y_label = ax.get_ylabel() # yラベルをゲット ax.set_xlabel(my_x_label, fontsize=my_fontsize) # フォントサイズを指定して再セット ax.set_ylabel(my_y_label, fontsize=my_fontsize) # フォントサイズを指定して再セット ax.set_title('MAE during training ' + target_Yname, fontsize=my_fontsize) # フォントサイズを指定して再セット ax.tick_params(labelsize=my_fontsize) # 目盛りのフォントサイズ plt.legend(loc='upper left', fontsize=my_fontsize) # 凡例 plt.tight_layout() plt.savefig(now + "_mean_absolute_error_" + target_Yname + '_' + title + '.png') plt.close() # 平均二乗誤差 ax = lgbm.plot_metric(evaluation_results, metric='l2', ylabel='mean_squared_error') my_x_label = ax.get_xlabel() # xラベルをゲット my_y_label = ax.get_ylabel() # yラベルをゲット ax.set_xlabel(my_x_label, fontsize=my_fontsize) # フォントサイズを指定して再セット ax.set_ylabel(my_y_label, fontsize=my_fontsize) # フォントサイズを指定して再セット ax.set_title('MSE during training ' + target_Yname, fontsize=my_fontsize) # フォントサイズを指定して再セット ax.tick_params(labelsize=my_fontsize) # 目盛りのフォントサイズ plt.legend(loc='upper left', fontsize=my_fontsize) # 凡例 plt.tight_layout() plt.savefig(now + "_mean_squared_error_" + target_Yname + '_' + title + '.png') plt.close() # 平均二乗偏差 ax = lgbm.plot_metric(evaluation_results, metric='rmse', ylabel='root_mean_squared_error') my_x_label = ax.get_xlabel() # xラベルをゲット my_y_label = ax.get_ylabel() # yラベルをゲット ax.set_xlabel(my_x_label, fontsize=my_fontsize) # フォントサイズを指定して再セット ax.set_ylabel(my_y_label, fontsize=my_fontsize) # フォントサイズを指定して再セット ax.set_title('RMSE during training ' + target_Yname, fontsize=my_fontsize) # フォントサイズを指定して再セット ax.tick_params(labelsize=my_fontsize) # 目盛りのフォントサイズ plt.legend(loc='upper left', fontsize=my_fontsize) # 凡例 plt.tight_layout() plt.savefig(now + "_root_mean_squared_error_" + target_Yname + '_' + title + '.png') plt.close() # ユーザー定義メトリックをリストで抽出 train_metric_r2 = evaluation_results['train']['R2'] test_metric_r2 = evaluation_results['test']['R2'] # R2の変化の過程をグラフ化 fig = plt.figure(figsize=(7, 5)) # (縦軸, 横軸) ax = fig.add_subplot(1,1,1) ax.set_xlabel('Iterations', fontsize=my_fontsize) ax.set_ylabel('R2', fontsize=my_fontsize) ax.set_title('R2 during training ' + target_Yname, fontsize=my_fontsize) ax.plot(train_metric_r2, label='train', c='royalblue', lw=2) ax.plot(test_metric_r2, label='test', c='orange', lw=2) ax.tick_params(labelsize=my_fontsize) plt.tight_layout() plt.grid() plt.legend(loc='upper left', numpoints=1, fontsize=my_fontsize) plt.savefig(now + "_R2_" + target_Yname + '_' + title + '.png') # 重要因子のグラフ化 ax = lgbm.plot_importance(regr, max_num_features=10, color='magenta', xlabel='Feature value', ylabel='', title='important factor') my_x_label = ax.get_xlabel() # xラベルをゲット my_title = ax.get_title() # タイトルをゲット ax.set_xlabel(my_x_label, fontsize=my_fontsize) # フォントサイズを指定して再セット ax.set_title(my_title + ' ' + target_Yname, fontsize=my_fontsize) # フォントサイズを指定して再セット ax.tick_params(labelsize=my_fontsize) # 目盛りのフォントサイズ plt.tight_layout() plt.savefig(now + "_importance_" + target_Yname + '_' + title + '.png') plt.close() # 重要因子のテキスト出力 importance_df = pd.DataFrame(regr.feature_importance(), index=my_feature_names, columns=['importance']) sorted_importance_df = importance_df.sort_values(by='importance', ascending=False) print(sorted_importance_df) if sorted_importance_df.shape[0] > 5: important_factor = list(sorted_importance_df.index[:5]) else: important_factor = list(sorted_importance_df) # 行列散布図 print(important_factor) my_scatter_matrix(X_train, target_Yname, important_factor, y=y_train, title='train_data') my_scatter_matrix(X_test, target_Yname, important_factor, y=y_test, title='test_data') my_scatter_matrix(X_df, target_Yname, important_factor, y=Y_df, title='all_data') ''' ax = lgbm.plot_split_value_histogram(regr, feature='RM', bins='auto') plt.savefig(now + "_split_value_histogram_" + target_Yname + '_' + title + '.png') plt.close() ax = lgbm.plot_tree(regr, tree_index=49, figsize=(12, 12), show_info=['split_gain']) plt.savefig(now + "_tree_" + target_Yname + '_' + title + '.png') plt.close() ''' # 決定木の分岐の可視化 graph = lgbm.create_tree_digraph(regr, tree_index=39, name='Tree40') graph.render(view=False) # テストデータで回帰予測して、決定係数R2を出力 y_pred = regr.predict(X_test, num_iteration=regr.best_iteration) R2 = r2_score(y_test, y_pred) print('R2={:.1%}'.format(R2)) # 全データで回帰予測（後ほどのグラフ化による可視化のため） pred_Y_np = regr.predict(X_df, num_iteration=regr.best_iteration) pred_Y_s = pd.Series(pred_Y_np) # 回帰モデルをオブジェクトファイルで保存 dump(regr, now + '_LightGBM_' + target_Yname + '_' + title \ + '_' + 'R2_{:.0}'.format(R2) + '.joblib') return pred_Y_s # 全Xデータに対する予測値YをpandasのSeries形式で返す ######################################## グラフ化 # 行列散布図 def my_scatter_matrix(X_df, target_Yname, important_factor, *, y=None, title='NA'): # *を挟むと、以降の変数をキーワード引数で強制する fig = plt.figure(figsize=(12,12)) ax = fig.add_subplot(1,1,1) x = X_df.loc[:,important_factor] # 抽出する列名を指定 if y is None: xy_df = x else: xy_df = pd.concat([x,y], axis=1) # 横方向に連結。axis=1 sns.set(style="ticks", font_scale=1.2) #, palette='gray', color_codes=True) g = sns.pairplot(xy_df, kind="reg", markers=".", diag_kind="kde", diag_kws=dict(shade=True), ) g.fig.suptitle(title + ' ' + target_Yname) #, fontsize=12) g.fig.subplots_adjust(top=0.9) plt.savefig(now + '_scatter_matrix_' + title + '_' + target_Yname + ".png") plt.close() # 相関性の確認1 オリジナル生データを青点で、回帰式によるデータを赤点でプロット def plot_scatter(Xname_for_graph, X, target_Yname, Y, pred_Y, title): save_fname = now + "_LightGBM_" + target_Yname + '_' + title #fig = plt.figure(num=2, figsize=(8, 6), dpi=100) fig = plt.figure(dpi=120) ax = fig.add_subplot(1,1,1) ax.set_title(title, fontsize=14) ax.set_xlabel(Xname_for_graph, fontsize=16) ax.set_ylabel(target_Yname, fontsize=16) ax.scatter(X[Xname_for_graph], Y, facecolors='none', edgecolors='b', label='Raw data') ax.scatter(X[Xname_for_graph], pred_Y, facecolors='none', edgecolors='r', label='LightGBM Regressor') ax.legend(loc='lower right', fontsize=14) ax.tick_params(labelsize=16) plt.grid(True) plt.tight_layout() #plt.show() plt.savefig(save_fname + '_scat.png') plt.close() # Close a figure window ### 近似式 y = ax 用の関数 def approximation_expression(x, a): return a * x # 相関性の確認2 横軸Y：縦軸Pred_Y def plot_reg(X, target_Yname, Y, pred_Y, title="original"): save_fname = target_Yname + "_LinearRegression_" + title # 欠損値を省いたことにより、連番でなくなっている可能性のあるインデックスをリセット Y = Y.reset_index(drop=True) # drop=Trueにより、元を削除して、新規インデックスを0から連番で作成 DF = pd.concat([Y.rename(target_Yname), pred_Y.rename(target_Yname + "_pred")], axis=1) DF.to_csv('PredY_vs_Y_' + title + ".csv", sep=',', index=False, header=True, encoding='utf-8') popt, pcov = optimize.curve_fit(approximation_expression, Y, pred_Y) ax = plt.figure(num=0, dpi=120).gca() #plt.rcParams["axes.labelsize"] = 15 ax.set_title("pred vs real " + title, fontsize=14) ax.set_xlabel(target_Yname, fontsize=16) ax.set_ylabel("Pred\n" + target_Yname, fontsize=16) rp = ax.scatter(x = target_Yname, y = target_Yname + "_pred", data = DF, facecolors="none", edgecolors='black') # purple # 生データの最小値と最大値 y_min = DF[target_Yname].min() y_max = DF[target_Yname].max() # 予測データの最小値と最大値 y_pred_min = DF[target_Yname + '_pred'].min() y_pred_max = DF[target_Yname + '_pred'].max() # プロットするためのxデータの範囲を決める x_min = min(y_min, y_pred_min) x_max = max(y_max, y_pred_max) # グラフのレンジを決める x_range = x_max - x_min # データの範囲が100以下か以上かで分ける print("x_range = x_max - x_min = " + str(x_range)) if x_max > 1: min_lim = 0 if x_range <= 10: max_lim = math.floor(x_max + 1) else: max_lim = math.floor(x_max + 10) rp.axes.set_xlim(min_lim, max_lim) rp.axes.set_ylim(min_lim, max_lim) else: max_lim = 0 if x_range <= 100: max_lim = math.floor(x_max - 1) else: max_lim = math.floor(x_max - 10) rp.axes.set_xlim(min_lim, max_lim) rp.axes.set_ylim(min_lim, max_lim) # 近似式プロットのためのデータを作成 x_approximation = np.linspace(min_lim, max_lim, 10) # numpyでxデータを作成 y_approximation = popt[0] * x_approximation # 近似式に代入してyデータを作成 line_approximation = ax.plot(x_approximation, y_approximation, linestyle = 'dashed', linewidth = 3, color='r') rp.axes.set_aspect('equal', adjustable='box') plt.grid(True) # 凡例を指定して記入する場合は、リストで指定する ax.legend([line_approximation[0]], ["y = {:.2f}x".format(popt[0])], loc='upper left', numpoints=1, fontsize=15) plt.tick_params(labelsize=15) plt.tight_layout() plt.savefig(save_fname + '_real_vs_pred.png') plt.close() ######################################## 指定ファイルを複数一括削除 def my_initialize(): for f in glob.glob('*.png'): os.remove(f) for f in glob.glob('data*'): os.remove(f) for f in glob.glob('20*'): os.remove(f) for f in glob.glob('Tree*'): os.remove(f) ######################################## メイン関数 def main(): # 初期化(前回実行時に生成されたファイルを削除) my_initialize() # 目的変数が複数ある場合、順番に実行する for target_Yname in target_Yname_list: # データセットからデータ抽出 X_df, Y_df, my_feature_names = my_load_csv(target_Yname) # 説明変数Xの前処理 if scalerjudge == 1: # 標準化する場合 X_df = my_StandardScaler(X_df) elif scalerjudge == 2: # 正規化する場合 X_df = my_MinMaxScaler(X_df) else: pass X_df.describe().to_csv('dataX_describe.csv', index=True, encoding='utf-8', mode='a', header=True) # 回帰分析の実行 pred_Y = my_lgbm(X_df, Y_df, my_feature_names, target_Yname, title) # グラフ作成 plot_scatter(Xname_for_graph, X_df, target_Yname, Y_df, pred_Y, title) plot_reg(X_df, target_Yname, Y_df, pred_Y, title) ######################################## if __name__ == "__main__": ##### parameter file_path = "boston_XYdata.csv" # 読み込む生データcsvファイル title = file_path[:-4] input_factors = 10 Xname_for_graph = 'RM' target_Yname_list = ['PRICE','LSTAT'] # 目的変数Y リストで複数指定できる ### 生データcsvファイルの読み込み設定 start_row = 1 # 読み込み開始行。1行目なら1, 3行目なら3, … target_value2ln = False # Yを対数変換するか？ check_flag = True # True:NaNやNoneなどの欠損値を除外する場合。False:指定文字列を残す場合 check_name = "ok" # 上記check_flagが「False」の場合に機能する。指定文字列を残す check_duplicates_columns=[] # ダブりを省く列名を指定。複数可。省かない場合は[]と空リストにする check_list = [check_flag, check_name, check_duplicates_columns] input_columns = (0, input_factors) # 設計因子の列番号(開始, 終了) out_start_column = 100 # この列以降のデータについて、欠損値の処理を行う scalerjudge = 0 # Xデータの前処理：生データそのままの場合は0、 標準化する場合は1、正規化する場合は2 ##### call function main() print("finished")
●参考資料
・ Parameters — LightGBM 2.3.2 documentation
・ Parameters Tuning — LightGBM 2.3.2 documentation
●関連リンク
ちなみに、決定木関連のひとつで、Random forest(ランダムフォレスト)があります。これは決定木を並列にブートストラップにより分析するデータをランダムに作成して分析した結果を多数決をとります。一方、本記事のLightGBMは決定木を直列に更新してゆく勾配ブースティング法のひとつであり、決定木を直列or並列に作成するかの違いです。
以上
<広告>
リンク
