化学者のためのベイズ最適化入門｜手元の実験データから「次に試すべき5条件」をPythonで出す

メーカーで化学の研究開発をしている者です。

実験は1回が高コストです。温度・圧力・触媒量…と振りたい変数が3つあるだけで、組み合わせは数百〜数千。全部試すのは不可能で、結局「勘と経験」で次の条件を決めがちです。

この記事では、過去の実験データ（20件程度）から「次に試すと一番おいしい条件」をAIに提案させるベイズ最適化（Bayesian Optimization, BO）の原理とコア実装を、コピペで動くPythonで示します。数式は最小限。ここまで読めば、BOが「なぜ効くか」と「最小実装」が分かります。

まず、この記事のゴール（出力）を先に見せます。

これまでの最高収率: 87.1%

 priority  temperature  pressure  catalyst_amount  predicted_yield  uncertainty    EI
        1         95.0       6.0              3.0            90.40         2.55  3.403
        2        100.0       6.0              3.0            90.33         2.31  3.301
        3         95.0       5.0              3.0            90.21         2.61  3.245
        4        100.0       6.0              3.5            90.25         1.89  3.173
        5         95.0       6.0              3.5            90.11         2.26  3.096

20回の探索的な実験では最高 87.1% 止まり。BOは**まだ誰も試していない最適域（約95℃・6atm・3.0mol%）**を指し示し、90%超を予測しました。

環境はGoogle Colab（無料）想定。numpy / pandas / scikit-learn / scipy だけで動きます。

ベイズ最適化の直感（30秒）

BOが狙うのは「予測収率が高い場所」だけではありません。狙うのは——

予測は高い × でもまだ誰も試していなくて不確かな場所

です。当たれば大きいからです。この「攻め（活用）と探索（未踏地）のバランス」を自動でとってくれるのがBOの肝。やみくもに振るより3〜5倍効率的に最適条件へ近づけます。

そのために必要な道具は2つだけ。

予測値 μ と不確実性 σ を両方返すモデル → ガウス過程回帰（GP）
「次に試す価値」を数値化する関数 → 獲得関数（ここではEI）

① データを用意する（自社CSVに差し替える場所）

ここはデモ用のダミーデータ生成です。実務では df_exp = pd.read_csv("your_data.csv") に置き換えるだけ。

import numpy as np, pandas as pd

rng = np.random.default_rng(42)

# 真の応答曲面（本来は未知。デモのために用意）。最適は temp95 / press6 / cat3 付近
def true_yield(T, P, C):
    return 90 - 0.010*(T-95)**2 - 0.45*(P-6)**2 - 1.8*(C-3.0)**2

# 初期20実験は探索的に散らばり、最適域はまだ踏めていない設定にする
def near_opt(T, P, C):
    return (abs(T-95) < 13) and (abs(P-6) < 2.2) and (abs(C-3.0) < 1.1)

n0 = 20
Ts, Ps, Cs = [], [], []
while len(Ts) < n0:
    t = rng.uniform(60,120); p = rng.uniform(1,10); c = rng.uniform(0.5,5.0)
    if near_opt(t,p,c):           # 最適域は初期サンプルから除外
        continue
    Ts.append(t); Ps.append(p); Cs.append(c)
T, P, C = np.array(Ts), np.array(Ps), np.array(Cs)
y = true_yield(T,P,C) + rng.normal(0, 1.5, n0)

df_exp = pd.DataFrame({"temperature":T, "pressure":P, "catalyst_amount":C, "yield":y})
print("これまでの最高収率: %.1f%%" % y.max())

② なぜ Random Forest ではなくガウス過程(GP)なのか

ここで多くの人がつまずきます。Random ForestやPyCaretの普通のモデルは「予測値」しか返しません。 でもBOでは「不確実性（自信のなさ）」も知りたい。まだ試していない＝不確実な場所こそ伸びしろがあるからです。

ガウス過程回帰(GP)は、予測値 μ と不確実性 σ の両方を返します。 だからBOの土台にはGPを使います。

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern, ConstantKernel
from sklearn.preprocessing import StandardScaler

feat = ["temperature", "pressure", "catalyst_amount"]
Xtr = df_exp[feat].values
ytr = df_exp["yield"].values

# スケールを揃える（温度60-120と触媒量0.5-5は桁が違うため必須）
scaler = StandardScaler().fit(Xtr)
Xtr_s = scaler.transform(Xtr)

kernel = ConstantKernel(1.0) * Matern(length_scale=1.0, nu=2.5)
gp = GaussianProcessRegressor(kernel=kernel, alpha=1e-2, normalize_y=True,
                              n_restarts_optimizer=5, random_state=42)
gp.fit(Xtr_s, ytr)

最大のつまずきポイント：StandardScaler を必ず噛ませること。スケールを揃えないと距離計算が温度に支配されてGPが壊れます。alpha は測定ノイズの大きさ。

③④ 候補を大量生成して、全部に μ と σ を出す

「振れる範囲」は人間が決めます。各変数を格子状に切って、その組み合わせ全部を候補プールにします。

T_grid = np.linspace(60, 120, 13)
P_grid = np.linspace(1, 10, 10)
C_grid = np.linspace(0.5, 5.0, 10)
TT, PP, CC = np.meshgrid(T_grid, P_grid, C_grid, indexing="ij")
X_cand = np.column_stack([TT.ravel(), PP.ravel(), CC.ravel()])   # 1300候補

mu, sigma = gp.predict(scaler.transform(X_cand), return_std=True)  # μ と σ を同時取得

これで1300個の候補それぞれに「予測収率μ」と「不確実性σ」が付きました。

⑤ 獲得関数 EI（Expected Improvement）— 心臓部

「次に試す価値」をどう数値化するか。最も使われるのが EI（期待改善量）：今の最高値をどれだけ超えられそうかの期待値です。μが高いほど、σが大きいほどEIは大きくなる＝攻めと探索を自動でバランスします。

from scipy.stats import norm

def expected_improvement(mu, sigma, y_best, xi=0.01):
    sigma = np.maximum(sigma, 1e-9)
    z = (mu - y_best - xi) / sigma
    return (mu - y_best - xi) * norm.cdf(z) + sigma * norm.pdf(z)

y_best = ytr.max()
ei = expected_improvement(mu, sigma, y_best)

⑥ EI上位5件を出力（メソッドの完成）

res = pd.DataFrame(X_cand, columns=feat)
res["predicted_yield"] = mu.round(2)
res["uncertainty"]     = sigma.round(2)
res["EI"]              = ei.round(3)
res = res.sort_values("EI", ascending=False).reset_index(drop=True)

top5 = res.head(5).copy()
top5.insert(0, "priority", range(1, 6))
print(top5.to_string(index=False))

実行すると冒頭の表が出ます。初期の最高87.1%に対し、提案5条件はいずれも最適域（真の最適は95℃/6atm/3.0mol%）を正確に指し、90%超を予測。np.random.default_rng(42) 固定なので手元でも同じ値が出ます。

ここまでで、BOの原理（GP→EI→上位提案）とコア実装は完成です。 汎用的なベイズ最適化としては、この記事だけで動かせます。

ここから先 —「実験室で本当に回す」には差が出る

トイデータでは動きました。でも自分の現場でMIを"成果"にするには、メソッドの周りに実務の層が必要です。具体的には：

自社の多変数CSVにそのまま差し替えるテンプレ（列名・前処理・欠損の扱い）
危険な組み合わせ（高温×高圧など）を候補から除外する安全制約の入れ方
提案を実験者・上司に"通す"ための根拠図（過去データと提案を重ねた1枚）
実験 → 結果を追記 → 再提案のループを5〜10サイクル回す運用
xi（探索の強さ）を煮詰まり具合で調整する勘所

このあたりは「汎用のベイズ最適化」では語られない、化学・材料の現場特有のノウハウです。コピペで自社データに差し替えられる完成Colabノートとあわせて、note有料記事 【Ch.7】ベイズ最適化 にまとめています（記事末リンク）。

つまずき・限界（最低限これだけは）

スケーリング必須：StandardScaler を入れないとGPが壊れる
設計空間がすべて：grid の範囲外に最適解があっても絶対に見つからない
出るのは相関であって因果ではない：「なぜ良いか」は人間がメカニズムで考察する

まとめ

BOは「予測が高い × 不確実で伸びしろがある」場所をEIで選び、次の実験を賢く絞る手法
GPを使うのは「予測値μ＋不確実性σ」の両方が要るから
原理と最小実装は GP → EI → 上位N の3ステップ。この記事だけで動かせます

関連（化学者のためのMI教材）

この記事は、現役の化学研究者が書いている学習教材「化学者のためのMI（マテリアルズインフォマティクス）」の汎用メソッド部分を抜粋・再構成したものです。

無料：MIで結局なにができるのか（3つの「そんなことできるの？」体験）→ note記事リンクhttps://note.com/mi_chemist_rino/n/n58b98f8e4440
主力（有料）：【Ch.7】ベイズ最適化 — 自社CSVテンプレ＋安全制約＋提案を通す根拠図＋運用ループ＋完成Colabノートまで（公開1ヶ月は980円→300円、全10章まとめ買いマガジンは2,980円）→ note記事リンクhttps://note.com/mi_chemist_rino/m/m9e7f5a072d09

統計→Python→PyCaret→記述子→本テーマと、Ch.1から順に自分の業務データで回せるよう設計しています。「自分のデータで詰まった」はコメントで歓迎します。