ベイズ最適化で「次の実験」を提案する — ダミーデータで体験する実験計画の自動化

導入

材料開発・触媒探索・プロセス条件探索・医工学系の実験条件最適化など、研究現場には「1回の実験（測定）が高い」問題がたくさんあります。

• 試料作製に時間がかかる
• 測定が高価／予約が必要
• 条件数が多く、全探索が現実的でない
• ノイズや個体差で結果がぶれる

こういうとき、闇雲に条件を振るのではなく、少ない実験回数で良い条件にたどり着きたい。
そのための代表的な手法がベイズ最適化（Bayesian Optimization, BO）です。

本記事では、Google Colab で動く最小コードを用意し、

• ダミーデータ（合成の真の関数 + ノイズ付き測定）で「高コスト実験」を模擬
• ランダム探索とベイズ最適化を比較
• ベイズ最適化が「次にどの条件を試すべきか」をどう提案するかを可視化

して、直感と理解を一致させます。

このデモは教育用です。現実の材料・化学・生体現象を忠実に再現するものではありません。
ただし「高コスト実験を少回数で最適化する」という方法論は、多くの○○インフォマティクス領域に共通します。

---

TL;DR

• ベイズ最適化は、少ない試行で「良い条件」を見つけるために、
  (1) 代理モデル（例：ガウス過程） と (2) 次の候補を決めるルール（獲得関数） を組み合わせる。
• 今回は 2 変数（温度 T と濃度 C）の条件探索をダミーで作り、
  ランダム探索 vs ベイズ最適化（Expected Improvement） を比較する。
• Colabでコードを回すと、図（探索点・ベスト値推移・獲得関数）が出る。
  予測で終わらず「次の実験提案」まで示した。

---

1. なぜ「次の実験提案」が重要か（○○インフォマティクス的な意義）

多くの○○インフォマティクス（マテリアルズ／ケモ／ナノ／プロセス等）で共通するのは、

• データ収集（実験・測定）にコストがかかる
• 全条件を試すのは無理
• でも意思決定（次に何を試すか）は毎回必要

という構造です。

ベイズ最適化は、単に「関係を当てる」だけでなく、
次に測るべき点（＝次の意思決定）を提案することで、研究の進み方を変えます。

---

2. ベイズ最適化についての簡単な説明

ベイズ最適化の基本は、次の2つです。

2.1 代理モデル（surrogate model）

本当の目的関数 f(x)（例：収率、性能、強度、毒性の低さ etc.）は未知で、測るのが高い。
そこで、すでに測ったデータから 「それっぽい関数」を学習して近似します。

代表例がガウス過程回帰（Gaussian Process Regression）です。

• 「この辺は高そう」「この辺は不確か」を同時に出せる（平均と不確実性）

2.2 獲得関数（acquisition function）

次に試す x_next を決めるルールです。
今回は Expected Improvement (EI) を使います。

EI(x) = E[\max(0, f(x) - f_{best} - \\xi)]

• 今までの最高値 f_best をどれくらい更新できそうか（期待値）
• 高そう（活用） と 不確か（探索） のバランスを自然に取る

---

3. 今回のデモ設定（ダミーデータ）

• 探索する実験条件：2変数
  • 温度 T：40–100 ℃
  • 濃度 C：0–1（単位は仮）
• 「真の関数」は非線形で、良い点が局所的に存在するように合成（物理モデルではない）
• 測定はノイズあり（現実の測定ブレを模擬）

---

4. Colabで動かす（コピペでOK）

以下のセルを 上から順に Colab に貼り付けて実行してください。
実行後、fig/ フォルダに画像が保存されます。

ConvergenceWarning が出ることがありますが、これは「ガウス過程のハイパーパラメータ最適化が境界に近い」等の警告です。
このデモでは 無視してOK（動作自体は継続します）。

---

セル1：準備（ライブラリ・設定）

import os
import numpy as np
import matplotlib.pyplot as plt

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern, WhiteKernel, ConstantKernel
from sklearn.preprocessing import StandardScaler
from scipy.stats import norm
from scipy.stats import qmc

# 再現性
SEED = 42
rng_global = np.random.default_rng(SEED)

# 探索空間 [T, C]
bounds = np.array([
    [40.0, 100.0],  # Temperature T (°C)
    [0.0, 1.0],     # Concentration C (arb.)
])

# 実験回数
N_INIT = 8       # 初期点（空間を広く見るための“最初の実験”）
N_ITER = 20      # 追加で回す提案回数（合計 N_INIT + N_ITER）
NOISE_SIGMA = 1.5

# 獲得関数を評価する候補点数（大きいほど提案は良くなりやすいが重くなる）
N_CANDIDATES = 8000

os.makedirs("fig", exist_ok=True)

---

セル2：真の関数（ダミー）と「高コスト実験」を定義

def _raw_true(T, C):
    """滑らかな2D関数：主ピーク + 局所ピーク + ゆるい尾根 + 小さな揺らぎ
    ※物理モデルではありません（デモ用の地形）。
    """
    main_peak  = 1.4*np.exp(-((T-82.0)/9.0)**2  - ((C-0.30)/0.10)**2)
    local_peak = 0.7*np.exp(-((T-60.0)/12.0)**2 - ((C-0.75)/0.08)**2)
    ridge      = 0.25*np.exp(-((T-95.0)/18.0)**2) * (1.0 - 0.8*np.abs(C-0.55))
    wiggle     = 0.05*np.sin(T/7.5) * np.cos(8.0*C)
    return main_peak + local_peak + ridge + wiggle

def build_true_function(bounds):
    # 0〜100にスケーリングして見やすくする（※現実の収率の意味は持たない）
    T = np.linspace(bounds[0,0], bounds[0,1], 501)
    C = np.linspace(bounds[1,0], bounds[1,1], 501)
    TT, CC = np.meshgrid(T, C)
    raw = _raw_true(TT, CC)
    mn, mx = raw.min(), raw.max()

    def f(X):
        X = np.asarray(X)
        T_, C_ = X[:,0], X[:,1]
        y = 100.0 * (_raw_true(T_, C_) - mn) / (mx - mn)
        return y

    return f

f_true = build_true_function(bounds)

def expensive_experiment(X, noise_sigma, rng):
    """ノイズ付き測定を返す（収率は0〜100にクリップ）"""
    y = f_true(np.asarray(X))
    y_noisy = y + rng.normal(0.0, noise_sigma, size=y.shape)
    return np.clip(y_noisy, 0.0, 100.0)

---

セル3：ベイズ最適化（ガウス過程 + Expected Improvement）

def fit_gp(X_train, y_train, seed):
    # 入力を標準化（GPのハイパーパラメータが安定しやすい）
    scaler = StandardScaler()
    Xs = scaler.fit_transform(X_train)

    # カーネル：定数 × Matern + ノイズ
    kernel = ConstantKernel(1.0, (1e-2, 1e3)) * Matern(
        length_scale=np.ones(X_train.shape[1]),
        length_scale_bounds=(1e-2, 1e2),
        nu=2.5,
    ) + WhiteKernel(noise_level=1.0, noise_level_bounds=(1e-5, 1e1))

    gp = GaussianProcessRegressor(
        kernel=kernel,
        normalize_y=True,
        n_restarts_optimizer=3,
        random_state=seed,
    )
    gp.fit(Xs, y_train)
    return gp, scaler

def expected_improvement(X_cand, X_train, y_train, gp, scaler, xi=0.01):
    Xc = scaler.transform(X_cand)
    mu, sigma = gp.predict(Xc, return_std=True)

    y_best = np.max(y_train)
    imp = mu - y_best - xi
    Z = imp / (sigma + 1e-12)
    ei = imp * norm.cdf(Z) + sigma * norm.pdf(Z)
    ei[sigma < 1e-12] = 0.0
    return ei

def propose_next_point(bounds, gp, scaler, X_train, y_train, rng, n_candidates):
    X_cand = rng.uniform(bounds[:,0], bounds[:,1], size=(n_candidates, bounds.shape[0]))
    ei = expected_improvement(X_cand, X_train, y_train, gp, scaler)
    return X_cand[np.argmax(ei)].reshape(1, -1)

---

セル4：ランダム探索 vs ベイズ最適化を実行して図を保存

def sobol_initial_points(bounds, n_points, seed):
    # 初期点は“空間をまんべんなく”見ると安定しやすい（Sobol）
    sampler = qmc.Sobol(d=bounds.shape[0], scramble=True, seed=seed)
    X01 = sampler.random_base2(m=int(np.log2(n_points)))  # n_pointsは2のべき乗が前提
    return qmc.scale(X01, bounds[:,0], bounds[:,1])

def run_random_search(bounds, X_init, y_init, n_iter, rng):
    X, y = X_init.copy(), y_init.copy()
    for _ in range(n_iter):
        x_new = rng.uniform(bounds[:,0], bounds[:,1], size=(1, bounds.shape[0]))
        y_new = expensive_experiment(x_new, noise_sigma=NOISE_SIGMA, rng=rng)
        X = np.vstack([X, x_new])
        y = np.concatenate([y, y_new])
    return X, y

def run_bayes_opt(bounds, X_init, y_init, n_iter, rng, gp_seed):
    X, y = X_init.copy(), y_init.copy()
    for _ in range(n_iter):
        gp, scaler = fit_gp(X, y, seed=gp_seed)
        x_new = propose_next_point(bounds, gp, scaler, X, y, rng=rng, n_candidates=N_CANDIDATES)
        y_new = expensive_experiment(x_new, noise_sigma=NOISE_SIGMA, rng=rng)
        X = np.vstack([X, x_new])
        y = np.concatenate([y, y_new])
    return X, y

# 初期点（N_INIT=8 なので2^3）
X_init = sobol_initial_points(bounds, N_INIT, seed=SEED)
y_init = expensive_experiment(X_init, noise_sigma=NOISE_SIGMA, rng=rng_global)

# 比較実行（初期点は同じにしてフェアに）
rng_rs = np.random.default_rng(SEED + 1)
rng_bo = np.random.default_rng(SEED + 2)

X_rs, y_rs = run_random_search(bounds, X_init, y_init, n_iter=N_ITER, rng=rng_rs)
X_bo, y_bo = run_bayes_opt(bounds, X_init, y_init, n_iter=N_ITER, rng=rng_bo, gp_seed=SEED)

best_rs = np.maximum.accumulate(y_rs)
best_bo = np.maximum.accumulate(y_bo)

print("Random search best:", float(best_rs[-1]))
print("BayesOpt best     :", float(best_bo[-1]))

# ---- plotting helpers ----
def eval_on_grid(bounds, f, n=300):
    T = np.linspace(bounds[0,0], bounds[0,1], n)
    C = np.linspace(bounds[1,0], bounds[1,1], n)
    TT, CC = np.meshgrid(T, C)
    Xg = np.column_stack([TT.ravel(), CC.ravel()])
    Z = f(Xg).reshape(TT.shape)
    return TT, CC, Xg, Z

def plot_surface(title, savepath):
    TT, CC, _, Z = eval_on_grid(bounds, f_true)
    plt.figure(figsize=(7,4.8))
    cf = plt.contourf(TT, CC, Z, levels=30)
    plt.colorbar(cf, label="Yield (0-100, synthetic)")
    plt.xlabel("Temperature T (°C)")
    plt.ylabel("Concentration C (arb.)")
    plt.title(title)
    plt.tight_layout()
    plt.savefig(savepath, dpi=150)
    plt.show()

def plot_points(X, title, savepath):
    TT, CC, _, Z = eval_on_grid(bounds, f_true)
    plt.figure(figsize=(7,4.8))
    cf = plt.contourf(TT, CC, Z, levels=30)
    plt.colorbar(cf, label="Yield (0-100, synthetic)")
    plt.scatter(X[:,0], X[:,1], s=26, edgecolor="k", linewidth=0.4)
    plt.xlabel("Temperature T (°C)")
    plt.ylabel("Concentration C (arb.)")
    plt.title(title)
    plt.tight_layout()
    plt.savefig(savepath, dpi=150)
    plt.show()

def plot_best(best_rs, best_bo, savepath):
    n = len(best_rs)
    x = np.arange(1, n+1)
    plt.figure(figsize=(7,4.5))
    plt.plot(x, best_rs, label="Random search")
    plt.plot(x, best_bo, label="Bayesian optimization (EI)")
    plt.xlabel("Number of experiments")
    plt.ylabel("Best observed yield so far")
    plt.title("Sample efficiency comparison")
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.savefig(savepath, dpi=150)
    plt.show()

plot_surface("Ground truth (only for this toy demo)", "fig/fig1_ground_truth.png")
plot_points(X_rs, "Random search: sampled points", "fig/fig2_random_points.png")
plot_points(X_bo, "Bayesian optimization: sampled points", "fig/fig3_bo_points.png")
plot_best(best_rs, best_bo, "fig/fig4_best_so_far.png")

print("Saved figures:")
for fn in ["fig1_ground_truth.png","fig2_random_points.png","fig3_bo_points.png","fig4_best_so_far.png"]:
    print(" -", os.path.join("fig", fn))

---

セル5：なぜそこを選んだ？（代理モデル平均・不確実性・EIを可視化）

def plot_gp_maps(X_train, y_train, save_prefix):
    gp, scaler = fit_gp(X_train, y_train, seed=SEED)

    TT, CC, Xg, _ = eval_on_grid(bounds, f_true, n=250)
    Xg_s = scaler.transform(Xg)
    mu, sigma = gp.predict(Xg_s, return_std=True)
    mu = mu.reshape(TT.shape)
    sigma = sigma.reshape(TT.shape)

    # GP mean
    plt.figure(figsize=(7,4.8))
    cf = plt.contourf(TT, CC, mu, levels=30)
    plt.colorbar(cf, label="GP mean (predicted yield)")
    plt.scatter(X_train[:,0], X_train[:,1], s=24, edgecolor="k", linewidth=0.4)
    plt.xlabel("Temperature T (°C)")
    plt.ylabel("Concentration C (arb.)")
    plt.title("Surrogate model (GP) mean")
    plt.tight_layout()
    plt.savefig(f"{save_prefix}_gp_mean.png", dpi=150)
    plt.show()

    # GP std
    plt.figure(figsize=(7,4.8))
    cf = plt.contourf(TT, CC, sigma, levels=30)
    plt.colorbar(cf, label="GP std (uncertainty)")
    plt.scatter(X_train[:,0], X_train[:,1], s=24, edgecolor="k", linewidth=0.4)
    plt.xlabel("Temperature T (°C)")
    plt.ylabel("Concentration C (arb.)")
    plt.title("Surrogate uncertainty (GP std)")
    plt.tight_layout()
    plt.savefig(f"{save_prefix}_gp_std.png", dpi=150)
    plt.show()

    # EI
    ei = expected_improvement(Xg, X_train, y_train, gp, scaler).reshape(TT.shape)
    plt.figure(figsize=(7,4.8))
    cf = plt.contourf(TT, CC, ei, levels=30)
    plt.colorbar(cf, label="Expected Improvement")
    plt.scatter(X_train[:,0], X_train[:,1], s=24, edgecolor="k", linewidth=0.4)
    plt.xlabel("Temperature T (°C)")
    plt.ylabel("Concentration C (arb.)")
    plt.title("Acquisition function (EI): where to try next?")
    plt.tight_layout()
    plt.savefig(f"{save_prefix}_ei.png", dpi=150)
    plt.show()

plot_gp_maps(X_bo, y_bo, "fig/fig5")

print("Saved extra figures:")
for fn in ["fig5_gp_mean.png","fig5_gp_std.png","fig5_ei.png"]:
    print(" -", os.path.join("fig", fn))

---

セル6：画像をまとめてダウンロード（Colab用）

import shutil
shutil.make_archive("figures", "zip", "fig")
print("Created figures.zip")

（Colab左のファイルツリーから figures.zip を右クリック → Download）

---

5. 本記事掲載コードを実行して得られる図

図1：真の関数

図2：ランダム探索が試した点

闇雲に散らばる（＝当たりを引くまで運）。

図3：ベイズ最適化が試した点

序盤は広く見て、途中から“良さそうな近辺”に寄る動きが出ます。

図4：ベスト値の推移（重要）

「同じ実験回数で、どっちが早く良い条件に到達するか」を示します。

図5：GP平均・不確実性・EI（仕組みの理解に効く）

• 平均：どこが高そうか
• 不確実性：どこが未知か
• EI：更新が期待できる場所（探索×活用）
  
  
  

上記の図5（GP mean / GP std / EI）は「この時点までに得られた実験データ」から計算した地図です。
今回の実行ではすでにベスト値が上限（100）付近まで到達しているため、EI（期待改善量）は全体に小さく（暗く）見えます。これは「これ以上の大幅な改善が見込みにくい＝収束に近い」ことを意味します。

---

6. 実務・研究で使うときの注意点（ここが“インフォマティクス”）

ベイズ最適化は便利ですが、現実では次の設計が重要です。

6.1 探索空間（範囲）の決め方が8割

• 物理的に不可能な条件は 最初から除外（安全・設備制約）
• 「広すぎる探索空間」は、少試行では学習できず失敗しやすい
  → まずは“妥当な範囲”に絞り、必要なら段階的に広げる

6.2 ノイズが大きいときは「再測定（リピート）」も設計する

• 1回の測定がぶれるなら、最良候補の近辺で 複数回測って平均を取る戦略が効く
• 目的は「運が良かった点」ではなく「再現する良い点」を見つけること

6.3 目的が複数あるなら「多目的最適化」

• 性能を上げたいがコストも下げたい、安全性も守りたい
  → パレート最適・制約付き最適化の枠組みが必要

6.4 カテゴリ変数が多いと難しくなる

• 触媒種類、材料系、製法の種類など
  → 連続値前提のGPだけでは扱いにくいことがある（工夫が必要）

---

7. 発展課題

• 実験コストが条件によって違う（高温は高い等）設定を入れて「コスト付きBO」をやる
• 目的を2つにして（例：性能↑・コスト↓）パレートフロントを可視化する
• 制約（危険領域）を入れて「安全制約付きBO」を作る
• 3〜6次元に拡張し、N_CANDIDATES と収束の関係を観察する
• ノイズを増やし、再測定（リピート）戦略の効果を比較する

---

まとめ

本記事では、ダミーデータで「高コスト実験」を模擬し、
ベイズ最適化が“次に試すべき実験条件”を提案できることを、Colabで再現しました。

○○インフォマティクスの共通テーマは、モデル精度だけでなく、

• データ収集（実験計画）
• 評価設計（何をもって良いとするか）
• 運用（次の意思決定に落とす）

までを、再現可能な手順でつなぐことです。

まずはこの最小例を動かし、あなたの分野の「次の一手」に置き換えてみてください。