いろんなモデルを勉強がてら使ってみよう
ということで、MLモデルを使って有機材料のPL波長を予測して行きましょう。前回は記述子の計算を行いデータセットの作成を行いました。
第1回：記述子計算とデータ可視化
https://qiita.com/Osarunokagoya/items/cc0f79e3a3d959de8635
↑記述子の計算が終わっていない人はまずこちら！
今回はこの図でいうところの③をやっていきます。
前回作ったデータセットを元に、どの予測モデルが優れているのか見ていきましょう。今回は、PLS, RandomForest, XGBoost, LightGBM, NNを使っていきます。AutoML(pycaret)使って、さっと一覧を表示させるのもありですが、今回は中身を理解するため実際に一つずつ実装していきます。
1. データ読み込み
必要なライブラリを読み込みます。
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import pandas as pd
from sklearn.cross_decomposition import PLSRegression
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split, KFold, cross_val_predict
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
from sklearn.preprocessing import StandardScaler
import xgboost as xgb
import lightgbm as lgb
import optuna
import torch
import torch.nn as nn
from torch.utils.data import Dataset, DataLoader
import joblib
import os
前回同様、データセットを読み込みましょう。今回はmordredの3次元記述子を読み込みますが、RDKitなどで実行したい場合はdescriptor_type = 'rdkit'としてください。第1回で問題なく保存できていれば、どの記述子でも実行できるようにしています。
# 入力：読み込みたい記述子のタイプを選択
descriptor_type = 'mordred_3d' # 'rdkit' or 'mordred_2d' or 'mordred_3d'
# ベースのデータ
dataset = pd.read_csv('material_data.csv', index_col=0)
# 条件に応じて記述子ファイルを読み込む
if descriptor_type == 'rdkit':
des = pd.read_csv('descriptor_rdkit.csv', index_col=0)
elif descriptor_type == 'mordred_2d':
des = pd.read_csv('descriptor_mordred_2d.csv', index_col=0)
elif descriptor_type == 'mordred_3d':
des = pd.read_csv('descriptor_mordred_3d.csv', index_col=0)
else:
raise ValueError(f"未知のdescriptor_type: {descriptor_type}")
# 結合
dataset_full = pd.concat([dataset.reset_index(), des.reset_index(drop=True)], axis=1)
dataset_full = dataset_full.set_index('Material')
# PLに欠損値が入っている行を消して、学習データのみにする
dataset_train = dataset_full.dropna(subset='PL')
# SMILESとTypeも使わないので消しておく
dataset_train = dataset_train.drop(['SMILES', 'Type'], axis=1)
予測モデルの構築なので、学習データのみの読み込みにしています。テストデータ(未知の材料)はまだ使いません。
# infをNaNに置き換え
dataset_train = dataset_train.replace(np.inf, np.nan).fillna(np.nan)
dataset_train = dataset_train.drop(dataset_train.columns[dataset_train.isnull().any()], axis=1)
# 標準偏差が0の記述を削除
dataset_train = dataset_train.drop(dataset_train.columns[dataset_train.std() == 0], axis=1)
# 学習データのstdが0の列を特定
zero_std_cols = dataset_train.columns[dataset_train.std() == 0]
# 学習データから特定した列を削除
dataset_train = dataset_train.drop(columns=zero_std_cols)
print(dataset_train.shape)
(251, 1220)
前処理を行ったところ、列数が1220個となりました。それでは、説明変数と目的変数に分けた後、学習データを7:3に分割します。7割が学習データ、3割が汎化性能を確認するための検証データです。
# 目的変数と説明変数に分ける
y = dataset_train['PL']
X = dataset_train.drop('PL', axis=1)
# Hold-out
X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.3, random_state=1234)
# 標準化
X_scaler = StandardScaler()
autoscaled_X_train = X_scaler.fit_transform(X_train)
autoscaled_X_val = X_scaler.transform(X_val)
y_scaler = StandardScaler()
autoscaled_y_train = y_scaler.fit_transform(y_train.values.reshape(-1,1))
autoscaled_y_val = y_scaler.transform(y_val.values.reshape(-1,1))
# pandas形式に
autoscaled_X_train = pd.DataFrame(autoscaled_X_train, columns=X_train.columns, index=X_train.index)
autoscaled_X_val = pd.DataFrame(autoscaled_X_val, columns=X_val.columns, index=X_val.index)
autoscaled_y_train = pd.DataFrame(autoscaled_y_train, index=y_train.index, columns=["y"])
autoscaled_y_val = pd.DataFrame(autoscaled_y_val, index=y_val.index, columns=["y"])
いわゆるホールドアウトというやつです。ささっと汎化性能を知りたいので、まずはこれで行きましょう。後で最終的な交差検証するので大丈夫です。
モデルに入力する前に標準化を実行しています。決定木モデルではあまり効果ないですが、PLSやNNなどのモデルでは効果ありです(むしろ実行しないとうまくいかないことあり)。目的変数yに関しては、pd.Series(1次元)をnumpyの列ベクトル(2次元)に変換してから標準化を行っています。
そしてXもyも標準化後はnp.arrayのままになっているので、再度、pd.DataFrameに変換しています。後程の結果の可視化のためです。この可視化用の関数も定義しておきましょう。実測値と予測値をプロットする、y-y plotと呼ばれるものです。
# 可視化の関数
def yyplot(train_df, test_df):
# 指標の計算
rmse_tr = np.sqrt(mean_squared_error(train_df['true'], train_df['pred']))
r2_tr = r2_score(train_df['true'], train_df['pred'])
mae_tr = mean_absolute_error(train_df['true'], train_df['pred'])
rmse_te = np.sqrt(mean_squared_error(test_df['true'], test_df['pred']))
r2_te = r2_score(test_df['true'], test_df['pred'])
mae_te = mean_absolute_error(test_df['true'], test_df['pred'])
# グラフ作成
plt.figure(figsize=(7, 7))
ax = plt.subplot(111)
# 散布図
ax.scatter(train_df['true'], train_df['pred'], c='red', label='Train', alpha=0.6)
ax.scatter(test_df['true'], test_df['pred'], c='blue', label='Val', alpha=0.6)
# 理想線
all_true = pd.concat([train_df['true'], test_df['true']])
min_val = all_true.min()
max_val = all_true.max()
ax.plot([min_val, max_val], [min_val, max_val], 'k--')
# 軸ラベル
ax.set_xlabel('True Value')
ax.set_ylabel('Predicted Value')
# 範囲調整
margin = (max_val - min_val) * 0.05
ax.set_xlim(min_val - margin, max_val + margin)
ax.set_ylim(min_val - margin, max_val + margin)
# 評価値表示（右上）
ax.text(0.05, 0.95, f'Train RMSE = {rmse_tr:.2f}\nTrain MAE = {mae_tr:.2f}\nTrain R² = {r2_tr:.2f}',
transform=ax.transAxes, fontsize=11, color='red', verticalalignment='top')
ax.text(0.05, 0.80, f'Val RMSE = {rmse_te:.2f}\nVal MAE = {mae_te:.2f}\nVal R² = {r2_te:.2f}',
transform=ax.transAxes, fontsize=11, color='blue', verticalalignment='top')
# 凡例
ax.legend(loc='lower right')
plt.tight_layout()
plt.show()
コピペでOKです。「学習データの実測値と予測値」、「検証データの実測値と予測値」をそれぞれDataFrameにして、その各DataFrameをこの関数に入力すれば、予測結果が出てくる仕様にしています。
今回は評価指標として、MAE, RMSE, R^2(決定係数)を使います。これらの指標を学習データと検証データそれぞれで算出し、グラフと共に表示する関数になります。
前準備が長くなりましたが、それではモデル構築していきましょう！！
2. モデル構築 ~Hold Out編~
■ 部分的最小二乗回帰(Partial Least Squares Regression, PLS)
PCAと似ていますが、PLSは説明変数だけでなく目的変数も同時に考慮して成分を作り、そこから線形回帰を行います。モデルが軽くて解釈性もあり、まず初手で使うにはちょうどいい手法だと思います。
↓詳しい理論は金子先生のサイトにお任せします。
https://datachemeng.com/partialleastsquares/
# CVによる成分数の最適化
# 使用する主成分の最大数。説明変数の数より小さい必要があります
max_number_of_principal_components = 30
# N-fold CV の N
fold_number = 5
# 空のリストを用意
components = []
r2_in_cv_all = []
# reshapeしておく（Series → 1D ndarray）
y_train_1d = autoscaled_y_train.values.ravel()
for component in range(1, max_number_of_principal_components+1):
# PLS
model = PLSRegression(n_components=component)
# CV
estimated_y_cv = cross_val_predict(model, autoscaled_X_train, y_train_1d,
cv=fold_number)
estimated_y_cv = pd.DataFrame(estimated_y_cv)
# r2を計算
r2_in_cv = r2_score(y_train_1d, estimated_y_cv)
# 成分数とr2を表示
print(component, r2_in_cv)
# r2とcomponentを追加
components.append(component)
r2_in_cv_all.append(r2_in_cv)
# 成分数毎のr^2をプロットし、CV後のr^2が最大の時を、最適成分数として抽出
plt.scatter(components, r2_in_cv_all, c='blue')
plt.xlabel('number of components')
plt.ylabel('cross_validated r^2')
plt.show()
# 最適な主成分数を取得
optimal_component_number = components[r2_in_cv_all.index(max(r2_in_cv_all))]
print('最適な主成分数 :', optimal_component_number)
最適な主成分数 : 6
何をしているかというと、PLSはPCAをするので主成分の数を決めないといけません(主成分の数がハイパーパラメータ)。なので学習データで交差検証を用いて、決定係数が最も良い主成分を探しています。金子先生の書籍から引用したコードです。上のプロットから主成分が「6」の時が最も決定係数が良く、主成分が増えるにつれて決定係数が下がっていくのが分かります。
それでは主成分を6として、再度PLSモデルを構築し、検証データにおける予測精度を算出してみましょう。
# PLSモデルの構築
pls_model = PLSRegression(n_components=optimal_component_number)
# モデル学習
pls_model.fit(autoscaled_X_train, autoscaled_y_train)
# 予測
y_train_pred_scaled = pls_model.predict(autoscaled_X_train)
y_val_pred_scaled = pls_model.predict(autoscaled_X_val)
# 元のスケールに逆変換
y_train_pred = y_scaler.inverse_transform(y_train_pred_scaled)
y_val_pred = y_scaler.inverse_transform(y_val_pred_scaled)
# 可視化用データフレーム作成、yは2D配列から1Dベクトルに戻す
train_df = pd.DataFrame({
'true': y_train.reset_index(drop=True),
'pred': y_train_pred.flatten()
})
val_df = pd.DataFrame({
'true': y_val.reset_index(drop=True),
'pred': y_val_pred.flatten()
})
# 可視化
yyplot(train_df, val_df)
x軸が実測のPL波長、y軸が予測したPL波長になっており、もし予測がうまくいっていれば対角線上にプロットされます。赤のプロットが学習データ、青のプロットが検証データです。
MAE, RMSEは小さいほど予測がうまくいっており、R^2は大きいほど予測と実測のずれが小さいことになります(最大値は1)。意外とうまくいっていることが分かります。
考察は最後にまとめてしますので、ひとまず一つの結果として頭の片隅に置いておきましょう。次行きます。
■ RandomForest
決定木をたくさん作って、それらを並列に学習させる手法です（バギングといいます）。各木が出す予測を平均して最終出力とするため、個々の木は弱くても、全体として強力なモデルになります。精度も高く、特徴量の重要度も見えるので、こちらも初手でよく使われます。まさに「3人寄れば文殊の知恵」な手法です。
# デフォルト
model_rf = RandomForestRegressor(random_state=1234)
# モデル構築、予測
model_rf.fit(autoscaled_X_train, autoscaled_y_train)
y_train_pred_scaled = model_rf.predict(autoscaled_X_train)
y_val_pred_scaled = model_rf.predict(autoscaled_X_val)
# 元のスケールに逆変換（reshapeが必要）
y_train_pred = y_scaler.inverse_transform(y_train_pred_scaled.reshape(-1, 1))
y_val_pred = y_scaler.inverse_transform(y_val_pred_scaled.reshape(-1, 1))
# 可視化用データフレーム作成
train_df = pd.DataFrame({
'true': y_train.reset_index(drop=True),
'pred': y_train_pred.flatten()
})
val_df = pd.DataFrame({
'true': y_val.reset_index(drop=True),
'pred': y_val_pred.flatten()
})
# 可視化
yyplot(train_df, val_df)
先ほどのPLSモデルに比べると過学習となっており(学習データにフィッティングしすぎて、ノイズなども拾ってしまうこと)、検証データにおける指標が落ちています。悪くはないのですが、ハイパーパラメータの調整は必要です。
ハイパーパラメータのチューニングにはoptunaと呼ばれる日本のPrefferdNetworks社が開発した、ハイパーパラメータの自動最適化フレームワークを使用します。optunaはTPE（Tree-structured Parzen Estimator）というベイズ最適化の一種を採用しており、パラメータの探索を効率的に行うことができます。
それでは使ってみます。まず、どのハイパーパラメータを調整するか、何を基準にハイパーパラメータを調整するかの設定をobjective関数の中で定義します。
def objective_rf(trial):
# ハイパーパラメータのサンプリング
n_estimators = trial.suggest_int('n_estimators', 50, 300)
max_depth = trial.suggest_int('max_depth', 3, 20)
min_samples_split = trial.suggest_int('min_samples_split', 2, 10)
min_samples_leaf = trial.suggest_int('min_samples_leaf', 1, 10)
max_features = trial.suggest_categorical('max_features', ['sqrt', 'log2', None])
# モデル定義
model = RandomForestRegressor(
n_estimators=n_estimators,
max_depth=max_depth,
min_samples_split=min_samples_split,
min_samples_leaf=min_samples_leaf,
max_features=max_features,
random_state=0,
n_jobs=-1
)
# モデル学習
model.fit(autoscaled_X_train, autoscaled_y_train)
# 検証データでRMSEを評価
y_pred_scaled = model.predict(autoscaled_X_val)
y_pred = y_scaler.inverse_transform(y_pred_scaled.reshape(-1, 1))
rmse = np.sqrt(mean_squared_error(y_val, y_pred))
return rmse
今回はRMSEを評価指標にしました。この指標を下げるように学習させていきます。
# 学習
study = optuna.create_study(direction='minimize')
study.optimize(objective_rf, n_trials=50)
RMSEは下げたい指標なので、direction='minimize'にして下げるように明示します。R^2の場合は、上げたい指標なので、direction='maximize'に変更しましょう。n_trialsを変更することで、何回試行を繰り返すかを設定することができ、今回は50回試行しています。
学習が終わると、最も良かった時のハイパーパラメータを見ることができます。
# ハイパーパラメータ・スコアの確認
print("Best trial:")
trial = study.best_trial
print(f" RMSE: {trial.value:.4f}")
print(" Params:")
for key, value in trial.params.items():
print(f" {key}: {value}")
Best trial:
RMSE: 0.4616
Params:
n_estimators: 251
max_depth: 16
min_samples_split: 5
min_samples_leaf: 1
max_features: sqrt
上記のハイパーパラメータが良いそうです。それでは、この最適化されたハイパーパラメータを選択して、再度RandomForestモデルを構築させて、検証データを予測してみます。
# 最良ハイパーパラメータで再学習
best_params = study.best_params
model_rf_optuna = RandomForestRegressor(**best_params, random_state=0, n_jobs=-1)
model_rf_optuna.fit(autoscaled_X_train, autoscaled_y_train)
# 予測
y_train_pred_scaled = model_rf_optuna.predict(autoscaled_X_train)
y_val_pred_scaled = model_rf_optuna.predict(autoscaled_X_val)
# 元のスケールに逆変換（reshapeが必要）
y_train_pred = y_scaler.inverse_transform(y_train_pred_scaled.reshape(-1, 1))
y_val_pred = y_scaler.inverse_transform(y_val_pred_scaled.reshape(-1, 1))
# 可視化用データフレーム作成
train_df = pd.DataFrame({
'true': y_train.reset_index(drop=True),
'pred': y_train_pred.flatten()
})
val_df = pd.DataFrame({
'true': y_val.reset_index(drop=True),
'pred': y_val_pred.flatten()
})
# 可視化
yyplot(train_df, val_df)
わずかではありますが検証データにおけるRMSEが下がり、R^2は上がっているので、ハイパーパラメータチューニングの効果はあったみたいです。RandomForestなどの決定木モデルは特徴量重要度なるものが見れます。どの説明変数が予測に役立ったかを教えてくれます。見てみましょう。
# 特徴量重要度を出す
feature_map_rf = pd.DataFrame([X_train.columns, model_rf_optuna.feature_importances_]).T
feature_map_rf.columns = ['feature', 'importances']
# 降順に
feature_map_rf = feature_map_rf.sort_values('importances', ascending=False)
# 規格化
feature_map_rf['importances'] = feature_map_rf['importances'] / feature_map_rf['importances'].max()
# 可視化
sns.barplot(y=feature_map_rf['feature'][:10],
x=feature_map_rf['importances'][:10],
palette='viridis')
plt.xlabel('Normalized Feature Importance')
plt.ylabel('Feature')
plt.title('Top 10 Important Features (RandomForest)')
plt.show()
nHetero(ヘテロ原子の数)が役に立ったみたいです。窒素や酸素といったヘテロ原子は孤立電子対を持つので、これがπ共役系に関与することで、電子の非局在化を強める場合があります。また、ヘテロ原子は炭素原子と電気陰性度が異なるので、電子を与える性質があるものはHOMOを押し上げ、電子を引っ張る性質があるものはLUMOを下げます。結果として、波長に影響を与えるHOMOとLUMOのGapが変化し、波長が変化するのです。
# 予測に役立った説明変数と目的関数の関係をplot
# 最も重要な特徴量の名前を取り出す（文字列）
top_feature = feature_map_rf['feature'].iloc[0]
sns.scatterplot(x=top_feature,
y='PL',
data=dataset_train)
plt.show()
相関があるかどうかは可視化するとわかりやすいです。ヘテロ原子の数が増えるとPL波長が長波になっていることが分かります。
ひとまずOptunaによる改善が見られたので、こんな感じで次のモデルも見ていきましょう。
■ XGBoost
こちらも決定木ベースのアンサンブルモデルですが、RandomForestが並列に処理するのに対し、XGBoostは1本ずつ直列に木を構築していきます（ブースティングといいます）。前のモデルがうまく予測できなかった部分を補うように次の木を学習するため、精度が出やすい一方で過学習には少し注意が必要です。
# デフォルトXGB
model_xgb = xgb.XGBRegressor(random_state=1234)
# テストデータの組
eval_set = [(autoscaled_X_val, autoscaled_y_val)]
# モデル学習
model_xgb.fit(autoscaled_X_train, autoscaled_y_train, eval_set=eval_set)
# 予測
y_train_pred_scaled = model_xgb.predict(autoscaled_X_train)
y_val_pred_scaled = model_xgb.predict(autoscaled_X_val)
# 元のスケールに逆変換（reshapeが必要）
y_train_pred = y_scaler.inverse_transform(y_train_pred_scaled.reshape(-1, 1))
y_val_pred = y_scaler.inverse_transform(y_val_pred_scaled.reshape(-1, 1))
# 可視化用データフレーム作成
train_df = pd.DataFrame({
'true': y_train.reset_index(drop=True),
'pred': y_train_pred.flatten()
})
val_df = pd.DataFrame({
'true': y_val.reset_index(drop=True),
'pred': y_val_pred.flatten()
})
# 可視化
yyplot(train_df, val_df)
学習データに関してはかなりオーバーフィッティングしていることが分かります。こちらもoptunaでチューニングしましょう。
def objective_xgb(trial):
# ハイパーパラメータのサンプリング
param = {
'n_estimators': trial.suggest_int('n_estimators', 50, 300),
'max_depth': trial.suggest_int('max_depth', 3, 10),
'learning_rate': trial.suggest_float('learning_rate', 1e-3, 0.3, log=True),
'subsample': trial.suggest_float('subsample', 0.5, 1.0),
'colsample_bytree': trial.suggest_float('colsample_bytree', 0.5, 1.0),
'reg_alpha': trial.suggest_float('reg_alpha', 1e-8, 10.0, log=True),
'reg_lambda': trial.suggest_float('reg_lambda', 1e-8, 10.0, log=True),
'random_state': 42,
'tree_method': 'hist', # CPUで高速にする場合
}
model = xgb.XGBRegressor(**param)
# 学習と予測
model.fit(autoscaled_X_train, autoscaled_y_train, eval_set = [(autoscaled_X_val, autoscaled_y_val)], verbose=0)
y_pred_scaled = model.predict(autoscaled_X_val)
y_pred = y_scaler.inverse_transform(y_pred_scaled.reshape(-1, 1))
# 検証データでRMSEを評価
rmse = np.sqrt(mean_squared_error(y_val, y_pred))
return rmse
# 学習
study = optuna.create_study(direction='minimize')
study.optimize(objective_xgb, n_trials=50, timeout=600) # 試行回数や時間は任意
こちらもRMSEを評価指標にして、50回試行して、最適なハイパーパラメータを選定します。
# ハイパーパラメータ・スコアの確認
print("Best trial:")
trial = study.best_trial
print(f" RMSE: {trial.value:.4f}")
print(" Params:")
for key, value in trial.params.items():
print(f" {key}: {value}")
Best trial:
RMSE: 0.4311
Params:
n_estimators: 170
max_depth: 9
learning_rate: 0.08883328502362148
subsample: 0.5964959689234826
colsample_bytree: 0.8465192376330929
reg_alpha: 0.0009159791865774875
reg_lambda: 0.6832941280882167
最適なハイパーパラメータが得られたので、これらをモデルにセットして再度評価します。
# 再学習
best_params = trial.params
model_xgb_optuna = xgb.XGBRegressor(**best_params)
model_xgb_optuna.fit(autoscaled_X_train, autoscaled_y_train, eval_set=eval_set, verbose=0)
# 予測
y_train_pred_scaled = model_xgb_optuna.predict(autoscaled_X_train)
y_val_pred_scaled = model_xgb_optuna.predict(autoscaled_X_val)
# 元のスケールに逆変換（reshapeが必要）
y_train_pred = y_scaler.inverse_transform(y_train_pred_scaled.reshape(-1, 1))
y_val_pred = y_scaler.inverse_transform(y_val_pred_scaled.reshape(-1, 1))
# 可視化用データフレーム作成
train_df = pd.DataFrame({
'true': y_train.reset_index(drop=True),
'pred': y_train_pred.flatten()
})
val_df = pd.DataFrame({
'true': y_val.reset_index(drop=True),
'pred': y_val_pred.flatten()
})
# 可視化
yyplot(train_df, val_df)
相変わらずオーバーフィッティング気味ですが、検証データにおける評価指標は改善しています。いい感じです。XGBoostも特徴量重要度が見れます。
# 特徴量重要度を取得
feature_map_xgb = pd.DataFrame([X_train.columns, model_xgb_optuna.feature_importances_]).T
feature_map_xgb.columns = ['feature', 'importances']
# 降順に並べ替え
feature_map_xgb = feature_map_xgb.sort_values('importances', ascending=False)
# 規格化（最大値で割る）
feature_map_xgb['importances'] = feature_map_xgb['importances'] / feature_map_xgb['importances'].max()
# 上位10個の特徴量を可視化
plt.figure(figsize=(8, 6))
sns.barplot(y=feature_map_xgb['feature'][:10],
x=feature_map_xgb['importances'][:10],
palette='viridis')
plt.xlabel('Normalized Feature Importance')
plt.ylabel('Feature')
plt.title('Top 10 Important Features (XGBoost)')
plt.tight_layout()
plt.show()
こちらもヘテロ原子の数が予測に役立ったようで、モデルの考えは似ているようです。次行きましょう。
■ LightGBM
XGBoostと同様、決定木ベースのブースティング手法です。XGBoostよりも軽量で計算が速いです。特徴的なのは、「木の育て方」に工夫があるところです。普通の決定木やXGBoostでは“深さ方向”に木を伸ばす（レベル単位で分岐）のに対して、LightGBMは葉ごとに分岐を進める「Leaf-wise成長」を採用しています。
それでは、実装してみます。コードはXGBoostと大きく変わりません。
# デフォルト
model_lgb = lgb.LGBMRegressor(random_state=1234)
eval_set = [(autoscaled_X_val, autoscaled_y_val)]
# モデル学習
model_lgb.fit(autoscaled_X_train, autoscaled_y_train, eval_set=eval_set)
# 予測
y_train_pred_scaled = model_lgb.predict(autoscaled_X_train)
y_val_pred_scaled = model_lgb.predict(autoscaled_X_val)
# 元のスケールに逆変換（reshapeが必要）
y_train_pred = y_scaler.inverse_transform(y_train_pred_scaled.reshape(-1, 1))
y_val_pred = y_scaler.inverse_transform(y_val_pred_scaled.reshape(-1, 1))
# 可視化用データフレーム作成
train_df = pd.DataFrame({
'true': y_train.reset_index(drop=True),
'pred': y_train_pred.flatten()
})
val_df = pd.DataFrame({
'true': y_val.reset_index(drop=True),
'pred': y_val_pred.flatten()
})
# 可視化
yyplot(train_df, val_df)
初手でRandomForestよりかは良い精度を示しました。ハイパーパラメータチューニングもやってみます。
def objective_lgb(trial):
param = {
'objective': 'regression',
'metric': 'rmse', # LightGBMが内部で使う評価指標
'n_estimators': trial.suggest_int('n_estimators', 50, 300),
'learning_rate': trial.suggest_float('learning_rate', 0.01, 0.3, log=True),
'max_depth': trial.suggest_int('max_depth', 3, 10),
'num_leaves': trial.suggest_int('num_leaves', 10, 100),
'subsample': trial.suggest_float('subsample', 0.5, 1.0),
'colsample_bytree': trial.suggest_float('colsample_bytree', 0.5, 1.0),
'reg_alpha': trial.suggest_float('reg_alpha', 1e-8, 1.0, log=True),
'reg_lambda': trial.suggest_float('reg_lambda', 1e-8, 1.0, log=True),
'min_data_in_leaf' : trial.suggest_int('min_data_in_leaf', 1, 10),
'min_child_samples' :trial.suggest_int('min_child_samples', 1, 10),
'random_state': 42,
'verbosity': -1,
}
model = lgb.LGBMRegressor(**param)
# 学習&予測
model.fit(autoscaled_X_train, autoscaled_y_train, eval_set=[(autoscaled_X_val, autoscaled_y_val)])
y_pred_scaled = model.predict(autoscaled_X_val)
y_pred = y_scaler.inverse_transform(y_pred_scaled.reshape(-1, 1))
# 汎化性能算出
rmse = np.sqrt(mean_squared_error(y_val, y_pred))
return rmse
# 最適化
study = optuna.create_study(direction='minimize')
study.optimize(objective_lgb, n_trials=50, timeout=600) # 50試行 or 600秒
# 結果表示
print("Best trial:")
print(f" RMSE: {study.best_trial.value:.4f}")
for key, value in study.best_trial.params.items():
print(f" {key}: {value}")
Best trial:
RMSE: 47.6359
n_estimators: 281
learning_rate: 0.08481241820280185
max_depth: 3
num_leaves: 55
subsample: 0.7051369888415676
colsample_bytree: 0.6189382024481481
reg_alpha: 6.386476273078106e-06
reg_lambda: 3.7597161363457143e-08
min_data_in_leaf: 6
min_child_samples: 9
最適化が終わりました。私のPCでは、
XGBoost : 3m50s
LightGBM : 26s
という実行時間となり、LightGBMが軽量モデルであることが改めて分かります。それでは、最後に最適化されたモデルの性能を見てみます。
# モデル再構築
best_params = study.best_trial.params
model_lgb_optuna = lgb.LGBMRegressor(**best_params)
model_lgb_optuna.fit(autoscaled_X_train, autoscaled_y_train)
# 予測
y_train_pred_scaled = model_lgb_optuna.predict(autoscaled_X_train)
y_val_pred_scaled = model_lgb_optuna.predict(autoscaled_X_val)
# 元のスケールに逆変換（reshapeが必要）
y_train_pred = y_scaler.inverse_transform(y_train_pred_scaled.reshape(-1, 1))
y_val_pred = y_scaler.inverse_transform(y_val_pred_scaled.reshape(-1, 1))
# 可視化用データフレーム作成
train_df = pd.DataFrame({
'true': y_train.reset_index(drop=True),
'pred': y_train_pred.flatten()
})
val_df = pd.DataFrame({
'true': y_val.reset_index(drop=True),
'pred': y_val_pred.flatten()
})
# 可視化
yyplot(train_df, val_df)
こちらは若干ですが検証データにおけるMAEが下がっています。もともとデフォルトの段階で最適化されている状態なのかもしれません。LightGBMももちろん特徴量重要度が可視化できます。
# 特徴量重要度の取得
feature_map_lgb = pd.DataFrame({
'feature': X_train.columns,
'importances': model_lgb_optuna.feature_importances_
})
# 降順に並べ替え
feature_map_lgb = feature_map_lgb.sort_values('importances', ascending=False)
# 規格化（最大値で割る）
feature_map_lgb['importances'] = feature_map_lgb['importances'] / feature_map_lgb['importances'].max()
# 上位10個の特徴量を可視化
plt.figure(figsize=(8, 6))
sns.barplot(y=feature_map_lgb['feature'][:10],
x=feature_map_lgb['importances'][:10],
palette='viridis')
plt.xlabel('Normalized Feature Importance')
plt.ylabel('Feature')
plt.title('Top 10 Important Features (LightGBM)')
plt.tight_layout()
plt.show()
RFやXGBとは異なる結果になりました。モデルの思考が変わっています。プロットしてみましょう。
# 予測に役立った説明変数と目的関数の関係をplot
# 最も重要な特徴量の名前を取り出す（文字列）
top_feature = feature_map_lgb['feature'].iloc[0]
sns.scatterplot(x=top_feature,
y='PL',
data=dataset_train)
plt.show()
SaaNはSum of solvent-accessible surface area on Nitrogen atomsのことで、窒素原子に対する溶媒接触面積を合計した値だそうです。分子内に露出している窒素原子が多く、その表面積が大きい場合、SaaNは大きくなります。はっきりと線形の関係にあるわけではないですが、LGBがこの特徴量を重視したみたいです。
長々と記載しましたが、ひとまず決定木モデル系は終わりです。一旦、ブレイクしましょう☕
■ Neural Network(NN)
ニューラルネットワーク（NN）は、人間の脳の神経回路を模したモデルで、入力→隠れ層→出力という構造でデータを処理します。入力された情報は、複数の「ニューロン（ノード）」を通って、非線形な変換を繰り返しながら出力に変換されます。
非線形な関係性や複雑な特徴の組み合わせを学習できるのが最大の強みで、線形モデルは表現しきれないような高度なパターンも扱えます。その反面、データ量が少ないと過学習しやすく、パラメータ調整や学習の安定化が難しいという一面もあります。
最近はPyTorchやTensorFlowといったフレームワークのおかげで実装の敷居も下がり、カスタマイズや自動最適化もかなり簡単になってきています。
それでは、実装行きます。私はPyTorchユーザーなのでこちらでコーディングさせていただきます。GPUがある場合はまず下記コードでGPUを認識しているか確認しましょう。
# GPU set
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
print(device)
cudaと返ってきたらGPUが使える証拠です。ちなみにCPUでも問題なく実装はできるので、CPUの場合は今後出てくるコードのdeviceという部分を消して実行してください。
まずNNの構成を定義します。
# NN定義
class MLPRegressor(nn.Module):
def __init__(self, input_dim):
super(MLPRegressor, self).__init__()
self.model = nn.Sequential(
nn.Linear(input_dim, 64),
nn.BatchNorm1d(64),
nn.ReLU(),
nn.Dropout(0.2),
nn.Linear(64, 32),
nn.BatchNorm1d(32),
nn.ReLU(),
nn.Dropout(0.2),
nn.Linear(32, 16),
nn.BatchNorm1d(16),
nn.ReLU(),
nn.Linear(16, 1)
)
def forward(self, x):
return self.model(x)
# モデル構築
model_nn = MLPRegressor(input_dim=X.shape[1]).to(device)
# 損失関数
criterion = nn.MSELoss()
# 最適化関数
optimizer = torch.optim.Adam(model_nn.parameters(), lr=1e-4)
細かい説明は省きますが、中間層が3つあり、ノード数は64→32→16となっている構成です。回帰問題(PL波長を当てる)なので、最後の出力層で1ノードの出力としています。
活性化関数を通すことで、非線形な関係を表現できるようになるのもNNの強みの一つです。今回はその関数にはReLU関数というものを使っています。
そして過学習防止のために、ドロップアウトやバッチ正規化も使用しています。
また、損失関数はMSE、最適化関数をAdamで定義しています。損失(誤差)を算出し、逆伝播して勾配計算を行い、得られた勾配をもとにパラメータを更新していくのがNNの基本の流れです。
次にDataset、DataLoaderを定義します。
# データセット定義
class RegressionDataset(Dataset):
def __init__(self, X, y):
self.X = torch.tensor(X, dtype=torch.float32)
self.y = torch.tensor(y, dtype=torch.float32).view(-1, 1)
def __len__(self):
return len(self.X)
def __getitem__(self, idx):
return self.X[idx], self.y[idx]
# Datasetに渡す（.valuesでnumpy化）
train_ds = RegressionDataset(autoscaled_X_train.values, autoscaled_y_train.values)
val_ds = RegressionDataset(autoscaled_X_val.values, autoscaled_y_val.values)
# DataLoaderに渡す
train_dl = DataLoader(train_ds, batch_size=16, shuffle=True)
val_dl = DataLoader(val_ds, batch_size=16)
DatasetはPyTorchにデータをどう渡すかを定義するクラスです。説明変数と目的変数をテンソルにしています。テンソルは3次元以上の数の配列です。あまり深く考えなくてよく、高速処理するためにはテンソルにしないといけない程度に覚えておけばいいです。
DataLoaderはDatasetからデータを取り出して、モデルに渡す「準備」をします。バッチ単位でデータを取り出すことができ(今回はバッチサイズは16)、シャッフルしてデータの偏りを防いだり、並列処理で高速化ができます。
最後に学習と評価の関数も定義しておきます。あらかじめ定義しておくと便利なので、先に関数化しています。
# 学習loop
def train(model, dataloader, optimizer, criterion, device):
model.train()
total_loss = 0
for X_batch, y_batch in dataloader:
X_batch, y_batch = X_batch.to(device), y_batch.to(device)
optimizer.zero_grad()
pred = model(X_batch)
loss = criterion(pred, y_batch)
loss.backward()
optimizer.step()
total_loss += loss.item()
return total_loss / len(dataloader)
# 評価
def evaluate(model, dataloader, criterion, device):
model.eval()
total_loss = 0.0
with torch.no_grad():
for X_batch, y_batch in dataloader:
X_batch, y_batch = X_batch.to(device), y_batch.to(device)
pred = model(X_batch)
loss = criterion(pred, y_batch)
total_loss += loss.item()
return total_loss / len(dataloader)
学習loopでは、先ほど定義したDataLoaderからデータをバッチとして取り出し、モデル予測→損失計算→誤差逆伝播(勾配計算)→パラメータ更新を行っています。そしてバッチ毎にlossを算出し、最後に平均化しています。評価でもほとんど同じ処理を行っていますが、一つ違うのは、パラメータの更新は行っていません。
長くなりましたが、それではNNの学習・評価を行います。試行回数(epochs)は300回でやってみます。
# 学習実行
n_epochs = 300
train_losses = []
val_losses = []
for epoch in range(n_epochs):
train_loss = train(model_nn, train_dl, optimizer, criterion, device)
val_loss = evaluate(model_nn, val_dl, criterion, device)
train_losses.append(train_loss)
val_losses.append(val_loss)
print(f"Epoch {epoch:3d} | Train Loss: {train_loss:.4f} | Val Loss: {val_loss:.4f}")
# 学習推移の可視化
plt.plot(train_losses, label='Train Loss')
plt.plot(val_losses, label='Val Loss')
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.legend()
plt.grid()
plt.show()
学習が問題なく進むと、lossが着々と減っていることが分かります。100epochsくらいでval_lossは下がりきっていることも分かりますね。
今までのモデルと同様に、yyplotを描写してみます。先のlossを返すevaluate関数とは別に、実際の予測値を返すpredict関数も新たに定義します。
# 推論関数
def predict(model, dataloader, device):
model.eval()
preds = []
with torch.no_grad():
for X_batch, _ in dataloader:
X_batch = X_batch.to(device)
y_pred = model(X_batch).cpu().numpy()
preds.append(y_pred)
return np.vstack(preds)
# 推論用（shuffleなし）DataLoaderを用意
train_dl_eval = DataLoader(train_ds, batch_size=32, shuffle=False)
val_dl_eval = DataLoader(val_ds, batch_size=32, shuffle=False)
# 推論
y_train_pred_scaled = predict(model_nn, train_dl_eval, device)
y_val_pred_scaled = predict(model_nn, val_dl_eval, device)
# 逆変換
y_train_pred = y_scaler.inverse_transform(y_train_pred_scaled)
y_val_pred = y_scaler.inverse_transform(y_val_pred_scaled)
# 可視化用データフレーム作成
train_df = pd.DataFrame({
'true': y_train.reset_index(drop=True),
'pred': y_train_pred.flatten()
})
val_df = pd.DataFrame({
'true': y_val.reset_index(drop=True),
'pred': y_val_pred.flatten()
})
# 可視化
yyplot(train_df, val_df)
改めて推論用のDataLoaderを定義して、モデルに入れて予測値を算出させています。そこそこの精度ですが、もう少し改善できそうです。
推論用の train_dl を shuffle=True にしてしまうとy_trainとpredの順序が一致しなくなって結果がズレます。shuffle=Falseで推論させてください。元のデータであるy_trainと予測値であるpredの順序を揃える必要があるためです。学習は既に終わっているので、shuffleする必要はありません。もちろん学習の際はshuffleして、偏った学習をするのを防ぎましょう。
最後にNNとoptunaを組み合わせます。少しコードが複雑ですが、ゆっくり解説していきます。まず条件を振るためのNNを定義します。
def define_model(trial, input_dim):
layers = []
n_layers = trial.suggest_int("n_layers", 2, 5)
hidden_dim = trial.suggest_int("hidden_dim", 64, 1024)
dropout_rate = trial.suggest_float("dropout_rate", 0.0, 0.5)
activation_name = trial.suggest_categorical("activation", ["relu", "leaky_relu"])
use_batchnorm = trial.suggest_categorical("use_batchnorm", [True, False])
in_dim = input_dim
for i in range(n_layers):
layers.append(nn.Linear(in_dim, hidden_dim))
# 活性化関数の前に入れる
if use_batchnorm:
layers.append(nn.BatchNorm1d(hidden_dim))
if activation_name == "relu":
layers.append(nn.ReLU())
else:
layers.append(nn.LeakyReLU())
layers.append(nn.Dropout(dropout_rate))
in_dim = hidden_dim
layers.append(nn.Linear(in_dim, 1))
return nn.Sequential(*layers)
layersという空リストにNNを構成する層を入れていきます。n_layersがNNの層の数、hidden_dimは隠れ層のノード数、dropout_rateはドロップアウト率(どれくらいノードを消去するか)、activation_nameは活性化関数にreluかleaky_reluのどちらを使うか、use_batchnormはバッチ正規化を使うかです。これらの要素をチューニングします。
forループの中身で入力層から中間層を設定し、最後に出力層(nn.Linear)を設定して1ノードの出力としています。
それでは、この条件振りNNをobjective関数に入れます。
def objective_nn(trial):
# データ準備、バッチサイズも振る
batch_size = trial.suggest_categorical("batch_size", [32, 64, 128])
train_dl = DataLoader(train_ds, batch_size=batch_size, shuffle=True)
val_dl = DataLoader(val_ds, batch_size=batch_size)
# 先ほど定義した条件振りNN
model = define_model(trial, input_dim=X.shape[1]).to(device)
# 学習率も振る
lr = trial.suggest_float("lr", 1e-5, 1e-2, log=True)
# 最適化関数も振る
optimizer_name = trial.suggest_categorical("optimizer", ["Adam", "SGD"])
if optimizer_name == "Adam":
optimizer = torch.optim.Adam(model.parameters(), lr=lr)
else:
optimizer = torch.optim.SGD(model.parameters(), lr=lr)
criterion = nn.MSELoss()
# 学習ループ
for epoch in range(30):
train(model, train_dl, optimizer, criterion, device)
# 検証
val_loss = evaluate(model, val_dl, criterion, device)
return val_loss
# 学習
study = optuna.create_study(direction="minimize")
study.optimize(objective_nn, n_trials=50)
# 結果表示
print("Best trial:")
print(f" RMSE: {study.best_trial.value:.4f}")
for key, value in study.best_trial.params.items():
print(f" {key}: {value}")
Best trial:
RMSE: 0.1685
batch_size: 128
n_layers: 2
hidden_dim: 450
dropout_rate: 0.43519210062923575
activation: leaky_relu
use_batchnorm: True
lr: 0.0007018617917129689
optimizer: Adam
検討した結果がでました。これらのチューニングされたハイパーパラメータを使ってNNを再定義し、学習させます。
# モデル再定義
class TunedMLP(nn.Module):
def __init__(self, input_dim):
super(TunedMLP, self).__init__()
layers = []
in_dim = input_dim
hidden_dim = study.best_params['hidden_dim']
n_layers = study.best_params['n_layers']
dropout_rate = study.best_params['dropout_rate']
activation = nn.LeakyReLU() if study.best_params['activation'] == "leaky_relu" else nn.ReLU()
use_batchnorm = study.best_params.get('use_batchnorm', False)
for _ in range(n_layers):
layers.append(nn.Linear(in_dim, hidden_dim))
if use_batchnorm:
layers.append(nn.BatchNorm1d(hidden_dim))
layers.append(activation)
layers.append(nn.Dropout(dropout_rate))
in_dim = hidden_dim
layers.append(nn.Linear(in_dim, 1))
self.model = nn.Sequential(*layers)
def forward(self, x):
return self.model(x)
# DataLoader（batch_sizeも反映）
train_dl = DataLoader(train_ds, batch_size=study.best_params['batch_size'], shuffle=True)
val_dl = DataLoader(val_ds, batch_size=study.best_params['batch_size'])
# モデル・損失関数・最適化手法
model_nn = TunedMLP(input_dim=X.shape[1]).to(device)
criterion = nn.MSELoss()
optimizer = torch.optim.Adam(model_nn.parameters(), lr=study.best_params['lr'])
# 学習
train_losses, val_losses = [], []
for epoch in range(1, 100):
train_loss = train(model_nn, train_dl, optimizer, criterion, device)
val_loss = evaluate(model_nn, val_dl, criterion, device)
train_losses.append(train_loss)
val_losses.append(val_loss)
print(f"Epoch {epoch:03d} | Train Loss: {train_loss:.4f} | Val Loss: {val_loss:.4f}")
100epochs学習させて、学習データと検証データにおける損失を出力させています。最後にy-y plotを描写しましょう。
# 推論用DataLoader（shuffle=Falseで固定）
train_dl = DataLoader(train_ds, batch_size=study.best_params['batch_size'], shuffle=False)
val_dl = DataLoader(val_ds, batch_size=study.best_params['batch_size'], shuffle=False)
# 推論 & 逆変換
y_train_pred = y_scaler.inverse_transform(predict(model_nn, train_dl_eval, device))
y_val_pred = y_scaler.inverse_transform(predict(model_nn, val_dl_eval, device))
# 可視化データ作成
train_df = pd.DataFrame({'true': y_train.reset_index(drop=True), 'pred': y_train_pred.flatten()})
val_df = pd.DataFrame({'true': y_val.reset_index(drop=True), 'pred': y_val_pred.flatten()})
# 可視化
yyplot(train_df, val_df)
精度向上が見られました。チューニングする箇所が多いので、どうしてもコードが煩雑になってしまいます。NNの場合はoptunaを使わずに、最初に定義したMLPRegressorのノード数やドロップアウト率を変更するだけでも精度は動くので、最初はこのような調整方法でも良いかもしれません。
3. モデル構築 ~交差検証編~
最後に交差検証して、モデルの最終確認します。ホールドアウトで十分やんと思いたいのですが、データセットが小さい場合はデータに偏りが出てしまい、予測精度が大きく変化してしまいます。
ハイパーパラメータチューニング(ハイチューと略したい)も一部のデータのみに適合してしまう恐れもあります。ですので、平均的な精度も見ておいた方がいいです。
今回はRF, LGB, NNの三つで交差検証やってみます。
データセットの読み込み・前処理まではホールドアウトの時と同じで、交差検証の場合は全部の説明変数と目的変数を用意しておきます。用意した後はそのままにしておきます。
# ↓ここが違う！！
# 目的変数と説明変数に分ける
y = dataset_train['PL']
X = dataset_train.drop('PL', axis=1)
■ RamdomForest
それでは、まずRFで交差検証します。下のコードを実行してみます。注意点としては、標準化は各CV毎に行う必要があります。必ず各CVにおける学習データを基準に平均と標準偏差を求めて、これらの値を検証データに適用しましょう。
# デフォルト
model_rf = RandomForestRegressor(random_state=1234)
# 5分割交差検証
kf = KFold(n_splits=5, shuffle=True, random_state=1234)
# スコア保存用
rmse_scores = []
mae_scores = []
r2_scores = []
for train_index, val_index in kf.split(X):
# 訓練と検証に分類
X_train, X_val = X.iloc[train_index, :], X.iloc[val_index, :]
y_train, y_val = y.iloc[train_index], y.iloc[val_index]
# 標準化、DataFrameに戻す
X_scaler = StandardScaler()
autoscaled_X_train = pd.DataFrame(X_scaler.fit_transform(X_train), columns=X_train.columns, index=X_train.index)
autoscaled_X_val = pd.DataFrame(X_scaler.transform(X_val), columns=X_val.columns, index=X_val.index)
y_scaler = StandardScaler()
autoscaled_y_train = pd.DataFrame(y_scaler.fit_transform(y_train.values.reshape(-1,1)), index=y_train.index, columns=['y'])
autoscaled_y_val = pd.DataFrame(y_scaler.transform(y_val.values.reshape(-1,1)), index=y_val.index, columns=['y'])
# 学習
model_rf.fit(autoscaled_X_train, autoscaled_y_train.values.ravel())
y_pred_scaled = model_rf.predict(autoscaled_X_val)
y_pred = y_scaler.inverse_transform(y_pred_scaled.reshape(-1,1))
# 性能チェック
rmse = np.sqrt(mean_squared_error(y_val, y_pred))
mae = mean_absolute_error(y_val, y_pred)
r2 = r2_score(y_val, y_pred)
rmse_scores.append(rmse)
mae_scores.append(mae)
r2_scores.append(r2)
print(f'Fold RMSE : {rmse:.4f}')
print(f'Fold MAE : {mae:.4f}')
print(f'Fold R2 : {r2:.4f}')
print()
print(f'平均RMSE : {np.mean(rmse_scores):.4f}')
print(f'平均MAE : {np.mean(mae_scores):.4f}')
print(f'平均R2 : {np.mean(r2_scores):.4f}')
Fold RMSE : 51.2594
Fold MAE : 36.0020
Fold R2 : 0.7491
Fold RMSE : 44.5128
Fold MAE : 31.6281
Fold R2 : 0.8050
Fold RMSE : 48.6313
Fold MAE : 31.3504
Fold R2 : 0.8328
Fold RMSE : 47.4534
Fold MAE : 32.8682
Fold R2 : 0.8079
Fold RMSE : 53.2576
Fold MAE : 32.8291
Fold R2 : 0.7229
平均RMSE : 49.0229
平均MAE : 32.9355
平均R2 : 0.7835
こんな感じで、平均的な出力が出てきます。ホールドアウトの時の比べて評価指標が変化していますね。やはり部分的なホールドアウトはリスキーなことが改めて分かります。次はOptunaでハイチューをしましょう。
def objective_rf(trial):
# 検証するパラメータ
params = {
"n_estimators" : trial.suggest_int('n_estimators', 5, 1000),
"max_depth" : trial.suggest_int('max_depth', 3, 50),
"min_samples_split" : trial.suggest_int('min_samples_split', 2, 10),
"min_samples_leaf" : trial.suggest_int('min_samples_leaf', 1, 10),
"max_features" : trial.suggest_categorical('max_features', ['sqrt', 'log2', None]),
'random_state' : 0,
'n_jobs' : -1,
}
# モデル定義
model_rf = RandomForestRegressor(**params)
# 分割、CV
kf = KFold(n_splits=5, shuffle=True, random_state=1234)
# 評価指標
rmse_scores = []
# mae_scores = []
# r2_scores = []
for train_index, val_index in kf.split(X):
# 訓練と検証に分類
X_train, X_val = X.iloc[train_index, :], X.iloc[val_index, :]
y_train, y_val = y.iloc[train_index], y.iloc[val_index]
# 標準化、DataFrameに戻す
X_scaler = StandardScaler()
autoscaled_X_train = pd.DataFrame(X_scaler.fit_transform(X_train), columns=X_train.columns, index=X_train.index)
autoscaled_X_val = pd.DataFrame(X_scaler.transform(X_val), columns=X_val.columns, index=X_val.index)
y_scaler = StandardScaler()
autoscaled_y_train = pd.DataFrame(y_scaler.fit_transform(y_train.values.reshape(-1,1)), index=y_train.index, columns=['y'])
autoscaled_y_val = pd.DataFrame(y_scaler.transform(y_val.values.reshape(-1,1)), index=y_val.index, columns=['y'])
# 学習&予測
model_rf.fit(autoscaled_X_train, autoscaled_y_train)
y_pred_scaled = model_rf.predict(autoscaled_X_val)
y_pred = y_scaler.inverse_transform(y_pred_scaled.reshape(-1,1))
# 性能チェック
rmse = np.sqrt(mean_squared_error(y_val, y_pred))
# mae = mean_absolute_error(y_val, y_pred)
# r2 = r2_score(y_val, y_pred)
# 結果格納
rmse_scores.append(rmse)
# mae_scores.append(mae)
# r2_scores.append(r2)
return np.mean(rmse_scores) # 最適化したい評価指標を選ぶ
Optunaの仕様上、最適化したい評価指標は一つしか選べません。ホールドアウトの時と同様、RMSEを最小化するようハイチューします。結果を見てみましょう。
# 最適化
study = optuna.create_study(direction='minimize', study_name='regression')
study.optimize(objective_rf, n_trials=50)
# ハイパーパラメータ・スコアの確認
print("Best trial:")
trial = study.best_trial
print(f" RMSE: {trial.value:.4f}")
print(" Params:")
for key, value in trial.params.items():
print(f" {key}: {value}")
Best trial:
RMSE: 47.2173
Params:
n_estimators: 980
max_depth: 27
min_samples_split: 4
min_samples_leaf: 2
max_features: sqrt
上記のハイパラの時に、交差検証において最もRMSEが下がる結果となりました。それではこのハイパラを用いて最後に確認してみます。
# optunaで最適化されたパラメータをセットし、予測
model_rf_op = RandomForestRegressor(**study.best_params)
# 5分割交差検証
kf = KFold(n_splits=5, shuffle=True, random_state=1234)
# スコア保存用
rmse_scores = []
mae_scores = []
r2_scores = []
for train_index, val_index in kf.split(X):
# 訓練と検証に分類
X_train, X_val = X.iloc[train_index, :], X.iloc[val_index, :]
y_train, y_val = y.iloc[train_index], y.iloc[val_index]
# 標準化、DataFrameに戻す
X_scaler = StandardScaler()
autoscaled_X_train = pd.DataFrame(X_scaler.fit_transform(X_train), columns=X_train.columns, index=X_train.index)
autoscaled_X_val = pd.DataFrame(X_scaler.transform(X_val), columns=X_val.columns, index=X_val.index)
y_scaler = StandardScaler()
autoscaled_y_train = pd.DataFrame(y_scaler.fit_transform(y_train.values.reshape(-1,1)), index=y_train.index, columns=['y'])
autoscaled_y_val = pd.DataFrame(y_scaler.transform(y_val.values.reshape(-1,1)), index=y_val.index, columns=['y'])
# 学習&予測
model_rf_op.fit(autoscaled_X_train, autoscaled_y_train.values.ravel())
y_pred_scaled = model_rf_op.predict(autoscaled_X_val)
y_pred = y_scaler.inverse_transform(y_pred_scaled.reshape(-1,1))
# 性能チェック
rmse = np.sqrt(mean_squared_error(y_val, y_pred))
mae = mean_absolute_error(y_val, y_pred)
r2 = r2_score(y_val, y_pred)
rmse_scores.append(rmse)
mae_scores.append(mae)
r2_scores.append(r2)
print(f'Fold RMSE : {rmse:.4f}')
print(f'Fold MAE : {mae:.4f}')
print(f'Fold R2 : {r2:.4f}')
print()
print(f'平均RMSE : {np.mean(rmse_scores):.4f}')
print(f'平均MAE : {np.mean(mae_scores):.4f}')
print(f'平均R2 : {np.mean(r2_scores):.4f}')
Fold RMSE : 47.2413
Fold MAE : 36.0150
Fold R2 : 0.7869
Fold RMSE : 45.1172
Fold MAE : 33.3644
Fold R2 : 0.7997
Fold RMSE : 50.3175
Fold MAE : 33.4597
Fold R2 : 0.8210
Fold RMSE : 52.3280
Fold MAE : 39.1730
Fold R2 : 0.7663
Fold RMSE : 42.5120
Fold MAE : 29.8906
Fold R2 : 0.8235
平均RMSE : 47.5032
平均MAE : 34.3805
平均R2 : 0.7995
若干ですが、改善が見られました。このハイパラ条件がRFの最終モデルとして良さそうです。この条件を使って、最後に全てのデータを再学習して、保存しておきましょう。次回以降、未知のデータ（本当のテストデータ）を予測するためです。
# 最終モデルの構築
model_rf_final = RandomForestRegressor(**study.best_params)
# すべての学習データを標準化して学習させる
X_scaler_final = StandardScaler()
y_scaler_final = StandardScaler()
autoscaled_X = X_scaler_final.fit_transform(X)
autoscaled_y = y_scaler_final.fit_transform(y.values.reshape(-1, 1))
# 学習
model_rf_final = model_rf_final.fit(autoscaled_X, autoscaled_y.ravel())
# ディレクトリ作成
os.makedirs('models/rf', exist_ok=True)
# モデルとスケーラーの保存（joblib使用）
joblib.dump(model_rf_final, 'models/rf/model_rf.pkl')
joblib.dump(X_scaler_final, 'models/rf/X_scaler.pkl')
joblib.dump(y_scaler_final, 'models/rf/y_scaler.pkl')
■ LightGBM
同様にしてLightGBMでも交差検証し、最後にOptunaでハイパラを最適化しましょう。まずはチューニング前の性能を見ます。
# デフォルト
model_lgb = lgb.LGBMRegressor(random_state=1234, verbose=-1)
# 5分割交差検証
kf = KFold(n_splits=5, shuffle=True, random_state=1234)
# スコア保存用
rmse_scores = []
mae_scores = []
r2_scores = []
for train_index, val_index in kf.split(X):
# 訓練と検証に分類
X_train, X_val = X.iloc[train_index, :], X.iloc[val_index, :]
y_train, y_val = y.iloc[train_index], y.iloc[val_index]
# 標準化、DataFrameに戻す
X_scaler = StandardScaler()
autoscaled_X_train = pd.DataFrame(X_scaler.fit_transform(X_train), columns=X_train.columns, index=X_train.index)
autoscaled_X_val = pd.DataFrame(X_scaler.transform(X_val), columns=X_val.columns, index=X_val.index)
y_scaler = StandardScaler()
autoscaled_y_train = pd.DataFrame(y_scaler.fit_transform(y_train.values.reshape(-1,1)), index=y_train.index, columns=['y'])
autoscaled_y_val = pd.DataFrame(y_scaler.transform(y_val.values.reshape(-1,1)), index=y_val.index, columns=['y'])
# 学習
model_lgb.fit(autoscaled_X_train, autoscaled_y_train, eval_set=[(autoscaled_X_val, autoscaled_y_val)])
y_pred_scaled = model_lgb.predict(autoscaled_X_val)
y_pred = y_scaler.inverse_transform(y_pred_scaled.reshape(-1,1))
# 性能チェック
rmse = np.sqrt(mean_squared_error(y_val, y_pred))
mae = mean_absolute_error(y_val, y_pred)
r2 = r2_score(y_val, y_pred)
rmse_scores.append(rmse)
mae_scores.append(mae)
r2_scores.append(r2)
print(f'Fold RMSE : {rmse:.4f}')
print(f'Fold MAE : {mae:.4f}')
print(f'Fold R2 : {r2:.4f}')
print()
print(f'平均RMSE : {np.mean(rmse_scores):.4f}')
print(f'平均MAE : {np.mean(mae_scores):.4f}')
print(f'平均R2 : {np.mean(r2_scores):.4f}')
Fold RMSE : 45.1902
Fold MAE : 31.6159
Fold R2 : 0.8050
Fold RMSE : 44.6490
Fold MAE : 32.9797
Fold R2 : 0.8038
Fold RMSE : 52.8957
Fold MAE : 34.8776
Fold R2 : 0.8022
Fold RMSE : 42.6535
Fold MAE : 30.5197
Fold R2 : 0.8448
Fold RMSE : 50.3521
Fold MAE : 34.6076
Fold R2 : 0.7523
平均RMSE : 47.1481
平均MAE : 32.9201
平均R2 : 0.8016
チューニング前でも十分な予測精度を示しています。実行時間も短くてありがたいです。改善がみられるかどうかハイチューしておきます。
def objective(trial):
params = {
'objective': 'regression',
'metric': 'rmse', # LightGBMが内部で使う評価指標
'n_estimators': trial.suggest_int('n_estimators', 50, 300),
'learning_rate': trial.suggest_float('learning_rate', 0.01, 0.3, log=True),
'max_depth': trial.suggest_int('max_depth', 3, 10),
'num_leaves': trial.suggest_int('num_leaves', 10, 100),
'subsample': trial.suggest_float('subsample', 0.5, 1.0),
'colsample_bytree': trial.suggest_float('colsample_bytree', 0.5, 1.0),
'reg_alpha': trial.suggest_float('reg_alpha', 1e-8, 1.0, log=True),
'reg_lambda': trial.suggest_float('reg_lambda', 1e-8, 1.0, log=True),
'random_state': 42,
'verbosity': -1,
'device' : 'cpu'
}
# モデル定義
model_lgb = lgb.LGBMRegressor(**params)
# 分割、CV
kf = KFold(n_splits=5, shuffle=True, random_state=1234)
# 評価指標
rmse_scores = []
# mae_scores = []
# r2_scores = []
for train_index, val_index in kf.split(X):
# 訓練と検証に分類
X_train, X_val = X.iloc[train_index, :], X.iloc[val_index, :]
y_train, y_val = y.iloc[train_index], y.iloc[val_index]
# 標準化、DataFrameに戻す
X_scaler = StandardScaler()
autoscaled_X_train = pd.DataFrame(X_scaler.fit_transform(X_train), columns=X_train.columns, index=X_train.index)
autoscaled_X_val = pd.DataFrame(X_scaler.transform(X_val), columns=X_val.columns, index=X_val.index)
y_scaler = StandardScaler()
autoscaled_y_train = pd.DataFrame(y_scaler.fit_transform(y_train.values.reshape(-1,1)), index=y_train.index, columns=['y'])
autoscaled_y_val = pd.DataFrame(y_scaler.transform(y_val.values.reshape(-1,1)), index=y_val.index, columns=['y'])
# 学習
model_lgb.fit(autoscaled_X_train, autoscaled_y_train, eval_set=[(autoscaled_X_val, autoscaled_y_val)])
y_pred_scaled = model_lgb.predict(autoscaled_X_val)
y_pred = y_scaler.inverse_transform(y_pred_scaled.reshape(-1,1))
# 性能チェック
rmse = np.sqrt(mean_squared_error(y_val, y_pred))
# mae = mean_absolute_error(y_val, y_pred)
# r2 = r2_score(y_val, y_pred)
# 結果格納
rmse_scores.append(rmse)
# mae_scores.append(mae)
# r2_scores.append(r2)
return np.mean(rmse_scores) # 最適化したい評価指標を選ぶ
# 最適化
study = optuna.create_study(direction='minimize', study_name='regression')
study.optimize(objective, n_trials=100)
# ハイパーパラメータ・スコアの確認
print("Best trial:")
trial = study.best_trial
print(f" RMSE: {trial.value:.4f}")
print(" Params:")
for key, value in trial.params.items():
print(f" {key}: {value}")
Best trial:
RMSE: 45.9234
Params:
n_estimators: 191
learning_rate: 0.042456292193131684
max_depth: 3
num_leaves: 72
subsample: 0.8129992448601999
colsample_bytree: 0.9852002818804191
reg_alpha: 4.660036691738851e-05
reg_lambda: 8.065184435318226e-05
RFの時と最適化する流れはほとんど同じです。こちらも交差検証の場合においてRMSEが小さくなるよう最適化を実行しています。Optunaによる最適なハイパラの条件が出ました。この条件で最後に予測してみて、精度を確認します。
# optunaで最適化されたパラメータをセットし、予測
model_lgb_op = lgb.LGBMRegressor(**study.best_params)
# 5分割交差検証
kf = KFold(n_splits=5, shuffle=True, random_state=1234)
# スコア保存用
rmse_scores = []
mae_scores = []
r2_scores = []
for train_index, val_index in kf.split(X):
# 訓練と検証に分類
X_train, X_val = X.iloc[train_index, :], X.iloc[val_index, :]
y_train, y_val = y.iloc[train_index], y.iloc[val_index]
# 標準化、DataFrameに戻す
X_scaler = StandardScaler()
autoscaled_X_train = pd.DataFrame(X_scaler.fit_transform(X_train), columns=X_train.columns, index=X_train.index)
autoscaled_X_val = pd.DataFrame(X_scaler.transform(X_val), columns=X_val.columns, index=X_val.index)
y_scaler = StandardScaler()
autoscaled_y_train = pd.DataFrame(y_scaler.fit_transform(y_train.values.reshape(-1,1)), index=y_train.index, columns=['y'])
autoscaled_y_val = pd.DataFrame(y_scaler.transform(y_val.values.reshape(-1,1)), index=y_val.index, columns=['y'])
# 学習
model_lgb_op.fit(autoscaled_X_train, autoscaled_y_train, eval_set=[(autoscaled_X_val, autoscaled_y_val)])
y_pred_scaled = model_lgb_op.predict(autoscaled_X_val)
y_pred = y_scaler.inverse_transform(y_pred_scaled.reshape(-1,1))
# 性能チェック
rmse = np.sqrt(mean_squared_error(y_val, y_pred))
mae = mean_absolute_error(y_val, y_pred)
r2 = r2_score(y_val, y_pred)
rmse_scores.append(rmse)
mae_scores.append(mae)
r2_scores.append(r2)
print(f'Fold RMSE : {rmse:.4f}')
print(f'Fold MAE : {mae:.4f}')
print(f'Fold R2 : {r2:.4f}')
print()
print(f'平均RMSE : {np.mean(rmse_scores):.4f}')
print(f'平均MAE : {np.mean(mae_scores):.4f}')
print(f'平均R2 : {np.mean(r2_scores):.4f}')
Fold RMSE : 46.1686
Fold MAE : 33.6581
Fold R2 : 0.7965
Fold RMSE : 45.0651
Fold MAE : 31.7902
Fold R2 : 0.8001
Fold RMSE : 50.3386
Fold MAE : 31.2969
Fold R2 : 0.8209
Fold RMSE : 42.1905
Fold MAE : 31.5493
Fold R2 : 0.8481
Fold RMSE : 50.2708
Fold MAE : 35.7656
Fold R2 : 0.7531
平均RMSE : 46.8067
平均MAE : 32.8120
平均R2 : 0.8037
あまりチューニングしても精度は大きく変わらないことが分かりました。平均的な予測精度は安定して優れていることが分かります。
最後に全てのデータを学習させたこのモデルも保存しておきましょう。
# 最終モデルの構築
model_lgb_final = lgb.LGBMRegressor(**study.best_params)
# すべての学習データを標準化して、学習させる
X_scaler_final = StandardScaler()
y_scaler_final = StandardScaler()
autoscaled_X = X_scaler_final.fit_transform(X)
autoscaled_y = y_scaler_final.fit_transform(y.values.reshape(-1, 1))
# 学習
model_lgb_final = model_lgb_final.fit(autoscaled_X, autoscaled_y.ravel())
# ディレクトリ作成（なければ）
os.makedirs('models/lgb', exist_ok=True)
# モデルとスケーラーの保存
joblib.dump(model_lgb_final, 'models/lgb/model_lgb.pkl')
joblib.dump(X_scaler_final, 'models/lgb/X_scaler.pkl')
joblib.dump(y_scaler_final, 'models/lgb/y_scaler.pkl')
■ Neural Network
最後にNNで交差検証します。だいぶ長くなってきたので、こやつはいきなりOptunaで最適化します。ホールドアウトの時と同様に、Datasetを準備して、便宜上、学習と評価の関数を定義しておきます。
# データセット定義
class RegressionDataset(Dataset):
def __init__(self, X, y):
self.X = torch.tensor(X, dtype=torch.float32)
self.y = torch.tensor(y, dtype=torch.float32).view(-1, 1)
def __len__(self):
return len(self.X)
def __getitem__(self, idx):
return self.X[idx], self.y[idx]
# 学習loop
def train(model, dataloader, optimizer, criterion, device):
model.train()
total_loss = 0
for X_batch, y_batch in dataloader:
X_batch, y_batch = X_batch.to(device), y_batch.to(device)
optimizer.zero_grad()
pred = model(X_batch)
loss = criterion(pred, y_batch)
loss.backward()
optimizer.step()
total_loss += loss.item()
return total_loss / len(dataloader)
# 評価
def evaluate(model, dataloader, criterion, device):
model.eval()
total_loss = 0.0
with torch.no_grad():
for X_batch, y_batch in dataloader:
X_batch, y_batch = X_batch.to(device), y_batch.to(device)
pred = model(X_batch)
loss = criterion(pred, y_batch)
total_loss += loss.item()
return total_loss / len(dataloader)
そしてOptunaの条件設定です。define_model関数の中で、条件検討したいハイパラの候補を挙げて、これをobjective関数の中で呼び出します。この辺もホールドアウトの時と大きく変わりません。
# Optuna設定
def define_model(trial, input_dim):
layers = []
n_layers = trial.suggest_int("n_layers", 2, 5)
hidden_dim = trial.suggest_int("hidden_dim", 32, 1024)
dropout_rate = trial.suggest_float("dropout_rate", 0.0, 0.5)
activation_name = trial.suggest_categorical("activation", ["relu", "leaky_relu"])
use_batchnorm = trial.suggest_categorical("use_batchnorm", [True, False])
in_dim = input_dim
for i in range(n_layers):
layers.append(nn.Linear(in_dim, hidden_dim))
# 活性化関数の前に入れる
if use_batchnorm:
layers.append(nn.BatchNorm1d(hidden_dim))
if activation_name == "relu":
layers.append(nn.ReLU())
else:
layers.append(nn.LeakyReLU())
layers.append(nn.Dropout(dropout_rate))
in_dim = hidden_dim
layers.append(nn.Linear(in_dim, 1))
return nn.Sequential(*layers)
def objective(trial):
# ハイパーパラメータ提案
batch_size = trial.suggest_categorical("batch_size", [16, 32, 64, 128])
lr = trial.suggest_float("lr", 1e-5, 1e-3, log=True)
optimizer_name = trial.suggest_categorical("optimizer", ["Adam", "SGD"])
epochs = trial.suggest_int("epochs", 30, 100)
# モデル定義
model = define_model(trial, X.shape[1]).to(device)
criterion = nn.MSELoss()
if optimizer_name == 'Adam':
optimizer = torch.optim.Adam(model.parameters(), lr=lr)
else:
optimizer = torch.optim.SGD(model.parameters(), lr=lr)
kf = KFold(n_splits=5, shuffle=True, random_state=1234)
val_losses = []
for train_idx, val_idx in kf.split(X):
X_train, X_val = X.iloc[train_idx], X.iloc[val_idx]
y_train, y_val = y.iloc[train_idx], y.iloc[val_idx]
# 標準化
X_scaler = StandardScaler()
autoscaled_X_train = X_scaler.fit_transform(X_train)
autoscaled_X_val = X_scaler.transform(X_val)
y_scaler = StandardScaler()
autoscaled_y_train = y_scaler.fit_transform(y_train.values.reshape(-1, 1))
autoscaled_y_val = y_scaler.transform(y_val.values.reshape(-1, 1))
# 各CV毎にDataset, DataLoaderを設定
train_ds = RegressionDataset(autoscaled_X_train, autoscaled_y_train)
val_ds = RegressionDataset(autoscaled_X_val, autoscaled_y_val)
train_dl = DataLoader(train_ds, batch_size=batch_size, shuffle=True)
val_dl = DataLoader(val_ds, batch_size=batch_size, shuffle=False)
# 学習、評価
for epoch in range(epochs):
train(model, train_dl, optimizer, criterion, device)
val_loss = evaluate(model, val_dl, criterion, device)
val_losses.append(val_loss)
return np.mean(val_losses)
ホールドアウトとの唯一の違いは、学習データと検証データに分けた後、分けたデータをDataset、DataLoaderに入れるという操作を各CV毎に行うということです。当たり前っちゃ当たり前なんですが、意外と混乱します。各CVにおいてミニバッチ学習をするため、このような流れになります。
それでは、学習スタートです。
# 実行
print('Optuna最適化開始...')
study = optuna.create_study(direction='minimize')
study.optimize(objective, n_trials=50) # 試行回数を増やす
print('Best trial :', study.best_trial.params)
Best trial : {'batch_size': 32, 'lr': 0.0001366811864427072, 'optimizer': 'Adam', 'epochs': 32, 'n_layers': 2, 'hidden_dim': 667, 'dropout_rate': 0.10335891509292101, 'activation': 'leaky_relu', 'use_batchnorm': False}
中間層のノードは多いですが、層数は少ない構成になりました。このチューンされたNNの実力を見てみます。
# optunaで最適化されたパラメータをセットし、予測
# 評価指標格納用
rmse_scores = []
mae_scores = []
r2_scores = []
kf = KFold(n_splits=5, shuffle=True, random_state=1234)
for train_index, val_index in kf.split(X):
# 分割
X_train = X.iloc[train_index]
X_val = X.iloc[val_index]
y_train = y.iloc[train_index]
y_val = y.iloc[val_index]
# 標準化
X_scaler = StandardScaler()
autoscaled_X_train = X_scaler.fit_transform(X_train)
autoscaled_X_val = X_scaler.transform(X_val)
y_scaler = StandardScaler()
autoscaled_y_train = y_scaler.fit_transform(y_train.values.reshape(-1, 1))
autoscaled_y_val = y_scaler.transform(y_val.values.reshape(-1, 1))
# DataLoader
train_ds = RegressionDataset(autoscaled_X_train, autoscaled_y_train)
val_ds = RegressionDataset(autoscaled_X_val, autoscaled_y_val)
train_dl = DataLoader(train_ds, batch_size=study.best_trial.params['batch_size'], shuffle=True)
val_dl = DataLoader(val_ds, batch_size=study.best_trial.params['batch_size'], shuffle=False)
# モデル構築
model_nn_op = define_model(study.best_trial, input_dim=autoscaled_X_train.shape[1]).to(device)
optimizer_name = study.best_trial.params['optimizer']
lr = study.best_trial.params['lr']
optimizer = torch.optim.Adam(model_nn_op.parameters(), lr=lr) if optimizer_name == "Adam" else torch.optim.SGD(model_nn_op.parameters(), lr=lr)
criterion = nn.MSELoss()
# 学習ループ
epochs = study.best_trial.params['epochs'] # Optuna で最適化されたエポック数を使う
for epoch in range(epochs):
train(model_nn_op, train_dl, optimizer, criterion, device)
# 推論と逆変換
model_nn_op.eval()
preds = []
with torch.no_grad():
for X_batch, _ in val_dl:
X_batch = X_batch.to(device)
y_pred_std = model_nn_op(X_batch).cpu().numpy()
preds.append(y_pred_std)
y_pred = y_scaler.inverse_transform(np.vstack(preds))
y_true = y_scaler.inverse_transform(autoscaled_y_val)
# 評価指標
rmse = np.sqrt(mean_squared_error(y_true, y_pred))
mae = mean_absolute_error(y_true, y_pred)
r2 = r2_score(y_true, y_pred)
print(f'Fold RMSE : {rmse:.4f}')
print(f'Fold MAE : {mae:.4f}')
print(f'Fold R2 : {r2:.4f}\n')
rmse_scores.append(rmse)
mae_scores.append(mae)
r2_scores.append(r2)
# 結果出力
print(f'平均RMSE: {np.mean(rmse_scores):.4f} ± {np.std(rmse_scores):.4f}')
print(f'平均MAE: {np.mean(mae_scores):.4f} ± {np.std(mae_scores):.4f}')
print(f'平均R2: {np.mean(r2_scores):.4f} ± {np.std(r2_scores):.4f}')
Fold RMSE : 45.0466
Fold MAE : 33.4744
Fold R2 : 0.8063
Fold RMSE : 43.9383
Fold MAE : 30.4366
Fold R2 : 0.8100
Fold RMSE : 69.4819
Fold MAE : 37.9011
Fold R2 : 0.6587
Fold RMSE : 42.3546
Fold MAE : 33.2047
Fold R2 : 0.8469
Fold RMSE : 58.4803
Fold MAE : 37.3230
Fold R2 : 0.6659
平均RMSE: 51.8604 ± 10.5251
平均MAE: 34.4680 ± 2.7847
平均R2: 0.7576 ± 0.0791
今回はばらつきが大きいことが分かります。もう少し改善はできそうですが、ひとまずは一結果として飲み込んでおきましょう。最後に全部の学習データを学習させて保存して終わります。
# すべてのデータを学習させて、保存
# 全データで標準化
X_scaler_final = StandardScaler()
y_scaler_final = StandardScaler()
autoscaled_X = X_scaler_final.fit_transform(X)
autoscaled_y = y_scaler_final.fit_transform(y.values.reshape(-1, 1))
# Dataset, DataLoader
final_ds = RegressionDataset(autoscaled_X, autoscaled_y)
final_dl = DataLoader(final_ds, batch_size=study.best_trial.params['batch_size'], shuffle=True)
# モデル定義
final_model = define_model(study.best_trial, input_dim=autoscaled_X.shape[1]).to(device)
optimizer_name = study.best_trial.params['optimizer']
lr = study.best_trial.params['lr']
optimizer = torch.optim.Adam(final_model.parameters(), lr=lr) if optimizer_name == "Adam" else torch.optim.SGD(final_model.parameters(), lr=lr)
criterion = nn.MSELoss()
epochs = study.best_trial.params['epochs']
# 学習ループ
for epoch in range(epochs):
train(final_model, final_dl, optimizer, criterion, device)
# 保存先ディレクトリを作成（なければ）
os.makedirs('models/nn', exist_ok=True)
# PyTorchモデルの保存（state_dict）、モデルの重み
torch.save(final_model.state_dict(), 'models/nn/final_model.pth')
# Optunaのbest_trial(ハイパラ)
joblib.dump(study.best_trial.params, 'models/nn/hparams.pkl')
# 特徴量名（列順の再現用)
joblib.dump(list(X.columns), 'models/nn/feature_names.pkl')
# スケーラーの保存（joblib）
joblib.dump(X_scaler_final, 'models/nn/X_scaler.pkl')
joblib.dump(y_scaler_final, 'models/nn/y_scaler.pkl')
結果まとめ
今までのモデルの精度まとめを記載しておきます。すべてOptunaで最適化した結果です。
表1. ホールドアウト結果まとめ
| 特徴量 | モデル | RMSE (Train/Test) | MAE (Train/Test) | R2 (Train/Test) |
|---|---|---|---|---|
| RDKit | PLS | 25.49 / 72.90 | 19.37 / 50.16 | 0.95 / 0.47 |
| RF | 24.09 / 51.29 | 17.59 / 38.04 | 0.95 / 0.74 | |
| XGB | 6.68 / 53.86 | 1.86 / 37.52 | 1.00 / 0.71 | |
| LGB | 7.76 / 52.56 | 3.44 / 36.37 | 1.00 / 0.72 | |
| NN | 19.96 / 49.77 | 16.63 / 31.78 | 0.97 / 0.75 | |
| Moredred 2D (3D) | PLS | 36.76 / 51.16 | 27.45 / 38.81 | 0.89 / 0.74 |
| RF | 25.19 / 52.19 | 18.04 / 39.03 | 0.95 / 0.73 | |
| XGB | 6.55 / 50.29 | 1.56 / 35.52 | 1.00 / 0.75 | |
| LGB | 6.54 / 50.86 | 1.63 / 34.64 | 1.00 / 0.74 | |
| NN | 9.96 / 48.49 | 6.95 / 33.61 | 0.99 / 0.76 |
表2. 交差検証結果まとめ(検証データにおける評価指標のみ記載)
| 特徴量 | モデル | RMSE (CV) | MAE (CV) | R2 (CV) |
|---|---|---|---|---|
| RDKit | PLS※ | 60.32 | 43.06 | 0.65 |
| RF | 47.17 | 33.57 | 0.80 | |
| XGB | 44.14 | 30.78 | 0.83 | |
| LGB | 46.38 | 31.50 | 0.81 | |
| NN※ | 52.65 | 38.67 | 0.74 | |
| Moredred 2D (3D) | PLS | 56.74 | 39.54 | 0.70 |
| RF | 47.48 | 34.29 | 0.80 | |
| XGB | 44.58 | 30.53 | 0.82 | |
| LGB | 46.38 | 32.57 | 0.80 | |
| NN | 54.72 | 35.34 | 0.73 |
※RDKitにおけるPLS, NNの結果は4-fold分の平均値(1つのfoldで評価指標が異常値であったため)
特徴量の数はrdkitが152個、mordred_2dで1219個、mordred_3dも欠損値が入っている列を消してみれば1219個で2dと同じでした。
まずホールドアウトと交差検証で結果が結構変わっていることが分かります。ホールドアウトの結果だけ見るとNNが優秀ですが、交差検証の結果だと決定木系のモデルが優秀です。
RF, XGB, LGBはデータセットの影響を受けにくいようで、説明変数の数が少ないRDKit記述子でもうまく動いています。一方、PLSやNNはRDKitを使うと、異常な評価値がでてしまうfoldがあり、うまく評価できない場合がありました。Mordredを使った場合だと、このような異常値が見られず、どのモデルでも安定した出力が得られます。
よって結論として
「Mordred記述子を使えばどのモデルも安定した予測値が得られる。その中でも決定木系は評価指標の値が優れている。NNも良い精度を示す時があるが、データセットに依存しやすい」
「最も良いモデルでも、MAEが30程度で±30nmは予測を外している。この辺をどう修正していくかがポイント」
データセットにおけるPL波長の分布が青から赤と結構広いので、青限定のデータセットとか赤限定のデータセットとして、サンプルを増やしてデータ拡張するのが良いと思われます。しかしデータを増やすのはなかなか時間がかかるので、計算値などで補填していく方が現実的かもしれません。
最後に
第2回お疲れさまでした。かなり長くて途中で断念した方もいるかと思います。なぜこんなことをしたかというと、MLの教科書など見るとホールドアウトの記載が多く、交差検証でハイチューする内容があまりありませんでした。ホールドアウトでハイチューしてしまうと、データセットが小さい場合はうまくいかないことが多く、もやもやしていました。
交差検証+GridSearchなどでハイチューする例は載っていたりするのですが、これもすごく時間がかかってしまい、実務向きではなく、なんとかoptunaと交差検証を組み合わすことができないかと色々試行錯誤した結果、今回記載したコードに落ち着いたわけです。最初に作ったコードを生成AIたちに修正してもらい、何とかここまで持っていけました(NNはまだコード修正が必要ですが、、、）
興味ある方ぜひ使ってみてください。データセットがそれなりに大きければ動くと思われます。
次回は未知の材料を予測して、その予測値が信頼できるかどうか、適用範囲(Application Domain)を考慮して眺めてみようと思います。
最後まで読んでいただきありがとうございました。
引用文献
・化学のためのPythonによるデータ解析・機械学習入門
https://datachemeng.com/post-4014/
前回同様、今回も引用です。いつもお世話になっています。PLSなどの理論が記載されています。
・機械学習はこれ一本！pythonインストール～機械学習実装まで完全理解講座
https://brain-market.com/u/rin-effort/a/bEjMwYjMgoTZsNWa0JXY
yyプロットの書き方など参考にさせていただきました。有料の記事ですが、初学者に読みやすい内容となっています。