2. Metropolis-Hastings算法：算法原理、提议分布的选择、接受率计算、随机游走Metropolis、实现步骤与代码示例

2.1 算法原理：从拒绝到接受的艺术

Metropolis-Hastings算法，说白了就是一套「有策略地碰运气」的方法。它的核心思想很简单：我们想从一个很难直接抽样的分布中采样，那就先随便猜一个点，然后根据某种规则决定要不要「接受」这个新点。

我个人习惯把MH算法理解成一个「挑剔的评委」：

• 你提出一个新样本（候选点）
• 评委计算一下这个新样本比当前样本好多少
• 如果更好，直接接受；如果更差，也有一定概率接受——但概率要小一些

这样走下来，最终留下的样本就会服从我们想要的目标分布。嗯，这里要注意：这个「一定概率接受更差点」的设计，恰恰是MH算法的精髓。如果没有这一步，算法就会卡在局部最优出不来。

2.2 提议分布的选择：这步走对了，后面就顺了

提议分布 q(x' | x) 决定了你下一步往哪儿跳。我在项目中遇到过最头疼的问题，十有八九都是提议分布没选好。

选提议分布，有两条黄金法则：

1. 对称性优先：如果 q(x' | x) = q(x | x')，公式可以简化，接受率计算也快。随机游走Metropolis就是典型例子。
2. 步长要适中：步长太小，样本之间高度相关，收敛慢；步长太大，接受率低，浪费计算资源。

你想想看，如果提议分布选得太「激进」，每次跳得很远，那大部分候选点都会被拒绝，算法几乎原地踏步。反过来，如果跳得太「保守」，虽然接受率高，但样本在参数空间里爬得比蜗牛还慢。

2.3 接受率计算：别被公式吓到

接受率的计算其实就三步：

1. 计算目标分布在新旧两个点的比值：π(θ*) / π(θ_t)
2. 计算提议分布在两个方向上的比值：q(θ_t | θ*) / q(θ* | θ_t)
3. 两个比值相乘，取 min(1, 结果)

我曾经犯过一个低级错误：忘了取对数。当目标分布的值非常小（比如10^-100），直接计算比值会下溢成0。后来我学乖了，一律用对数空间计算，最后再exp回来。

2.4 随机游走Metropolis：最常用的变体

随机游走Metropolis（RWM）是MH算法最经典的变体。它的提议分布是：

θ* = θ_t + ε，其中 ε ~ N(0, σ²)

因为正态分布是对称的，所以 q(θ* | θ_t) = q(θ_t | θ*)，接受率简化为：

α = min(1, π(θ*) / π(θ_t))

说白了，就是只看目标分布的比值。这大大简化了计算。我在做期权定价模型校准的时候，用的就是RWM，配合自适应步长调整，效果相当不错。

RWM的步长σ是关键参数。太小了样本像喝醉了酒一样原地打转，太大了又容易「跳崖」。我一般用以下策略调参：

• 先跑500步，统计接受率
• 如果接受率 < 20%，σ乘以0.8
• 如果接受率 > 50%，σ乘以1.2
• 重复直到接受率稳定在30%左右

2.5 实现步骤：手把手教你写代码

下面是我常用的MH算法实现模板。这个代码我用了好几年，改来改去，现在算是比较成熟了。

import numpy as np
import matplotlib.pyplot as plt

def metropolis_hastings(target_log_pdf, initial, n_samples, proposal_std):
    """
    Metropolis-Hastings采样器
    
    参数:
        target_log_pdf: 目标分布的对数概率密度函数
        initial: 初始值
        n_samples: 采样数量
        proposal_std: 提议分布的标准差（步长）
    
    返回:
        samples: 采样结果数组
        acceptance_rate: 接受率
    """
    samples = np.zeros(n_samples)
    samples[0] = initial
    current = initial
    n_accepted = 0
    
    for i in range(1, n_samples):
        # 1. 从提议分布生成候选点
        candidate = current + np.random.normal(0, proposal_std)
        
        # 2. 计算接受率（对数形式）
        log_accept_ratio = target_log_pdf(candidate) - target_log_pdf(current)
        
        # 3. 判断是否接受
        if np.log(np.random.uniform()) < log_accept_ratio:
            current = candidate
            n_accepted += 1
        
        samples[i] = current
    
    acceptance_rate = n_accepted / (n_samples - 1)
    return samples, acceptance_rate

# 示例：从混合高斯分布采样
def target_log_pdf(x):
    # 混合高斯：0.3*N(-2, 0.5) + 0.7*N(3, 1)
    log_p1 = np.log(0.3) - 0.5*np.log(2*np.pi*0.5) - (x+2)**2/(2*0.5)
    log_p2 = np.log(0.7) - 0.5*np.log(2*np.pi*1) - (x-3)**2/(2*1)
    return np.logaddexp(log_p1, log_p2)

# 运行采样
samples, acc_rate = metropolis_hastings(
    target_log_pdf, 
    initial=0.0, 
    n_samples=10000, 
    proposal_std=2.0
)

print(f"接受率: {acc_rate:.2%}")

2.6 代码解读与实战技巧

这段代码有几个关键点，我挨个说一下：

• 对数空间计算：所有概率都用对数形式，避免数值下溢。np.logaddexp是个好函数，专门处理两个对数概率相加的场景。
• 接受判断技巧：用 np.log(np.random.uniform()) < log_accept_ratio 代替直接比较概率，既避免了除法，又保持了数值稳定性。
• 初始值选择：别太离谱。我建议先用优化方法找个大概位置，再跑MCMC。否则前几百步都在「热身」，浪费计算资源。

2.7 知识体系总览

下面这张图总结了MH算法的核心脉络，我画了很久才满意：

2.8 总结与实用建议

MH算法是MCMC家族里最灵活、最通用的成员。它不要求目标分布是标准形式，只要你能写出概率密度函数（或者正比于它的函数），就能用。

最后分享几个我踩过的坑：

• 别用太差的初始值：我曾经从一个离目标分布很远的点开始跑，结果前2000步全在「爬山」，浪费了大量计算。现在我会先用优化器找一下大概位置。
• 监控接受率：每次跑完都看一眼接受率。如果低于5%，基本可以断定步长设置有问题。
• 多链并行：如果计算资源允许，同时跑3-5条链，从不同初始值出发。这样能判断算法是否收敛到了正确的分布。

---

公众号:蓝海资料掘金营，微信deep3321