1. はじめに

連続空間でのベイズ最適化を理解するために、マーケティング担当者が直面する最適化問題であるデジタル広告の入札単価設定を取り上げる。入札単価を上げれば他の入札に勝ちやすくなり獲得数は増えるが、1件あたりの粗利は減少する。このトレードオフのなかで純利益を最大化する単価をどう見つけるか、という問題である。

ナイーブなアプローチとして、グリッドサーチや A/B テストがある。しかしどちらも「まず均等に試し、後から判断する」という受動的な戦略であり、観測コスト（実際に入札する）が高い場面では非効率である。

ベイズ最適化（Bayesian Optimization, BO）は、これまでの観測から確率モデルを構築し、「次にどこを試すべきか」を能動的に決定する適応的な戦略である。本稿では、ガウス過程（Gaussian Process, GP）を用いて入札単価と純利益の関係をモデル化する。次点選択戦略として、獲得関数を用いる GP-UCB とサンプリングベースの GP-TS の2手法を実装し、探索と活用のトレードオフを可視化する。

2. データ：広告入札シミュレーター

実際の広告配信データを使う代わりに、市場反応を模したシミュレーターを構築する。

入力変数：入札単価 $x \in [100, 2000]$（JPY）

目的変数：純利益（Net Profit）

$$f(x) = \underbrace{\frac{L}{1 + e^{-k(x - x_0)}}}_{\text{コンバージョン数}} \times (V - x)$$

コンバージョン数はロジスティック曲線でモデル化し、1件あたりの粗利は $(V - x)$（価値 $V$ から入札単価を差し引いたもの）とする。

パラメータ	値	意味
$L$	20	最大コンバージョン数（上限）
$k$	0.01	ロジスティック曲線の傾き
$x_0$	1000 JPY	中点（感度が最も高い入札単価）
$V$	3000 JPY	1コンバージョンあたりの価値
$\sigma_{\text{noise}}$	3000 JPY	観測ノイズ（日次変動を模倣）

観測値には日次変動を模したガウスノイズ $\varepsilon \sim \mathcal{N}(0, \sigma_{\text{noise}}^2)$ を加える。

$$y = f(x) + \varepsilon$$

この関数は $x \approx 1274$ JPY 付近に内点最適解を持ち、それより高い単価では粗利の減少が獲得増加を上回って純利益が下落する。真の最大純利益は約 32,427 JPY である。

3. 手法：ガウス過程によるベイズ最適化

3.1 ガウス過程（GP）

GP は関数全体に対する確率分布であり、任意の有限点集合上の同時分布がガウス分布になるという性質を持つ。カーネル関数 $K(x, x')$ が関数の滑らかさを規定する。

本稿では Matérn 3/2 カーネルを用いる。（gaussianカーネルも試したが、ロジスティック曲線の急変化を捉えきれなかったため）

$$K(x, x') = \left(1 + \frac{\sqrt{3},r}{l}\right)\exp!\left(-\frac{\sqrt{3},r}{l}\right), \qquad r = |\tilde{x} - \tilde{x}'|$$

ここで $\tilde{x} = (x - x_{\min}) / (x_{\max} - x_{\min})$ は入札単価を $[0,1]$ に正規化したもの、$l$ は長さスケール（length scale）である。RBF（ガウス）カーネルは無限回微分可能で過度に滑らかな関数を仮定するが、Matérn 3/2 は1回微分可能な関数クラスを対象とし、局所的な急変化を捉えやすい。入札単価と純利益のような実務的な関係は一般に RBF ほど滑らかではないため、Matérn 3/2 の方が適している。

観測データ $\mathbf{X}_{\text{obs}} = {x_1, \ldots, x_n}$ 、

観測値 $\mathbf{r} = {r_1, \ldots, r_n}$ が与えられたとき、未観測点 $X_*$ における GP 事後分布は以下で与えられる。

$$\mu_* = K_*^\top (K + \sigma_n^2 I)^{-1} \mathbf{r}$$

$$\Sigma_{*} = K_{**} - K_{*}^\top (K + \sigma_n^2 I)^{-1} K_*$$

ここで $K$ は観測点間のカーネル行列、$K_*$ は観測点と予測点のクロスカーネル、$K_{**}$ は予測点間のカーネル行列、$\sigma_n^2$ は GP の正則化パラメータ（観測ノイズの分散に対応）である。シミュレーターの観測ノイズ $\sigma_{\text{noise}} = 3000$ JPY はデータ生成用の設定値であり、$\sigma_n^2$ は GP 回帰においてモデルが推定する別のパラメータである点に注意する。

GP は予測平均 $\mu_{*}$ （活用の根拠）と予測分散 $\text{diag}(\Sigma_*)$（不確実性の指標）を同時に提供する。

3.2 GP-UCB（Upper Confidence Bound）

UCB 獲得関数は平均と不確実性を明示的に組み合わせる。

$$a_{\text{UCB}}(x) = \mu_{*}(x) + \alpha \ \sigma_{*}(x)$$

ここで $\sigma_{*}(x) = \sqrt{\Sigma_{*xx}}$ は点 $x$ における GP 予測標準偏差である。$\alpha$ は探索（不確実性重視）と活用（平均重視）のバランスを制御するハイパーパラメータであり（本稿では $\alpha = 2.0$、すなわち予測平均＋2σ の上限を採用）、各ステップで $a_{\text{UCB}}(x)$ を最大化する入札単価を次のサンプル点として選択する [1]。

3.3 GP-TS（Thompson Sampling）

トンプソンサンプリングは、確率的な探索戦略である。GP 事後分布から関数をまるごとサンプリングし、そのサンプル関数の最大点を選択する。

$$\hat{f} \sim \mathcal{N}(\mu_{*}, \Sigma_{*})$$

$$x_{\text{next}} = \arg\max_x \hat{f}(x)$$

UCB が決定論的であるのに対し、TS は同じ観測から異なる点を選ぶことがあり、より多様な探索パターンを示す [2]。また、各サンプリングは独立に実行できるため、複数点を同時に選ぶバッチ BO に自然に拡張できる（詳細は §5 を参照）。

3.4 ハイパーパラメータの MAP 推定

$l$（長さスケール）と $\sigma_n^2$（GP 正則化パラメータ）は事前には未知であり、観測データから推定する。観測ごとに対数事後分布（log posterior）を最大化する MAP 推定を行う [3]。

対数周辺尤度（MLE 目的関数）は以下で与えられる。ここで $\mathbf{r}n = \mathbf{r} / Y{scale}}$ は正規化済み観測値、$K_n$ は観測点間のカーネル行列である。

$$\log p(\mathbf{r} \mid \mathbf{X}, l, \sigma_n^2) = -\frac{1}{2}\mathbf{r}_n^\top (K_n + \sigma_n^2 I)^{-1}\mathbf{r}_n - \frac{1}{2}\log|K_n + \sigma_n^2 I| - \frac{n}{2}\log 2\pi$$

第1項はデータへの当てはまり、第2項はモデル複雑さへのペナルティである。観測数が少ない BO の初期段階では最尤推定（MLE）だと $l \to 0$ に縮退するリスクがあるため、対数正規事前分布 $\pi(l, \sigma_n^2)$ を加えた MAP 推定を採用する。

$$\hat{l}, \hat{\sigma}n^2 = \arg\max{l,, \sigma_n^2} \left[ \log p(\mathbf{r} \mid \mathbf{X}, l, \sigma_n^2) + \log \pi(l, \sigma_n^2) \right]$$

$$\log \pi(l, \sigma_n^2) = -\frac{(\log l - \log l_0)^2}{2\tau_l^2} - \frac{(\log \sigma_n^2 - \log \sigma_{n,0}^2)^2}{2\tau_\sigma^2}$$

初期値は「100 JPY 離れるとほぼ無相関」というドメイン知識から $l_0 = 100$ JPY とし、最適化は対数空間で L-BFGS-B を用いる。

3.5 スケール正規化の重要性

y の正規化スケール $Y_{\text{scale}}$ は、初期3点の観測値の標準偏差から自動設定する。

$$Y_{\text{scale}} = \max!\left(\text{std}(r_{\text{init}}),\ 1000\right) \text{ JPY}$$

GP の内部計算は正規化値 $r_n = r / Y_{\text{scale}}$ で行う。これにより GP 予測平均 $\tilde{\mu}$ と予測標準偏差 $\tilde{\sigma}_{\ast}$ がともに O(1) に収まる。UCB 獲得関数を元の JPY スケールに戻すと

$$a_{\text{UCB}}(x) = \tilde{\mu}(x) \cdot Y_{\text{scale}} + \alpha, \tilde{\sigma}\ast(x) \cdot Y{\text{scale}}$$

となり、探索ボーナス $\alpha \tilde{\sigma}\ast Y{\text{scale}} \approx \alpha \cdot Y_{\text{scale}}$ と予測平均 $\tilde{\mu} Y_{\text{scale}}$ が同一スケールで競合できる。$Y_{\text{scale}}$ を固定値（例: 15,000 JPY）ではなく初期観測の分散から導くことで、ドメインに依存しない自動設定が実現する。

4. 結果

初期観測を探索空間の等間隔3点（575 / 1,050 / 1,525 JPY）でサンプリングし、10 ステップの BO ループを実行した（乱数シード 42）。

初期観測の標準偏差から $Y_{\text{scale}} = 12{,}353$ JPY が自動設定された。

4.1 GP の事後分布の更新

各ステップで可視化した GP 予測分布（平均 ± 2σ）から、以下の挙動が確認された。

初期ステップ：不確実性が大きく（±2σ 幅が広い）、エージェントは真の関数から離れた領域も積極的に探索する。
ステップ進行：観測点が増えるにつれて ±2σ の幅が縮小し、GP 予測平均が真の関数に近づく。特に観測点の密度が高い領域での不確実性の減少が顕著である。
収束：後半ステップでは最適付近への集中が見られ、「活用」フェーズへの移行が確認できる。

GP-UCGの結果

GP-TSの結果

4.2 2手法の探索挙動の比較

探索した入札単価の推移（Step 1〜10）を比較すると、2手法に異なる挙動が見られた。

GP-UCB：UCB 値が最大の点を決定論的に選ぶため、各ステップの選択が安定している。不確実性の高い領域を体系的に消去しながら最適点に近づく。
GP-TS：事後分布からのサンプリングに確率性があるため、ステップ間の選択にばらつきが生じる。同じ事後分布からでも異なる点を選ぶことがあり、探索の多様性が高い。

観測純利益の推移では、両手法とも序盤は低報酬点も観測しながら（探索）、後半に向けて高報酬点への集中（活用）が確認された。

5. 考察

探索と活用のトレードオフ

BO の本質的な優位性は、観測コストが高いときに際立つ。本シミュレーションで 13 点の観測から真の最適値の 98% 以上の純利益を達成したことは、グリッドサーチ（200 点評価）や A/B テスト（均等配分）と比較して観測効率が大幅に高いことを示している。

手法	観測数（目安）	アプローチ
グリッドサーチ	200 点	均等に試す
A/B テスト	数十〜数百点	事前設定した候補のみ
ベイズ最適化	13 点で高精度	不確実な領域を優先して試す

GP-UCB vs GP-TS の使い分け

2手法の本質的な差は2点に集約される。

観点	GP-UCB	GP-TS
選択の決定性	決定論的（同じ事後分布なら同じ点を選ぶ）	確率的（毎回異なる点を選びうる）
探索制御	$\alpha$ で明示的に調整	事後分布の不確実性に自動追従

GP-UCB は $\alpha$ を小さくすることで活用を重視し、失敗コストが高い場面でのリスクを抑えられる。また選択が決定論的であるため、「なぜこの単価を試したのか」という根拠をレポートで説明しやすい。

GP-TS は探索の強度をチューニングするパラメータを持たず、不確実性の大きさに応じた探索が自動的に行われる。確率的な選択により多峰性を持つ目的関数でも局所最適に陥りにくい。

実務への適用

この枠組みは実際の広告配信システムにも適用可能である。

初期の数キャンペーンで複数の入札単価を試し、GPの事前観測を構築する。
各日（またはキャンペーン単位）を1ステップとして観測を追加し、GPを更新する。
次の入札単価を UCB または TS で決定する。

観測コスト（実際の広告費）が高いほど、BO による効率的な探索は経済的合理性を持つ。

今後の課題：カーネル関数の改善

本稿では Matérn 3/2 カーネルを採用したが、13 点の観測では真の関数形（ロジスティック曲線×線形）を完全には再現できていない。より適切なカーネルの候補として、スペクトラルミクスチャーカーネル（複数の周波数成分を重ね合わせて複雑な形状を表現）や、目的関数の構造知識を組み込んだカスタムカーネルが挙げられる。観測数が増えるにつれてカーネル選択の影響は大きくなるため、実運用では複数のカーネルを周辺尤度で比較することが推奨される。

6. まとめ

本稿では、広告入札単価の最適化問題をベイズ最適化で解くシミュレーションを実装した。

ガウス過程を代理モデルとすることで、関数全体の不確実性を確率分布として表現し、探索と活用を統一的な枠組みで扱うことができる。
長さスケール $l$ と GP 正則化パラメータ $\sigma_n^2$ を MAP 推定で観測ごとに更新することで、ドメイン知識が足りない場合でも、ドメインに合わせたモデルの構築が可能となる。
GP-UCB と GP-TS の両手法とも、わずか 13 点の観測から真の最大純利益の 98% 以上を達成した。A/B テストやグリッドサーチが「均等に試す」受動的戦略であるのに対し、BO は「不確実な領域を優先して試す」能動的戦略であり、観測コストが高い実務環境での優位性は大きい。

参考文献

[1] Srinivas, N., Krause, A., Kakade, S. M., & Seeger, M. (2010). Gaussian Process Optimization in the Bandit Setting: No Regret and Experimental Design. ICML. — GP-UCB 獲得関数と累積後悔の理論的上界を提案した原論文。
[2] Thompson, W. R. (1933). On the likelihood that one unknown probability exceeds another in view of the evidence of two samples. Biometrika, 25(3/4), 285–294. — Thompson Sampling の原論文。GP-TS はこの手法を連続ガウス過程に拡張したもの。
[3] Shahriari, B., Swersky, K., Wang, Z., Adams, R. P., & de Freitas, N. (2016). Taking the Human Out of the Loop: A Review of Bayesian Optimization. Proceedings of the IEEE, 104(1), 148–175. — BO の包括的サーベイ。ハイパーパラメータ推定（MAP/MLE）・獲得関数・応用事例を網羅的にレビュー。
[4] 飯塚修平 (2020). ウェブ最適化ではじめる機械学習 . オライリー・ジャパン. ISBN: 978-4-87311-916-8. — 本稿の実装はこの書籍の構成に沿ったシミュレーション設計を参考にしている。

ソースコード

pyproject.toml

[project]
name = "web-optimization"
version = "0.1.0"
requires-python = ">=3.11"
dependencies = [
    "numpy",
    "matplotlib",
    "scipy",
]

[dependency-groups]
dev = [
    "ipykernel>=7.1.0",
]

main.py

import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import minimize

plt.style.use("grayscale")
np.random.seed(42)

# ---------------------------------------------------------------------------
# パラメータ設定
# ---------------------------------------------------------------------------
X_MIN, X_MAX = 100, 2000
N_GRID = 200
L = 20.0
K_SLOPE = 0.01
X0 = 1000
V = 3000.0
SIGMA_NOISE = 3000.0

L_SCALE_JPY = 100.0
S_INIT = 0.1
ALPHA = 2.0
N_INIT = 3
N_ITER = 10

X_STAR = np.linspace(X_MIN, X_MAX, N_GRID)

# ---------------------------------------------------------------------------
# シミュレーター
# ---------------------------------------------------------------------------

def true_net_profit(x):
    conversions = L / (1 + np.exp(-K_SLOPE * (x - X0)))
    return conversions * (V - x)


def observe(x):
    return true_net_profit(x) + np.random.normal(0, SIGMA_NOISE)

# ---------------------------------------------------------------------------
# カーネル / GP
# ---------------------------------------------------------------------------

def _norm_x(x):
    return (np.asarray(x) - X_MIN) / (X_MAX - X_MIN)


def _l_norm(l_jpy):
    return l_jpy / (X_MAX - X_MIN)


def matern32_kernel(x1, x2, l):
    r = np.abs(x1 - x2)
    t = np.sqrt(3) * r / l
    return (1.0 + t) * np.exp(-t)


def build_K(X, l):
    n = len(X)
    K = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            K[i, j] = matern32_kernel(X[i], X[j], l)
    return K


def build_K_star(X, X_star, l):
    K_s = np.zeros((len(X), len(X_star)))
    for i, xi in enumerate(X):
        for j, xj in enumerate(X_star):
            K_s[i, j] = matern32_kernel(xi, xj, l)
    return K_s


def gp_posterior(X_obs, r_obs, X_star, l_norm, s, y_scale):
    xn_obs = _norm_x(X_obs)
    rn_obs = np.asarray(r_obs) / y_scale
    xn_star = _norm_x(X_star)

    K = build_K(xn_obs, l_norm)
    K_star = build_K_star(xn_obs, xn_star, l_norm)
    K_ss = build_K(xn_star, l_norm)
    A = np.linalg.inv(K + s * np.eye(len(xn_obs)))

    mu = (K_star.T @ A @ rn_obs) * y_scale
    sigma = (K_ss - K_star.T @ A @ K_star) * y_scale**2
    return mu, sigma

# ---------------------------------------------------------------------------
# MAP ハイパーパラメータ推定
# ---------------------------------------------------------------------------

def _log_marginal_likelihood(log_l_norm, log_s, x_norm, r_norm):
    l = np.exp(log_l_norm)
    s = np.exp(log_s)
    n = len(x_norm)
    K = build_K(x_norm, l)
    try:
        L_chol = np.linalg.cholesky(K + s * np.eye(n))
        alpha = np.linalg.solve(L_chol.T, np.linalg.solve(L_chol, r_norm))
        log_det = 2 * np.sum(np.log(np.diag(L_chol)))
        return -0.5 * (r_norm @ alpha + log_det + n * np.log(2 * np.pi))
    except np.linalg.LinAlgError:
        return -1e10


def _log_prior(log_l_norm, log_s):
    l_norm_init = _l_norm(L_SCALE_JPY)
    lp_l = -0.5 * ((log_l_norm - np.log(l_norm_init)) / 1.0) ** 2
    lp_s = -0.5 * ((log_s - np.log(S_INIT)) / 1.0) ** 2
    return lp_l + lp_s


def optimize_hyperparams(X_obs, r_obs, l_jpy_current, s_current, y_scale):
    x_norm = _norm_x(X_obs)
    r_norm = np.asarray(r_obs) / y_scale

    def neg_map(params):
        log_l, log_s = params
        return -(
            _log_marginal_likelihood(log_l, log_s, x_norm, r_norm)
            + _log_prior(log_l, log_s)
        )

    res = minimize(
        neg_map,
        x0=[np.log(_l_norm(l_jpy_current)), np.log(s_current)],
        method="L-BFGS-B",
        options={"maxiter": 100},
    )
    return np.exp(res.x[0]) * (X_MAX - X_MIN), np.exp(res.x[1])

# ---------------------------------------------------------------------------
# エージェント
# ---------------------------------------------------------------------------

class GPUCBAgent:
    def __init__(self, X_star, y_scale, alpha=ALPHA):
        self.X_star = X_star
        self.y_scale = y_scale
        self.alpha = alpha
        self.xs, self.rs = [], []
        self.l_jpy = L_SCALE_JPY
        self.s = S_INIT
        self.mu = np.zeros(len(X_star))
        self.sigma = build_K(_norm_x(X_star), _l_norm(self.l_jpy)) * y_scale**2

    def get_arm(self):
        ucb = self.mu + self.alpha * np.sqrt(np.maximum(np.diag(self.sigma), 0))
        return self.X_star[np.argmax(ucb)]

    def sample(self, x, r):
        self.xs.append(x)
        self.rs.append(r)
        self.l_jpy, self.s = optimize_hyperparams(
            np.array(self.xs), np.array(self.rs), self.l_jpy, self.s, self.y_scale
        )
        self.mu, self.sigma = gp_posterior(
            np.array(self.xs), np.array(self.rs), self.X_star,
            _l_norm(self.l_jpy), self.s, self.y_scale,
        )


class GPTSAgent:
    def __init__(self, X_star, y_scale):
        self.X_star = X_star
        self.y_scale = y_scale
        self.xs, self.rs = [], []
        self.l_jpy = L_SCALE_JPY
        self.s = S_INIT
        self.mu = np.zeros(len(X_star))
        self.sigma = build_K(_norm_x(X_star), _l_norm(self.l_jpy)) * y_scale**2

    def get_arm(self):
        sigma_safe = self.sigma + 1e-6 * np.eye(len(self.sigma))
        f_sample = np.random.multivariate_normal(self.mu, sigma_safe)
        return self.X_star[np.argmax(f_sample)]

    def sample(self, x, r):
        self.xs.append(x)
        self.rs.append(r)
        self.l_jpy, self.s = optimize_hyperparams(
            np.array(self.xs), np.array(self.rs), self.l_jpy, self.s, self.y_scale
        )
        self.mu, self.sigma = gp_posterior(
            np.array(self.xs), np.array(self.rs), self.X_star,
            _l_norm(self.l_jpy), self.s, self.y_scale,
        )

# ---------------------------------------------------------------------------
# 可視化
# ---------------------------------------------------------------------------

def plot_gp_1d(agent, step, x_next, method_name, ax=None):
    if ax is None:
        fig, ax = plt.subplots(figsize=(8, 4))

    std = np.sqrt(np.maximum(np.diag(agent.sigma), 0))

    ax.plot(X_STAR, true_net_profit(X_STAR), linestyle="dotted", color="black", label="True Net Profit")
    ax.plot(X_STAR, agent.mu, color="black", label="GP mean")
    ax.fill_between(X_STAR, agent.mu - 2 * std, agent.mu + 2 * std, alpha=0.3, color="gray", label="±2σ")
    if agent.xs:
        ax.scatter(agent.xs, agent.rs, marker="x", color="black", zorder=5, label="Observations")
    ax.axvline(x_next, linestyle="--", color="black", linewidth=0.8)
    ax.scatter([x_next], [true_net_profit(x_next)], marker="*", s=150, color="black", zorder=6, label="Next sample")

    ax.set_xlabel("Bid Price (JPY)")
    ax.set_ylabel("Net Profit (JPY)")
    ax.set_title(f"Step {step}: {method_name}")
    ax.legend(fontsize=7, loc="upper left")
    ax.set_xlim(X_MIN, X_MAX)

# ---------------------------------------------------------------------------
# メイン
# ---------------------------------------------------------------------------

def main():
    # 初期観測
    x_init = np.linspace(X_MIN, X_MAX, N_INIT + 2)[1:-1]
    r_init = np.array([observe(x) for x in x_init])
    y_scale = float(max(np.std(r_init), 1000.0))
    print(f"x_init : {x_init.astype(int)}")
    print(f"r_init : {r_init.astype(int)}")
    print(f"Y_SCALE: {y_scale:.0f} JPY")

    n_cols, n_rows = 5, N_ITER // 5

    # GP-UCB
    ucb_agent = GPUCBAgent(X_STAR, y_scale)
    for x, r in zip(x_init, r_init):
        ucb_agent.sample(x, r)

    fig, axes = plt.subplots(n_rows, n_cols, figsize=(20, n_rows * 4))
    for step, ax in enumerate(axes.flatten()):
        x_next = ucb_agent.get_arm()
        plot_gp_1d(ucb_agent, step + 1, x_next, "GP-UCB", ax=ax)
        ucb_agent.sample(x_next, observe(x_next))
    plt.tight_layout()
    plt.savefig("gp_ucb.png", dpi=100)
    plt.close()

    best_idx = np.argmax(ucb_agent.mu)
    print(f"GP-UCB optimal bid: {X_STAR[best_idx]:.0f} JPY  "
          f"(net profit: {ucb_agent.mu[best_idx]:.0f} JPY)")

    # GP-TS
    np.random.seed(42)
    ts_agent = GPTSAgent(X_STAR, y_scale)
    for x, r in zip(x_init, r_init):
        ts_agent.sample(x, r)

    fig, axes = plt.subplots(n_rows, n_cols, figsize=(20, n_rows * 4))
    for step, ax in enumerate(axes.flatten()):
        x_next = ts_agent.get_arm()
        plot_gp_1d(ts_agent, step + 1, x_next, "GP-TS", ax=ax)
        ts_agent.sample(x_next, observe(x_next))
    plt.tight_layout()
    plt.savefig("gp_ts.png", dpi=100)
    plt.close()

    best_idx = np.argmax(ts_agent.mu)
    print(f"GP-TS  optimal bid: {X_STAR[best_idx]:.0f} JPY  "
          f"(net profit: {ts_agent.mu[best_idx]:.0f} JPY)")

    # 比較プロット
    steps = np.arange(1, N_ITER + 1)
    fig, axes = plt.subplots(1, 2, figsize=(12, 4))

    axes[0].plot(steps, ucb_agent.xs[N_INIT:], marker="o", color="black", label="GP-UCB")
    axes[0].plot(steps, ts_agent.xs[N_INIT:], marker="x", linestyle="--", color="gray", label="GP-TS")
    axes[0].axhline(X_STAR[np.argmax(true_net_profit(X_STAR))], linestyle="dotted", color="black", label="True optimum")
    axes[0].set_xlabel("Step")
    axes[0].set_ylabel("Bid Price (JPY)")
    axes[0].set_title("Sampled Bid Price per Step")
    axes[0].legend()

    axes[1].plot(steps, ucb_agent.rs[N_INIT:], marker="o", color="black", label="GP-UCB")
    axes[1].plot(steps, ts_agent.rs[N_INIT:], marker="x", linestyle="--", color="gray", label="GP-TS")
    axes[1].axhline(true_net_profit(X_STAR[np.argmax(true_net_profit(X_STAR))]),
                    linestyle="dotted", color="black", label="True max")
    axes[1].set_xlabel("Step")
    axes[1].set_ylabel("Observed Net Profit (JPY)")
    axes[1].set_title("Observed Net Profit per Step")
    axes[1].legend()

    plt.tight_layout()
    plt.savefig("comparison.png", dpi=100)
    plt.close()
    print("Saved: gp_ucb.png, gp_ts.png, comparison.png")


if __name__ == "__main__":
    main()