マルコフ連鎖モンテカルロ（MCMC）の基礎：メトロポリス法とギブスサンプリング

import numpy as np
import matplotlib.pyplot as plt

def metropolis_hastings(log_target, initial, n_samples, proposal_std=1.0):
    """メトロポリス・ヘイスティングス法"""
    samples = [initial]
    current = initial
    accepted = 0

    for _ in range(n_samples):
        # 提案分布（ガウスランダムウォーク）
        proposal = current + np.random.normal(0, proposal_std)

        # 対数受容確率
        log_alpha = log_target(proposal) - log_target(current)

        # 受容/棄却
        if np.log(np.random.random()) < log_alpha:
            current = proposal
            accepted += 1

        samples.append(current)

    acceptance_rate = accepted / n_samples
    return np.array(samples), acceptance_rate

# --- 目標分布: 混合ガウス分布 ---
def log_target(x):
    """対数目標密度: 0.3*N(-2,1) + 0.7*N(3,0.5)"""
    from scipy.special import logsumexp
    log_p1 = np.log(0.3) - 0.5 * (x + 2)**2
    log_p2 = np.log(0.7) - (x - 3)**2  # sigma=0.5 → 1/(2*0.25)=2
    return logsumexp([log_p1, log_p2])

# --- 実行 ---
np.random.seed(42)
samples, acc_rate = metropolis_hastings(log_target, initial=0.0,
                                        n_samples=50000, proposal_std=1.5)
print(f"受容率: {acc_rate:.3f}")

# バーンイン期間を除去
burn_in = 5000
samples = samples[burn_in:]

# --- 可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# トレースプロット
axes[0].plot(samples[:2000], alpha=0.5, linewidth=0.5)
axes[0].set_xlabel('Iteration')
axes[0].set_ylabel('Sample value')
axes[0].set_title('Trace Plot')
axes[0].grid(True, alpha=0.3)

# ヒストグラム vs 真の分布
x = np.linspace(-6, 6, 200)
true_pdf = 0.3 * np.exp(-0.5 * (x + 2)**2) / np.sqrt(2*np.pi) + \
           0.7 * np.exp(-0.5 * ((x - 3)/0.5)**2) / (0.5 * np.sqrt(2*np.pi))
axes[1].hist(samples, bins=100, density=True, alpha=0.5, label='MCMC samples')
axes[1].plot(x, true_pdf, 'r-', linewidth=2, label='True density')
axes[1].set_xlabel('x')
axes[1].set_ylabel('Density')
axes[1].set_title('Histogram vs True Distribution')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()