🤖
シンボリック回帰PySRで外挿予測性能評価
はじめに
前回はシンボリック回帰のインストールから簡単なチュートリアルによる動作確認までを行いました。
この記事では、裏に方程式が隠れている
環境
- Ubuntu 20.04.1
- Python 3.9.15
- japanize-matplotlib 1.1.3
- matplotlib 3.6.2
- numpy 1.23.5
- optuna 2.10.1
- pandas 1.5.2
- pysr 0.12.3
- scikit-learn 1.1.3
データ
関数を定義します。ここでは、説明変数に相当する
外挿予測性能を評価するため、
ベンチマーク
学習データについて5-fold cvでMSEが最小になるようにoptunaでハイパーパラメータを決定し、テストデータに対する予測性能をMSEで評価しました。以下コードです。
import pandas as pd
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import mean_squared_error # MSEで比較
from sklearn.linear_model import ElasticNet
from sklearn.svm import SVR
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import cross_val_score, KFold
import optuna
def plot_results(y_tr, y_tr_pred, y_te, y_te_pred):
"""
Display a scatter plot comparing the actual and predicted values for training and testing data.
Parameters:
-----------
y_tr : array-like
Array containing the actual target values for training data.
y_tr_pred : array-like
Array containing the predicted target values for training data.
y_te : array-like
Array containing the actual target values for testing data.
y_te_pred : array-like
Array containing the predicted target values for testing data.
Returns:
--------
None
"""
__min = min([y_tr.min(), y_tr_pred.min(), y_te.min(), y_te_pred.min()])
__max = max([y_tr.max(), y_tr_pred.max(), y_te.max(), y_te_pred.max()])
min1 = __min - (__max - __min)*0.1
max1 = __max + (__max - __min)*0.1
mse_te = mean_squared_error(y_te, y_te_pred)
plt.figure(figsize=(5,5))
plt.title(f"MSE for test={round(mse_te,2)}")
plt.scatter(y_tr, y_tr_pred, zorder=2, label="trained data")
plt.scatter(y_te, y_te_pred, zorder=3, label="tested data")
plt.xlim(min1, max1)
plt.ylim(min1, max1)
plt.grid()
plt.plot([min1, max1],[min1, max1], color="gray", zorder=1)
plt.legend()
# データの生成
np.random.seed(1)
X = 2 * np.random.randn(100, 5)
noise = np.random.normal(0,1,y.size)
y = 2.5382 * np.cos(X[:, 3]) + X[:, 0] ** 2 - 0.5 + noise
# 可視化
df = pd.DataFrame(
np.c_[y,X], columns=["y"]+[f"x{i}" for i in range(5)]
)
df["label"] = "train"
df.loc[df["y"]>=15, "label"]="test"
sns.pairplot(df, hue = "label")
# テストデータと学習データに分割
df_te = df.query("y>=15")
df_tr = df.query("y<15")
X_te = df_te.drop("y",axis=1)
X_tr = df_tr.drop("y",axis=1)
y_te = df_te["y"]
y_tr = df_tr["y"]
Random Forest
def objective(trial):
ne = trial.suggest_int("n_estimators", 10,200)
regr = RandomForestRegressor(
n_estimators=ne,
random_state=123
)
score = cross_val_score(regr, X_tr, y_tr, cv=5, scoring="neg_mean_squared_error")
mse_mean = score.mean()
return mse_mean
study = optuna.create_study(direction='maximize')
study.optimize(objective, n_trials=100)
clf_rf_best = RandomForestRegressor(
**study.best_params, random_state=123
)
clf_rf_best.fit(X_tr, y_tr)
y_te_pred = clf_rf_best.predict(X_te)
y_tr_pred = clf_rf_best.predict(X_tr)
print("MSE :", mean_squared_error(y_te, y_te_pred))
plot_results(y_tr, y_tr_pred, y_te, y_te_pred)
plt.show()
Random Forestは原理上学習データの外側を予測値として算出不可のため当然の結果です。
SVR
def objective(trial):
g = trial.suggest_loguniform("gamma", 2**-10, 2**10)
c = trial.suggest_loguniform("C", 2**-10, 2**10)
e = trial.suggest_loguniform("epsilon", 2**-10, 2**10)
regr = SVR(
kernel="rbf",
gamma=g,
C=c,
epsilon=e
)
score = cross_val_score(regr, X_tr, y_tr, cv=5, scoring="neg_mean_squared_error")
mse_mean = score.mean()
return mse_mean
study = optuna.create_study(direction='maximize')
study.optimize(objective, n_trials=200)
clf_svr_best = SVR(
**study.best_params
)
clf_svr_best.fit(X_tr, y_tr)
y_te_pred = clf_svr_best.predict(X_te)
y_tr_pred = clf_svr_best.predict(X_tr)
print("MSE :", mean_squared_error(y_te, y_te_pred))
plot_results(y_tr, y_tr_pred, y_te, y_te_pred)
plt.show()
rbfカーネルを使用したSVRであればそれなりに外挿予測できます。リーズナブルですね。ただし、SVRの場合は
ElasticNet
def objective(trial):
a = trial.suggest_loguniform("alpha", 2**-12, 2**10)
l = trial.suggest_uniform("l1_ratio", 0, 1)
regr = ElasticNet(
alpha=a, l1_ratio=l
)
score = cross_val_score(regr, X_tr, y_tr, cv=5, scoring="neg_mean_squared_error")
mse_mean = score.mean()
return mse_mean
study = optuna.create_study(direction='maximize')
study.optimize(objective, n_trials=200)
clf_en_best = ElasticNet(
**study.best_params
)
clf_en_best.fit(X_tr, y_tr)
y_te_pred = clf_en_best.predict(X_te)
y_tr_pred = clf_en_best.predict(X_tr)
print("MSE :", mean_squared_error(y_te, y_te_pred))
plot_results(y_tr, y_tr_pred, y_te, y_te_pred)
plt.show()
Xからyを生成した関数に線形の要素が存在しないので、線形回帰の手法では学習データに関してもモデル化できないのは当然の結果です。
シンボリック回帰
さて、本題のPySRを使ったシンボリック回帰です。
from pysr import PySRRegressor
model = PySRRegressor(
niterations=40, # < Increase me for better results
binary_operators=["+", "*", "/"],
unary_operators=[
"cos",
"exp",
"sin",
"inv(x) = 1/x",
# ^ Custom operator (julia syntax)
],
extra_sympy_mappings={"inv": lambda x: 1 / x},
# ^ Define operator for SymPy as well
loss="loss(prediction, target) = (prediction - target)^2",
# ^ Custom loss function (julia syntax)
)
# 学習データでfit
model.fit(X_tr,y_tr)
# bestモデルで予測
y_tr_pred = model.predict(X_tr) # best_modelで予測値を返す
y_te_pred = model.predict(X_te) # best_modelで予測値を返す
print("MSE :", mean_squared_error(y_te, y_te_pred))
plot_results(y_tr, y_tr_pred, y_te, y_te_pred)
plt.show()
テストデータに対するMSE=1.49という良い性能が得られることを確認できました。
シンボリック回帰の良いところは、得られたモデルが式として自明になっていることです。見てみましょう。
model.get_best()
equationの行をきれいに書くと、
となります。
真の関数は
でしたので、概ね真の関数を再現できていることがわかります。切片項が消えてしまいましたが、これは切片項を-0.5、
さいごに
無闇に機械学習手法を適用すると学習データに過学習したモデルしか得られず、外挿予測など汎用性に欠如してしまうことがあります。データの裏側に何かしらの方程式が想定される場合には、シンボリック回帰を適用してみることで外挿予測性も良く、非線形であっても解釈性の良いモデルを得ることができまる場合があります。ビジネス課題解決に従事するデータサイエンティストであれば一手法として手札に持っておきたいですね。
以上
Discussion