第1章：ベイズ最適化の基礎

1.1 ブラックボックス最適化問題

プロセス産業では、入力と出力の関係が複雑で解析的に記述できない「ブラックボックス」な問題が多く存在します。例えば、化学反応器の最適操作条件を探索する場合、温度・圧力・反応時間などのパラメータと収率の関係は、実験やシミュレーションでしか評価できません。

Example 1: ブラックボックス問題の定式化（化学反応器）

重合反応プロセスを想定し、温度・圧力・触媒濃度の3パラメータを最適化する問題を定義します。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

# ===================================
# Example 1: 化学反応器のブラックボックス最適化問題
# ===================================

class ChemicalReactor:
    """化学反応器のブラックボックスシミュレータ

    パラメータ:
        - 温度 (T): 300-400 K
        - 圧力 (P): 1-5 bar
        - 触媒濃度 (C): 0.1-1.0 mol/L

    目的: 収率 (Yield) を最大化
    """

    def __init__(self, noise_level=0.02):
        self.noise_level = noise_level
        self.bounds = np.array([[300, 400], [1, 5], [0.1, 1.0]])
        self.dim_names = ['Temperature (K)', 'Pressure (bar)', 'Catalyst (mol/L)']
        self.optimal_x = np.array([350, 3.0, 0.5])  # 真の最適解

    def evaluate(self, x):
        """収率を計算（実際の実験/シミュレーションに相当）

        Args:
            x: [温度, 圧力, 触媒濃度]

        Returns:
            yield: 収率 [0-1] + ノイズ
        """
        T, P, C = x

        # 温度依存性（Arrhenius型）
        T_opt = 350
        temp_factor = np.exp(-((T - T_opt) / 30)**2)

        # 圧力依存性（放物線型）
        P_opt = 3.0
        pressure_factor = 1 - 0.3 * ((P - P_opt) / 2)**2

        # 触媒濃度依存性（Langmuir型）
        catalyst_factor = C / (0.2 + C)

        # 相互作用項（温度と圧力の協調効果）
        interaction = 0.1 * np.sin((T - 300) / 50 * np.pi) * (P - 1) / 4

        # 収率計算
        yield_val = 0.85 * temp_factor * pressure_factor * catalyst_factor + interaction

        # ノイズ追加（測定誤差）
        noise = np.random.normal(0, self.noise_level)

        return float(np.clip(yield_val + noise, 0, 1))

    def plot_landscape(self, fixed_catalyst=0.5):
        """目的関数の可視化（触媒濃度固定）"""
        T_range = np.linspace(300, 400, 50)
        P_range = np.linspace(1, 5, 50)
        T_grid, P_grid = np.meshgrid(T_range, P_range)

        Y_grid = np.zeros_like(T_grid)
        for i in range(len(T_range)):
            for j in range(len(P_range)):
                Y_grid[j, i] = self.evaluate([T_grid[j, i], P_grid[j, i], fixed_catalyst])

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

        # 3D surface
        ax1 = fig.add_subplot(121, projection='3d')
        surf = ax1.plot_surface(T_grid, P_grid, Y_grid, cmap=cm.viridis, alpha=0.8)
        ax1.set_xlabel('Temperature (K)')
        ax1.set_ylabel('Pressure (bar)')
        ax1.set_zlabel('Yield')
        ax1.set_title('Chemical Reactor Response Surface')

        # Contour
        ax2 = fig.add_subplot(122)
        contour = ax2.contourf(T_grid, P_grid, Y_grid, levels=20, cmap='viridis')
        ax2.set_xlabel('Temperature (K)')
        ax2.set_ylabel('Pressure (bar)')
        ax2.set_title('Yield Contour (Catalyst=0.5 mol/L)')
        plt.colorbar(contour, ax=ax2, label='Yield')

        plt.tight_layout()
        return fig

# 使用例
reactor = ChemicalReactor(noise_level=0.02)

# 初期条件で評価
x_initial = np.array([320, 2.0, 0.3])
yield_initial = reactor.evaluate(x_initial)
print(f"Initial conditions: T={x_initial[0]}K, P={x_initial[1]}bar, C={x_initial[2]}mol/L")
print(f"Initial yield: {yield_initial:.3f}")

# 最適条件付近で評価
x_optimal = np.array([350, 3.0, 0.5])
yield_optimal = reactor.evaluate(x_optimal)
print(f"\nOptimal conditions: T={x_optimal[0]}K, P={x_optimal[1]}bar, C={x_optimal[2]}mol/L")
print(f"Optimal yield: {yield_optimal:.3f}")

# 可視化
fig = reactor.plot_landscape()
plt.show()

1.2 逐次設計戦略

ベイズ最適化の核心は「逐次設計（Sequential Design）」にあります。前回までの観測結果を活用して、次に評価すべき点を賢く選択することで、効率的な最適化を実現します。

Example 2: 逐次設計 vs ランダムサンプリング

同じ評価回数でも、戦略の違いで最適化性能が大きく変わることを示します。

import numpy as np
import matplotlib.pyplot as plt

# ===================================
# Example 2: 逐次設計戦略のデモンストレーション
# ===================================

def simple_objective(x):
    """1次元テスト関数（解析解あり）"""
    return -(x - 3)**2 * np.sin(5 * x) + 2

class SequentialDesigner:
    """逐次設計による最適化"""

    def __init__(self, objective_func, bounds, n_initial=3):
        self.objective = objective_func
        self.bounds = bounds
        self.X_observed = []
        self.Y_observed = []

        # 初期点（Latin Hypercube Sampling）
        np.random.seed(42)
        for _ in range(n_initial):
            x = np.random.uniform(bounds[0], bounds[1])
            y = objective_func(x)
            self.X_observed.append(x)
            self.Y_observed.append(y)

    def select_next_point(self):
        """次の評価点を選択（簡易版: 探索と活用のバランス）"""
        # 候補点を生成
        candidates = np.linspace(self.bounds[0], self.bounds[1], 100)

        # 既観測点からの距離（探索）
        min_distances = []
        for c in candidates:
            distances = [abs(c - x) for x in self.X_observed]
            min_distances.append(min(distances))

        # 現在の最良値からの期待改善（活用）
        best_y = max(self.Y_observed)
        improvements = [max(0, self.objective(c) - best_y) for c in candidates]

        # スコア計算（探索60% + 活用40%）
        scores = 0.6 * np.array(min_distances) + 0.4 * np.array(improvements)

        return candidates[np.argmax(scores)]

    def optimize(self, n_iterations=10):
        """最適化実行"""
        for i in range(n_iterations):
            x_next = self.select_next_point()
            y_next = self.objective(x_next)
            self.X_observed.append(x_next)
            self.Y_observed.append(y_next)

            current_best = max(self.Y_observed)
            print(f"Iteration {i+1}: x={x_next:.2f}, y={y_next:.3f}, best={current_best:.3f}")

        return self.X_observed, self.Y_observed

# 比較実験
bounds = [0, 5]
n_total = 13  # 初期3点 + 追加10点

# 1. 逐次設計
print("=== Sequential Design ===")
seq_designer = SequentialDesigner(simple_objective, bounds, n_initial=3)
X_seq, Y_seq = seq_designer.optimize(n_iterations=10)

# 2. ランダムサンプリング
print("\n=== Random Sampling ===")
np.random.seed(42)
X_random = np.random.uniform(bounds[0], bounds[1], n_total)
Y_random = [simple_objective(x) for x in X_random]
for i, (x, y) in enumerate(zip(X_random, Y_random), 1):
    current_best = max(Y_random[:i])
    print(f"Sample {i}: x={x:.2f}, y={y:.3f}, best={current_best:.3f}")

# 可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 真の関数
x_true = np.linspace(0, 5, 200)
y_true = simple_objective(x_true)

# 逐次設計
axes[0].plot(x_true, y_true, 'k-', linewidth=2, label='True function', alpha=0.7)
axes[0].scatter(X_seq, Y_seq, c=range(len(X_seq)), cmap='viridis',
                s=100, edgecolor='black', linewidth=1.5, label='Sequential samples', zorder=5)
axes[0].scatter(X_seq[np.argmax(Y_seq)], max(Y_seq), color='red', s=300,
                marker='*', edgecolor='black', linewidth=2, label='Best found', zorder=6)
axes[0].set_xlabel('x')
axes[0].set_ylabel('f(x)')
axes[0].set_title('Sequential Design (Best: {:.3f})'.format(max(Y_seq)))
axes[0].legend()
axes[0].grid(alpha=0.3)

# ランダムサンプリング
axes[1].plot(x_true, y_true, 'k-', linewidth=2, label='True function', alpha=0.7)
axes[1].scatter(X_random, Y_random, c=range(len(X_random)), cmap='plasma',
                s=100, edgecolor='black', linewidth=1.5, label='Random samples', zorder=5)
axes[1].scatter(X_random[np.argmax(Y_random)], max(Y_random), color='red', s=300,
                marker='*', edgecolor='black', linewidth=2, label='Best found', zorder=6)
axes[1].set_xlabel('x')
axes[1].set_ylabel('f(x)')
axes[1].set_title('Random Sampling (Best: {:.3f})'.format(max(Y_random)))
axes[1].legend()
axes[1].grid(alpha=0.3)

plt.tight_layout()
plt.show()

# 収束比較
best_seq = [max(Y_seq[:i+1]) for i in range(len(Y_seq))]
best_random = [max(Y_random[:i+1]) for i in range(len(Y_random))]

plt.figure(figsize=(8, 5))
plt.plot(best_seq, 'o-', linewidth=2, markersize=8, label='Sequential Design')
plt.plot(best_random, 's-', linewidth=2, markersize=8, label='Random Sampling')
plt.xlabel('Number of Evaluations')
plt.ylabel('Best Objective Value')
plt.title('Convergence Comparison')
plt.legend()
plt.grid(alpha=0.3)
plt.show()

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

ベイズ最適化の最も重要な概念が「探索（Exploration）」と「活用（Exploitation）」のバランスです。

• 活用（Exploitation）: 現在の最良点付近を集中的に調査し、局所的な改善を追求
• 探索（Exploration）: 未知の領域を調査し、より良い大域的最適解を発見

Example 3: 探索 vs 活用の可視化

異なるバランスが最適化に与える影響を視覚的に理解します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# ===================================
# Example 3: 探索と活用のトレードオフ
# ===================================

class ExplorationExploitationDemo:
    """探索・活用のバランスを可視化"""

    def __init__(self):
        # サンプル関数: 複数の局所最適解を持つ
        self.x_range = np.linspace(0, 10, 200)
        self.true_func = lambda x: np.sin(x) + 0.3 * np.sin(3*x) + 0.5 * np.cos(5*x)

        # 既観測点（5点）
        self.X_obs = np.array([1.0, 3.0, 4.5, 7.0, 9.0])
        self.Y_obs = self.true_func(self.X_obs) + np.random.normal(0, 0.1, len(self.X_obs))

    def predict_with_uncertainty(self, x):
        """簡易的な予測平均と不確実性（ガウス過程の近似）"""
        # 距離ベースの重み付け平均
        weights = np.exp(-((self.X_obs[:, None] - x) / 1.0)**2)
        weights = weights / (weights.sum(axis=0) + 1e-10)

        mean = (weights.T @ self.Y_obs)

        # 不確実性（既観測点からの距離に依存）
        min_dist = np.min(np.abs(self.X_obs[:, None] - x), axis=0)
        uncertainty = 0.5 * (1 - np.exp(-min_dist / 2.0))

        return mean, uncertainty

    def exploitation_strategy(self):
        """活用戦略: 予測平均が最大の点を選択"""
        mean, _ = self.predict_with_uncertainty(self.x_range)
        return self.x_range[np.argmax(mean)]

    def exploration_strategy(self):
        """探索戦略: 不確実性が最大の点を選択"""
        _, uncertainty = self.predict_with_uncertainty(self.x_range)
        return self.x_range[np.argmax(uncertainty)]

    def balanced_strategy(self, alpha=0.5):
        """バランス戦略: UCB (Upper Confidence Bound)"""
        mean, uncertainty = self.predict_with_uncertainty(self.x_range)
        ucb = mean + alpha * uncertainty
        return self.x_range[np.argmax(ucb)]

    def visualize(self):
        """可視化"""
        mean, uncertainty = self.predict_with_uncertainty(self.x_range)

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

        # (1) 真の関数と予測
        ax = axes[0, 0]
        ax.plot(self.x_range, self.true_func(self.x_range), 'k--',
                linewidth=2, label='True function', alpha=0.7)
        ax.plot(self.x_range, mean, 'b-', linewidth=2, label='Predicted mean')
        ax.fill_between(self.x_range, mean - uncertainty, mean + uncertainty,
                        alpha=0.3, label='Uncertainty (±1σ)')
        ax.scatter(self.X_obs, self.Y_obs, c='red', s=100, zorder=5,
                  edgecolor='black', linewidth=1.5, label='Observations')
        ax.set_xlabel('x')
        ax.set_ylabel('f(x)')
        ax.set_title('Model Prediction with Uncertainty')
        ax.legend()
        ax.grid(alpha=0.3)

        # (2) 活用戦略
        ax = axes[0, 1]
        x_exploit = self.exploitation_strategy()
        ax.plot(self.x_range, mean, 'b-', linewidth=2, label='Predicted mean')
        ax.scatter(self.X_obs, self.Y_obs, c='gray', s=80, zorder=4, alpha=0.5)
        ax.axvline(x_exploit, color='red', linestyle='--', linewidth=2,
                  label=f'Next point (Exploitation)\nx={x_exploit:.2f}')
        ax.set_xlabel('x')
        ax.set_ylabel('Predicted mean')
        ax.set_title('Exploitation Strategy: Maximize Mean')
        ax.legend()
        ax.grid(alpha=0.3)

        # (3) 探索戦略
        ax = axes[1, 0]
        x_explore = self.exploration_strategy()
        ax.plot(self.x_range, uncertainty, 'g-', linewidth=2, label='Uncertainty')
        ax.scatter(self.X_obs, np.zeros_like(self.X_obs), c='gray', s=80,
                  zorder=4, alpha=0.5, label='Observations')
        ax.axvline(x_explore, color='blue', linestyle='--', linewidth=2,
                  label=f'Next point (Exploration)\nx={x_explore:.2f}')
        ax.set_xlabel('x')
        ax.set_ylabel('Uncertainty')
        ax.set_title('Exploration Strategy: Maximize Uncertainty')
        ax.legend()
        ax.grid(alpha=0.3)

        # (4) バランス戦略（UCB）
        ax = axes[1, 1]
        alpha = 1.5
        x_balanced = self.balanced_strategy(alpha=alpha)
        ucb = mean + alpha * uncertainty
        ax.plot(self.x_range, mean, 'b-', linewidth=1.5, label='Mean', alpha=0.7)
        ax.plot(self.x_range, ucb, 'purple', linewidth=2, label=f'UCB (α={alpha})')
        ax.fill_between(self.x_range, mean, ucb, alpha=0.2, color='purple')
        ax.scatter(self.X_obs, self.Y_obs, c='gray', s=80, zorder=4, alpha=0.5)
        ax.axvline(x_balanced, color='purple', linestyle='--', linewidth=2,
                  label=f'Next point (Balanced)\nx={x_balanced:.2f}')
        ax.set_xlabel('x')
        ax.set_ylabel('Value')
        ax.set_title('Balanced Strategy: UCB = Mean + α × Uncertainty')
        ax.legend()
        ax.grid(alpha=0.3)

        plt.tight_layout()
        return fig

# 実行
demo = ExplorationExploitationDemo()

print("=== Exploration vs Exploitation ===")
print(f"Exploitation (最良予測点): x = {demo.exploitation_strategy():.2f}")
print(f"Exploration (最大不確実性): x = {demo.exploration_strategy():.2f}")
print(f"Balanced (UCB, α=0.5): x = {demo.balanced_strategy(alpha=0.5):.2f}")
print(f"Balanced (UCB, α=1.5): x = {demo.balanced_strategy(alpha=1.5):.2f}")

fig = demo.visualize()
plt.show()

1.4 ベイズ最適化の基本ループ

ベイズ最適化のアルゴリズムは、以下の4ステップを反復します：

1. サロゲートモデルの構築: 観測データから目的関数を近似
2. 獲得関数の計算: 次に評価すべき点の有望度を定量化
3. 次点の選択: 獲得関数を最大化する点を選択
4. 評価と更新: 実際に評価し、観測データに追加

Example 4: シンプルなベイズ最適化の実装

最小限の実装で全体の流れを理解します。

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

# ===================================
# Example 4: Simple Bayesian Optimization Loop
# ===================================

class SimpleBayesianOptimization:
    """シンプルなベイズ最適化（1次元）"""

    def __init__(self, objective_func, bounds, noise_std=0.01):
        self.objective = objective_func
        self.bounds = bounds
        self.noise_std = noise_std

        self.X_obs = []
        self.Y_obs = []

        # 初期サンプリング（3点）
        np.random.seed(42)
        for _ in range(3):
            x = np.random.uniform(bounds[0], bounds[1])
            y = self.evaluate(x)
            self.X_obs.append(x)
            self.Y_obs.append(y)

    def evaluate(self, x):
        """目的関数の評価（ノイズ付き）"""
        return self.objective(x) + np.random.normal(0, self.noise_std)

    def gaussian_kernel(self, x1, x2, length_scale=0.5):
        """RBFカーネル"""
        return np.exp(-0.5 * ((x1 - x2) / length_scale)**2)

    def predict(self, x_test):
        """ガウス過程による予測（簡易版）"""
        X_obs_array = np.array(self.X_obs).reshape(-1, 1)
        Y_obs_array = np.array(self.Y_obs).reshape(-1, 1)
        x_test_array = np.array(x_test).reshape(-1, 1)

        # カーネル行列
        K = np.zeros((len(self.X_obs), len(self.X_obs)))
        for i in range(len(self.X_obs)):
            for j in range(len(self.X_obs)):
                K[i, j] = self.gaussian_kernel(self.X_obs[i], self.X_obs[j])

        # ノイズ項追加
        K += self.noise_std**2 * np.eye(len(self.X_obs))

        # テスト点とのカーネル
        k_star = np.array([self.gaussian_kernel(self.X_obs[i], x_test)
                          for i in range(len(self.X_obs))])

        # 予測平均
        K_inv = np.linalg.inv(K)
        mean = k_star.T @ K_inv @ Y_obs_array

        # 予測分散
        k_star_star = self.gaussian_kernel(x_test, x_test)
        variance = k_star_star - k_star.T @ K_inv @ k_star
        std = np.sqrt(np.maximum(variance, 0))

        return mean.flatten(), std.flatten()

    def expected_improvement(self, x):
        """Expected Improvement 獲得関数"""
        mean, std = self.predict(x)
        best_y = max(self.Y_obs)

        # EI計算
        with np.errstate(divide='ignore', invalid='ignore'):
            z = (mean - best_y) / (std + 1e-9)
            ei = (mean - best_y) * norm.cdf(z) + std * norm.pdf(z)
            ei[std == 0] = 0.0

        return -ei  # 最小化問題に変換

    def select_next_point(self):
        """次の評価点を選択（EI最大化）"""
        result = minimize(
            lambda x: self.expected_improvement(x),
            x0=np.random.uniform(self.bounds[0], self.bounds[1]),
            bounds=[self.bounds],
            method='L-BFGS-B'
        )
        return result.x[0]

    def optimize(self, n_iterations=10, verbose=True):
        """最適化実行"""
        for i in range(n_iterations):
            # 次点選択
            x_next = self.select_next_point()
            y_next = self.evaluate(x_next)

            # データ更新
            self.X_obs.append(x_next)
            self.Y_obs.append(y_next)

            current_best = max(self.Y_obs)
            if verbose:
                print(f"Iter {i+1}: x={x_next:.3f}, y={y_next:.3f}, best={current_best:.3f}")

        best_idx = np.argmax(self.Y_obs)
        return self.X_obs[best_idx], self.Y_obs[best_idx]

    def plot_progress(self):
        """最適化の進捗を可視化"""
        x_plot = np.linspace(self.bounds[0], self.bounds[1], 200)
        y_true = [self.objective(x) for x in x_plot]
        mean, std = self.predict(x_plot)

        fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))

        # サロゲートモデル
        ax1.plot(x_plot, y_true, 'k--', linewidth=2, label='True function', alpha=0.7)
        ax1.plot(x_plot, mean, 'b-', linewidth=2, label='GP mean')
        ax1.fill_between(x_plot, mean - 2*std, mean + 2*std, alpha=0.3, label='95% CI')
        ax1.scatter(self.X_obs, self.Y_obs, c='red', s=100, zorder=5,
                   edgecolor='black', linewidth=1.5, label='Observations')
        best_idx = np.argmax(self.Y_obs)
        ax1.scatter(self.X_obs[best_idx], self.Y_obs[best_idx],
                   c='gold', s=300, marker='*', zorder=6,
                   edgecolor='black', linewidth=2, label='Best')
        ax1.set_xlabel('x')
        ax1.set_ylabel('f(x)')
        ax1.set_title('Gaussian Process Surrogate Model')
        ax1.legend()
        ax1.grid(alpha=0.3)

        # 獲得関数
        ei_values = [-self.expected_improvement(x) for x in x_plot]
        ax2.plot(x_plot, ei_values, 'g-', linewidth=2, label='Expected Improvement')
        ax2.fill_between(x_plot, 0, ei_values, alpha=0.3, color='green')
        ax2.axvline(self.X_obs[-1], color='red', linestyle='--',
                   linewidth=2, label=f'Last selected: x={self.X_obs[-1]:.3f}')
        ax2.set_xlabel('x')
        ax2.set_ylabel('EI(x)')
        ax2.set_title('Acquisition Function (Expected Improvement)')
        ax2.legend()
        ax2.grid(alpha=0.3)

        plt.tight_layout()
        return fig

# テスト関数
def test_function(x):
    return -(x - 2.5)**2 * np.sin(10 * x) + 3

# 実行
print("=== Simple Bayesian Optimization ===\n")
bo = SimpleBayesianOptimization(test_function, bounds=[0, 5], noise_std=0.05)
x_best, y_best = bo.optimize(n_iterations=12, verbose=True)

print(f"\n=== Final Result ===")
print(f"Best x: {x_best:.4f}")
print(f"Best y: {y_best:.4f}")

fig = bo.plot_progress()
plt.show()

1.5 比較: BO vs グリッドサーチ vs ランダムサーチ

ベイズ最適化の優位性を定量的に評価します。

Example 5: 3手法の性能比較

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

# ===================================
# Example 5: BO vs Grid Search vs Random Search
# ===================================

# テスト関数（2次元）
def branin_function(x):
    """Branin関数（グローバル最適化のベンチマーク）"""
    x1, x2 = x[0], x[1]
    a, b, c = 1, 5.1/(4*np.pi**2), 5/np.pi
    r, s, t = 6, 10, 1/(8*np.pi)

    term1 = a * (x2 - b*x1**2 + c*x1 - r)**2
    term2 = s * (1 - t) * np.cos(x1)
    term3 = s

    return -(term1 + term2 + term3)  # 最大化問題に変換

class OptimizationComparison:
    """3手法の比較実験"""

    def __init__(self, objective, bounds, budget=30):
        self.objective = objective
        self.bounds = np.array(bounds)  # [[x1_min, x1_max], [x2_min, x2_max]]
        self.budget = budget
        self.dim = len(bounds)

    def grid_search(self):
        """グリッドサーチ"""
        n_per_dim = int(np.ceil(self.budget ** (1/self.dim)))

        grid_1d = [np.linspace(b[0], b[1], n_per_dim) for b in self.bounds]
        grid = np.meshgrid(*grid_1d)

        X_grid = np.column_stack([g.ravel() for g in grid])[:self.budget]
        Y_grid = [self.objective(x) for x in X_grid]

        return X_grid, Y_grid

    def random_search(self, seed=42):
        """ランダムサーチ"""
        np.random.seed(seed)
        X_random = np.random.uniform(
            self.bounds[:, 0], self.bounds[:, 1],
            size=(self.budget, self.dim)
        )
        Y_random = [self.objective(x) for x in X_random]

        return X_random, Y_random

    def bayesian_optimization(self, seed=42):
        """ベイズ最適化（簡易版）"""
        np.random.seed(seed)

        # 初期サンプル（5点）
        n_init = 5
        X = np.random.uniform(self.bounds[:, 0], self.bounds[:, 1],
                             size=(n_init, self.dim))
        Y = [self.objective(x) for x in X]

        # 逐次最適化
        for _ in range(self.budget - n_init):
            # 簡易GP予測
            def gp_mean_std(x_test):
                distances = np.linalg.norm(X - x_test, axis=1)
                weights = np.exp(-distances**2 / 2.0)
                weights = weights / (weights.sum() + 1e-10)

                mean = weights @ Y
                std = 1.0 * np.exp(-np.min(distances) / 1.5)

                return mean, std

            # EI獲得関数
            def neg_ei(x):
                mean, std = gp_mean_std(x)
                best_y = max(Y)
                z = (mean - best_y) / (std + 1e-9)
                ei = (mean - best_y) * norm.cdf(z) + std * norm.pdf(z)
                return -ei

            # 次点選択
            result = minimize(
                neg_ei,
                x0=np.random.uniform(self.bounds[:, 0], self.bounds[:, 1]),
                bounds=self.bounds,
                method='L-BFGS-B'
            )

            x_next = result.x
            y_next = self.objective(x_next)

            X = np.vstack([X, x_next])
            Y.append(y_next)

        return X, Y

    def compare(self, n_trials=5):
        """複数回試行して平均性能を比較"""
        results = {
            'Grid Search': [],
            'Random Search': [],
            'Bayesian Optimization': []
        }

        for trial in range(n_trials):
            print(f"\n=== Trial {trial + 1}/{n_trials} ===")

            # Grid Search
            X_grid, Y_grid = self.grid_search()
            best_grid = [max(Y_grid[:i+1]) for i in range(len(Y_grid))]
            results['Grid Search'].append(best_grid)
            print(f"Grid Search best: {max(Y_grid):.4f}")

            # Random Search
            X_rand, Y_rand = self.random_search(seed=trial)
            best_rand = [max(Y_rand[:i+1]) for i in range(len(Y_rand))]
            results['Random Search'].append(best_rand)
            print(f"Random Search best: {max(Y_rand):.4f}")

            # Bayesian Optimization
            X_bo, Y_bo = self.bayesian_optimization(seed=trial)
            best_bo = [max(Y_bo[:i+1]) for i in range(len(Y_bo))]
            results['Bayesian Optimization'].append(best_bo)
            print(f"Bayesian Optimization best: {max(Y_bo):.4f}")

        return results

    def plot_comparison(self, results):
        """比較結果の可視化"""
        fig, ax = plt.subplots(figsize=(10, 6))

        colors = {'Grid Search': 'blue', 'Random Search': 'orange',
                 'Bayesian Optimization': 'green'}

        for method, trials in results.items():
            trials_array = np.array(trials)
            mean_curve = trials_array.mean(axis=0)
            std_curve = trials_array.std(axis=0)

            x_axis = np.arange(1, len(mean_curve) + 1)
            ax.plot(x_axis, mean_curve, linewidth=2.5, label=method,
                   color=colors[method], marker='o', markersize=4)
            ax.fill_between(x_axis, mean_curve - std_curve, mean_curve + std_curve,
                           alpha=0.2, color=colors[method])

        ax.set_xlabel('Number of Evaluations', fontsize=12)
        ax.set_ylabel('Best Objective Value Found', fontsize=12)
        ax.set_title('Optimization Performance Comparison (Mean ± Std over 5 trials)',
                    fontsize=13, fontweight='bold')
        ax.legend(fontsize=11)
        ax.grid(alpha=0.3)

        return fig

# 実験実行
bounds = [[-5, 10], [0, 15]]  # Branin関数の定義域
comparison = OptimizationComparison(branin_function, bounds, budget=30)

print("Running optimization comparison...")
results = comparison.compare(n_trials=5)

# 最終性能のサマリー
print("\n=== Final Performance Summary ===")
for method, trials in results.items():
    final_values = [trial[-1] for trial in trials]
    print(f"{method:25s}: {np.mean(final_values):.4f} ± {np.std(final_values):.4f}")

fig = comparison.plot_comparison(results)
plt.show()

1.6 収束解析とイテレーション追跡

最適化の進捗を定量的に評価し、収束判定の方法を学びます。

Example 6: 収束診断ツール

import numpy as np
import matplotlib.pyplot as plt

# ===================================
# Example 6: 収束解析とイテレーション追跡
# ===================================

class ConvergenceAnalyzer:
    """ベイズ最適化の収束性を分析"""

    def __init__(self, X_history, Y_history, true_optimum=None):
        self.X_history = np.array(X_history)
        self.Y_history = np.array(Y_history)
        self.true_optimum = true_optimum
        self.n_iter = len(Y_history)

    def compute_metrics(self):
        """収束メトリクスを計算"""
        # 累積最良値
        cumulative_best = [max(self.Y_history[:i+1]) for i in range(self.n_iter)]

        # 改善量（各イテレーションでの向上）
        improvements = [0]
        for i in range(1, self.n_iter):
            improvements.append(max(0, cumulative_best[i] - cumulative_best[i-1]))

        # 最適性ギャップ（真の最適値が既知の場合）
        if self.true_optimum is not None:
            optimality_gap = [self.true_optimum - cb for cb in cumulative_best]
        else:
            optimality_gap = None

        # 収束率（直近5点の改善の標準偏差）
        convergence_rate = []
        window = 5
        for i in range(self.n_iter):
            if i < window:
                convergence_rate.append(np.nan)
            else:
                recent_improvements = improvements[i-window+1:i+1]
                convergence_rate.append(np.std(recent_improvements))

        return {
            'cumulative_best': cumulative_best,
            'improvements': improvements,
            'optimality_gap': optimality_gap,
            'convergence_rate': convergence_rate
        }

    def is_converged(self, tolerance=1e-3, patience=5):
        """収束判定"""
        metrics = self.compute_metrics()
        improvements = metrics['improvements']

        # 直近patience回のイテレーションで改善がtolerance以下
        if len(improvements) < patience:
            return False

        recent_improvements = improvements[-patience:]
        return all(imp < tolerance for imp in recent_improvements)

    def plot_diagnostics(self):
        """診断プロット"""
        metrics = self.compute_metrics()

        fig, axes = plt.subplots(2, 2, figsize=(14, 10))
        iterations = np.arange(1, self.n_iter + 1)

        # (1) 累積最良値の推移
        ax = axes[0, 0]
        ax.plot(iterations, metrics['cumulative_best'], 'b-', linewidth=2, marker='o')
        if self.true_optimum is not None:
            ax.axhline(self.true_optimum, color='red', linestyle='--',
                      linewidth=2, label=f'True optimum: {self.true_optimum:.3f}')
            ax.legend()
        ax.set_xlabel('Iteration')
        ax.set_ylabel('Best Value Found')
        ax.set_title('Convergence: Cumulative Best')
        ax.grid(alpha=0.3)

        # (2) 各イテレーションでの改善量
        ax = axes[0, 1]
        ax.bar(iterations, metrics['improvements'], color='green', alpha=0.7)
        ax.set_xlabel('Iteration')
        ax.set_ylabel('Improvement')
        ax.set_title('Improvement per Iteration')
        ax.grid(alpha=0.3, axis='y')

        # (3) 最適性ギャップ（対数スケール）
        ax = axes[1, 0]
        if metrics['optimality_gap'] is not None:
            gap = np.array(metrics['optimality_gap'])
            gap[gap <= 0] = 1e-10  # 負値を処理
            ax.semilogy(iterations, gap, 'r-', linewidth=2, marker='s')
            ax.set_xlabel('Iteration')
            ax.set_ylabel('Optimality Gap (log scale)')
            ax.set_title('Distance to True Optimum')
            ax.grid(alpha=0.3, which='both')
        else:
            ax.text(0.5, 0.5, 'True optimum unknown',
                   ha='center', va='center', fontsize=14)
            ax.set_xticks([])
            ax.set_yticks([])

        # (4) 収束率（改善の変動性）
        ax = axes[1, 1]
        valid_idx = ~np.isnan(metrics['convergence_rate'])
        ax.plot(iterations[valid_idx], np.array(metrics['convergence_rate'])[valid_idx],
               'purple', linewidth=2, marker='d')
        ax.axhline(1e-3, color='orange', linestyle='--',
                  linewidth=2, label='Convergence threshold')
        ax.set_xlabel('Iteration')
        ax.set_ylabel('Convergence Rate (Std of recent improvements)')
        ax.set_title('Convergence Rate Indicator')
        ax.legend()
        ax.grid(alpha=0.3)

        plt.tight_layout()
        return fig

    def print_summary(self):
        """サマリーレポート"""
        metrics = self.compute_metrics()

        print("=== Convergence Analysis Summary ===")
        print(f"Total iterations: {self.n_iter}")
        print(f"Best value found: {max(self.Y_history):.6f}")
        print(f"Final improvement: {metrics['improvements'][-1]:.6f}")

        if self.true_optimum is not None:
            final_gap = self.true_optimum - max(self.Y_history)
            print(f"True optimum: {self.true_optimum:.6f}")
            print(f"Optimality gap: {final_gap:.6f} ({final_gap/self.true_optimum*100:.2f}%)")

        converged = self.is_converged()
        print(f"Converged: {'Yes' if converged else 'No'}")

        # 最大改善が起きたイテレーション
        max_imp_iter = np.argmax(metrics['improvements'])
        print(f"Largest improvement at iteration: {max_imp_iter + 1} "
              f"(Δy = {metrics['improvements'][max_imp_iter]:.6f})")

# デモ実行（Example 4の結果を使用）
np.random.seed(42)

def test_func(x):
    return -(x - 2.5)**2 * np.sin(10 * x) + 3

# BOを実行してhistoryを取得
from scipy.stats import norm
from scipy.optimize import minimize

X_hist, Y_hist = [], []
bounds = [0, 5]

# 初期サンプリング
for _ in range(3):
    x = np.random.uniform(bounds[0], bounds[1])
    X_hist.append(x)
    Y_hist.append(test_func(x) + np.random.normal(0, 0.02))

# 逐次最適化（15イテレーション）
for iteration in range(15):
    # 簡易GP予測
    def gp_predict(x_test):
        dists = np.abs(np.array(X_hist) - x_test)
        weights = np.exp(-dists**2 / 0.5)
        mean = weights @ Y_hist / (weights.sum() + 1e-10)
        std = 0.5 * (1 - np.exp(-np.min(dists) / 1.0))
        return mean, std

    # EI獲得関数
    def neg_ei(x):
        mean, std = gp_predict(x)
        z = (mean - max(Y_hist)) / (std + 1e-9)
        ei = (mean - max(Y_hist)) * norm.cdf(z) + std * norm.pdf(z)
        return -ei

    # 次点選択
    res = minimize(neg_ei, x0=np.random.uniform(bounds[0], bounds[1]),
                  bounds=[bounds], method='L-BFGS-B')

    x_next = res.x[0]
    y_next = test_func(x_next) + np.random.normal(0, 0.02)

    X_hist.append(x_next)
    Y_hist.append(y_next)

# 収束解析
analyzer = ConvergenceAnalyzer(X_hist, Y_hist, true_optimum=3.62)
analyzer.print_summary()
fig = analyzer.plot_diagnostics()
plt.show()

1.7 ハイパーパラメータチューニング実践例

ベイズ最適化の代表的な応用例として、機械学習モデルのハイパーパラメータチューニングを実装します。

Example 7: Random Forestのハイパーパラメータ最適化

import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import cross_val_score
from sklearn.datasets import make_regression
from scipy.stats import norm
from scipy.optimize import minimize

# ===================================
# Example 7: ハイパーパラメータチューニング
# ===================================

# サンプルデータセット生成（化学プロセスのセンサーデータを想定）
np.random.seed(42)
X_data, y_data = make_regression(n_samples=200, n_features=10, noise=10, random_state=42)

class HyperparameterOptimizer:
    """Random Forestのハイパーパラメータをベイズ最適化"""

    def __init__(self, X_train, y_train):
        self.X_train = X_train
        self.y_train = y_train

        # 最適化する3つのハイパーパラメータ
        self.param_bounds = {
            'n_estimators': [10, 200],      # 決定木の数
            'max_depth': [3, 20],            # 最大深さ
            'min_samples_split': [2, 20]     # 分割最小サンプル数
        }

        self.X_obs = []
        self.Y_obs = []

        # 初期ランダムサンプリング（5点）
        for _ in range(5):
            params = self._sample_random_params()
            score = self._evaluate(params)
            self.X_obs.append(params)
            self.Y_obs.append(score)

    def _sample_random_params(self):
        """ランダムにハイパーパラメータをサンプリング"""
        return [
            np.random.randint(self.param_bounds['n_estimators'][0],
                            self.param_bounds['n_estimators'][1] + 1),
            np.random.randint(self.param_bounds['max_depth'][0],
                            self.param_bounds['max_depth'][1] + 1),
            np.random.randint(self.param_bounds['min_samples_split'][0],
                            self.param_bounds['min_samples_split'][1] + 1)
        ]

    def _evaluate(self, params):
        """ハイパーパラメータの性能評価（5-fold CV）"""
        n_est, max_d, min_split = [int(p) for p in params]

        model = RandomForestRegressor(
            n_estimators=n_est,
            max_depth=max_d,
            min_samples_split=min_split,
            random_state=42,
            n_jobs=-1
        )

        # Cross-validation R² score
        scores = cross_val_score(model, self.X_train, self.y_train,
                                cv=5, scoring='r2')
        return scores.mean()

    def _gp_predict(self, params_test):
        """簡易ガウス過程予測"""
        X_obs_array = np.array(self.X_obs)
        params_array = np.array(params_test).reshape(1, -1)

        # 正規化
        X_obs_norm = (X_obs_array - X_obs_array.mean(axis=0)) / (X_obs_array.std(axis=0) + 1e-10)
        params_norm = (params_array - X_obs_array.mean(axis=0)) / (X_obs_array.std(axis=0) + 1e-10)

        # 距離ベースの重み
        dists = np.linalg.norm(X_obs_norm - params_norm, axis=1)
        weights = np.exp(-dists**2 / 2.0)

        mean = weights @ self.Y_obs / (weights.sum() + 1e-10)
        std = 0.2 * (1 - np.exp(-np.min(dists) / 1.5))

        return mean, std

    def _expected_improvement(self, params):
        """EI獲得関数"""
        mean, std = self._gp_predict(params)
        best_y = max(self.Y_obs)

        z = (mean - best_y) / (std + 1e-9)
        ei = (mean - best_y) * norm.cdf(z) + std * norm.pdf(z)

        return -ei  # 最小化問題

    def optimize(self, n_iterations=15):
        """ベイズ最適化実行"""
        print("=== Hyperparameter Optimization ===\n")

        for i in range(n_iterations):
            # 次のハイパーパラメータを選択
            bounds_array = [
                self.param_bounds['n_estimators'],
                self.param_bounds['max_depth'],
                self.param_bounds['min_samples_split']
            ]

            result = minimize(
                self._expected_improvement,
                x0=self._sample_random_params(),
                bounds=bounds_array,
                method='L-BFGS-B'
            )

            params_next = [int(p) for p in result.x]
            score_next = self._evaluate(params_next)

            self.X_obs.append(params_next)
            self.Y_obs.append(score_next)

            current_best = max(self.Y_obs)
            print(f"Iter {i+1}: n_est={params_next[0]}, max_depth={params_next[1]}, "
                  f"min_split={params_next[2]} → R²={score_next:.4f} (best={current_best:.4f})")

        # 最良パラメータ
        best_idx = np.argmax(self.Y_obs)
        best_params = self.X_obs[best_idx]
        best_score = self.Y_obs[best_idx]

        print(f"\n=== Best Hyperparameters ===")
        print(f"n_estimators: {best_params[0]}")
        print(f"max_depth: {best_params[1]}")
        print(f"min_samples_split: {best_params[2]}")
        print(f"Best R² score: {best_score:.4f}")

        return best_params, best_score

    def plot_optimization_history(self):
        """最適化履歴の可視化"""
        fig, axes = plt.subplots(2, 2, figsize=(14, 10))

        iterations = np.arange(1, len(self.Y_obs) + 1)
        cumulative_best = [max(self.Y_obs[:i+1]) for i in range(len(self.Y_obs))]

        # (1) R² scoreの推移
        ax = axes[0, 0]
        ax.plot(iterations, self.Y_obs, 'o-', linewidth=2, markersize=8,
               label='Observed R²', alpha=0.7)
        ax.plot(iterations, cumulative_best, 's-', linewidth=2.5, markersize=8,
               color='red', label='Best R²')
        ax.set_xlabel('Iteration')
        ax.set_ylabel('R² Score')
        ax.set_title('Optimization Progress')
        ax.legend()
        ax.grid(alpha=0.3)

        # (2) n_estimatorsの探索軌跡
        ax = axes[0, 1]
        n_estimators = [x[0] for x in self.X_obs]
        scatter = ax.scatter(iterations, n_estimators, c=self.Y_obs,
                           cmap='viridis', s=100, edgecolor='black', linewidth=1.5)
        ax.set_xlabel('Iteration')
        ax.set_ylabel('n_estimators')
        ax.set_title('Parameter Exploration: n_estimators')
        plt.colorbar(scatter, ax=ax, label='R² Score')
        ax.grid(alpha=0.3)

        # (3) max_depthの探索軌跡
        ax = axes[1, 0]
        max_depth = [x[1] for x in self.X_obs]
        scatter = ax.scatter(iterations, max_depth, c=self.Y_obs,
                           cmap='viridis', s=100, edgecolor='black', linewidth=1.5)
        ax.set_xlabel('Iteration')
        ax.set_ylabel('max_depth')
        ax.set_title('Parameter Exploration: max_depth')
        plt.colorbar(scatter, ax=ax, label='R² Score')
        ax.grid(alpha=0.3)

        # (4) パラメータ空間の探索（2D投影: n_estimators vs max_depth）
        ax = axes[1, 1]
        scatter = ax.scatter(n_estimators, max_depth, c=self.Y_obs,
                           s=150, cmap='viridis', edgecolor='black', linewidth=1.5)
        best_idx = np.argmax(self.Y_obs)
        ax.scatter(n_estimators[best_idx], max_depth[best_idx],
                  s=500, marker='*', color='red', edgecolor='black', linewidth=2,
                  label='Best', zorder=5)
        ax.set_xlabel('n_estimators')
        ax.set_ylabel('max_depth')
        ax.set_title('Parameter Space Exploration')
        plt.colorbar(scatter, ax=ax, label='R² Score')
        ax.legend()
        ax.grid(alpha=0.3)

        plt.tight_layout()
        return fig

# 実行
optimizer = HyperparameterOptimizer(X_data, y_data)
best_params, best_score = optimizer.optimize(n_iterations=15)

fig = optimizer.plot_optimization_history()
plt.show()

# ベースラインとの比較
print("\n=== Comparison with Default Parameters ===")
default_model = RandomForestRegressor(random_state=42, n_jobs=-1)
default_score = cross_val_score(default_model, X_data, y_data, cv=5, scoring='r2').mean()
print(f"Default R² score: {default_score:.4f}")
print(f"Optimized R² score: {best_score:.4f}")
print(f"Improvement: {(best_score - default_score) / default_score * 100:.2f}%")

学習目標の確認

この章を完了すると、以下を説明・実装できるようになります：

基本理解

• ✅ ブラックボックス最適化問題の特徴と課題を説明できる
• ✅ 逐次設計戦略の利点をランダムサーチと比較して述べられる
• ✅ 探索と活用のトレードオフの概念を理解している

実践スキル

• ✅ 化学プロセスのブラックボックス問題を定式化できる
• ✅ シンプルなベイズ最適化ループを実装できる
• ✅ 3手法（BO/Grid/Random）の性能を比較評価できる
• ✅ 収束診断ツールを使って最適化の進捗を分析できる

応用力

• ✅ 機械学習モデルのハイパーパラメータチューニングに応用できる
• ✅ 実務問題に対して適切な最適化手法を選択できる
• ✅ 最適化結果の信頼性を評価できる

次のステップ

第1章では、ベイズ最適化の基本概念と実装を学びました。次章では、ベイズ最適化の中核技術である「ガウス過程」について詳しく学びます。