import numpy as np
import matplotlib.pyplot as plt

# 目標分布の定義（正規分布の確率密度関数）
def target_func(x):
    return np.exp(-0.5 * x**2) / np.sqrt(2 * np.pi)

# 対数確率密度の勾配（スコア関数）
def grad_log_target(x):
    return -x  # 正規分布 N(0,1) の対数微分

# ランジュバン・モンテカルロ法
def langevin_monte_carlo(num_samples=10000, num_steps=100, step_size=0.01):
    samples = []
    x = np.random.randn()  # 初期値
    for _ in range(num_samples):
        for _ in range(num_steps):
            grad = grad_log_target(x)
            noise = np.random.randn()
            x += step_size * grad + np.sqrt(2 * step_size) * noise
        samples.append(x)
    return np.array(samples)

# サンプル生成
samples = langevin_monte_carlo()

# 可視化
x_vals = np.linspace(-4, 4, 200)
plt.figure(figsize=(8, 4))
plt.hist(samples, bins=50, density=True, alpha=0.6, label="Langevin Samples")
plt.plot(x_vals, target_func(x_vals), 'r-', label="Target Distribution")
plt.title("Langevin Monte Carlo vs Target Distribution")
plt.legend()
plt.grid(True)
plt.show()