MCMC 采样

一、直观理解：用“随机游走”访问难以直接抽样的分布

MCMC（Markov Chain Monte Carlo，马尔可夫链蒙特卡洛）用于解决一类常见困难：目标分布 $\pi(x)$ 的密度可以计算（或可计算到差一个未知归一化常数），但无法直接独立抽样。

直观上，算法在状态空间中构造一条马尔可夫链，使链的平稳分布恰好是 $\pi$。链运行足够长时间并舍弃预热期后，轨迹上的样本可以近似看作来自 $\pi$，再用来估计期望、后验均值、可信区间等。

可以把 MCMC 想象成在一座起伏的山地（对应概率密度）上按规则随机漫步：高密度区域更容易停留，低密度区域也会偶尔访问。时间足够长后，停留时间比例近似等于概率质量比例。

在数学建模与贝叶斯统计中，MCMC 常用于：

• 后验分布没有闭式解时的参数推断；
• 高维积分与归一化常数未知的模型；
• 层级模型、状态空间模型中的隐变量采样。

二、适用场景与优缺点

类型	说明
适用场景	目标密度已知至多差常数；需要后验均值、分位数或边缘分布；直接独立抽样不可行；维数中等且能设计合理提议分布。
核心优点	不要求归一化常数；框架统一（MH、Gibbs、HMC 等）；可估计完整后验而非点估计；易于扩展到层级模型。
局限与缺点	样本相关，有效样本量小于名义样本量；需要诊断收敛与混合；预热与步长选择有技巧；高维强相关时朴素 MH 效率差。

使用前应确认：

1. 目标密度 $\pi(x)$（或 $\tilde\pi(x)\propto\pi(x)$）能否稳定计算；
2. 参数空间是否有约束（正值、概率单纯形等）；
3. 是否有共轭结构可优先用 Gibbs；
4. 需要估计哪些后验泛函；
5. 计算预算允许的链长与链数。

三、核心原理：构造以目标分布为平稳分布的马尔可夫链

1. Metropolis-Hastings（MH）算法

设当前状态为 $x$，从提议分布 $q(\cdot\mid x)$ 生成候选 $y$，接受概率为

$$\alpha(x,y)=\min\left(1,\frac{\pi(y)q(x\mid y)}{\pi(x)q(y\mid x)}\right)$$

以概率 $\alpha$ 转移到 $y$，否则停留在 $x$。当 $q$ 对称（如随机游走正态提议）时，比率简化为 $\pi(y)/\pi(x)$。

2. 细致平衡与平稳分布

在正则条件下，MH 核满足细致平衡，从而 $\pi$ 是平稳分布。因此长时间运行后 $X_t\sim\pi$（近似）。

3. 用 MCMC 样本估计期望

对可积函数 $g$，

$$\hat\mu=\frac{1}{M}\sum_{t=B+1}^{B+M}g(X_t)$$

其中 $B$ 为预热（burn-in）长度，$M$ 为保留样本数。由于样本自相关，方差估计应使用有效样本量（ESS）或批次均值法，而不是简单的 $s/\sqrt{M}$。

4. Gibbs 采样（简述）

若满条件分布 $\pi(x_j\mid x_{-j})$ 易于抽样，可依次更新各坐标。Gibbs 是 MH 的特例，接受率恒为 1。

四、关键实现决策

• 提议步长：过小则链爬行缓慢，过大则拒绝率过高；随机游走 MH 常将接受率调到大致 20%–40%（维数相关）；
• 预热 $B$：观察迹图，丢弃明显未进入高概率区的前期样本；
• 链长与多链：使用多条不同初值的链，比较链间/链内方差（如 $\hat R$ 思想）；
• 收敛诊断：迹图、自相关图、有效样本量；不要只看“好像平稳”；
• 再参数化：强相关参数可变换坐标以改善混合。

五、MATLAB 示例

目标：用随机游走 Metropolis 从标准正态与一个双峰分布采样，并估计期望。

%% 示例 1：标准正态的随机游走 Metropolis
rng(1);
nIter = 8000;
burnin = 2000;
x = zeros(nIter, 1);
x(1) = 0;
sigmaProp = 0.8;          % 提议标准差
logTarget = @(z) -0.5*z.^2;  % 标准正态（差常数）

nAccept = 0;
for t = 2:nIter
    y = x(t-1) + sigmaProp * randn;
    logAlpha = logTarget(y) - logTarget(x(t-1));
    if log(rand) < logAlpha
        x(t) = y;
        nAccept = nAccept + 1;
    else
        x(t) = x(t-1);
    end
end

samples = x(burnin+1:end);
accRate = nAccept / (nIter - 1);
muHat = mean(samples);
varHat = var(samples);

fprintf("接受率 = %.3f\n", accRate);
fprintf("后验均值估计 = %.4f（真值 0）\n", muHat);
fprintf("后验方差估计 = %.4f（真值 1）\n", varHat);

figure;
subplot(2,1,1);
plot(samples);
xlabel("迭代"); ylabel("x"); title("迹图（预热后）"); grid on;
subplot(2,1,2);
histogram(samples, 40, "Normalization", "pdf"); hold on;
xx = linspace(-4, 4, 200);
plot(xx, normpdf(xx), "r", "LineWidth", 1.5);
xlabel("x"); ylabel("密度"); title("样本直方图 vs 标准正态");
legend("MCMC", "N(0,1)"); grid on;

%% 示例 2：双峰目标密度 \propto exp(-x^4/4 + x^2/2)
rng(2);
nIter = 12000;
burnin = 3000;
x = zeros(nIter, 1);
x(1) = 0;
sigmaProp = 1.0;
logTarget = @(z) -0.25*z.^4 + 0.5*z.^2;

for t = 2:nIter
    y = x(t-1) + sigmaProp * randn;
    if log(rand) < logTarget(y) - logTarget(x(t-1))
        x(t) = y;
    else
        x(t) = x(t-1);
    end
end

samples = x(burnin+1:end);
fprintf("双峰分布样本均值 = %.4f\n", mean(samples));
fprintf("样本落在正半轴比例 = %.3f\n", mean(samples > 0));

figure;
histogram(samples, 50, "Normalization", "pdf");
xlabel("x"); ylabel("密度");
title("双峰目标的 MCMC 样本");
grid on;

六、如何解释结果

报告 MCMC 结果时应包括：

1. 目标分布/后验的定义与先验选择；
2. 算法类型（MH/Gibbs 等）与提议分布；
3. 链数、链长、预热、稀释（thinning，如使用）；
4. 接受率、迹图与收敛诊断结论；
5. 参数的后验均值、标准差、分位数或可信区间；
6. 有效样本量或自相关程度的说明；
7. 对实质性结论的解释（哪一参数显著偏离 0、哪一情景概率更高等）。

注意：MCMC 给出的是分布层面的信息，点估计只是后验的一个摘要。

七、常见误区

1. 把相关样本当成独立样本

名义上有 1 万个点，有效信息可能远少于此。方差与置信区间必须考虑自相关。

2. 不做收敛诊断

链尚未混合就报告“后验均值”，结论可能完全取决于初值。

3. 接受率极端仍强行使用

接受率接近 0 或接近 1 通常说明步长不当，应调整提议尺度或改用更好算法。

4. 忽略多峰

单链可能长期困在一个峰。应尝试多初值、并行回火或更强混合的采样器。

5. 把 MCMC 神化成“自动贝叶斯”

先验、似然与识别性仍需认真设定；采样器只是计算工具。

八、小结

MCMC 通过构造以目标分布为平稳分布的马尔可夫链，在难以直接抽样时仍能估计后验与复杂积分。Metropolis-Hastings 是最基础的通用框架；诊断收敛、报告不确定性、检查有效样本量，与写出正确的目标密度同等重要。

在建模论文中，MCMC 的价值不在于“运行了采样器”，而在于把参数不确定性完整地传递到决策相关的后验结论中。