Metropolis 算法详解：从原理到 Python 实现

import numpy as np
import matplotlib.pyplot as plt

# --- 1. 参数设置 ---
# 目标分布 pi (岛A: 0.2, 岛B: 0.5, 岛C: 0.3)
pi = np.array([0.2, 0.5, 0.3])
states_map = {0: 'Island A (0.2)', 1: 'Island B (0.5)', 2: 'Island C (0.3)'}

# 提议矩阵 Q (对称，均匀跳跃)
# Q[i][j] 代表从 i 提议去 j 的概率
Q = np.array([
    [1/3, 1/3, 1/3],
    [1/3, 1/3, 1/3],
    [1/3, 1/3, 1/3]
])

def discrete_metropolis(n_iter, initial_state):
    samples = []
    current_state = initial_state
    
    # 预计算 Q 的累积分布，用于 U1 的分位数查找
    Q_cumsum = np.cumsum(Q, axis=1)
    
    for t in range(n_iter):
        # --- 第一步：生成提议 (U1) ---
        u1 = np.random.uniform(0, 1)
        # 根据 Q 的当前行累积概率确定候选状态 j
        # 对应数学公式：sum(q_ik) < u1 <= sum(q_ik)
        proposed_state = np.searchsorted(Q_cumsum[current_state], u1) # 逆变换采样 (Inverse Transform Sampling) 的离散实现
        
        # --- 第二步：计算接受率 alpha ---
        # 这里的 pi_j / pi_i
        alpha = min(1, pi[proposed_state] / pi[current_state])
        
        # --- 第三步：判定是否转移 (U2) ---
        u2 = np.random.uniform(0, 1)
        if u2 <= alpha:
            # 接受提议
            current_state = proposed_state
        else:
            # 拒绝提议，状态保持不变
            pass
            
        samples.append(current_state)
        
    return np.array(samples)

# --- 2. 运行模拟 ---
N = 5000  # 迭代次数
samples = discrete_metropolis(N, initial_state=0)

# --- 3. 可视化展示 ---
plt.figure(figsize=(15, 10))

# 图 1: 轨迹图 (Trace Plot) - 只看前 100 步
# 这让我们看到具体的“跳跃”动作
plt.subplot(2, 2, 1)
plt.step(range(100), samples[:100], where='mid', color='blue', linewidth=1.5)
plt.yticks([0, 1, 2], [states_map[0], states_map[1], states_map[2]])
plt.title("1. Trace Plot (First 100 Steps)\nWatch the jumper move!", fontsize=12)
plt.grid(True, axis='y', linestyle='--', alpha=0.5)
plt.xlabel("Step")

# 图 2: 收敛图 (Convergence of Frequencies)
# 展示频率是如何随时间逼近真实概率的
plt.subplot(2, 2, 2)
# 计算动态频率
iterations = np.arange(1, N + 1)
prob_0 = np.cumsum(samples == 0) / iterations
prob_1 = np.cumsum(samples == 1) / iterations
prob_2 = np.cumsum(samples == 2) / iterations

plt.plot(prob_0, label='Simulated A', color='red', alpha=0.6)
plt.plot(prob_1, label='Simulated B', color='green', alpha=0.6)
plt.plot(prob_2, label='Simulated C', color='blue', alpha=0.6)
# 画出理论横线
plt.axhline(pi[0], color='red', linestyle='--', label='True A (0.2)')
plt.axhline(pi[1], color='green', linestyle='--', label='True B (0.5)')
plt.axhline(pi[2], color='blue', linestyle='--', label='True C (0.3)')

plt.title("2. Convergence Plot\nLaw of Large Numbers in action", fontsize=12)
plt.xlabel("Iterations")
plt.ylabel("Estimated Probability")
plt.legend()
plt.grid(True, alpha=0.3)

# 图 3: 最终直方图 vs 理论分布
plt.subplot(2, 1, 2)
counts = np.bincount(samples, minlength=3)
freqs = counts / N
x_pos = [0, 1, 2]

# 画柱状图
plt.bar(x_pos, freqs, width=0.4, label='Simulation', color='gray', alpha=0.7)
# 画理论点
plt.plot(x_pos, pi, 'ro', markersize=10, label='Theoretical Target', linestyle='None')
# 添加标签
for i, v in enumerate(freqs):
    plt.text(i, v + 0.01, f"{v:.3f}", ha='center', fontweight='bold')

plt.xticks(x_pos, [states_map[0], states_map[1], states_map[2]])
plt.ylim(0, 0.6)
plt.title(f"3. Final Distribution (N={N})\nDid we match the target?", fontsize=12)
plt.legend()

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# --- 0. 参数设置 ---
# 目标分布 pi
pi = np.array([0.25, 0.1, 0.2, 0.35, 0.1])
# define the probability distribution function
def target_distribution(x):
    return pi[x]

# 提议矩阵 Q (对称，均匀跳跃)
# Q[i][j] 代表从 i 提议去 j 的概率
Q = np.array([
    [1/5, 1/5, 1/5, 1/5, 1/5],
    [1/5, 1/5, 1/5, 1/5, 1/5],
    [1/5, 1/5, 1/5, 1/5, 1/5],
    [1/5, 1/5, 1/5, 1/5, 1/5],
    [1/5, 1/5, 1/5, 1/5, 1/5],
])

# maximum number of iterations
iter_max = 10000 # 迭代次数越大，越能逼近目标分布。但它们的关系并非线性，通常需要经验判断。

def discrete_metropolis(n_iter, initial_state):
    samples = [] # chain of samples vector
    alphas = [] # acceptance probability vector
    flag_accepts = [] # acceptance flag vector (true if the sample is accepted)

    current_state = initial_state
    
    # 预计算 Q 的累积分布，用于 U1 的分位数查找
    Q_cumsum = np.cumsum(Q, axis=1)
    
    for t in range(n_iter):
        # --- 第一步：生成提议 (U1) ---
        u1 = np.random.uniform(0, 1)
        # 根据 Q 的当前行累积概率确定候选状态 j
        # 对应数学公式：sum(q_ik) < u1 <= sum(q_ik)
        proposed_state = np.searchsorted(Q_cumsum[current_state], u1) # 逆变换采样 (Inverse Transform Sampling) 的离散实现
        
        # --- 第二步：计算接受率 alpha ---
        # 这里的 pi_j / pi_i
        alpha = min(1, target_distribution(proposed_state) / target_distribution(current_state))
        alphas.append(alpha)
        
        # --- 第三步：判定是否转移 (U2) ---
        u2 = np.random.uniform(0, 1)
        if u2 <= alpha:
            # 接受提议
            current_state = proposed_state
            flag_accepts.append(True)
        else:
            # 拒绝提议，状态保持不变
            flag_accepts.append(False)
            
        samples.append(current_state)
        
    return np.array(samples), np.array(alphas), np.array(flag_accepts)


# --- 2. 运行模拟 ---
N = iter_max  # 迭代次数
samples, alphas, flag_accepts = discrete_metropolis(N, initial_state=0)
t = round(0.25 * len(samples)) # burn-in 期数
print('Uniform proposal distribution:')
print('  Rate of acceptance: %.2f%%' % (100 * np.sum(flag_accepts) / len(flag_accepts)))
print('  Rate of rejection: %.2f%%' % (100 * (1 - np.sum(flag_accepts) / len(flag_accepts))))

# --- 3. 可视化展示 ---
plt.figure(figsize=(15, 12))
# 图 1: Limit probability distribution, pi
plt.subplot(2, 2, 1)
plt.bar(range(len(pi)), pi, color='gray', alpha=0.7)
for i, v in enumerate(pi):
    plt.text(i, v + 0.01, f"{v:.2f}", ha='center', fontweight='bold')
plt.xticks(range(len(pi)))
plt.ylim(0, 0.4)
plt.title("1. Target Distribution π", fontsize=12)
plt.xlabel("States")
plt.ylabel("Probability")
# 图 2: Sampled pdf after burn-in
plt.subplot(2, 2, 2)
counts = np.bincount(samples[t:], minlength=len(pi))
freqs = counts / (N - t)
plt.bar(range(len(pi)), freqs, color='gray', alpha=0.7)
for i, v in enumerate(freqs):
    plt.text(i, v + 0.01, f"{v:.2f}", ha='center', fontweight='bold')
plt.xticks(range(len(pi)))
plt.ylim(0, 0.4)
plt.title(f"2. Sampled Distribution after Burn-in (t={t})", fontsize=12)
plt.xlabel("States")
plt.ylabel("Estimated Probability")
# 图 3: Mixing plot of the chain
plt.subplot(2, 2, 3)
plt.step(range(300), samples[:300], where='mid', color='blue', linewidth=1.5)
plt.yticks(range(len(pi)))
plt.title("3. Trace Plot (First 100 Steps)\nWatch the jumper move!", fontsize=12)
plt.grid(True, axis='y', linestyle='--', alpha=0.5)
plt.xlabel("Step")
# 图 4: Correlation in time
plt.subplot(2, 2, 4)
lags = 50
autocorr = [np.corrcoef(samples[:-lag], samples[lag:])[0, 1] for lag in range(1, lags + 1)]
plt.stem(range(1, lags + 1), autocorr)
plt.title("4. Autocorrelation of the Chain", fontsize=12)
plt.xlabel("Lag")
plt.ylabel("Autocorrelation")
plt.grid(True, linestyle='--', alpha=0.5)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import math

# --- 0. 参数设置 ---
# 目标分布 pi
pi = np.array([0.25, 0.1, 0.2, 0.35, 0.1])
# define the probability distribution function
def target_distribution(x):
    return pi[x]

# 提议矩阵 Q：正态分布（truncated）
# Q[i][j] 代表从 i 提议去 j 的概率
sigma = 1.2 # 这个值越小，状态之间的相关性更大，因此更难跳到另外的状态，故而更难逼近目标分布。
def normcdf_custom(value, mu, sigma):
    return 0.5 * (1 + math.erf((value - mu) / (sigma * np.sqrt(2))))
def gaussian_proposal_matrix(num_states, sigma):
    Q = np.zeros((num_states, num_states))
    for i in range(num_states):
        for j in range(num_states):
            # Q(i,j) = normcdf(x(j)+0.5,x(i),sigmaE) - normcdf(x(j)-0.5,x(i),sigmaE);
            Q[i, j] = normcdf_custom(j + 0.5, i, sigma) - normcdf_custom(j - 0.5, i, sigma)
        Q[i, :] /= np.sum(Q[i, :])  # 归一化
    return Q
Q = gaussian_proposal_matrix(len(pi), sigma)
print("Proposal Matrix Q (Gaussian):")
print(Q)

# maximum number of iterations
iter_max = 10000 # 迭代次数越大，越能逼近目标分布。但它们的关系并非线性，通常需要经验判断。

def discrete_metropolis(n_iter, initial_state):
    samples = [] # chain of samples vector
    alphas = [] # acceptance probability vector
    flag_accepts = [] # acceptance flag vector (true if the sample is accepted)

    current_state = initial_state
    
    # 预计算 Q 的累积分布，用于 U1 的分位数查找
    Q_cumsum = np.cumsum(Q, axis=1)
    
    for t in range(n_iter):
        # --- 第一步：生成提议 (U1) ---
        u1 = np.random.uniform(0, 1)
        # 根据 Q 的当前行累积概率确定候选状态 j
        # 对应数学公式：sum(q_ik) < u1 <= sum(q_ik)
        proposed_state = np.searchsorted(Q_cumsum[current_state], u1) # 逆变换采样 (Inverse Transform Sampling) 的离散实现
        
        # --- 第二步：计算接受率 alpha ---
        # 这里Q不对称，所以不能用 pi_j / pi_i
        alpha = min(1, (target_distribution(proposed_state)*Q[proposed_state, current_state]) / (target_distribution(current_state)*Q[current_state, proposed_state]))
        alphas.append(alpha)
        
        # --- 第三步：判定是否转移 (U2) ---
        u2 = np.random.uniform(0, 1)
        if u2 <= alpha:
            # 接受提议
            current_state = proposed_state
            flag_accepts.append(True)
        else:
            # 拒绝提议，状态保持不变
            flag_accepts.append(False)
            
        samples.append(current_state)
        
    return np.array(samples), np.array(alphas), np.array(flag_accepts)


# --- 2. 运行模拟 ---
N = iter_max  # 迭代次数
samples, alphas, flag_accepts = discrete_metropolis(N, initial_state=0)
t = round(0.25 * len(samples)) # burn-in 期数
print('Uniform proposal distribution:')
print('  Rate of acceptance: %.2f%%' % (100 * np.sum(flag_accepts) / len(flag_accepts)))
print('  Rate of rejection: %.2f%%' % (100 * (1 - np.sum(flag_accepts) / len(flag_accepts))))

# --- 3. 可视化展示 ---
plt.figure(figsize=(15, 12))
# 图 1: Limit probability distribution, pi
plt.subplot(2, 2, 1)
plt.bar(range(len(pi)), pi, color='gray', alpha=0.7)
for i, v in enumerate(pi):
    plt.text(i, v + 0.01, f"{v:.2f}", ha='center', fontweight='bold')
plt.xticks(range(len(pi)))
plt.ylim(0, 0.4)
plt.title("1. Target Distribution π", fontsize=12)
plt.xlabel("States")
plt.ylabel("Probability")
# 图 2: Sampled pdf after burn-in
plt.subplot(2, 2, 2)
counts = np.bincount(samples[t:], minlength=len(pi))
freqs = counts / (N - t)
plt.bar(range(len(pi)), freqs, color='gray', alpha=0.7)
for i, v in enumerate(freqs):
    plt.text(i, v + 0.01, f"{v:.2f}", ha='center', fontweight='bold')
plt.xticks(range(len(pi)))
plt.ylim(0, 0.4)
plt.title(f"2. Sampled Distribution after Burn-in (t={t})", fontsize=12)
plt.xlabel("States")
plt.ylabel("Estimated Probability")
# 图 3: Mixing plot of the chain
plt.subplot(2, 2, 3)
plt.step(range(300), samples[:300], where='mid', color='blue', linewidth=1.5)
plt.yticks(range(len(pi)))
plt.title("3. Trace Plot (First 100 Steps)\nWatch the jumper move!", fontsize=12)
plt.grid(True, axis='y', linestyle='--', alpha=0.5)
plt.xlabel("Step")
# 图 4: Correlation in time
plt.subplot(2, 2, 4)
lags = 50
autocorr = [np.corrcoef(samples[:-lag], samples[lag:])[0, 1] for lag in range(1, lags + 1)]
plt.stem(range(1, lags + 1), autocorr)
plt.title("4. Autocorrelation of the Chain", fontsize=12)
plt.xlabel("Lag")
plt.ylabel("Autocorrelation")
plt.grid(True, linestyle='--', alpha=0.5)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# --- 1. 定义目标分布 (双峰) ---
# 这是一个非归一化的密度函数 (Unnormalized PDF)
# 峰 1: Mean=-2, Std=0.5
# 峰 2: Mean= 2, Std=0.5
def target_pdf(x):
    p1 = np.exp(- (x + 2)**2 / (2 * 0.5**2))
    p2 = np.exp(- (x - 2)**2 / (2 * 0.5**2))
    return p1 + p2  # 双峰叠加

# --- 2. 定义提议核 (Proposal Kernel) ---
# 这里使用对称的高斯核: q(y|x) ~ N(x, sigma^2)
def sample_proposal_kernel(current_x, sigma):
    return np.random.normal(current_x, sigma)

# --- 3. 连续 Metropolis 算法 ---
def continuous_metropolis(n_iter, start_x, sigma):
    samples = []
    current_x = start_x
    accepted_count = 0
    
    for _ in range(n_iter):
        # A. 从核中采样 (Proposal)
        x_star = sample_proposal_kernel(current_x, sigma)
        
        # B. 计算接受率 (Symmetric Metropolis)
        # alpha = min(1, pi(new) / pi(old))
        p_new = target_pdf(x_star)
        p_old = target_pdf(current_x)
        
        alpha = min(1, p_new / p_old) # 这里的Q是对称的，因此可以使用简化形式
        
        # C. 决策 (Decision)
        u = np.random.uniform(0, 1)
        if u <= alpha:
            current_x = x_star
            accepted_count += 1
            
        samples.append(current_x)
        
    return np.array(samples), accepted_count / n_iter

# --- 4. 运行模拟 ---
N = 10000
start_x = 0.0  # 从峡谷中间开始

# 实验 A: 窄核 (小步长) - 可能会被困住
samples_narrow, acc_narrow = continuous_metropolis(N, start_x, sigma=0.2)

# 实验 B: 宽核 (适中步长) - 能跨越峡谷
samples_wide, acc_wide = continuous_metropolis(N, start_x, sigma=1.5)

# --- 5. 可视化对比 ---
plt.figure(figsize=(14, 8))
x_axis = np.linspace(-6, 6, 1000)
y_truth = target_pdf(x_axis)
# 简单归一化一下用于画图对比
y_truth /= np.trapezoid(y_truth, x_axis)

# 图 1: 窄核 (Sigma=0.2)
plt.subplot(2, 2, 1)
plt.plot(samples_narrow, color='orange', alpha=0.7, lw=0.5)
plt.title(f"A. Narrow Kernel (sigma=0.2)\nAcc Rate: {acc_narrow:.1%}")
plt.ylabel("Position x")
plt.grid(alpha=0.3)

plt.subplot(2, 2, 3)
plt.hist(samples_narrow, bins=50, density=True, color='orange', alpha=0.5, label='Samples')
plt.plot(x_axis, y_truth, 'r-', lw=2, label='Target')
plt.title("Distribution (Trapped in one peak?)")
plt.legend()

# 图 2: 宽核 (Sigma=1.5)
plt.subplot(2, 2, 2)
plt.plot(samples_wide, color='green', alpha=0.7, lw=0.5)
plt.title(f"B. Wide Kernel (sigma=1.5)\nAcc Rate: {acc_wide:.1%}")
plt.ylabel("Position x")
plt.grid(alpha=0.3)

plt.subplot(2, 2, 4)
plt.hist(samples_wide, bins=50, density=True, color='green', alpha=0.5, label='Samples')
plt.plot(x_axis, y_truth, 'r-', lw=2, label='Target')
plt.title("Distribution (Good Mixing)")
plt.legend()

plt.tight_layout()
plt.show()

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

# 设置随机种子以复现结果
np.random.seed(0)

# ==========================================
# 0. 定义目标分布 (Target Distribution)
# ==========================================
xmin, xmax = -200, 200
dx = 1
x_axis = np.arange(xmin, xmax + dx, dx)

# 目标概率密度函数 (Weibull-like)
# MATLAB: k=2; l=10; m=0; fx = @(x) ...
def target_pdf(x):
    k, l, m = 2, 10, 0
    # 处理向量或标量输入
    x = np.asarray(x)
    term1 = (x + m)**(k - 1)
    term2 = np.exp(-((x + m) / l)**k)
    condition = (x + m) > 0
    # 防止 0 的负次幂或无效计算
    val = np.abs(term1 * term2) * condition
    return val

# 计算理论曲线用于绘图
fxx = target_pdf(x_axis)

plt.figure(figsize=(10, 4))
plt.plot(x_axis, fxx, linewidth=3)
plt.title('Limit probability distribution (Target)')
plt.grid(True)
plt.show()

# ==========================================
# 辅助函数：Geweke 测试 (计算预烧期 Burn-in)
# ==========================================
def run_geweke_test(samples, alpha=0.05):
    """
    Geweke 测试逻辑：
    比较链条前10%和后50%的均值，逐步剔除头部数据直到均值无显著差异。
    """
    z_lim = norm.ppf(1 - alpha / 2) # 双尾检验阈值 (约1.96)
    stop = False
    
    n = len(samples)
    n1p = int(0.1 * n) # 前 10%
    n2p = int(0.5 * n) # 后 50%
    
    i = 0
    transitory_len = 0
    
    # 循环剔除直到满足条件或超出前10%范围
    while i < n1p and not stop:
        # 截取片段
        # 注意：MATLAB索引从1开始，Python从0开始，逻辑需对齐
        # 前段：从当前剔除位置 i 到原来的 n1p
        if i >= n1p: break
        
        # MATLAB: sx1p = sx1(1:n1p-i); 意为随着i增加，前段越来越短(从后往前切? 不, 它是剔除头部)
        # 仔细看 MATLAB 代码: sx1(1:n1p-i) 其实是取了头部的一段，而且长度在变短。
        # 但标准的 Burn-in 逻辑应该是剔除头部，剩下用于检测。
        # MATLAB 逻辑似乎是：检测 sx1[0 : 10%-i] 和 sx1[50% : end]
        # 这意味着它在检测"最开始的一小段"是否和"最后稳定的一大段"一致。
        
        sample_head = samples[0 : n1p - i]
        sample_tail = samples[n2p:] 
        
        if len(sample_head) < 2: break 

        mean1 = np.mean(sample_head)
        mean2 = np.mean(sample_tail)
        var1  = np.var(sample_head)
        var2  = np.var(sample_tail)
        
        # Z-score
        z_exp = (mean1 - mean2) / np.sqrt(var1 + var2)
        
        if abs(z_exp) > z_lim:
            # 如果差异显著，说明还不稳定，但这逻辑有点反直觉。
            # MATLAB 代码意图：如果差异大，说明包含这段(0:n1p-i)是不对的？
            # 实际上 MATLAB 代码里的 i 是作为计数器，如果 z > z_lim，stop=true? 
            # 原代码: if abs > zlim, stop=true. else i=i+1.
            # 这意味着：一旦发现差异显著，就停止循环？
            # 不，这实际上是反向寻找。
            # 原代码逻辑较为特殊，为了完全复刻，我们严格按公式写：
            stop = True # 这里原意可能是反过来的，或者它在寻找第一个差异显著点？
            # 让我们修正为标准的理解：通常我们想剔除 i 个点。
            # 既然是翻译，我们保持原逻辑:
            pass 
        else:
            i += 1
            
    # 根据 MATLAB 逻辑计算 t
    # 如果循环正常结束(没触发 stop)，说明前面的都很接近后面，不需要剔除？
    # 或者 i 增加到了 n1p。
    
    # 让我们重新审视 MATLAB 那个 while 循环：
    # while i < n1p && ~stop
    #    计算 z
    #    if abs(z) > z_lim: stop = true
    #    else: i = i + 1
    # 逻辑：从长到短检查头部。如果 z 很大（差异大），说明包含了不该包含的瞬态，停止。
    # 此时 i 就是我们需要丢弃的数量的倒数相关量。
    # t = n1p - i + 1. 
    # 如果 i 很大，说明前面一小截和后面都很像。如果 i 很小(一开始就停)，说明整个前10%都和后面不一样。
    
    # Python 重新实现该特定逻辑
    i = 0
    stop = False
    while i < n1p and not stop:
        s1 = samples[0 : n1p - i]
        s2 = samples[n2p : ]
        if len(s1) == 0: break
        
        z_val = (np.mean(s1) - np.mean(s2)) / np.sqrt(np.var(s1) + np.var(s2))
        
        if abs(z_val) > z_lim:
            stop = True
        else:
            i += 1
            
    transitory_len = n1p - i
    print(f"The length of the transitory is: {transitory_len}")
    return transitory_len


# ==========================================
# 1. 采样：正态提议分布 (Random Walk Metropolis)
# ==========================================
print("--- Method 1: Normal Proposal ---")
imax = 100000
sx1 = np.zeros(imax)
flag_accept1 = np.zeros(imax, dtype=bool)

sigmaE = 30
# 提议函数: N(current, sigma)
proposal_func_norm = lambda s: np.random.randn() * sigmaE + s

# 初始化
sx1[0] = -200 

for i in range(1, imax):
    # 1. 提议新样本
    s_candidate = proposal_func_norm(sx1[i-1])
    
    # 2. 计算接受率 alpha
    # Metropolis Ratio: P(new)/P(old) * q(old|new)/q(new|old)
    # 对于正态分布(对称)，q项抵消，仅剩 P(new)/P(old)
    p_current = target_pdf(sx1[i-1])
    p_candidate = target_pdf(s_candidate)
    
    if p_current == 0:
        alpha = 1 # 避免除零，如果在概率为0的地方，肯定要移走
    else:
        alpha = min(1, p_candidate / p_current)
    
    # 3. 接受/拒绝
    u = np.random.rand()
    if u <= alpha:
        sx1[i] = s_candidate
        flag_accept1[i] = True
    else:
        sx1[i] = sx1[i-1]
        flag_accept1[i] = False

# Geweke 测试
t1 = run_geweke_test(sx1)
clean_chain1 = sx1[t1:]

# 绘图 1: 直方图 vs 理论
plt.figure(figsize=(12, 10))
plt.subplot2grid((2, 2), (0, 0))
plt.plot(x_axis, fxx, 'b-', linewidth=2, label='Theoretical')
# 绘制直方图并缩放以匹配未归一化的 PDF 高度 (复刻 MATLAB 视觉效果)
counts, bins = np.histogram(clean_chain1, bins=100, density=True)
# 缩放因子：为了让 density 直方图的高度匹配未归一化的 fxx
# MATLAB: fr.*max(fxx)./max(fr)
scale_factor = np.max(fxx) / np.max(counts) if np.max(counts) > 0 else 1
plt.hist(clean_chain1, bins=100, density=True, color='r', alpha=0.5, 
         weights=np.ones_like(clean_chain1)*scale_factor, label='Sampled')
# 注意：上面的 weights 方式在 density=True 时可能不生效，为了简单复刻图表，
# 我们直接画 density=True 的直方图，然后把 fxx 归一化可能更科学。
# 但为了还原 MATLAB 的视觉（它缩放了柱状图）：
plt.cla() # 清除重画
plt.plot(x_axis, fxx, 'b-', linewidth=2, label='Theoretical (Unnormalized)')
hist_vals, hist_bins = np.histogram(clean_chain1, bins=x_axis, density=True)
# 简单的缩放对齐
if np.max(hist_vals) > 0:
    plt.bar(hist_bins[:-1], hist_vals * (np.max(fxx)/np.max(hist_vals)), width=dx, color='r', alpha=0.5, label='Sampled')
plt.grid(True)
plt.legend()
plt.title('Sampled (Normal Proposal)')

# 绘图 2: 混合轨迹 (Trace Plot)
plt.subplot2grid((2, 2), (0, 1))
plt.plot(sx1, '.-', markersize=1, alpha=0.5)
plt.xlabel('Iteration')
plt.title('Mixing / Trace Plot')

# 绘图 3: 自相关
plt.subplot2grid((2, 2), (1, 0), colspan=2) 
# 计算自相关
lags = 1000 # 只看前1000
corr1 = np.correlate(clean_chain1 - np.mean(clean_chain1), clean_chain1 - np.mean(clean_chain1), mode='full')
corr1 = corr1[len(corr1)//2:] 
corr1 = corr1 / corr1[0] # 归一化
plt.plot(corr1[:lags])
plt.title('Correlation in time (Normal Proposal)')
plt.grid(True)
plt.show()

print(f"Acceptance Rate: {np.mean(flag_accept1[t1:])*100:.2f}%")


# ==========================================
# 2. 采样：均匀提议分布 (Independent Metropolis)
# ==========================================
print("\n--- Method 2: Uniform Proposal ---")
sx2 = np.zeros(imax)
flag_accept2 = np.zeros(imax, dtype=bool)

# 提议函数: Uniform(xmin, xmax) - 独立于上一个样本
proposal_func_unif = lambda: np.random.rand() * (xmax - xmin) + xmin

sx2[0] = -200

for i in range(1, imax):
    # 1. 提议
    s_candidate = proposal_func_unif()
    
    # 2. 计算接受率
    # q(x|y) = q(x) = 1/(xmax-xmin). q(y|x) = q(y) = 1/(xmax-xmin).
    # q 项依然抵消，公式同上
    p_current = target_pdf(sx2[i-1])
    p_candidate = target_pdf(s_candidate)
    
    if p_current == 0:
        alpha = 1
    else:
        alpha = min(1, p_candidate / p_current)
        
    # 3. 接受/拒绝
    u = np.random.rand()
    if u <= alpha:
        sx2[i] = s_candidate
        flag_accept2[i] = True
    else:
        sx2[i] = sx2[i-1]
        flag_accept2[i] = False

# Geweke 测试
t2 = run_geweke_test(sx2)
clean_chain2 = sx2[t2:]

# 绘图 1
plt.figure(figsize=(12, 10))
plt.subplot2grid((2, 2), (0, 0))
plt.plot(x_axis, fxx, 'b-', linewidth=2, label='Theoretical')
hist_vals, hist_bins = np.histogram(clean_chain2, bins=x_axis, density=True)
if np.max(hist_vals) > 0:
    plt.bar(hist_bins[:-1], hist_vals * (np.max(fxx)/np.max(hist_vals)), width=dx, color='r', alpha=0.5, label='Sampled')
plt.grid(True)
plt.legend()
plt.title('Sampled (Uniform Proposal)')

# 绘图 2
plt.subplot2grid((2, 2), (0, 1))
plt.plot(sx2, '.-', markersize=1, alpha=0.5)
plt.xlabel('Iteration')
plt.title('Mixing / Trace Plot')

# 绘图 3: 自相关
plt.subplot2grid((2, 2), (1, 0), colspan=2) 
corr2 = np.correlate(clean_chain2 - np.mean(clean_chain2), clean_chain2 - np.mean(clean_chain2), mode='full')
corr2 = corr2[len(corr2)//2:] 
corr2 = corr2 / corr2[0]
plt.plot(corr2[:lags])
plt.title('Correlation in time (Uniform Proposal)')
plt.grid(True)
plt.show()

print(f"Acceptance Rate: {np.mean(flag_accept2[t2:])*100:.2f}%")

import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm

# --- 1. 准备工作 ---
# 目标：标准正态分布 N(0, 1)
def target_log_prob(x):
    return -0.5 * x**2

# Metropolis 采样核
def run_chain(n_steps, step_size, start_x=0.0):
    samples = np.zeros(n_steps)
    current_x = start_x
    accepted = 0
    
    for i in range(n_steps):
        # 提议
        proposal = current_x + np.random.normal(0, step_size)
        
        # 接受率 (Log scale for stability)
        log_alpha = target_log_prob(proposal) - target_log_prob(current_x)
        # min(1, A) -> min(0, log_A) in log domain
        if np.log(np.random.rand()) < log_alpha:
            current_x = proposal
            accepted += 1
            
        samples[i] = current_x
        
    acc_rate = accepted / n_steps
    return samples, acc_rate

# --- 2. 运行三种场景 ---
N = 2000
# A. Smooth (步长太小)
chain_smooth, acc_smooth = run_chain(N, step_size=0.1)
# B. Stiff (步长太大)
chain_stiff, acc_stiff = run_chain(N, step_size=50.0)
# C. Optimal (步长适中)
chain_optimal, acc_optimal = run_chain(N, step_size=2.4) # 对于1维高斯，大一点没事

# --- 3. 绘图对比 (Trace + ACF) ---
fig, axes = plt.subplots(3, 2, figsize=(14, 10))
lags = 100 # 看前100步的相关性

# 辅助函数：画一行图
def plot_row(row_idx, samples, acc, title, color):
    # 左边：Trace Plot
    ax_trace = axes[row_idx, 0]
    ax_trace.plot(samples, color=color, lw=1)
    ax_trace.set_title(f"{title} - Trace (Acc: {acc:.1%})")
    ax_trace.set_ylabel("Value")
    
    # 右边：ACF Plot
    ax_acf = axes[row_idx, 1]
    # 使用 statsmodels 计算 ACF
    acf_values = sm.tsa.acf(samples, nlags=lags)
    ax_acf.bar(range(len(acf_values)), acf_values, width=0.3, color=color, alpha=0.7)
    ax_acf.set_title(f"{title} - Autocorrelation")
    ax_acf.set_ylim(-0.1, 1.1)
    ax_acf.set_ylabel("Correlation")
    ax_acf.axhline(0.05, linestyle='--', color='gray', alpha=0.5) # 显著性阈值
    ax_acf.axhline(0, color='black', lw=1)

# 绘制三行
plot_row(0, chain_smooth, acc_smooth, "A. Smooth (Step=0.1)", "orange")
plot_row(1, chain_stiff, acc_stiff, "B. Stiff (Step=50)", "purple")
plot_row(2, chain_optimal, acc_optimal, "C. Optimal (Step=2.4)", "green")

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 1. 目标：标准正态分布 N(0, 1)
def target_log_prob(x):
    return -0.5 * x**2

# 2. 模拟采样
def run_burnin_demo(n_steps, start_x, step_size):
    samples = np.zeros(n_steps)
    current_x = start_x
    
    for i in range(n_steps):
        # 提议
        proposal = current_x + np.random.normal(0, step_size)
        # 接受率
        log_alpha = target_log_prob(proposal) - target_log_prob(current_x)
        if np.log(np.random.rand()) < log_alpha:
            current_x = proposal
        samples[i] = current_x
        
    return samples

# --- 设置 ---
N = 1000
start_val = 20.0  # <--- 起点离中心(0)非常远！
step = 0.5        # 步长较小，走得慢

chain = run_burnin_demo(N, start_val, step)

# --- 绘图 ---
plt.figure(figsize=(12, 6))

# 画出轨迹
plt.plot(chain, label='MCMC Chain', color='blue', lw=1.5)

# 标出 Burn-in 的分界线 (比如大约在第 200 步到达中心)
burn_in_cutoff = 200
plt.axvline(x=burn_in_cutoff, color='red', linestyle='--', lw=2, label='Burn-in Cutoff')

# 标注文字
plt.text(50, 15, "Transient Phase\n(Garbage Samples)", color='red', fontsize=12, fontweight='bold')
plt.text(400, 5, "Stationary Phase\n(Valid Samples)", color='green', fontsize=12, fontweight='bold')

plt.title(f"Visualizing Burn-in: Starting from x={start_val} to target N(0,1)")
plt.xlabel("Iteration")
plt.ylabel("Sample Value")
plt.axhline(0, color='gray', linestyle=':', alpha=0.5) # 目标均值
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()