'20/07/12更新：graphvizを使用しない（インストール作業が不要な）verです。
分析データはcsvファイルであれば良く、お試し用の作成例は右記参照Python scikit-learn 付属のボストン市の住宅価格データ（Boston house prices dataset）をcsvファイル化する - HK29’s blog
▼本プログラム実行後のフォルダ内には、下図のようなファイルを出力する仕様です。
1. 行列散布図：サンプルに使用した分布と、指標(目的変数)に対する設計因子(説明変数)の傾向がわかる。
2. 個別散布図：詳細な傾向がわかる
3. 特徴量の順：指標(目的変数)を決める重要因子(設計因子)とその感度がわかる
4. 適合度の変化の様子をグラフ化：決定係数R2のあてはまりの良さの変化を示す。
5. 適合度90％を超えた時点の回帰モデルをオブジェクトファイルで保存する(.joblib)
■本プログラム
#!/usr/bin/env python3 # -*- coding: utf-8 -*- import os, sys, glob import shutil import itertools import pandas as pd import numpy as np import scipy as sp from scipy import optimize import math from sklearn import linear_model from sklearn.ensemble import RandomForestRegressor from sklearn.preprocessing import StandardScaler from sklearn.preprocessing import MinMaxScaler import matplotlib matplotlib.use('Agg') # バックグラウンドで画像生成などする宣言 import matplotlib.pyplot as plt plt.rcParams.update({'figure.max_open_warning': 0}) import matplotlib.ticker as tick 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(path, check_column, check_duplicate_columns, start_row, input_columns, out_start_column, target_Yname, target_value2ln): df = pd.read_csv(path, skiprows=start_row-1, encoding="utf-8") 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へ戻す file_path = os.path.join(target_Yname, "data_nan.csv") df.to_csv(file_path, sep=",", header=True, index=False) print('check NaN -> ' + str(df.isnull().any())) # ひとつでもNaNがある列はTrue、そうでなければFalseと表示される。 print(df.dtypes) if check_duplicate_columns: #指定した列名があった場合、重複行を削除する df2 = df.drop_duplicates(check_duplicate_columns, keep='first') #最初の行が重複とみなされずに残す else: df2 = df print("df2 -> " + str(df2)) if check_column[0]: df3 = df2.dropna(subset=[check_column[1]]) # NaNやNoneなどの欠損値を除外する場合 else: df3 = df2[df2[check_column[1]].str.contains(check_column[2], 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)) file_path = os.path.join(target_Yname, "data_df4.csv") df4.to_csv(file_path, 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] DF_columns = x_DF.columns print("DF_columns -> " + str(DF_columns)) file_path = os.path.join(target_Yname, "data_DF.csv") DF.to_csv(file_path, sep=",", header=True, index=False) print('check2 NaN -> ' + str(DF.isnull().any())) return x_DF, y, DF_columns, target_Yname ######################################## 説明変数Xを標準化 ######################################## 各因子のデータを平均0,標準偏差を1へ変換 def my_scaler(X_df, my_feature_names): scaler = StandardScaler() x_sc = scaler.fit_transform(X_df) # sklearnの上記のようなfit関連の実行後、データの書式がnumpy形式となりカラム(列)名が無くなる。 # そのため、再度、pandas形式に変換する(次の2行)。これはグラフの軸名で使用するため。 x_sc = pd.DataFrame(x_sc) # pandasデータフレーム形式へ変換 x_sc.columns = my_feature_names # カラム名を追記 file_path = os.path.join(target_Yname, 'data_scalerX.csv') x_sc.to_csv(file_path, sep=',', index=False, encoding='utf-8', header=True) x_sc.describe().to_csv('data_describe.csv', index=True, encoding='utf-8', mode='w', header=True) return x_sc ######################################## 説明変数Xを正規化 ######################################## 各因子のデータ範囲を0～1へ変換 def my_mscaler(X_df, my_feature_names): mscaler = MinMaxScaler() x_msc = mscaler.fit_transform(X_df) # sklearnの上記のようなfit関連の実行後、データの書式がnumpy形式となりカラム(列)名が無くなる。 # そのため、再度、pandas形式に変換する(次の2行)。これはグラフの軸名で使用するため。 x_msc = pd.DataFrame(x_msc) # pandasデータフレーム形式へ変換 x_msc.columns = my_feature_names # カラム名を追記 file_path = os.path.join(target_Yname, 'data_scalerX.csv') x_msc.to_csv(file_path, sep=',', index=False, encoding='utf-8', header=True) file_path = os.path.join(target_Yname, 'data_describe.csv') x_msc.describe().to_csv(file_path, index=True, encoding='utf-8', mode='w', header=True) return x_msc ######################################## 線形重回帰分析 ######################################## 暫定的に重要因子の抽出(可視化のため横軸にしてグラフ化する) def my_linear_regression(X_df, Y, target_Yname): file_name = target_Yname + '_LinearRegression.png' file_path = os.path.join(target_Yname, file_name) clf = linear_model.LinearRegression(fit_intercept=True, normalize=False, copy_X=True, n_jobs=1) clf.fit(X_df, Y) print("LinearRegression") df = pd.DataFrame({"Name":X_df.columns, "Coefficients":clf.coef_}).sort_values(by='Coefficients') df = round(df, 5) print(df) print(clf.intercept_) # 切片 print(clf.score(X_df,Y)) # R^2 plt.figure(dpi=120) #df2 = df.sort_values(by='Coefficients', ascending=False)[:6] # 降順にソートして6つ抽出 df2 = df.sort_values(by='Coefficients', ascending=False)[:] df2 = df2.sort_values(by='Coefficients', ascending=True) # 昇順にソート x=df2.Name y=df2.Coefficients plt.barh(range(len(x)), y, tick_label=x, align="center", color="magenta", height=0.7) # 軸が文字列の場合、ソートされてしまう対策 plt.title('LinearRegression', fontsize=14) plt.xlabel('Coefficients', fontsize = 14) plt.tick_params(labelsize=16) plt.tight_layout() #plt.show() plt.savefig(file_path) print("principal_feature") max_value = max(df.Coefficients) print(max_value) print(df.dtypes) df=df.astype(str) print(df.dtypes) df3 = df[df['Coefficients'].str.contains(str(max_value))] print(df3) principal_feature = df3.Name.values.tolist()[0] return principal_feature ######################################## ランダムフォレストによる回帰分析 ######################################## 決定係数R^2を指標として、回帰関数を選定する def my_random_forest(X_df, Y, target_Yname, trees, depth, step, R2judge, pc, inter_factors): cnt=1 cntlist=[] scorelist=[] judge=0 n_depth = (depth[1] - depth[0] + step)//step n_trees = (trees[1] - trees[0] + step)//step for j in range(depth[0], depth[1], step): for i in range(trees[0], trees[1], step): print("cnt -> " + str(cnt)) # ランダムフォレストの実行（回帰式の作成） regr = RandomForestRegressor(n_estimators=i, max_depth=j, random_state=0) # 設定条件と決定係数の表示 regr = regr.fit(X_df,Y) # 設定条件全てを出力（default引数はは上記4つの他にもあるため） print("### regr") print(regr) score = regr.score(X_df,Y) # 決定係数R^2 print("### score") print(score) scorelist.append(score) cntlist.append(cnt) if (judge==0) & (score>=R2judge): # 重要因子の可視化 # 散布図 print("### regr.estimators_") estimators = regr.estimators_ print(estimators) print("### regr.feature_importances_" ) importances = regr.feature_importances_ print(importances) print("### sort") num = np.sort(importances) print(num) print("### element of sorted") indices = np.argsort(importances) print(indices) #fig = plt.figure(num=1, figsize=(8, 6), dpi=100) fig = plt.figure(dpi=120) ax = fig.add_subplot(1,1,1) ax.set_title('Important Factor Num' + str(cnt), fontsize=14) ax.barh(range(len(indices)), importances[indices], color='magenta', align='center') print(range(len(indices))) print(X_df.columns[indices]) loc = list(range(len(indices))) labels = X_df.columns[indices] ax.set_yticks(loc) ax.set_yticklabels(labels) ax.set_xlim([0, 1]) ax.tick_params(labelsize=16) plt.tight_layout() file_name = target_Yname + "_RandomForest_important_factor_" + str('{:02g}'.format(cnt)) + '.png' file_path = os.path.join(target_Yname, file_name) plt.savefig(file_path) plt.close() # Close a figure window # 横棒グラフ select_condition=[score,cnt,i,j] if input_factors >= 5: important_factor=X_df.columns[indices][-5:] # 優位上から5因子を抽出 else: important_factor=X_df.columns[-1*input_factors:] judge=1 print("judge No -> " + str(cnt)) # 回帰モデルをオブジェクトファイルで保存 file_name = now + '_' + target_Yname + "_RandomForestRegressor_" + str('{:02g}'.format(cnt)) + '.joblib' file_path = os.path.join(target_Yname, file_name) dump(regr, file_path) # 回帰式にXを代入してYを生成 pred_Y_np = regr.predict(X_df) print(pred_Y_np) print(type(pred_Y_np)) pred_Y = pd.Series(pred_Y_np) # 再現性を散布図で確認 plot_scatter(X_df, Y, target_Yname, pc, pred_Y, cnt) plot_reg(Y, target_Yname, pc, pred_Y, cnt) cnt += 1 # 決定係数の変化を可視化 #my_fontsize = 12 ax = plt.figure(num=0, dpi=120).gca() #, figsize=(8,8)) ax.set_title('the evolution of the coefficient of determination') #, fontsize=my_fontsize) ax.scatter(cntlist, scorelist, c="red") #, s=20*5) ax.plot(cntlist, scorelist, color="red") ax.set_xlabel('Number of Cycle') #, fontsize=my_fontsize) ax.set_ylabel('R^2') #, fontsize=my_fontsize) for i, j in zip(cntlist, scorelist): ax.text(i, j, '{:.2f}'.format(j), ha='center', va='bottom', color='b') #, fontsize=my_fontsize) ax.set_ylim([0, 1]) ax.xaxis.set_major_locator(tick.MultipleLocator(1)) ax.grid(which='both') #ax.tick_params(labelsize=my_fontsize) # 目盛りのサイズ plt.tight_layout() #plt.show() file_name = target_Yname + '_R2.png' file_path = os.path.join(target_Yname, file_name) plt.savefig(file_path) plt.close() # Close a figure window if judge==0: select_condition='N/A: MaxR2=' + str(score) + '<' + str(R2judge) print(select_condition) print('### scorelist') print(scorelist) sys.exit() return select_condition, important_factor ######################################## 散布図 # 相関性の確認1 オリジナル生データを青点で、回帰式によるデータを赤点でプロット def plot_scatter(X_df, Y, target_Yname, pc, pred_Y, cnt): file_name = target_Yname + "scatter_" + str('{:02g}'.format(cnt)) + '.png' file_path = os.path.join(target_Yname, file_name) #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('Scatter Num' + str(cnt), fontsize=14) ax.set_xlabel(pc, fontsize=16) ax.set_ylabel(target_Yname, fontsize=16) ax.scatter(X_df.loc[:,pc], Y, facecolors='none', edgecolors='b', label='Raw data') ax.scatter(X_df.loc[:,pc], pred_Y, facecolors='none', edgecolors='r', label='RandomForest') ax.legend(loc='lower right', fontsize=14) ax.tick_params(labelsize=16) plt.grid() plt.tight_layout() #plt.show() plt.savefig(file_path) plt.close() # Close a figure window ### 近似式 y = ax 用の関数 def approximation_expression(x, a): return a * x ### 近似式 y = ax + b 用の関数 def approximation_expression_with_y_intercept(x, a, b): return a * x + b # 相関性の確認2 横軸Y：縦軸Pred_Y def plot_reg(Y, target_Yname, pc, pred_Y, cnt): file_name = target_Yname + "_RandomForest_" + "{:02g}".format(cnt) + '_real_vs_pred' file_path = os.path.join(target_Yname, file_name) # 欠損値を省いたことにより、連番でなくなっている可能性のあるインデックスをリセット 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(file_path + ".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 " + "{:02g}".format(cnt), 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) print(x_min) print(x_max) # グラフのレンジを決める x_range = x_max - x_min # データの範囲が100以下か以上かで分ける print("x_range = x_max - x_min = " + str(x_range)) if x_max > 1: print("range A") min_lim = 0 if x_range <= 100: max_lim = math.floor(x_max + 1) else: max_lim = math.floor(x_max + 10) elif x_max >= 0: print("range B") min_lim = -0.1 max_lim = 1.1 elif x_min >= -1: print("range C") min_lim = -1.1 max_lim = 0.1 else: print("range D") max_lim = 0.1 if x_range <= 100: min_lim = math.floor(x_min - 1) else: min_lim = math.floor(x_min - 10) print(min_lim, max_lim) 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(file_path + '.png') plt.close() ######################################## 重要因子の上位6つ(important_factor)で作成する。 ### 行列散布図 def my_scatter_matrix(X_df, important_factor, *, y=None, n=1): # *を挟むと、以降の変数をキーワード引数で強制する file_name = target_Yname + '_2d_scatter_matrix.png' file_path = os.path.join(target_Yname, file_name) fig = plt.figure(figsize=(8,8)) 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.pairplot(xy_df, kind="reg", markers=".", diag_kind="kde", diag_kws=dict(shade=True), ) plt.savefig(file_path) plt.close() # Close a figure window ### 特徴量2因子の散布図xy pairplot def my_scatter_xy_pairplot(X_df, tick_setup, important_factor, *, y=None, n=1): # *を挟むと、以降の変数をキーワード引数で強制する file_name = target_Yname + '_2d_scatter_xy_pairplot.png' file_path = os.path.join(target_Yname, file_name) mySpecie = important_factor[-1] fig = plt.figure(figsize=(14,14)) 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 Yname = xy_df.columns[-1] Xnamelist = xy_df.columns[:-1].values.tolist() # 自分自身を含む場合は-1,含まない場合は-2 sns.set_style("whitegrid") # グラフに格子状を入れる g = sns.pairplot(xy_df, hue = mySpecie, palette='bwr', # 'gray' x_vars=Xnamelist, y_vars=[Yname], ) if tick_setup: g.set(ylim=(tick_setup[0], tick_setup[1])) #目盛り区間 plt.savefig(file_path) plt.close() # Close a figure window ### 特徴量2因子の散布図xy lmplot def my_scatter_xy_lmplot(X_df, target_Yname, tick_setup, important_factor, *, y=None, n=1): # *を挟むと、以降の変数をキーワード引数で強制する Yname = target_Yname Xname = important_factor[-1] mySpecie = important_factor[-2] file_name = target_Yname + '_2d_scatter_xy_lmplot_' + Xname +"_"+ mySpecie + ".png" file_path = os.path.join(target_Yname, file_name) fig, ax = plt.subplots(figsize=(12,10)) 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('notebook', 'whitegrid', 'dark', font_scale=1.4) # グラフに格子状を入れる g = sns.lmplot(Xname, Yname, xy_df, hue = mySpecie, palette='bwr', # 'Blues_d', # 'Reds_d', fit_reg = False, ) if tick_setup: g.set(ylim=(tick_setup[0], tick_setup[1])) #目盛り区間 plt.yticks(np.arange(tick_setup[0], tick_setup[1] + 0.01, tick_setup[2])) #目盛り間隔。+0.01としてる理由は0とすると0が表記されないため。 plt.axhline(y=tick_setup[3], lw=2, color='red', linestyle='dashed') # 判断基準のライン #ax.axhspan(ymin=, ymax=, color="lightcyan", alpha=0.3) # 指定範囲を塗りつぶす。縦方向の場合はaxvspan plt.savefig(file_path) plt.close() # Close a figure window ######################################## main関数 def main(): ### 初期化(前回実行時に生成されたファイルを削除) if os.path.isdir(target_Yname): shutil.rmtree(target_Yname) # 生成ファイルを保存するフォルダを作成 os.mkdir(target_Yname) ### データセットからデータ抽出 data = my_load_csv(path, check_column, check_duplicates_columns, start_row, input_columns, out_start_column, target_Yname, target_value2ln) X = data[0] Y = data[1] feature_names = data[2] print('### extraction data') print(X) print(Y) print(feature_names) ### 説明変数Xを標準化,正規化する？ if scalerjudge == 1: # 標準化する場合 X = my_scaler(X, feature_names) elif scalerjudge == 2: # 正規化する場合 X = my_mscaler(X, feature_names) else: pass print('### scalerX data') print(X) ### 線形重回帰分析により暫定重要因子を特定(グラフ化のために使用するだけ) pc = my_linear_regression(X, Y, target_Yname) print("### principal_feature ###") print(pc) ### ランダムフォレストにより回帰関数の作成(指定の決定係数を超えた関数を選定する) result = my_random_forest(X, Y, target_Yname, trees, depth, step, R2judge, pc, input_factors) print('### select_condition: 0score, 1num, 2trees, 3depth -> '+str(result[0])) print('### importnat factor -> '+str(result[1])) # 散布図 my_scatter_matrix(X, result[1], y=Y, n=1) my_scatter_xy_pairplot(X, tick_setup, result[1], y=Y, n=1) for contour_x1x2 in list(itertools.permutations(result[1], 2)): print("contour_x1x2 -> " + str(contour_x1x2)) print(tick_setup) my_scatter_xy_lmplot(X, target_Yname, tick_setup, contour_x1x2, y=Y, n=1) ###################################################################### メイン関数 if __name__ == "__main__": ### 基本パラメータ ### path = "boston.csv" # 読み込む生データcsvファイル input_factors = 13 target_Yname = 'PRICE' # 目的変数Y ### 生データcsvファイルの読み込み設定 start_row = 1 # 読み込み開始行。1行目なら1, 2行目なら2, … target_value2ln = False # Yを対数変換するか？ check_flag = True # True:NaNやNoneなどの欠損値を除外する場合。False:指定文字列を残す場合 check_column = (check_flag, target_Yname, 'ok') # 要素3は、Falseの場合に機能する。指定文字列を残す check_duplicates_columns=[] # ダブりを省く列名を指定。複数可。省かない場合は[]と空リストにする input_columns = (0, input_factors) # 設計因子の列番号(開始, 終了) out_start_column = 13 # この列以降のデータについて、欠損値の処理を行う ### データ分析の設定 scalerjudge = 0 # データXを標準化する場合は1, 正規化する場合は2, しない場合はその他任意 R2judge = 0.9 # 決定係数の判断基準。この値を初めて超えた場合の回帰を決定木で視覚化する。 contour_method = 'linear' # 'nearest' or 'linear' or 'cubic' contour_range = 6 #tick_setup = (0, 60, 10, 25) # 縦軸Yの目盛り設定で、(最小, 最大, 間隔, 判断基準)の順に入力 tick_setup = () #オートレンジにする場合は空() ### ランダムフォレストの設定 trees = [1, 8] # num_of_trees [min, max] depth = [1, 8] # max_depth [min, max] step = 3 # num of step ###### call function ###################### main() print("finished")
■ランダムフォレストに関する記事
下記リンクは、短いコードですがグラフ化は簡素にしています。単にランダムフォレストのコードを見たい人向け。
下記リンクは、ランダムフォレストに関する総まとめです。枝葉に分岐した結果も可視化します。dtreevizとgraphvizを使用します。
以上
<広告>
リンク
