用Python手搓MCMC采样器：逃离数学公式的实战指南

当概率分布的解析解遥不可及时，MCMC就像一台概率世界的3D打印机——通过随机游走的轨迹，将抽象分布具象化为可触摸的样本点。本文将以Metropolis-Hastings算法为例，带你用Python从零构建一个会"学习"的采样器，它能在未知的概率地形中自主探索，最终绘制出完整的分布地图。

1. 准备知识工具箱

在开始编码之前，我们需要理解几个核心概念：

• 马尔可夫链 ：下一个状态只依赖当前状态的记忆丧失型随机过程，就像醉酒者的行走路径
• 细致平衡条件 ：保证采样器最终收敛到目标分布的关键数学性质
• 接受率 ：决定是否跳转到新采样点的概率阀门

安装必要的Python环境：

pip install numpy matplotlib scipy

2. 构建目标分布：城堡房间的隐喻

假设我们要采样的目标分布模拟城堡房间的人气值，定义为一个双峰分布：

import numpy as np
from scipy.stats import norm

def target_distribution(x):
    """双峰高斯混合分布"""
    return 0.3*norm.pdf(x, loc=-2, scale=1) + 0.7*norm.pdf(x, loc=3, scale=0.5)

可视化这个分布：

import matplotlib.pyplot as plt

x = np.linspace(-5, 6, 1000)
plt.plot(x, target_distribution(x))
plt.title("城堡房间人气分布")
plt.xlabel("房间位置")
plt.ylabel("相对人气值")

3. Metropolis-Hastings算法实现

算法核心流程可分为四个步骤：

1. 初始化 ：从任意起点开始（这里选择x=0）
2. 提议移动 ：根据建议分布生成新位置（常用正态分布）
3. 接受判断 ：计算接受概率决定是否跳转
4. 记录样本 ：保存当前状态到样本链

完整实现代码：

def metropolis_hastings(target, n_samples, initial=0, proposal_width=1):
    """
    Metropolis-Hastings采样器实现
    
    参数：
        target: 目标分布函数
        n_samples: 采样数量
        initial: 初始位置
        proposal_width: 建议分布的标准差
        
    返回：
        samples: 采样链
        accepted: 接受率
    """
    samples = [initial]
    current = initial
    accepted = 0
    
    for _ in range(n_samples):
        # 提议新位置（对称建议分布）
        proposed = np.random.normal(current, proposal_width)
        
        # 计算接受概率
        prob_accept = min(1, target(proposed) / target(current))
        
        # 决定是否接受
        if np.random.rand() < prob_accept:
            current = proposed
            accepted += 1
        
        samples.append(current)
    
    return np.array(samples), accepted / n_samples

4. 采样过程可视化与分析

运行采样器并观察收敛过程：

samples, acceptance_rate = metropolis_hastings(target_distribution, 10000)

print(f"接受率: {acceptance_rate:.2%}")  # 理想值通常在20-50%之间

绘制采样轨迹和分布对比：

plt.figure(figsize=(12, 6))

# 采样轨迹
plt.subplot(1, 2, 1)
plt.plot(samples[:500], alpha=0.6)
plt.title("前500次采样轨迹")
plt.xlabel("迭代次数")
plt.ylabel("采样位置")

# 分布对比
plt.subplot(1, 2, 2)
plt.hist(samples[500:], bins=50, density=True, alpha=0.6, label="采样分布")
plt.plot(x, target_distribution(x), 'r', label="真实分布")
plt.title("采样分布 vs 真实分布")
plt.legend()

5. 调参实战：建议分布的影响

建议分布的宽度（proposal_width）显著影响采样效率：

宽度值	接受率	收敛速度	样本相关性
0.1	>80%	慢	高
1.0	40-60%	适中	中
10.0	<20%	快	低

测试不同参数的效果：

for width in [0.1, 1.0, 10.0]:
    samples, rate = metropolis_hastings(target_distribution, 5000, proposal_width=width)
    print(f"宽度 {width}: 接受率 {rate:.2%}")

6. 高级技巧：自适应MCMC

基础算法的一个改进是让建议分布自动调整宽度：

def adaptive_mh(target, n_samples, initial_width=1, target_rate=0.4):
    samples = [0]
    current = 0
    accepted = 0
    width = initial_width
    
    for i in range(n_samples):
        proposed = np.random.normal(current, width)
        prob_accept = min(1, target(proposed) / target(current))
        
        if np.random.rand() < prob_accept:
            current = proposed
            accepted += 1
        
        samples.append(current)
        
        # 每100步调整一次建议分布宽度
        if i % 100 == 0 and i > 0:
            current_rate = accepted / 100
            width *= 1 + 0.5*(current_rate - target_rate)
            accepted = 0
    
    return np.array(samples)

7. 诊断采样质量

评估MCMC结果可靠性的三种方法：

1. 轨迹图 ：观察采样链是否充分混合

   plt.plot(samples[::10])  # 降采样显示
   

2. 自相关函数 ：检测样本独立性

   from statsmodels.tsa.stattools import acf
   plt.plot(acf(samples, nlags=50))
   

3. 多链比较 ：运行多个独立链检查收敛

   chains = [metropolis_hastings(target_distribution, 2000)[0] for _ in range(4)]
   plt.plot(np.array(chains).T, alpha=0.6)
   

提示：当R-hat统计量接近1时，表示各链已收敛到相同分布