モンテカルロ最適化（CEM/SA/GA/MPPI/PSO）の基礎・比較・Python実装

import numpy as np

def monte_carlo_optimize(
    f, sampler, updater, init_state, n_iter=50, n_samples=100, seed=0
):
    """
    汎用モンテカルロ最適化ループ。

    Parameters
    ----------
    f          : callable  目的関数 f(x) を最小化
    sampler    : callable  state -> shape (n_samples, dim) のサンプル配列
    updater    : callable  (state, samples, scores) -> 新しい state
    init_state : 初期状態（手法ごとに分布パラメータ・個体集団など）
    n_iter     : 反復回数
    n_samples  : 1反復あたりのサンプル数
    """
    rng = np.random.default_rng(seed)
    state = init_state
    history = []
    for t in range(n_iter):
        samples = sampler(state, n_samples, rng)
        scores = np.array([f(x) for x in samples])
        state = updater(state, samples, scores)
        best = scores.min()
        history.append(best)
    return state, np.array(history)


# === 例1: CEMライクな更新（エリート分位点で正規分布をフィット） ===
def cem_sampler(state, n, rng):
    mu, sigma = state
    return rng.normal(mu, sigma, size=(n, mu.shape[0]))


def cem_updater(state, samples, scores, elite_frac=0.2):
    mu, sigma = state
    k = max(1, int(len(scores) * elite_frac))
    elite_idx = np.argsort(scores)[:k]
    elite = samples[elite_idx]
    return elite.mean(axis=0), elite.std(axis=0) + 1e-6


# === 例2: PSOライクな更新（速度更新） ===
def make_pso(dim, n_particles, w=0.7, c1=1.5, c2=1.5):
    def sampler(state, n, rng):
        x, v, pbest, pbest_val, gbest = state
        return x  # PSO は state そのものを評価する

    def updater(state, samples, scores):
        x, v, pbest, pbest_val, gbest = state
        # 個体ベスト更新
        improved = scores < pbest_val
        pbest = np.where(improved[:, None], x, pbest)
        pbest_val = np.where(improved, scores, pbest_val)
        # 大域ベスト更新
        g_idx = pbest_val.argmin()
        gbest = pbest[g_idx]
        # 速度・位置更新
        rng = np.random.default_rng()
        r1, r2 = rng.random(x.shape), rng.random(x.shape)
        v = w * v + c1 * r1 * (pbest - x) + c2 * r2 * (gbest - x)
        x = x + v
        return x, v, pbest, pbest_val, gbest

    return sampler, updater


# === 目的関数: Rastrigin (多峰性ベンチマーク) ===
def rastrigin(x, A=10.0):
    n = x.shape[0]
    return A * n + np.sum(x ** 2 - A * np.cos(2 * np.pi * x))


# === CEM で最適化 ===
dim = 5
init_state = (np.zeros(dim), np.ones(dim) * 2.0)
state, hist = monte_carlo_optimize(
    rastrigin, cem_sampler, cem_updater, init_state,
    n_iter=40, n_samples=200
)
print(f"CEM best: {hist.min():.4f}")

import numpy as np

# --- 共通ベンチマーク: シフト付き Rastrigin (dim=5, 大域最小 x*=shift で 0) ---
dim = 5
shift = np.array([3.0, -2.0, 1.0, -2.5, 4.0])


def rastrigin(x, A=10.0):
    z = x - shift
    n = z.shape[0]
    return A * n + np.sum(z**2 - A * np.cos(2 * np.pi * z))


x0 = np.array([-4.0, 4.0, -4.0, 4.0, -4.0])  # 全手法共通の開始点


# --- 汎用ループ（updater にも rng を渡す拡張版） ---
def monte_carlo_optimize(f, sampler, updater, init_state, n_iter=40, n_samples=200, seed=0):
    rng = np.random.default_rng(seed)
    state = init_state
    history = []
    best_so_far = np.inf
    for t in range(n_iter):
        samples = sampler(state, n_samples, rng)
        scores = np.array([f(x) for x in samples])
        state = updater(state, samples, scores, rng)
        best_so_far = min(best_so_far, scores.min())
        history.append(best_so_far)
    return state, np.array(history)


# --- CEM: エリート分位点で正規分布をフィット ---
def cem_sampler(state, n, rng):
    mu, sigma = state
    return rng.normal(mu, sigma, size=(n, mu.shape[0]))


def cem_updater(state, samples, scores, rng, elite_frac=0.2):
    mu, sigma = state
    k = max(1, int(len(scores) * elite_frac))
    elite = samples[np.argsort(scores)[:k]]
    return elite.mean(axis=0), elite.std(axis=0) + 1e-6


# --- SA: 近傍から n_samples 個の提案を作り、逐次 Metropolis 受理 ---
def sa_sampler(state, n, rng, step=0.6):
    x, fx, T = state
    return x + rng.normal(0, step, size=(n, x.shape[0]))


def sa_updater(state, samples, scores, rng, cooling=0.93):
    x, fx, T = state
    for x_prop, f_prop in zip(samples, scores):
        if f_prop < fx or rng.random() < np.exp(-(f_prop - fx) / T):
            x, fx = x_prop, f_prop
    return x, fx, T * cooling


# --- GA: トーナメント選択 + 一様交叉 + 突然変異、エリート保存1個体 ---
def ga_sampler(state, n, rng):
    return state  # 現世代そのものを評価対象にする


def ga_updater(state, samples, scores, rng, mutation_rate=0.2, mutation_scale=0.5, k=3):
    pop = samples
    n, d = pop.shape
    new_pop = np.empty_like(pop)
    for i in range(n):
        idx1 = rng.integers(0, n, size=k)
        p1 = pop[idx1[np.argmin(scores[idx1])]]
        idx2 = rng.integers(0, n, size=k)
        p2 = pop[idx2[np.argmin(scores[idx2])]]
        # 一様交叉（遺伝子ごとにどちらかの親をコピー。blend交叉は分散が縮小し
        # 多峰性関数で早期収束しやすいため、離散的な遺伝継承を採用）
        cross_mask = rng.random(d) < 0.5
        child = np.where(cross_mask, p1, p2)
        mask = rng.random(d) < mutation_rate
        child = np.where(mask, child + rng.normal(0, mutation_scale, d), child)
        new_pop[i] = child
    new_pop[0] = pop[scores.argmin()]  # エリート保存
    return new_pop


# --- MPPI: 静的最適化への単純化版（軌道=パラメータそのもの） ---
def mppi_sampler(state, n, rng):
    u, sigma, lam = state
    return u + rng.normal(0, sigma, size=(n, u.shape[0]))


def mppi_updater(state, samples, scores, rng, lam=1.0):
    u, sigma, _ = state
    w = np.exp(-(scores - scores.min()) / lam)
    w /= w.sum()
    return (w[:, None] * samples).sum(axis=0), sigma, lam


# --- PSO: 慣性 + 認知 + 社会項による速度更新 ---
def pso_sampler(state, n, rng):
    x, v, pbest, pbest_val, gbest = state
    return x


def pso_updater(state, samples, scores, rng, w=0.7, c1=1.5, c2=1.5):
    x, v, pbest, pbest_val, gbest = state
    improved = scores < pbest_val
    pbest = np.where(improved[:, None], x, pbest)
    pbest_val = np.where(improved, scores, pbest_val)
    gbest = pbest[pbest_val.argmin()]
    r1, r2 = rng.random(x.shape), rng.random(x.shape)
    v = w * v + c1 * r1 * (pbest - x) + c2 * r2 * (gbest - x)
    return x + v, v, pbest, pbest_val, gbest


# --- 初期状態を共通の予算・開始点でそろえる ---
n_iter, n_samples = 40, 200
init_rng = np.random.default_rng(1)  # 初期集団生成専用

cem_init = (x0.copy(), np.ones(dim) * 2.0)
sa_init = (x0.copy(), rastrigin(x0), 10.0)
ga_init = x0 + init_rng.normal(0, 2.0, size=(n_samples, dim))
mppi_init = (x0.copy(), np.ones(dim) * 1.5, 1.0)

pso_x0 = x0 + init_rng.normal(0, 2.0, size=(n_samples, dim))
pso_v0 = np.zeros((n_samples, dim))
pso_pbest_val0 = np.array([rastrigin(xx) for xx in pso_x0])
pso_init = (pso_x0, pso_v0, pso_x0.copy(), pso_pbest_val0, pso_x0[pso_pbest_val0.argmin()])

methods = {
    "CEM": (cem_sampler, cem_updater, cem_init),
    "SA": (sa_sampler, sa_updater, sa_init),
    "GA": (ga_sampler, ga_updater, ga_init),
    "MPPI": (mppi_sampler, mppi_updater, mppi_init),
    "PSO": (pso_sampler, pso_updater, pso_init),
}

results = {}
for name, (sampler, updater, init_state) in methods.items():
    _, hist = monte_carlo_optimize(rastrigin, sampler, updater, init_state, n_iter, n_samples, seed=0)
    results[name] = hist
    print(f"{name:5s}: best@iter40 = {hist[-1]:.4f}")

CEM  : best@iter40 = 6.9648
SA   : best@iter40 = 10.1668
GA   : best@iter40 = 0.5372
MPPI : best@iter40 = 9.9631
PSO  : best@iter40 = 3.0682

目的関数の評価コストは？
├─ 超高コスト ──> ベイズ最適化
└─ 中〜低コスト
    │
    解は連続 or 離散？
    ├─ 離散 ──> SA / GA
    └─ 連続
        │
        リアルタイム性が必要？
        ├─ Yes ──> MPPI
        └─ No
            │
            個体ベストの履歴情報を活かしたい？
            ├─ Yes ──> PSO
            └─ No  ──> CEM