从猫的随机游走到MCMC：用生活场景理解马尔可夫链蒙特卡罗-平芜编程栈

从猫的随机游走到MCMC：用生活场景理解马尔可夫链蒙特卡罗

1. 当猫咪成为概率大师

我家那只橘猫每天在家里的活动轨迹，简直是一部活生生的随机游走教科书。早晨从猫窝（状态A）出发，有60%概率溜达到阳台晒太阳（状态B），30%概率跳到书桌捣乱（状态C），还有10%概率继续赖床。这种"下一步去哪只取决于现在位置"的特性，恰好揭示了马尔可夫链的核心思想——无记忆性。

猫的转移概率矩阵可以这样表示：

当前状态 \ 下一状态	猫窝(A)	阳台(B)	书桌(C)
猫窝(A)	0.1	0.6	0.3
阳台(B)	0.4	0.3	0.3
书桌(C)	0.2	0.2	0.6

观察几周后发现有趣的现象：无论猫咪当天从哪个位置开始，长期来看：

在猫窝出现的概率稳定在约23%
在阳台出现的概率约42%
在书桌的概率约35%

这个稳定分布就是马尔可夫链的平稳分布。类似地，MCMC方法正是通过构造特殊的"猫咪游走规则"，使得最终的稳定状态恰好是我们想研究的复杂概率分布。

2. 蒙特卡洛的魔法：用随机性解决确定性难题

想象我们要计算蝙蝠侠标志形状的客厅面积，但用尺子测量其不规则边缘太困难。于是我们：

在客厅地板上画一个规则的大矩形（覆盖整个蝙蝠侠标志）
向矩形区域随机抛撒猫粮颗粒
统计落在蝙蝠侠形状内的颗粒比例
用比例乘以矩形面积得到近似值

这就是蒙特卡洛方法的典型应用。当需要计算复杂形状的面积或积分时，公式推导可能极其困难，但通过随机采样却能获得令人满意的近似解。

import numpy as np def estimate_area(n_samples): in_shape = 0 for _ in range(n_samples): x, y = np.random.uniform(0,10), np.random.uniform(0,5) if (y < 2.5*np.sin(x) + 3) and (y > abs(x-5)**0.5): # 蝙蝠侠形状的简化条件 in_shape += 1 return (in_shape/n_samples) * 50 # 假设矩形面积50平米 print(estimate_area(100000)) # 输出近似面积

但当涉及复杂概率分布时，直接采样可能效率极低。就像如果蝙蝠侠标志只占客厅的1%，绝大多数猫粮都会浪费在无效区域。这时就需要更聪明的采样策略——MCMC。

3. 构建智能游走规则：Metropolis-Hastings算法

回到猫咪的例子，假设我们想让它按照特殊偏好分布游走：

厨房概率最高（因为有猫粮）
客厅次之（有阳光）
卧室最低（主人禁止上床）

但猫咪只会凭本能随机走动。这时可以设计这样的接受-拒绝规则：

猫咪当前在客厅，按本能想移动到厨房
计算接受概率 α = min(1, 厨房偏好/客厅偏好)
如果α=0.8，就掷一枚不均匀硬币（80%正面）
出现正面则接受移动，否则留在原地

这种策略就是Metropolis-Hastings算法的核心。通过精心设计的接受概率，最终会使猫咪的停留比例完全符合我们的目标分布。

算法步骤可视化：

从任意初始位置x₀开始（比如卧室）
生成候选位置x*（按简单提议分布，如均匀随机）
计算接受概率：
```
α = min(1, p(x*)/p(xₜ))
```
以概率α接受x*，否则保持xₜ
重复步骤2-4，记录所有访问位置

经过足够多次迭代后，位置序列的分布就会收敛到目标分布p(x)。关键在于提议分布可以很简单，但通过接受概率的调节，最终能采样复杂分布。

4. Gibbs采样：高维空间的优雅探索

当处理多变量分布时（比如同时追踪猫咪位置和姿势），Gibbs采样展现出独特优势。假设我们想研究：

位置：{客厅, 厨房, 卧室}
姿势：{趴着, 坐着, 站着}

Gibbs采样的步骤就像猫咪轮换调整每个变量：

固定当前姿势"趴着"，按条件概率更新位置
固定新位置"厨房"，按条件概率更新姿势
重复这个过程

# 伪代码示例 def gibbs_sampling(iterations): pos = "客厅" # 初始值 pose = "趴着" samples = [] for _ in range(iterations): # 固定pose更新pos pos = sample_from(p(pos | pose)) # 固定pos更新pose pose = sample_from(p(pose | pos)) samples.append((pos, pose)) return samples

这种方法特别适合各变量间存在依赖关系的场景。就像猫咪的姿势常与位置相关（厨房多站着，客厅多趴着），Gibbs采样能有效捕捉这种联合分布特征。

5. 实践中的技巧与陷阱

实际应用MCMC时，有几个关键注意事项：

燃烧期处理：

前1000次迭代可能还未收敛（就像猫咪刚开始在陌生环境探索）
应丢弃这部分"热身"样本

自相关性：

相邻样本高度相关（猫咪不会瞬间跨房间）
解决方案：
- 每k个样本保留一个
- 增加提议分布的方差

收敛诊断：

运行多条链观察是否趋于相同分布
计算Gelman-Rubin统计量

经验法则：当变异系数小于1.1时，可认为收敛

效率优化技巧：

对连续变量，建议接受率在20-50%之间
高维空间考虑分块更新
结合梯度信息（如HMC算法）

6. 从理论到实践：一个完整案例

假设我们要估计新型猫粮对猫咪活动的影响，建立如下贝叶斯模型：

参数θ表示猫粮效果
先验p(θ)为N(0,1)
似然p(data|θ)为二项分布

后验分布p(θ|data)难以直接计算，采用MCMC：

import numpy as np from scipy.stats import norm, binom def metropolis(data, n_iterations): theta_current = 0.5 # 初始值 samples = [] for i in range(n_iterations): theta_proposed = norm(theta_current, 0.1).rvs() p_current = norm(0,1).pdf(theta_current) * binom.pmf(sum(data), len(data), theta_current) p_proposed = norm(0,1).pdf(theta_proposed) * binom.pmf(sum(data), len(data), theta_proposed) alpha = min(1, p_proposed/p_current) if np.random.rand() < alpha: theta_current = theta_proposed samples.append(theta_current) return samples # 假设观察数据：30次观察中20次选择新猫粮 samples = metropolis([1]*20 + [0]*10, 10000) print(f"后验均值：{np.mean(samples[1000:]):.3f}") # 丢弃前1000次

这个简单实现已经能给出θ的可靠估计。在实际数据分析中，我们可能会使用PyMC3等专业库，但它们核心原理与此一致。