ベイズ最適化の基礎とPython実装

import numpy as np
from scipy.stats import norm

class GaussianProcess:
    def __init__(self, length_scale=1.0, signal_var=1.0, noise_var=1e-6):
        self.l = length_scale
        self.sf2 = signal_var
        self.sn2 = noise_var
        self.X_train = None
        self.y_train = None
        self.K_inv = None

    def kernel(self, X1, X2):
        """RBFカーネル（式3）"""
        dist_sq = np.sum(X1**2, axis=1, keepdims=True) \
                  - 2.0 * X1 @ X2.T \
                  + np.sum(X2**2, axis=1)
        return self.sf2 * np.exp(-0.5 * dist_sq / self.l**2)

    def fit(self, X, y):
        """観測データでGPを学習"""
        self.X_train = X.copy()
        self.y_train = y.copy()
        K = self.kernel(X, X) + self.sn2 * np.eye(len(X))
        self.K_inv = np.linalg.inv(K)

    def predict(self, X_new):
        """新しい点での予測平均と予測分散（式4, 5）"""
        k_star = self.kernel(self.X_train, X_new)
        mu = k_star.T @ self.K_inv @ self.y_train
        var = self.sf2 - np.diag(k_star.T @ self.K_inv @ k_star)
        var = np.maximum(var, 1e-10)
        return mu, var

def expected_improvement(mu, var, y_best, xi=0.01):
    """Expected Improvement（式7, 8）"""
    sigma = np.sqrt(var)
    with np.errstate(divide="ignore", invalid="ignore"):
        Z = (y_best - mu - xi) / sigma
        ei = (y_best - mu - xi) * norm.cdf(Z) + sigma * norm.pdf(Z)
        ei[sigma < 1e-10] = 0.0
    return ei

def upper_confidence_bound(mu, var, kappa=2.0):
    """UCB（最小化用：LCB）（式9）"""
    sigma = np.sqrt(var)
    return mu - kappa * sigma

def probability_of_improvement(mu, var, y_best, xi=0.01):
    """Probability of Improvement（式10）"""
    sigma = np.sqrt(var)
    with np.errstate(divide="ignore", invalid="ignore"):
        Z = (y_best - mu - xi) / sigma
        pi = norm.cdf(Z)
        pi[sigma < 1e-10] = 0.0
    return pi

def bayesian_optimization(f, bounds, n_init=3, n_iter=15,
                          acq_func="ei"):
    """ベイズ最適化のメインループ"""
    # 初期点のランダムサンプリング
    X_init = np.random.uniform(bounds[0], bounds[1],
                               size=(n_init, 1))
    y_init = np.array([f(x) for x in X_init]).reshape(-1)

    gp = GaussianProcess(length_scale=0.5, signal_var=1.0,
                         noise_var=1e-6)
    X_sample = X_init.copy()
    y_sample = y_init.copy()

    # 最適化ループ
    for i in range(n_iter):
        gp.fit(X_sample, y_sample)

        # 候補点での獲得関数の評価
        X_cand = np.linspace(bounds[0], bounds[1], 500).reshape(-1, 1)
        mu, var = gp.predict(X_cand)
        y_best = np.min(y_sample)

        if acq_func == "ei":
            acq = expected_improvement(mu, var, y_best)
            x_next = X_cand[np.argmax(acq)]
        elif acq_func == "ucb":
            acq = upper_confidence_bound(mu, var)
            x_next = X_cand[np.argmin(acq)]
        elif acq_func == "pi":
            acq = probability_of_improvement(mu, var, y_best)
            x_next = X_cand[np.argmax(acq)]

        # 新しい点を評価して追加
        y_next = f(x_next.reshape(1, -1))
        X_sample = np.vstack([X_sample, x_next.reshape(1, -1)])
        y_sample = np.append(y_sample, y_next)

    return X_sample, y_sample, gp

import matplotlib.pyplot as plt

def test_function(x):
    """複数の極小値を持つテスト関数"""
    return np.sin(3 * x) + x**2 * 0.1 - 0.5 * np.cos(7 * x)

np.random.seed(42)
bounds = (-3.0, 3.0)
X_sample, y_sample, gp = bayesian_optimization(
    test_function, bounds, n_init=3, n_iter=12, acq_func="ei"
)

# 最適化結果のプロット
X_plot = np.linspace(bounds[0], bounds[1], 300).reshape(-1, 1)
y_true = np.array([test_function(x) for x in X_plot])
mu, var = gp.predict(X_plot)
sigma = np.sqrt(var)

fig, axes = plt.subplots(2, 1, figsize=(10, 8))

# GP代理モデル
axes[0].plot(X_plot, y_true, "k--", label="True function")
axes[0].plot(X_plot, mu, "b-", label="GP mean")
axes[0].fill_between(X_plot.ravel(), mu - 2*sigma, mu + 2*sigma,
                     alpha=0.2, color="blue", label="95% CI")
axes[0].scatter(X_sample[:3], y_sample[:3],
                c="green", s=80, zorder=5, label="Initial points")
axes[0].scatter(X_sample[3:], y_sample[3:],
                c="red", s=80, zorder=5, label="BO samples")
axes[0].set_xlabel("x")
axes[0].set_ylabel("f(x)")
axes[0].set_title("Bayesian Optimization with GP Surrogate")
axes[0].legend()
axes[0].grid(True)

# EI獲得関数
acq_ei = expected_improvement(mu, var, np.min(y_sample))
axes[1].plot(X_plot, acq_ei, "r-", label="EI")
axes[1].set_xlabel("x")
axes[1].set_ylabel("EI(x)")
axes[1].set_title("Expected Improvement")
axes[1].legend()
axes[1].grid(True)

plt.tight_layout()
plt.savefig("bayesian_optimization_result.png", dpi=150)
plt.show()

best_idx = np.argmin(y_sample)
print(f"Best x: {X_sample[best_idx].item():.4f}")
print(f"Best f(x): {y_sample[best_idx]:.4f}")
print(f"Total evaluations: {len(y_sample)}")