第十九章：马尔可夫链蒙特卡罗法
19.1 引言
马尔可夫链蒙特卡罗法（Markov Chain Monte Carlo，简称MCMC）是一类基于马尔可夫链的随机模拟方法，主要用于从复杂概率分布中进行数值采样和积分计算。MCMC方法在统计学、机器学习、计算物理、生物信息学等领域有着广泛的应用，是现代统计计算的核心技术之一。
MCMC方法的兴起与发展与统计学从频率学派向贝叶斯学派的演进密切相关。贝叶斯统计的核心任务是根据观测数据推断未知参数的后验分布，而后验分布的计算往往涉及高维积分，这些积分在解析上通常无法求解。传统的数值积分方法在高维情况下面临"维数灾难"，计算复杂度随维数指数增长。MCMC方法通过构造一个马尔可夫链，使其平稳分布恰好等于目标后验分布，然后从该链上采集样本，用样本均值近似积分，从而有效克服了高维积分的困难。
本章系统介绍MCMC方法的基本原理、核心算法和收敛诊断技术。主要内容包括：蒙特卡罗积分的基本思想、马尔可夫链的基础理论、Metropolis-Hastings算法、 Gibbs采样算法以及收敛诊断方法。通过对这些内容的深入学习，读者将能够理解和掌握MCMC方法的核心思想，并具备将其应用于实际问题的能力。
19.2 蒙特卡罗积分
19.2.1 随机投点法
蒙特卡罗法（Monte Carlo Method）起源于20世纪40年代，由冯·诺依曼（John von Neumann）和乌拉姆（Stanislaw Ulam）在研究核物理问题时首次提出。蒙特卡罗法的基本思想是利用随机数进行数值计算，其名称来源于摩纳哥的著名赌城，暗示了方法的随机性特征。
考虑计算一个定积分：\(I = \int_a^b f(x) \, dx\)。假设\(f(x)\)在区间\([a,b]\)上有定义，我们可以将其转化为一概率期望问题。设随机变量\(X\)服从区间\([a,b]\)上的均匀分布，即\(X \sim U(a,b)\)，则\(f(X)\)的数学期望为：
由此可得：\(I = (b-a) \cdot E[f(X)]\)。当我们从均匀分布\(U(a,b)\)中独立抽取\(n\)个样本\(x_1, x_2, \ldots, x_n\)时，可以用样本均值近似期望：
这就是蒙特卡罗积分的基本估计量。根据大数定律，当\(n \to \infty\)时，\(\hat{I}_n\)依概率收敛到\(I\)。蒙特卡罗积分的误差来源于样本均值对期望的偏差，其标准误差约为\(O(n^{-1/2})\)，与积分维数无关，这是蒙特卡罗法相比传统数值积分方法的主要优势。
19.2.2 期望估计与重要性采样
在实际问题中，我们经常需要计算复杂分布的数学期望。设\(g(X)\)是某个随机变量\(X\)的函数，我们希望计算\(\theta = E_g[g(X)] = \int g(x) p(x) \, dx\)，其中\(p(x)\)是\(X\)的概率密度函数。如果直接从\(p(x)\)中采样困难，我们可以引入一个已知的提议分布（proposal distribution）\(q(x)\)，从中容易进行采样，从而得到重要性采样（Importance Sampling）估计。
假设\(q(x)\)与\(p(x)\)有相同的支撑集，且\(q(x) > 0\) whenever \(p(x) > 0\)，则：
设权重\(w(x) = \frac{p(x)}{q(x)}\)，称为重要性权重。从\(q(x)\)中抽取样本\(x_1, x_2, \ldots, x_n\)，则\(\theta\)的重要性采样估计为：
为了得到归一化的估计，可以令\(\tilde{w}_i = \frac{w(x_i)}{\sum_{j=1}^n w(x_j)}\)，则归一化估计为：
重要性采样的有效性依赖于提议分布\(q(x)\)的选择。一个好的\(q(x)\)应该与\(p(x)\)有相似的形状，并且使得权重\(w(x)\)的方差尽可能小。如果\(q(x)\)与\(p(x)\)差异太大，重要性采样可能会产生较大的估计误差，甚至导致估计量失效。
19.2.3 蒙特卡罗法的收敛性
蒙特卡罗积分具有坚实的统计理论基础。根据中心极限定理，对于独立同分布的样本，样本均值渐近服从正态分布：
其中方差\(\sigma^2 = \text{Var}(f(X)) = \int f^2(x) p(x) \, dx - \left(\int f(x) p(x) \, dx\right)^2\)。由此可以得到蒙特卡罗估计的置信区间：
蒙特卡罗法的收敛速度是\(O(n^{-1/2})\)，这意味着要提高一位精度（将标准误差减半），需要增加4倍的样本量。相比之下，高维数值积分的收敛速度通常为\(O(n^{-1/d})\)，其中\(d\)是维数。当\(d\)很大时，蒙特卡罗法的优势就非常明显了。
19.3 马尔可夫链基础
19.3.1 马尔可夫链的定义与性质
马尔可夫链（Markov Chain）是一种特殊的随机过程，具有"无记忆性"（memorylessness）的特征。设\(\{X_t\}_{t=0}^{\infty}\)是一个离散时间的随机序列，如果对任意时刻\(t\)和任意状态\(x_0, x_1, \ldots, x_{t+1}\)，满足：
则称\(\{X_t\}\)为一个马尔可夫链。通俗地说，马尔可夫链的下一个状态只与当前状态有关，与过去的状态无关，这一特性称为马尔可夫性质（Markov property）。
马尔可夫链的行为完全由两部分决定：状态空间（state space）和转移核（transition kernel）。设状态空间为\(\mathcal{X}\)，对于离散状态空间，转移核由转移概率矩阵\(P = (p_{ij})\)定义，其中\(p_{ij} = P(X_{t+1} = j | X_t = i)\)。对于连续状态空间，转移核由转移概率密度\(p(x'|x)\)定义，表示从状态\(x\)转移到状态\(x'\)的概率密度。
马尔可夫链的有限维分布由初始分布和转移核完全确定。设初始时刻的分布为\(\pi^{(0)}(x)\)，则在时刻\(t\)的边缘分布为：
或等价地，\(\pi^{(t)} = \pi^{(t-1)} \cdot P\)（离散情形）。
19.3.2 平稳分布与遍历定理
平稳分布（stationary distribution）是马尔可夫链理论中的核心概念。设\(\pi(x)\)是状态空间\(\mathcal{X}\)上的一个概率分布，如果对任意\(x \in \mathcal{X}\)满足：
或离散形式\(\pi_j = \sum_i \pi_i p_{ij}\)，则称\(\pi(x)\)为马尔可夫链的平稳分布。平稳分布的几何意义是：如果链在某个时刻服从平稳分布，那么在之后的任何时刻都保持这个分布不变。
平稳分布的存在性和唯一性取决于马尔可夫链的遍历性质。不可约（irreducible）和非周期（aperiodic）是保证平稳分布存在且唯一的两个关键条件。不可约意味着从任意状态出发，经过有限的转移步骤都可以到达任意其他状态；非周期意味着链不会陷入周期性的循环。
马尔可夫链的一个重要性质是遍历定理（Ergodic Theorem）：对于满足遍历条件的马尔可夫链，时间平均趋向于空间平均。具体而言，设\(\{X_t\}\)是一个遍历的马尔可夫链，平稳分布为\(\pi(x)\)，\(f\)是任意关于\(\pi\)可积的函数，则：
这一定理是MCMC方法的理论基础。它保证了我们可以通过在马尔可夫链上进行足够长时间的模拟，用样本均值来近似目标积分。
19.3.3细致平衡条件
细致平衡条件（Detailed Balance Condition）是构造满足给定平稳分布的马尔可夫链的重要工具。一个马尔可夫链如果满足细致平衡条件，则必定以某分布为平稳分布。
细致平衡条件要求对任意状态\(x, x' \in \mathcal{X}\)：
对上式两边对\(x'\)积分：
这正是平稳分布的定义。细致平衡条件给出了构造转移核的一个充分条件，许多MCMC算法（如Metropolis-Hastings算法）正是基于这一条件来设计转移机制的。
19.4 Metropolis-Hastings算法
19.4.1 算法原理
Metropolis-Hastings算法（简称M-H算法）是MCMC方法中最基本、最重要的算法之一，由Metropolis等人于1953年首次提出，后由Hastings于1970年推广到更一般的情形。
M-H算法的核心思想是构造一个满足细致平衡条件的马尔可夫链，使其平稳分布等于目标分布\(\pi(x)\)。设目标分布为\(\pi(x)\)（可以是归一化或未归一化的），提议分布为\(q(x'|x)\)，表示在当前状态\(x\)下提议转移到\(x'\)的概率。M-H算法的转移机制如下：
对于当前状态\(x\)：
1。 从提议分布\(q(x'|x)\)中抽取候选状态\(x^*\)； 2。 计算接受概率（acceptance probability）：
3。 以概率\(\alpha(x, x^*)\)接受候选状态，即\(x_{t+1} = x^*\)；以概率\(1 - \alpha(x, x^*)\)拒绝候选状态，即\(x_{t+1} = x\)。
可以验证，这样构造的马尔可夫链满足细致平衡条件：
因此\(\pi(x)\)是该马尔可夫链的平稳分布。
19.4.2 提议分布的选择
提议分布\(q(x'|x)\)的选择对M-H算法的效率有重要影响。常见的提议分布包括：
对称提议分布：如果\(q(x'|x) = q(x|x')\)，则接受概率简化为\(\alpha = \min\{1, \pi(x^*)/\pi(x)\}\)。常用的对称提议分布有： - 均匀随机游走：\(q(x'|x) \sim U(x - \delta, x + \delta)\)，在多维情况下通常选择球形或椭球形邻域； - 正态随机游走：\(q(x'|x) \sim N(x, \Sigma)\)，协方差矩阵\(\Sigma\)的选取影响采样效率。
独立提议分布：\(q(x'|x) = q(x')\)，即提议分布与当前状态无关。此时接受概率为：
独立提议分布要求与目标分布\(\pi(x)\)有相似的形状，否则接受率会很低。
自适应提议分布：在迭代过程中根据样本信息动态调整提议分布的参数，如协方差矩阵，以改进采样效率。
19.4.3 多维情况与随机游走
在高维参数空间中进行MCMC采样时，随机游走提议（Random Walk Proposal）是最常用的策略。设参数向量为\(\boldsymbol{\theta} \in \mathbb{R}^d\)，随机游走提议为：
其中\(\boldsymbol{\epsilon}_t\)通常服从均值为零、协方差矩阵为\(\Sigma\)的正态分布，即\(\boldsymbol{\epsilon}_t \sim N(\mathbf{0}, \Sigma)\)。
接受概率为：
由于正态分布是对称的，\(q(\boldsymbol{\theta}|\boldsymbol{\theta}^*) = q(\boldsymbol{\theta}^*|\boldsymbol{\theta})\)，因此提议分布在分子分母中相互抵消。
随机游走的步长选择对采样效率有重要影响。步长太大会导致接受率过低，大量候选状态被拒绝；步长太小会导致样本之间高度相关，有效样本量降低。一种实用的做法是调节步长使接受率维持在\(0.2\)至\(0.5\)之间，在高维情况下接受率可能更低。
19.5 Gibbs采样
19.5.1 条件分布采样
Gibbs采样是另一种重要的MCMC算法，由Geman和Geman于1984年提出，以物理学中的Gibbs分布命名。Gibbs采样特别适用于多维参数空间的情况，其基本思想是依次对每个参数在给定其他参数条件下的条件分布进行采样。
设\(\boldsymbol{\theta} = (\theta_1, \theta_2, \ldots, \theta_d)\)是\(d\)维参数向量，目标联合分布为\(p(\boldsymbol{\theta}) = p(\theta_1, \theta_2, \ldots, \theta_d)\)。Gibbs采样的一个迭代步骤如下：
对\(i = 1, 2, \ldots, d\)依次进行：
1。 从条件分布\(p(\theta_i | \theta_{-i})\)中抽取新的\(\theta_i\)，其中\(\theta_{-i}\)表示除\(\theta_i\)以外的所有其他参数当前值。
完整的Gibbs采样算法可以表示为：从某个初始值\(\boldsymbol{\theta}^{(0)}\)开始，重复以下步骤直到收敛：
这里的一个关键假设是条件分布\(p(\theta_i | \theta_{-i})\)具有解析形式或可以用其他MCMC方法容易地从中采样。
19.5.2 详细平衡与收敛性
Gibbs采样实际上是Metropolis-Hastings算法的一种特殊形式。考虑Gibbs采样的一个完整循环（更新所有\(d\)个参数），从状态\(\boldsymbol{\theta} = (\theta_1, \ldots, \theta_d)\)转移到\(\boldsymbol{\theta}' = (\theta_1', \ldots, \theta_d')\)的转移概率为：
可以验证，Gibbs采样满足细致平衡条件。考虑相邻的两个状态\(\boldsymbol{\theta}\)和\(\boldsymbol{\theta}'\)，它们只在第\(i\)个分量上不同。不失一般性，设\(\boldsymbol{\theta}'\)在第\(i\)个分量上变为\(\theta_i'\)，其他分量相同，则：
这是因为联合分布可以分解为\(p(\boldsymbol{\theta}) = p(\theta_i | \boldsymbol{\theta}_{-i}) p(\boldsymbol{\theta}_{-i})\)，且\(p(\boldsymbol{\theta}_{-i}) = p(\boldsymbol{\theta}_{-i}')\)。
因此，Gibbs采样构造的马尔可夫链以目标分布\(p(\boldsymbol{\theta})\)为平稳分布。根据遍历定理，样本轨道的时间平均收敛到空间平均，Gibbs采样的样本均值可以用于估计感兴趣的期望。
19.5.3 Gibbs采样的优势与局限
Gibbs采样相比一般的M-H算法具有以下优势：
第一，接受率始终为1，不存在候选被拒绝的问题，因此采样效率较高。当条件分布具有标准形式（如共轭分布）时，采样特别方便。
第二，条件分布采样通常比联合分布采样更容易实现。在许多层次化贝叶斯模型中，条件分布具有简洁的解析形式。
第三，样本之间具有较强的相关性但具有明确的结构，这使得诊断收敛性相对容易。
然而，Gibbs采样也有其局限性：首先，需要能够从所有条件分布中进行有效采样，如果某些条件分布形式复杂，可能需要嵌套其他MCMC方法；其次，当参数之间相关性较强时，Gibbs采样收敛可能很慢，存在"维度灾难"问题；第三，条件分布的形式取决于先验分布和似然函数的选择，不是所有模型都能得到方便的条件分布。
19.6 收敛诊断
19.6.1 收敛判断的困难
MCMC采样的一个核心挑战是如何判断马尔可夫链已经收敛到平稳分布，从而样本可以用于统计推断。在实际应用中，我们只能观察到有限长度的样本序列，无法直接验证马尔可夫链是否已经达到平稳分布。判断收敛是MCMC应用中最困难的问题之一。
收敛诊断的本质是检测链的" burn-in"阶段（预烧期），即链从初始值到达平稳分布所需的迭代次数。在burn-in期间，样本分布可能与目标分布有显著差异，不应纳入最终分析。
造成收敛困难的原因是多方面的：马尔可夫链可能存在多个平稳分布（不可约性问题）；从某些初始值出发，链可能需要很长时间才能到达高概率区域；参数空间可能存在"挡光"（blocking）结构，使得不同区域之间的转移困难。
19.6.2 Geweke检验
Geweke检验是一种基于频谱方法的收敛诊断技术，由Geweke于1991年提出。其基本思想是比较链开始部分和结尾部分样本的均值是否存在显著差异。
设从马尔可夫链中抽取了\(N\)个样本，取前\(n_1\)个样本和后\(n_2\)个样本分别计算均值：
Geweke统计量为：
其中\(\hat{S}_{\theta}^A(0)\)和\(\hat{S}_{\theta}^B(0)\)分别是两段样本的谱密度估计在频率零处的值。如果链已收敛，两段样本均值应大致相等，\(Z\)统计量应渐近服从标准正态分布。若\(|Z|\)值较大（通常以\(|Z| > 1.96\)为判据），则表明两段样本存在显著差异，链可能未收敛。
Geweke检验适用于检验单个参数或参数线性组合的收敛性，但不能用于检测多维参数的整体收敛性。
19.6.3 Gelman-Rubin诊断
Gelman-Rubin诊断（又称R-hat统计量）是一种基于多个马尔可夫链的收敛诊断方法，由Gelman和Rubin于1992年提出。该方法通过比较链内方差和链间方差来判断收敛性。
设运行了\(m\)条独立的马尔可夫链，每条链有\(n\)个样本（去除burn-in后）。对某个参数\(\theta\)，定义：
- 链内方差：\(W = \frac{1}{m} \sum_{j=1}^m s_j^2\)，其中\(s_j^2 = \frac{1}{n-1} \sum_{t=1}^n (\theta_{jt} - \bar{\theta}_j。)^2\)是第\(j\)条链的样本方差；
- 链间方差：\(B = \frac{n}{m-1} \sum_{j=1}^m (\bar{\theta}_j。 - \bar{\theta}_{...})^2\)，其中\(\bar{\theta}_j。\)是第\(j\)条链的均值，\(\bar{\theta}_{...}\)是所有链的总均值；
- 目标分布方差的无偏估计：\(\hat{\text{Var}}(\theta) = \frac{n-1}{n} W + \frac{1}{n} B\)。
R-hat统计量定义为：
当所有链都已收敛到同一平稳分布时，\(\hat{R}\)应接近1。通常认为\(\hat{R} < 1.1\)或\(\hat{R} < 1.2\)时链已收敛，可以开始采样。\(\hat{R}\)越大，表明链间方差越大，链可能尚未收敛或存在混合不良的问题。
19.7 公式汇总表
为了便于读者查阅和复习，本章涉及的主要公式汇总如下：
| 序号 | 公式名称 | 公式表达 | 说明 |
|---|---|---|---|
| 1 | 蒙特卡罗积分估计 | 均匀分布下定积分的蒙特卡罗估计 | |
| 2 | 重要性采样估计 | 基于提议分布\(q(x)\)的重要性采样 | |
| 3 | 细致平衡条件 | 马尔可夫链以\(\pi\)为平稳分布的充分条件 | |
| 4 | 平稳分布定义 | 马尔可夫链平稳分布的积分方程形式 | |
| 5 | M-H接受概率 | Metropolis-Hastings算法的接受准则 | |
| 6 | 对称提议的M-H接受率 | 当$q(x^* | |
| 7 | Gibbs条件分布采样 | Gibbs采样中第\(i\)个分量的更新 | |
| 8 | Geweke统计量 | 比较链首尾均值的收敛诊断 | |
| 9 | Gelman-Rubin R-hat | 基于多条链方差的收敛诊断统计量 | |
| 10 | 遍历定理 | 遍历马尔可夫链的时间平均收敛性 |
参考文献
1。 Metropolis, N。, Rosenbluth, A。 W。, Rosenbluth, M。 N。, Teller, A。 H。, & Teller, E。 (1953)。 Equations of State Calculations by Fast Computing Machines。 Journal of Chemical Physics, 21(6), 1087-1092。
2。 Hastings, W。 K。 (1970)。 Monte Carlo Sampling Methods Using Markov Chains and Their Applications。 Biometrika, 57(1), 97-109。
3。 Geman, S。, & Geman, D。 (1984)。 Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images。 IEEE Transactions on Pattern Analysis and Machine Intelligence, 6(6), 721-741。
4。 Gilks, W。 R。, Richardson, S。, & Spiegelhalter, D。 J。 (1996)。 Markov Chain Monte Carlo in Practice。 Chapman and Hall/CRC。
5。 Gelman, A。, Carlin, J。 B。, Stern, H。 S。, Dunson, D。 B。, Vehtari, A。, & Rubin, D。 B。 (2013)。 Bayesian Data Analysis (3rd ed。)。 CRC Press。
6。 Robert, C。 P。, & Casella, G。 (2004)。 Monte Carlo Statistical Methods (2nd ed。)。 Springer。