Python 「PyCaret」ベストな回帰モデルを自動選定するautoml()

'20/12/13更新：PyCaretのverup2.2に伴いコードや図を編集、差し替えしました。

　本記事では、機械学習でお馴染みのボストン住宅価格のデータセットを用いて、多数ある回帰モデルからベスト回帰モデルを抽出してオブジェクトファイルに保存する雛形コードを載せました。データ分析する例題のcsvファイルは、次のリンク先を例に作成できます。Python scikit-learn付属のボストン市の住宅価格データ（Boston house prices dataset）をcsvファイル化する - PythonとVBAで世の中を便利にする　

 　本コードを実行すると、はじめに下図のように複数の回帰モデルを実行して各々5つの指標の結果を出力します。そして、指定した指標で並び替えによるランキング付けを行い、トップ5の回帰モデルをオブジェクトに格納します。

次に、そのトップ5を元に下図のようにたった一行で各々の回帰モデルのハイパーパラメータを自動調整します。

調整されたハイパーパラメータを下図のようにテキストファイルに保存して確認できるようにもしています。参考までに、回帰モデルの実使用はオブジェクト（.pklファイル）を扱うため、このテキストファイルを直接操作することはありません。

最終的に、下図のようにautoml()でベストな回帰モデルを自動で選定します。そして、predict_model()で回帰モデルの確認評価をすると、今回のベストな回帰モデルは、CatBoost Regressorということがわかります。

そして、interpret_model(best)と僅か一行のコードで、下図のようにシャープレイ値を可視化します。これにより、目的変数に対する感度の高い説明変数を視覚的に確認することが出来ます。横軸は目的変数に対する影響度（シャープレイ値）、縦軸は各説明変数です。赤色から青色は説明変数(Feature)が大きい、小さい場合に対応します。例えば、RMは部屋数を表し、その赤点が右側に多くプロットされています。このことより、部屋数が多い程、住宅価格が高くなる傾向であると理解することが出来ます。

そして、上記結果より感度の高い因子を元に、下図のように行列散布図で傾向を視覚的に確認することもできます。

下図は今回のベスト回帰モデルCatBoostの適合度を視覚的に確認したグラフです。横軸は生データで縦軸は回帰モデルによる予測結果です。比例y=x線上に点が乗ることが適合度が高いことを意味します。実機データを分析する場合は過剰適合は避けます。理由はばらつきによる飛び値があることが普通のためです。一方、シミュレーション結果をデータ分析する場合は過剰適合するくらいが良いです。理由はシミュレーション結果は理屈による数値計算結果であるため、飛び値などはあり得ないためです。但し、シミュレータのバグや数値計算誤差などによって、そういうのが発生するかもしれないことは頭に入れておき、行列散布図等で計算結果を確認することは必要です。

参考までに、下図は線形重回帰分析による適合度を可視化したグラフです。傾向は再現していますが、絶対値の当てはまりがよくないことが見た目にもわかります。

　　

　本コードを実行すると下図のように様々なデータを作成します。分析するデータのサマリー、行列散布図、基本統計量。他には、PyCaretによるベスト5の回帰モデル、ベストに選定された回帰モデルの特徴量の感度を可視化した図等です。PyCaretによるautoml()で選定されたベスト回帰モデルは.pklファイルで保存します。回帰モデルの利用は、そのオブジェクトファイルを読み出して使用します。

最後に、PyCaretが優れているのは実行が簡単なだけでなく、分析結果も簡単に行えます。plot_mode(best, plot='ほげほげ')といった僅か一行のコードで、次のように相関図、残渣図、特徴量のランキング図等を作成することが出来ます。

■本プログラム
　Jupyter Labのコードを.pyファイルへ変換して一部編集したものです。Anacona PromptなどのCUIで実行できます。

#!/usr/bin/env python
# coding: utf-8

# In[1]:


# データセットを読み込む
import pandas as pd
file_path = './boston_dataset.csv'
df = pd.read_csv(file_path, # csvファイルをpandasのDataFrame型で読み込む
                 sep=',',
                 skiprows=0,
                 header=0,
                 encoding='utf-8')
print(df)


# In[2]:


# 欠損値データの確認
print(df.isna().sum())


# In[3]:


# 基本統計量の確認
import datetime
now = datetime.datetime.now()
now = now.strftime("%y%m%d")
df_describe = df.describe()
df_describe.to_csv(now + '_describe.csv', # csvファイルに保存する
                   sep=',',
                   index=True,
                   encoding='utf-8')
print(df_describe)


# In[4]:


# ヒストグラムを作成する
import matplotlib
matplotlib.use('Agg') # バックグラウンドで画像生成などする宣言
import matplotlib.pyplot as plt
import warnings # わかっている警告（表示しなくてよい警告）は非表示にする
warnings.resetwarnings()
warnings.simplefilter('ignore', UserWarning)
fig = plt.figure(figsize=(8, 8), dpi=120)
ax = fig.add_subplot(1, 1, 1)
plt.subplots_adjust(left=0.1, right=0.95, bottom=0.1, top=0.95)
df.hist(ax=ax) # ヒストグラムを表示する
plt.tight_layout() # グラフ表示を整える
plt.savefig(now + '_hist.png')
plt.close()


# In[5]:


# PyCaret用のデータセット生成（前処理）
from pycaret.regression import *
set_data = setup(data = df,
    normalize = False, # 標準化するかどうか
    #normalize_method = 'zscore', # 標準化'zscore' 正規化'minmax'
                                # 絶対値を1'maxabs' ロバスト'robust
    #ignore_features = ['ZN', 'CHAS', B'],
    train_size = 0.75, # 訓練データの割合
    session_id = 1, # ランダムシード
    target = 'PRICE', # 目的変数
    #categorical_features=['列名'], # カテゴリ型を指定する場合               
    #numeric_features = ['列名'], # 数値型を指定する場合
    #remove_outliers = False, # 外れ値の除去
    #outliers_threshold = 0.05, # 0.05の場合、分布の両端裾の0.025%を除去
    #remove_multicollinearity = False, # マルチコ（多重共線性）除去
    #multicollinearity_threshold = 0.9, # 強い相関のある説明変数の片方を削除する
    #create_clusters = False, # クラスタリングにより新規特徴量を作成するかどうか
    #cluster_iter = 20, # default 20
    #polynomial_features = False, # 交互作用項を作成するか（新規特徴量）
    #trigonometry_features = False # 三角関数で作成するか
    #polynomial_degree = 2, # ↑次数で作成するか。[1, a, b, a^2, ab, b^2]
    #polynomial_threshold = 0.1, # 特徴量を残す判定しきい値
    #feature_interaction = False, # 交互作用（a * b）
    #feature_ratio = False, # 交互作用（a / b）
    #interaction_threshold = 0.01, # 特徴量を残す判定しきい値
    #pca = False, # 主成分分析による次元削減
    #pca_method = 'linear', # 'linear' 'kernel' 'incremental'
    #pca_components = 0.99, # 残す特徴量数。int型では数を指定。float型は割合
    silent=True, # 型チェックの手動確認のため一時ストップするか
)


# In[8]:


# モデルを比較してトップ5のモデルを抽出
top5 = compare_models(n_select = 5,
                      sort = 'MAE',
                      verbose = False) #default is 'R2'
                             #MAE, MSE, RMSE, R2, RMSLE, MAPE 
with open(now + '_top5.txt', 'a') as f:
    for i, my_model in enumerate(top5, start=1):
        f.write(str(i) + '\n')
        f.write(str(my_model) + '\n')


# In[9]:


# 各モデルのハイパーパラメータをデフォルト範囲で自動調整する
tuned_top5 = [tune_model(i, verbose = False) for i in top5]
with open(now + '_tuned_top5.txt', 'a') as f:
    for i, my_model in enumerate(tuned_top5, start=1):
        f.write(str(i) + '\n')
        f.write(str(my_model) + '\n')


# In[10]:


# ベスト回帰モデルを自動で選定する
best = automl(optimize = 'MAE')


# In[11]:


# モデルを評価する
predict_model(best)


# In[12]:


# 結果の可視化
plot_model(best, plot = 'error')
plt.savefig(now + '_error.png')
plt.close()
plot_model(best, plot = 'residuals')
plt.savefig(now + '_residuals.png')
plt.close()
plot_model(best, plot = 'feature')
plt.savefig(now + '_feature.png')
plt.close()
#plot_model(best, plot = 'cooks')
#plot_model(best, plot = 'learning')
plot_model(best, plot = 'manifold')
plt.savefig(now + '_manifold.png')
plt.close()
#evaluate_model(best)
plot_model(best, plot = 'parameter')
plt.savefig(now + '_parameter.png')
plt.close()

# 特徴量のサマリー
try: 
    interpret_model(best)
    plt.savefig(now + '_interpret_model.png')
    plt.close()
    #interpret_model(best, plot = 'correlation')
    #interpret_model(best, plot = 'reason', observation = 0)
except: 
    pass
# 構築した回帰モデルをデプロイする
#deploy_model(model = best, model_name = 'tuned_reg_aws', platform = 'aws', 
#             authentication =  {'bucket'  : 'pycaret-test'})


# In[13]:


# 目的変数とそれに対する感度の高い説明変数の上位4つを指定して行列散布図で視覚的に確認する
import seaborn as sns
# 列名を指定してデータ抽出
df2 = df.loc[:, ['LSTAT', 'RM', 'PTRATIO', 'AGE', 'PRICE']]
# 行列散布図
sns.set_context('talk')
mypairplot = sns.pairplot(df2,
                          kind="reg",
                          markers=".",
                          diag_kind="kde",
                          diag_kws=dict(shade=True),
                         )
plt.savefig(now + '_pairplot.png')
plt.close()


# In[14]:


# 回帰モデルを保存する
save_model(best, model_name = 'best_model_' + now)


# In[15]:


# 保存した回帰モデルをロードする
load_reg_model = load_model('best_model_' + now)
print(load_reg_model)


# In[16]:


# 新しいデータを読み込んで予測する。ここでは便宜上、既存データ全てを読み込む。本来は未知データを読み込む。
predict_df = predict_model(load_reg_model, data=df)
print(predict_df)


# In[17]:


# 評価指標を確認する。ここではPRICEとLabelの相関性
from pycaret.utils import check_metric
metric_list = ['MAE', 'MSE', 'RMSE', 'R2', 'RMSLE']
for metric in metric_list:
    print(metric, check_metric(predict_df.PRICE, predict_df.Label, metric))


# In[18]:


# 予測結果をファイルに保存する
predict_df.to_csv(now + '_predict.csv', index=False, encoding='utf-8')


# In[19]:


# 適合度（相関性）をグラフで確認する。
# catboost以外のモデルの場合は、上記のplot_model(best, plot = 'error')のたった一行できる。
import numpy as np
from scipy import optimize
import math

### 近似式 y = ax 用の関数 
def approximation_expression(x, a):
    return a * x

# 相関性の確認 横軸Y：縦軸Pred_Y
def plot_reg(DF, target_Yname, pred_Yname):
    Y = DF[target_Yname].values.tolist() # 生データをリストへ格納
    pred_Y = DF[pred_Yname].values.tolist() # 予測結果をリストへ格納
    
    # 近似式の作成
    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 ", fontsize=14)
    ax.set_xlabel(target_Yname, fontsize=14)
    ax.set_ylabel("Pred\n" + target_Yname, fontsize=14)
    rp = ax.scatter(x = target_Yname,
                    y = pred_Y,
                    data = DF,
                    facecolors="none",
                    edgecolors='black') # purple

    # 生データの最小値と最大値
    y_min = DF[target_Yname].min()
    y_max = DF[target_Yname].max()
    # 予測データの最小値と最大値
    y_pred_min = DF[pred_Yname].min()
    y_pred_max = DF[pred_Yname].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:
        print("range A")
        min_lim = 0
        if x_range <= 10:
            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 = {:.3f}x".format(popt[0])],
               loc='upper left',
               numpoints=1,
               fontsize=14)
    plt.tick_params(labelsize=14)
    plt.tight_layout() 
    #plt.show()
    plt.savefig(now + '_real_vs_pred.png') 
    plt.close()
    
plot_reg(predict_df, 'PRICE', 'Label')

本記事の雛形コードをJupyter Labで順番に実行して、機械学習の手順がわかる動画を下記リンク先に作成しました。

youtu.be

 以上

<広告>

リンク