Python 「scikit-learn」による様々な回帰分析（線形/勾配ブースティング/ランダムフォレスト/サポートベクトルマシン/ニューラルネットワーク）

'19/12/22更新
・汎用性を持たせるためにcsvファイルを指定して読み込むようにした
・NN(Neural Network)の隠れ層をリストとして、複数パターンを指定できるようにした
・得られた回帰モデル(回帰式)をシリアライズにより「.joblib」ファイルで保存する
　本記事では、scikit-learn(sklearn, サイキットラーン)で回帰分析できる次の7つの方法と結果例を載せました。
1. LinearRegression：線形重回帰
2. GradientBoostingRegressor：勾配ブースティング回帰
3. RandomForestRegressor：ランダムフォレスト回帰
4. Support Vector Regression(SVR)：サポートベクトルマシン(SVM)の回帰版
▼Neural Network「Multilayer perceptron」：ニューラルネットワーク
　5. Logistic Neurons：活性化関数にSigmoid(シグモイド)関数
　6. TanH Neurons：活性化関数にtanh(ハイパボリックタンジェント)関数
　7. Relu：活性化関数にRelu(Rectified Linear Unit)関数

　Anacondaをインストールした環境であれば本コードをコピペで実行可能です。使用したデータはsklearnに含まれているボストンデータセットです。一旦、csvファイル化する方法はこちら→

Python scikit-learn付属のボストン市の住宅価格データ（Boston house prices dataset）をcsvファイル化する - HK29’s blog

■本プログラム

回帰分析の結果例は本コードの下に記載しています。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import os, sys, glob, datetime
import numpy as np
import pandas as pd
import scipy as sp
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import seaborn as sns

from sklearn.linear_model import LinearRegression
from sklearn.svm import SVR
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn import neural_network 

from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import MinMaxScaler
from sklearn.model_selection import train_test_split

import joblib
######################################## 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="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_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))
        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] 
    DF_columns = x_DF.columns
    print("DF_columns -> " + str(DF_columns))
    DF.to_csv("data_DF.csv", 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 # カラム名を追記
    x_sc.to_csv('data_scalerX.csv', sep=',', index=False, encoding='utf-8', header=True)
    x_sc.describe().to_csv('data_describe.csv', index=True, encoding='utf-8',
                 mode='a', 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 # カラム名を追記
    x_msc.to_csv('data_scalerX.csv', sep=',', index=False, encoding='utf-8', header=True)
    x_msc.describe().to_csv('data_describe.csv', index=True, encoding='utf-8',
                 mode='a', header=True)

    return x_msc

######################################## 散布図
# 相関性の確認1 オリジナル生データを青点で、回帰式によるデータを赤点でプロット
def plot_scatter(X_df, Y, target_Yname, Xname_for_graph, pred_Y, cnt, name, R2):
    #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(Xname_for_graph, fontsize=16)
    ax.set_ylabel(target_Yname + '(target)', fontsize=16)
    ax.scatter(X_df.loc[:,Xname_for_graph], Y, facecolors='none', edgecolors='b', label='Raw data')
    ax.scatter(X_df.loc[:,Xname_for_graph], pred_Y, facecolors='none', edgecolors='r', label=name)
    ax.legend(loc='lower right', fontsize=14)
    ax.tick_params(labelsize=16)
    plt.text(4, 40, 'R2={:.1%}'.format(R2), size = 16, color = "red")
    plt.tight_layout()
    plt.grid()
    plt.savefig("scatter_" + target_Yname + '_' + name + '.png')
    plt.cla()   # Clear axis
    plt.clf()   # Clear figure
    plt.close() # Close a figure window

# 相関性の確認2 横軸Y：縦軸Pred_Y
def plot_reg(Y, target_Yname, Xname_for_graph, pred_Y, cnt, name):
    DF = pd.concat([Y, pred_Y], axis=1)
    DF.columns = [target_Yname, target_Yname + "_pred"]
    slope, intercept, r_value, p_value, std_err = sp.stats.linregress(Y, pred_Y)
    r_pearsonr, p_pearsonr = sp.stats.pearsonr(Y, pred_Y)
    ax = plt.figure(num=0, dpi=120).gca() 
    plt.rcParams["axes.labelsize"] = 15
    ax.set_title(name, fontsize=14)
    sns.regplot(x = target_Yname, y = target_Yname + "_pred", data = DF, color="purple",
                line_kws={'label':"y={0:.2f}x+{1:.1f}".format(slope,intercept)})
    ax.set_xlim(0, 60)
    ax.set_ylim(0, 60)
    ax.set_aspect('equal', adjustable='box')
    plt.grid(True)
    plt.legend(loc='upper left', fontsize=15)
    plt.tick_params(labelsize=15) 
    plt.savefig(target_Yname + '_' + name + "_real_vs_pred.png")
    plt.close()

    buf1_str = ("slope, intercept, r_value, p_value, std_err, r_pearsonr, p_pearsonr")
    buf2_list = [slope, intercept, r_value, p_value, std_err, r_pearsonr, p_pearsonr]
    buf2_list_str = [str(n) for n in buf2_list]
    buf2_list_str_with_csv = ','.join(buf2_list_str)
    with open("data_" + target_Yname + '_' + name + "_r_p.csv", 'a') as f:
        f.write(buf1_str + '\n')
        f.writelines(buf2_list_str_with_csv)

######################################## 回帰分析の実行
def run_regressor(reg_list, X_df, Y, Xname_for_graph, target_Yname):
    x_train, x_test, y_train, y_test= train_test_split(X_df, Y, test_size=0.2, random_state=42)
    '''
    print("train data")
    print(x_train, y_train)
    print("test data")
    print(x_test, y_test)
    '''
    now = datetime.datetime.now()
    now = now.strftime("%y%m%d")
    i = 1
    for name, regr in reg_list:
        regr.fit(x_train, y_train)
        joblib.dump(regr, now + '_' + target_Yname + '_' + name + '.joblib') 
        R2 = regr.score(x_test, y_test)
        print('{0:}. {1:} -> R2:{2:.1%}'.format(i, name, R2))
        
        pred_Y_np = regr.predict(X_df)
        pred_Y = pd.Series(pred_Y_np)
        plot_scatter(X_df, Y, target_Yname, Xname_for_graph, pred_Y, i, name, R2)
        plot_reg(Y, target_Yname, Xname_for_graph, pred_Y, i, name)
        i += 1

######################################## 指定ファイルを複数一括削除
def my_del_files():
    for f1 in glob.glob('*.png'):
        os.remove(f1)
    for f2 in glob.glob('data_*'):
        os.remove(f2)

######################################## メイン
if __name__ == '__main__':
    ### parameter ###
    path = "zzz_boston_XYdata.csv" # 読み込む生データcsvファイル
    input_factors = 13
    target_Yname = 'PRICE' # 目的変数Y
    Xname_for_graph = 'RM'
    
    ### 生データ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 = 14 # この列以降のデータについて、欠損値の処理を行う
    scalerjudge = 0 # データXを標準化する場合は1, 正規化する場合は2, しない場合はその他任意

    ##### call function
    ### 初期化(前回実行時に生成されたファイルを削除)
    my_del_files()
    
    ### データセットからデータ抽出
    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]
    target_Yname = data[3]
    print(X)
    print(Y)
    print(feature_names)
    print(target_Yname)

    ### 説明変数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)
    
    ### ニューラルネットワーク以外の回帰アルゴリズム
    regr_algorithm_list = [
        ('LinearRegression', LinearRegression()),
        ('GradientBoostingRegressor', GradientBoostingRegressor(random_state=1, n_estimators=100, max_depth=5)),
        ('RandomForestRegressor', RandomForestRegressor(n_estimators=9, max_depth=9, random_state=1)),
        ('Support Vector Regression', SVR(gamma='scale', C=18.0, epsilon=1.4)),
        ]
    run_regressor(regr_algorithm_list, X, Y, Xname_for_graph, target_Yname)

    ### ニューラルネットワーク
    hidden_layers_list = [(200,100,50,),
                          (200,100,50,25,),
                          (200,100,50,25,5,)]
    for hidden_layers in hidden_layers_list:
        s = len(hidden_layers)
        nn_layer_name = '_{}layer'.format(s)
        regr_algorithm_list = [
            #('Logistic Neurons' + nn_layer_name, neural_network.MLPRegressor(activation='logistic', hidden_layer_sizes=hidden_layers)),
            ('NN_TanH Neurons'+ nn_layer_name, neural_network.MLPRegressor(activation='tanh', hidden_layer_sizes=hidden_layers ,solver='lbfgs')),
            ('NN_Relu' + nn_layer_name, neural_network.MLPRegressor(activation='relu', hidden_layer_sizes=hidden_layers, solver='lbfgs'))
            ]
        run_regressor(regr_algorithm_list, X, Y, Xname_for_graph, target_Yname)

■解析結果

　ニューラルネットワークの結果は活性化関数別に4,5,6であり、各々について隠れ層が3層、4層、5層の場合の結果を載せた。しかし、それらの設定(ハイパーパラメータ)は最適化されたものではないため、適合度がよろしくない様子がわかる。つまり、ニューラルネットワークは適当に設定して実行すれば良い結果が得られるわけではなく、解析しようとしている各々のデータ群に対して、適切な活性化関数の選択と隠れ層やノード数などの設定値を考慮や探索して設定する必要がある。

　一方、ランダムフォレストの設定は比較的楽であり、決定木の枝葉の数を増やせばモデルの適合度は増す。注意点は全ての回帰分析に言えることで過剰適合である。

1. LinearRegression

f:id:HK29:20191102175440p:plain