AI2L × Symbolic Regression：Q10 則（温度依存の指数則）をダミーデータから引き当てる

ここで扱うデータは すべてダミー であり、実際の実験データの話ではありません。
「AI2L 的に Symbolic Regression をどう使うか？」という 道具箱の紹介 を行います。

この記事でやること

化学・生物系でよく出てくる経験則のひとつに Q10 則 があります。

• 「温度を 10 ℃ 上げると反応速度がだいたい 2 倍になる」

というものです。これを少しきちんと書くと、例えば基準温度 20 ℃ での速度を 1 と決めるとき、
温度 (T) における反応速度 (v(T)) は

$$
v(T) = 2^{(T-20)/10} = \exp\left(\alpha (T-20)\right),
\quad \alpha = \frac{\ln 2}{10}
$$

と指数関数の形で表せます。

この記事では、この Q10 型の指数則から

1. ダミーデータを生成し
2. 線形回帰・多項式回帰・ランダムフォレストでフィットして
3. Symbolic Regression（シンボリック回帰） で式そのものを引き当てる

という流れを、AI2L（AI to Learn）の視点で整理します。

ポイントは、

• 学習範囲（0–40 ℃）ではどのモデルもそこそこ良く見えるのに、
• 外挿範囲（40–80 ℃）では
  式の形を理解している Symbolic Regression だけが指数関数に追従できる

という構図がはっきり出るように設計してあることです。

---

0. 本記事の結果を出したコードについて

• この記事のコードはすべて本文末尾に掲載しています（Google Colab でそのまま動作）。
• 画像は Colab から保存して Qiita に貼り付けたものです。

---

1. 題材：Q10 則を数式として書き直す

Q10 則を少し整理して、ダミーデータの「真の法則」として使える形にします。

• 基準温度：T_ref = 20 ℃
• Q10：Q10 = 2（10 ℃ 上昇で速度 2 倍）
• 基準温度での速度：v(T_ref) = 1 とする

このとき、温度 T における速度 v(T) は

$$
v(T) = 2^{(T-T_\text{ref})/10}
= \exp\left( \alpha (T-T_\text{ref}) \right),
\quad \alpha = \frac{\ln 2}{10}
$$

となります。

AI2L 的には、まず

• どの法則を「真」とみなしているのか
• それをどう数式で書いているのか

を明示しておくことが重要です。

---

2. ダミーデータの設計

2.1 温度範囲と外挿範囲

今回は次のように範囲を分けました。

• 学習用温度範囲：0–40 ℃
• 外挿の評価範囲：40–80 ℃

Q10 則に従って真の値を計算し、観測ノイズ（標準偏差 0.01）を足して観測値 v_obs を作ります。

T_ref = 20.0
Q10 = 2.0
alpha = np.log(Q10) / 10.0

def true_rate(T, alpha=alpha, T_ref=T_ref):
    return np.exp(alpha * (T - T_ref))

def to_feature(T):
    # 後で学習用に使う、中心化・スケーリング済み特徴量
    return (T - T_ref) / 20.0

# 学習範囲とテスト範囲
T_train = np.linspace(0.0, 40.0, 50)
T_test_in  = np.linspace(0.0, 40.0, 200)
T_test_out = np.linspace(40.0, 80.0, 400)

# ノイズ付き観測
noise_std = 0.01
v_train_true = true_rate(T_train)
v_train_obs  = v_train_true + np.random.normal(scale=noise_std,
                                               size=T_train.shape)

df_train = pd.DataFrame({
    "T": T_train,
    "v_obs": v_train_obs,
    "v_true": v_train_true
})
df_train.head()

先頭 5 行はこんな感じです：

T	v_obs	v_true
0.000000	0.254967	0.250000
0.816327	0.263171	0.264554
1.632653	0.286432	0.279955
2.448980	0.311483	0.296252
3.265306	0.311157	0.313499

2.2 データの可視化

plt.figure(figsize=(6, 4))
plt.scatter(df_train["T"], df_train["v_obs"],
            label="Observed (noisy)", alpha=0.7)
plt.plot(df_train["T"], df_train["v_true"],
         label="True law v(T)=exp(alpha*(T-T_ref))", linewidth=2)
plt.xlabel("Temperature T (°C)")
plt.ylabel("Reaction rate v(T)")
plt.title("Training data vs true Q10-type law (T=0-40 °C)")
plt.legend()
plt.grid(True)
plt.show()

0–40 ℃ の範囲で、指数関数としての曲がりがそれなりに見えているのが分かります。

---

3. ベースライン：線形／3 次多項式／ランダムフォレスト

AI2L 的には、いきなり高度な手法に飛びつくのではなく、

まず素朴なベースラインを用意し、
どこで限界が来るかを可視化する

ところから始めます。

3.1 特徴量の設計

指数関数を扱う都合上、温度そのもの T ではなく

$$
X = \frac{T - T_\text{ref}}{20}
$$

という 中心化＋スケーリングした特徴量 を全モデルで共通に使いました。
こうしておくと、Symbolic Regression でも「指数の形」が見つかりやすくなります。

3.2 モデルと評価関数

def evaluate_model(name, model, T_train, v_train,
                   T_test_in, T_test_out):
    X_train = to_feature(T_train).reshape(-1, 1)
    X_in  = to_feature(T_test_in).reshape(-1, 1)
    X_out = to_feature(T_test_out).reshape(-1, 1)

    model.fit(X_train, v_train)

    v_pred_in  = model.predict(X_in)
    v_pred_out = model.predict(X_out)

    rmse_in = np.sqrt(mean_squared_error(true_rate(T_test_in),  v_pred_in))
    rmse_out = np.sqrt(mean_squared_error(true_rate(T_test_out), v_pred_out))

    return {
        "name": name,
        "rmse_in": rmse_in,
        "rmse_out": rmse_out,
        "v_pred_in": v_pred_in,
        "v_pred_out": v_pred_out,
    }

lin_reg = LinearRegression()
poly_reg = Pipeline([
    ("poly", PolynomialFeatures(degree=3)),
    ("lin", LinearRegression())
])
rf_reg = RandomForestRegressor(
    n_estimators=200,
    max_depth=4,
    random_state=42
)

results_baseline = []
results_baseline.append(
    evaluate_model("LinearRegression", lin_reg,
                   T_train, v_train_obs, T_test_in, T_test_out)
)
results_baseline.append(
    evaluate_model("Poly3+Linear", poly_reg,
                   T_train, v_train_obs, T_test_in, T_test_out)
)
results_baseline.append(
    evaluate_model("RandomForest", rf_reg,
                   T_train, v_train_obs, T_test_in, T_test_out)
)

pd.DataFrame([
    {"model": res["name"],
     "RMSE_in(0-40)": res["rmse_in"],
     "RMSE_out(40-80)": res["rmse_out"]}
    for res in results_baseline
])

実際に得られた RMSE は次の通りです：

model	RMSE_in(0–40 ℃)	RMSE_out(40–80 ℃)
LinearRegression	0.340666	22.963070
Poly3+Linear	0.013763	12.652544
RandomForest	0.023305	24.224077

• 学習範囲内では 多項式回帰と RandomForest はかなり良い（RMSE ≈ 0.02 以下）
• しかし外挿範囲では どれも指数則から大きく外れる ことが分かります。

3.3 グラフで挙動を確認する

plt.figure(figsize=(10, 4))

# Training range
plt.subplot(1, 2, 1)
plt.scatter(T_train, v_train_obs,
            label="Training data (noisy)", alpha=0.6)
plt.plot(T_test_in, true_rate(T_test_in),
         label="True law", linewidth=2)
for res in results_baseline:
    plt.plot(T_test_in, res["v_pred_in"], label=res["name"])
plt.xlabel("Temperature T (°C)")
plt.ylabel("Reaction rate v(T)")
plt.title("Training range T=0-40 °C")
plt.legend()
plt.grid(True)

# Extrapolation range
plt.subplot(1, 2, 2)
plt.plot(T_test_out, true_rate(T_test_out),
         label="True law", linewidth=2)
for res in results_baseline:
    plt.plot(T_test_out, res["v_pred_out"], label=res["name"])
plt.xlabel("Temperature T (°C)")
plt.ylabel("Reaction rate v(T)")
plt.title("Extrapolation range T=40-80 °C")
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

左図だけを見ると、

• Poly3+Linear と RandomForest は「よく当たっているように見える」
• 線形回帰だけが少し怪しい

という印象です。

しかし右図（40–80 ℃）では

• 線形回帰：指数の立ち上がりにまったく追いつかず、大きく下振れ
• 多項式回帰：曲率をがんばって真似しているものの、80 ℃ ではまだ指数に追いつかない
• RandomForest：学習範囲内の値の補間しかできないため、6 付近でほぼ飽和してしまう

といった 外挿の弱さ がはっきり見えます。

---

4. Symbolic Regression で Q10 則を引き当てる

4.1 関数セットの設計

指数則を当ててもらいたいので、SymbolicRegressor には

• add（足し算）
• mul（掛け算）
• exp（指数関数）

だけを渡します。
また、複雑さへのペナルティ parsimony_coefficient を 0 にして、

「多少長くても MSE が良い式」を優先

させました（今回は 式の短さよりも、物理的に正しい形を優先）。

from gplearn.functions import make_function

def _exp(x):
    x = np.clip(x, -10, 10)  # 数値爆発の防止
    return np.exp(x)

exp = make_function(function=_exp, name="exp", arity=1)
function_set = ("add", "mul", exp)

sr = SymbolicRegressor(
    population_size=5000,
    generations=40,
    function_set=function_set,
    const_range=(-2.0, 2.0),
    p_crossover=0.7,
    p_subtree_mutation=0.1,
    p_hoist_mutation=0.05,
    p_point_mutation=0.1,
    max_samples=0.9,
    parsimony_coefficient=0.0,
    verbose=1,
    random_state=42
)

X_train = to_feature(T_train).reshape(-1, 1)
sr.fit(X_train, v_train_obs)

print("Best individual (LISP-style):")
print(sr)

学習ログの最後で得られた 最良個体（LISP 風表記） は：

Best individual (LISP-style):
exp(add(X0, mul(X0, 0.386)))

ここで X0 は

$$
X_0 = \frac{T - T_\text{ref}}{20}
$$

ですから、得られた式は

$$
\hat{v}(T)
= \exp\left( X_0 + 0.386 X_0 \right)
= \exp\left( 1.386 X_0 \right)
= \exp\left( \frac{1.386}{20} (T - T_\text{ref}) \right)
$$

となります。

一方、真の Q10 則は

$$
v(T) = \exp\left( \alpha (T - T_\text{ref}) \right),
\quad \alpha = \frac{\ln 2}{10}
$$

でした。1.386 は実は ln(4) にほぼ等しいので、

$$
\frac{1.386}{20} \approx \frac{\ln 4}{20}
= \frac{\ln 2}{10} = \alpha
$$

つまり

SymbolicRegressor は、Q10=2 の指数則そのものにかなり近い式を再発見した

と言えます。

4.2 性能比較（数値）

X_in  = to_feature(T_test_in).reshape(-1, 1)
X_out = to_feature(T_test_out).reshape(-1, 1)

v_sr_in  = sr.predict(X_in)
v_sr_out = sr.predict(X_out)

rmse_sr_in = np.sqrt(mean_squared_error(true_rate(T_test_in),  v_sr_in))
rmse_sr_out = np.sqrt(mean_squared_error(true_rate(T_test_out), v_sr_out))

rows = [
    {"model": res["name"],
     "RMSE_in(0-40)": res["rmse_in"],
     "RMSE_out(40-80)": res["rmse_out"]}
    for res in results_baseline
]
rows.append({"model": "SymbolicRegressor",
             "RMSE_in(0-40)": rmse_sr_in,
             "RMSE_out(40-80)": rmse_sr_out})

pd.DataFrame(rows)

結果は次の通りです：

model	RMSE_in(0–40 ℃)	RMSE_out(40–80 ℃)
LinearRegression	0.340666	22.963070
Poly3+Linear	0.013763	12.652544
RandomForest	0.023305	24.224077
SymbolicRegressor	0.000876	0.050834

学習範囲内でも外挿範囲でも、SymbolicRegressor が圧倒的に小さい誤差になっています。

4.3 グラフで見る Symbolic Regression の挙動

plt.figure(figsize=(10, 4))

# Training range
plt.subplot(1, 2, 1)
plt.scatter(T_train, v_train_obs,
            label="Training data (noisy)", alpha=0.6)
plt.plot(T_test_in, true_rate(T_test_in),
         label="True law", linewidth=2)
plt.plot(T_test_in, v_sr_in,
         label="SymbolicRegressor", linewidth=2)
plt.xlabel("Temperature T (°C)")
plt.ylabel("Reaction rate v(T)")
plt.title("Training range T=0-40 °C")
plt.legend()
plt.grid(True)

# Extrapolation range
plt.subplot(1, 2, 2)
plt.plot(T_test_out, true_rate(T_test_out),
         label="True law", linewidth=2)
plt.plot(T_test_out, v_sr_out,
         label="SymbolicRegressor", linewidth=2)
for res in results_baseline:
    plt.plot(T_test_out, res["v_pred_out"],
             label=res["name"], linestyle="--")
plt.xlabel("Temperature T (°C)")
plt.ylabel("Reaction rate v(T)")
plt.title("Extrapolation range T=40-80 °C")
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

外挿範囲（右図）で

• 青（真の指数則）とオレンジ（SymbolicRegressor）がほぼ完全に重なり、
• 緑（Poly3）、黄（Linear）、紫（RandomForest）はいずれも大きく下に外れている

ことが一目で分かります。

---

5. AI2L 的に見たこのデモのポイント

5.1 ブラックボックスを残さない

• 「真の法則」（Q10 型指数則）を最初に宣言
• ダミーデータ生成コードを公開
• SymbolicRegressor が見つけた式をそのまま表示し、
  どの項が $(T-T_\text{ref})$に対応しているかを人が確認できる

→ モデルの中身を 人間が読める形 で残せています。

5.2 説明責任（アカウンタビリティ）

• 線形／多項式／RandomForest が外挿でどれだけ外れるかを
  RMSE と図で説明できます。
• なぜ Symbolic Regression だけが外挿に強いのかも、
  「指数関数という構造を式として持っているから」と説明できます。

5.3 情報保護

• 完全にシミュレーションされたダミーデータのみを使用
• 実データや研究計画はこの記事には登場しません。

AI2L の観点では、こうした 安全なおもちゃ問題でワークフローを検証する のは
とても大事なステップだと思っています。

5.4 Green AI（省エネ・低コスト）

• 1 次元の小さな問題に対して、過剰に大きなモデルは使っていません。
• 一度「$(v(T) \approx \exp(\beta (T-T_\text{ref}))$」という式が得られれば、
  以後はこの式を評価するだけでよく、重いモデルの再学習は不要です。

---

6. まとめと今後の使い方のイメージ

• Q10 型温度依存則を題材に、
  • ダミーデータ → ベースライン ML → Symbolic Regression
    という流れを AI2L 的に整理しました。
• 結果として、
  • 汎用 ML モデルは学習範囲内ではよく当てはまるものの、外挿範囲では指数の伸びに追いつかない
  • Symbolic Regression は指数関数そのものを再発見し、外挿でも真の法則に追従する
    ことが確認できました。

本記事はあくまで ダミーデータでの道具箱テスト ですが、

• 既知の法則で Symbolic Regression の挙動を確認しておく
• 特徴量設計・関数セット・ペナルティの効き方の感触をつかむ
• ノートブックと図を丸ごと残しておく

といった AI2L 的なスタイルは、
どこかで式探索を試すときの テンプレート として役立つはずです。

---

7. おまけ：Google Colab 用フルコード

最後に、この記事の図と結果を出力した Colab 用のコードをまとめておきます。
（そのまま 1 つのノートブックに貼れば動きます。）

# 依存ライブラリのインストール（Colab / 新規環境向け）
try:
    import gplearn  # noqa: F401
    print("gplearn is already installed.")
except ImportError:
    print("Installing gplearn and related libraries...")
    !pip install -q gplearn scikit-learn matplotlib pandas
    print("Done.")

import numpy as np
import pandas as pd

import matplotlib
import matplotlib.pyplot as plt

from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error

from gplearn.genetic import SymbolicRegressor
from gplearn.functions import make_function

RANDOM_SEED = 42
np.random.seed(RANDOM_SEED)

matplotlib.rcParams["axes.unicode_minus"] = False
%matplotlib inline

# --- 真の Q10 型法則 ---

T_ref = 20.0
Q10 = 2.0
alpha = np.log(Q10) / 10.0

def true_rate(T, alpha=alpha, T_ref=T_ref):
    return np.exp(alpha * (T - T_ref))

def to_feature(T):
    return (T - T_ref) / 20.0

# training and test ranges
T_train = np.linspace(0.0, 40.0, 50)
T_test_in  = np.linspace(0.0, 40.0, 200)
T_test_out = np.linspace(40.0, 80.0, 400)

noise_std = 0.01
v_train_true = true_rate(T_train)
v_train_obs  = v_train_true + np.random.normal(scale=noise_std,
                                               size=T_train.shape)

df_train = pd.DataFrame({
    "T": T_train,
    "v_obs": v_train_obs,
    "v_true": v_train_true
})
df_train.head()

# --- トレーニングデータと真の曲線 ---

plt.figure(figsize=(6, 4))
plt.scatter(df_train["T"], df_train["v_obs"],
            label="Observed (noisy)", alpha=0.7)
plt.plot(df_train["T"], df_train["v_true"],
         label="True law v(T)=exp(alpha*(T-T_ref))", linewidth=2)
plt.xlabel("Temperature T (°C)")
plt.ylabel("Reaction rate v(T)")
plt.title("Training data vs true Q10-type law (T=0-40 °C)")
plt.legend()
plt.grid(True)
plt.show()

# --- ベースラインモデル ---

def evaluate_model(name, model, T_train, v_train,
                   T_test_in, T_test_out):
    X_train = to_feature(T_train).reshape(-1, 1)
    X_in  = to_feature(T_test_in).reshape(-1, 1)
    X_out = to_feature(T_test_out).reshape(-1, 1)

    model.fit(X_train, v_train)

    v_pred_in  = model.predict(X_in)
    v_pred_out = model.predict(X_out)

    rmse_in = np.sqrt(mean_squared_error(true_rate(T_test_in),  v_pred_in))
    rmse_out = np.sqrt(mean_squared_error(true_rate(T_test_out), v_pred_out))

    return {
        "name": name,
        "rmse_in": rmse_in,
        "rmse_out": rmse_out,
        "v_pred_in": v_pred_in,
        "v_pred_out": v_pred_out,
    }

lin_reg = LinearRegression()
poly_reg = Pipeline([
    ("poly", PolynomialFeatures(degree=3)),
    ("lin", LinearRegression())
])
rf_reg = RandomForestRegressor(
    n_estimators=200,
    max_depth=4,
    random_state=RANDOM_SEED
)

results_baseline = []
results_baseline.append(
    evaluate_model("LinearRegression", lin_reg,
                   T_train, v_train_obs, T_test_in, T_test_out)
)
results_baseline.append(
    evaluate_model("Poly3+Linear", poly_reg,
                   T_train, v_train_obs, T_test_in, T_test_out)
)
results_baseline.append(
    evaluate_model("RandomForest", rf_reg,
                   T_train, v_train_obs, T_test_in, T_test_out)
)

pd.DataFrame([
    {"model": res["name"],
     "RMSE_in(0-40)": res["rmse_in"],
     "RMSE_out(40-80)": res["rmse_out"]}
    for res in results_baseline
])

# --- ベースラインの予測曲線 ---

plt.figure(figsize=(10, 4))

plt.subplot(1, 2, 1)
plt.scatter(T_train, v_train_obs,
            label="Training data (noisy)", alpha=0.6)
plt.plot(T_test_in, true_rate(T_test_in),
         label="True law", linewidth=2)
for res in results_baseline:
    plt.plot(T_test_in, res["v_pred_in"], label=res["name"])
plt.xlabel("Temperature T (°C)")
plt.ylabel("Reaction rate v(T)")
plt.title("Training range T=0-40 °C")
plt.legend()
plt.grid(True)

plt.subplot(1, 2, 2)
plt.plot(T_test_out, true_rate(T_test_out),
         label="True law", linewidth=2)
for res in results_baseline:
    plt.plot(T_test_out, res["v_pred_out"], label=res["name"])
plt.xlabel("Temperature T (°C)")
plt.ylabel("Reaction rate v(T)")
plt.title("Extrapolation range T=40-80 °C")
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

# --- Symbolic Regression ---

def _exp(x):
    x = np.clip(x, -10, 10)
    return np.exp(x)

exp = make_function(function=_exp, name="exp", arity=1)
function_set = ("add", "mul", exp)

sr = SymbolicRegressor(
    population_size=5000,
    generations=40,
    function_set=function_set,
    const_range=(-2.0, 2.0),
    p_crossover=0.7,
    p_subtree_mutation=0.1,
    p_hoist_mutation=0.05,
    p_point_mutation=0.1,
    max_samples=0.9,
    parsimony_coefficient=0.0,
    verbose=1,
    random_state=RANDOM_SEED
)

X_train = to_feature(T_train).reshape(-1, 1)
sr.fit(X_train, v_train_obs)

print("Best individual (LISP-style):")
print(sr)

# --- SR の評価と比較 ---

X_in  = to_feature(T_test_in).reshape(-1, 1)
X_out = to_feature(T_test_out).reshape(-1, 1)

v_sr_in  = sr.predict(X_in)
v_sr_out = sr.predict(X_out)

rmse_sr_in = np.sqrt(mean_squared_error(true_rate(T_test_in),  v_sr_in))
rmse_sr_out = np.sqrt(mean_squared_error(true_rate(T_test_out), v_sr_out))

rows = [
    {"model": res["name"],
     "RMSE_in(0-40)": res["rmse_in"],
     "RMSE_out(40-80)": res["rmse_out"]}
    for res in results_baseline
]
rows.append({"model": "SymbolicRegressor",
             "RMSE_in(0-40)": rmse_sr_in,
             "RMSE_out(40-80)": rmse_sr_out})

df_scores = pd.DataFrame(rows)
df_scores

# --- 最終比較図 ---

plt.figure(figsize=(10, 4))

plt.subplot(1, 2, 1)
plt.scatter(T_train, v_train_obs,
            label="Training data (noisy)", alpha=0.6)
plt.plot(T_test_in, true_rate(T_test_in),
         label="True law", linewidth=2)
plt.plot(T_test_in, v_sr_in,
         label="SymbolicRegressor", linewidth=2)
plt.xlabel("Temperature T (°C)")
plt.ylabel("Reaction rate v(T)")
plt.title("Training range T=0-40 °C")
plt.legend()
plt.grid(True)

plt.subplot(1, 2, 2)
plt.plot(T_test_out, true_rate(T_test_out),
         label="True law", linewidth=2)
plt.plot(T_test_out, v_sr_out,
         label="SymbolicRegressor", linewidth=2)
for res in results_baseline:
    plt.plot(T_test_out, res["v_pred_out"],
             label=res["name"], linestyle="--")
plt.xlabel("Temperature T (°C)")
plt.ylabel("Reaction rate v(T)")
plt.title("Extrapolation range T=40-80 °C")
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()