概率采样：基本采样与马尔可夫链—蒙特卡洛采样

结构化分布与序列模型讨论了联合分布的分解、条件独立和隐马尔可夫模型。给定这些概率模型后，推断仍可能卡在一个具体的计算问题上：后验分布已经写出，却无法解析求出积分；联合分布可以逐点计算，却难以列举高维空间中的全部状态。

概率采样用一组随机样本近似这些积分。若目标分布容易抽样，可以生成相互独立的样本；若只能计算未归一化密度，可以借助拒绝采样、重要性采样或马尔可夫链。不同方法产生的样本形式不同，误差也不能用同一种方式判断。本文围绕同一个期望估计问题，说明各类方法成立所需的条件，以及采样结果应怎样检查。

蒙特卡洛估计#

设目标分布的密度或概率质量函数为 $\pi(x)$，希望计算函数 $h(X)$ 的期望

$$I =\mathbb{E}_{\pi}[h(X)] =\int h(x)\pi(x)\,\mathrm{d}x.$$

离散变量时把积分换成求和。许多概率推断任务都能写成这种形式。令 $h(x)=x$ 可以求后验均值；令 $h(x)=\mathbb{1}(x\in A)$ 可以求事件 $A$ 的概率；令 $h(x)=\ell(x)$ 可以计算某个决策在后验分布下的期望损失。

样本均值#

若能从 $\pi$ 独立同分布地取得 $N$ 个样本 $X^{(1)},\ldots,X^{(N)}$，蒙特卡洛估计量为

$$\widehat{I}_N =\frac{1}{N}\sum_{n=1}^{N}h\!\left(X^{(n)}\right).$$

只要 $\mathbb{E}_{\pi}[|h(X)|]$ 有限，大数定律就保证 $\widehat{I}_N$ 随样本数增加而收敛到 $I$。独立样本下，这个估计量无偏；若 $\sigma_h^2=\operatorname{Var}_{\pi}[h(X)]$ 有限，则

$$\operatorname{Var}(\widehat{I}_N)=\frac{\sigma_h^2}{N}.$$

误差的典型尺度因此按 $N^{-1/2}$ 下降。想把标准误缩小一半，通常需要把样本数增加到四倍。这个收敛速度不显式依赖积分维数，是蒙特卡洛方法能够处理高维积分的重要原因；维数仍会通过抽样难度、函数方差和样本相关性影响实际成本。

蒙特卡洛标准误#

独立样本下，可以用样本方差

$$s_h^2 =\frac{1}{N-1}\sum_{n=1}^{N} \left(h(X^{(n)})-\widehat{I}_N\right)^2$$

估计 $\sigma_h^2$，对应的蒙特卡洛标准误（Monte Carlo standard error，MCSE）为

$$\widehat{\operatorname{MCSE}}(\widehat{I}_N) =\sqrt{\frac{s_h^2}{N}}.$$

在有限方差和常见正则条件下，中心极限定理给出近似正态误差。MCSE 衡量的是有限次随机计算带来的不确定性，与模型是否正确、数据是否有偏和后验是否集中属于不同问题。增加采样次数可以降低 MCSE，无法修复错误的似然、先验或数据生成假设。

基本采样#

基本采样指直接利用已知分布的抽样规则生成独立样本。计算机通常先产生近似独立的均匀伪随机数，再通过变换、分解或专门算法得到目标分布样本。每次用新的随机数重复完整过程时，所得样本之间不需要烧入期，也没有接受率和链混合问题。

逆变换采样#

设一维连续分布的累积分布函数为

$$F(x)=P(X\leq x).$$

若 $U\sim\operatorname{Uniform}(0,1)$，则广义逆函数

$$X=F^{-1}(U)$$

服从目标分布。以指数分布为例，$F(x)=1-\exp(-\lambda x)$，因此

$$X=-\frac{1}{\lambda}\log(1-U).$$

由于 $1-U$ 仍服从 $(0,1)$ 上的均匀分布，实现中也常写成 $X=-\log U/\lambda$。离散分布可以把 $U$ 落入的累积概率区间映射到相应类别。类别很多且概率固定时，还可以预先构造别名表，将单次抽样的查找成本降到常数级。

逆变换需要能够计算或数值求解 $F^{-1}$。多元分布通常没有一个可直接使用的联合逆函数，实际算法会利用条件分布、矩阵分解或分布之间的构造关系。

变换采样#

考虑一个硬币参数的贝叶斯更新。令正面概率为 $\theta$，先验为 $\operatorname{Beta}(1,1)$。观察到 7 次正面和 3 次反面后，后验为

$$\pi(\theta) =\operatorname{Beta}(\theta\mid 8,4), \qquad 0<\theta<1,$$

其未归一化密度为

$$\widetilde{\pi}(\theta)=\theta^7(1-\theta)^3.$$

这个后验的期望有解析值

$$\mathbb{E}[\theta\mid\mathcal{D}] =\frac{8}{8+4} =\frac{2}{3}.$$

已知 Beta 分布与 Gamma 分布的关系时，可以先独立生成

$$G_1\sim\operatorname{Gamma}(8,1), \qquad G_2\sim\operatorname{Gamma}(4,1),$$

再令

$$\theta=\frac{G_1}{G_1+G_2}.$$

每次重复这组操作都会得到一个独立的 Beta 样本，样本均值应在 MCSE 所描述的随机范围内接近 $2/3$。这个例子有解析答案，适合检查算法；实际推断中，目标密度往往不再属于可以直接抽样的标准分布。

拒绝采样#

拒绝采样（rejection sampling）用一个容易抽样的提议分布 $q(x)$ 包住目标分布。设目标只知道未归一化形式 $\widetilde{\pi}(x)$，归一化常数为

$$Z=\int\widetilde{\pi}(x)\,\mathrm{d}x, \qquad \pi(x)=\frac{\widetilde{\pi}(x)}{Z}.$$

若存在有限常数 $M$，使目标支持集上的所有 $x$ 都满足

$$\widetilde{\pi}(x)\leq Mq(x),$$

则一次候选抽样包含两步：

1. 生成 $X\sim q$ 和 $U\sim\operatorname{Uniform}(0,1)$；
2. 当 $U\leq\widetilde{\pi}(X)/(Mq(X))$ 时接受 $X$，否则丢弃并重新生成。

被接受样本的密度正比于

$$q(x)\frac{\widetilde{\pi}(x)}{Mq(x)} =\frac{\widetilde{\pi}(x)}{M},$$

归一化后恰好等于 $\pi(x)$。每次从头独立提议，因此接受的样本仍相互独立。平均接受率为 $Z/M$；提议分布越贴近目标形状，包络空出的面积越少，所需候选数也越少。

在 Beta 后验例子中，可取 $q(\theta)=1$ 的均匀分布。$\widetilde{\pi}(\theta)$ 在 $\theta=0.7$ 处达到最大值，因此可令

$$M=0.7^7\times0.3^3.$$

又因为 $Z=B(8,4)=1/1320$，理论接受率约为 $34\%$。若把同样做法直接推广到高维，提议分布需要同时覆盖每个方向，$Z/M$ 往往随维数迅速下降。目标呈重尾、包含多个相距较远的峰，或者无法给出全局上界时，拒绝采样也难以实施。

重要性采样#

重要性采样（importance sampling）同样从容易抽样的 $q(x)$ 生成样本，但不丢弃候选。若在 $h(x)\pi(x)\neq 0$ 的区域都有 $q(x)>0$，则

$$I =\int h(x)\pi(x)\,\mathrm{d}x =\int h(x)\frac{\pi(x)}{q(x)}q(x)\,\mathrm{d}x =\mathbb{E}_q[h(X)w(X)],$$

其中

$$w(x)=\frac{\pi(x)}{q(x)}$$

称为重要性权重。对 $X^{(n)}\overset{\mathrm{i.i.d.}}{\sim}q$，估计量

$$\widehat{I}_{\mathrm{IS}} =\frac{1}{N}\sum_{n=1}^{N} h(X^{(n)})w(X^{(n)})$$

在相应期望存在时是无偏的。它的方差取决于 $h(x)\pi(x)$ 与 $q(x)$ 的匹配程度。某些区域在目标下很重要、在 $q$ 下却极少被抽到时，少数巨大权重会控制整个结果，样本数很大也可能不稳定。

自归一化重要性采样#

贝叶斯后验常只知道 $\widetilde{\pi}(x)$，无法计算 $Z$。此时使用未归一化权重

$$\widetilde{w}_n =\frac{\widetilde{\pi}(X^{(n)})}{q(X^{(n)})}$$

构造自归一化估计量

$$\widehat{I}_{\mathrm{SNIS}} =\frac{\sum_{n=1}^{N}\widetilde{w}_n h(X^{(n)})} {\sum_{n=1}^{N}\widetilde{w}_n}.$$

分子与分母中的 $Z$ 相互抵消。这个比值在有限样本下一般有偏，但在常见矩条件下是一致的。对前面的 Beta 后验，若仍取 $q(\theta)=1$，只需令 $\widetilde{w}_n=\theta_n^7(1-\theta_n)^3$，再计算加权平均。与拒绝采样相比，所有候选都进入计算；这不表示每个候选贡献相同，远离后验高密度区域的样本权重会很小。

有效样本量#

归一化权重为

$$\overline{w}_n =\frac{\widetilde{w}_n}{\sum_{j=1}^{N}\widetilde{w}_j}.$$

常用的权重有效样本量为

$$N_{\mathrm{eff}}^{\mathrm{IS}} =\frac{1}{\sum_{n=1}^{N}\overline{w}_n^2} =\frac{\left(\sum_n\widetilde{w}_n\right)^2} {\sum_n\widetilde{w}_n^2}.$$

所有权重相等时它等于 $N$；一个样本几乎占据全部权重时，它接近 1。这个量是权重集中程度的诊断，不是对任意 $h$ 都精确成立的独立样本等价数。选择提议分布时，除了覆盖目标支持集，还应让尾部足够重，并让概率质量集中在 $|h(x)|\pi(x)$ 对积分贡献较大的区域。

马尔可夫链#

直接采样、拒绝采样和重要性采样都以独立提议为起点。马尔可夫链—蒙特卡洛（Markov chain Monte Carlo，MCMC）改用一条具有目标平稳分布的随机序列。样本之间允许相关，只要链经过足够长的演化后按 $\pi$ 访问状态空间，沿链计算的平均值仍可在适当条件下收敛到目标期望。

转移核#

马尔可夫链 $X^{(0)},X^{(1)},\ldots$ 满足

$$P\!\left(X^{(t+1)}\in A \mid X^{(0)},\ldots,X^{(t)}\right) =P\!\left(X^{(t+1)}\in A\mid X^{(t)}\right).$$

给定当前状态后，下一步的条件分布不再依赖更早的完整路径。连续状态空间中，这个条件分布写成转移核 $K(x,A)$；有限状态空间中可以写成转移矩阵 $P$，其中 $P_{ij}=P(X^{(t+1)}=j\mid X^{(t)}=i)$。

例如两状态链的转移矩阵为

$$P= \begin{bmatrix} 0.9 & 0.1\\ 0.2 & 0.8 \end{bmatrix}.$$

链在状态 1 时有 $0.1$ 的概率转到状态 2，在状态 2 时有 $0.2$ 的概率转回状态 1。较大的自转移概率会形成连续重复状态，这正是相关样本的一个简单来源。

平稳分布#

若分布 $\pi$ 经过一步转移后保持不变，它就是链的平稳分布（stationary distribution）：

$$\pi(A)=\int K(x,A)\pi(x)\,\mathrm{d}x.$$

有限状态下写成行向量形式为

$$\boldsymbol{\pi}=\boldsymbol{\pi}P.$$

上面的两状态链满足 $\boldsymbol{\pi}=(2/3,1/3)$。若初始状态本来就按这个分布生成，之后每个时刻的边缘分布都相同；相邻状态仍可能相关。MCMC 的设计目标是构造一个以所需 $\pi$ 为平稳分布、又能在状态空间中有效移动的转移核。

详细平衡#

详细平衡（detailed balance）要求任意两个状态之间的平稳概率流相等：

$$\pi(x)K(x,\mathrm{d}x') =\pi(x')K(x',\mathrm{d}x).$$

对所有 $x'$ 积分后即可得到平稳性，因此详细平衡是保证平稳分布的一项充分条件。满足它的链称为可逆链。平稳链不一定可逆，详细平衡也不是构造 MCMC 的唯一途径；Metropolis–Hastings 算法采用它，是因为局部接受概率可以直接让正反方向的概率流配平。

遍历性#

拥有目标平稳分布还不够。若链被困在互不连通的区域，从某个初始值出发就无法访问整个目标支持集。有限状态链中，不可约性保证任意状态之间存在可达路径，非周期性排除只能按固定周期返回的振荡；两者成立时，链具有唯一平稳分布，并从任意初始分布收敛到它。一般状态空间需要相应的不可约性、再生性或 Harris 遍历条件，不能只凭转移公式外观判断收敛。

从任意初始值开始时，链的早期分布通常不等于 $\pi$。实践中会丢弃一段预热或烧入（warm-up/burn-in）样本，但固定长度的烧入无法证明链已经收敛。链若没有跨过相隔很远的众数，继续丢弃前缀只会隐藏问题。

自相关#

平稳链上令 $Y_t=h(X^{(t)})$，滞后 $k$ 的自相关为

$$\rho_k =\operatorname{Corr}(Y_t,Y_{t+k}).$$

在自相关可求和等条件下，样本均值的渐近方差可写成

$$\operatorname{Var}(\widehat{I}_N) \approx \frac{\sigma_h^2}{N} \left(1+2\sum_{k=1}^{\infty}\rho_k\right).$$

括号中的量称为积分自相关时间，对应的 MCMC 有效样本量为

$$N_{\mathrm{eff}}^{\mathrm{MCMC}} =\frac{N}{1+2\sum_{k=1}^{\infty}\rho_k}.$$

它针对具体函数 $h$ 和链的相关结构，与重要性采样中由权重集中程度定义的有效样本量含义不同。直接把相关的 MCMC 状态代入独立样本标准误，会系统性低估随机误差。

Metropolis–Hastings 算法#

Metropolis–Hastings（MH）算法只要求能够计算目标密度的相对值。设当前状态为 $x$，从提议核 $q(x'\mid x)$ 生成候选 $x'$，再以概率

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

接受候选。若拒绝，则令下一状态仍为 $x$。目标归一化常数在比值中抵消，这使 MH 适合只知道后验“先验乘似然”形式的情形。

接受概率把两个方向的概率流都压到较小的一侧：

$$\pi(x)q(x'\mid x)\alpha(x,x') =\min\{\pi(x)q(x'\mid x),\pi(x')q(x\mid x')\}.$$

交换 $x$ 与 $x'$ 后右侧不变，因此转移满足详细平衡。拒绝时保留当前状态也属于转移核的一部分；若只保存接受的候选，通常会改变目标分布，因为不同位置的拒绝概率不同。

随机游走 Metropolis#

随机游走提议常写成

$$x'=x+\varepsilon, \qquad \varepsilon\sim\mathcal{N}(0,\Sigma_q).$$

对称提议满足 $q(x'\mid x)=q(x\mid x')$，接受概率化为

$$\alpha(x,x') =\min\left\{1, \frac{\widetilde{\pi}(x')}{\widetilde{\pi}(x)} \right\}.$$

$\Sigma_q$ 太小时，候选几乎都能接受，但链只能缓慢移动；尺度太大时，候选经常落到低密度区而被拒绝。接受率本身不是目标，关键是单位计算成本能产生多少有效样本。文献中的某些最优接受率来自高维独立同分布目标和特定随机游走极限，不能当作所有模型的固定阈值。

参数受正值、区间或单纯形约束时，直接加高斯噪声可能频繁越界。常见做法是在无约束变换后的空间采样，并把变换的 Jacobian 行列式计入目标密度；也可以设计反射、截断或问题专用提议，但必须在 MH 比值中使用实际的正反提议概率。

独立 Metropolis#

独立 Metropolis 使用与当前状态无关的提议 $q(x')$，接受概率为

$$\alpha(x,x') =\min\left\{1, \frac{\widetilde{\pi}(x')q(x)} {\widetilde{\pi}(x)q(x')} \right\}.$$

它与重要性采样都依赖一个全局提议分布。$q$ 覆盖目标且形状接近时，链可以跨越较远区域；$q$ 在目标尾部过轻时，链进入某些高权重区域后可能很难离开。与拒绝采样不同，MH 不要求已知全局包络常数，但生成的是相关样本。

采样诊断#

单条链的轨迹看起来有波动，不足以证明它已经代表目标分布。较稳妥的检查同时使用多条从分散初值开始的链，观察轨迹是否访问相同区域，并报告有效样本量、MCSE 和链间诊断量。$\widehat{R}$ 接近 1 只说明链间与链内变异较为一致；所有链若都困在同一个局部众数，它也可能给出过于乐观的结果。

诊断应针对实际报告的量。参数均值混合较快，不代表尾部概率、极端分位数或两个众数之间的相对质量也估计准确。增加迭代次数只在链能够遍历目标区域时有效；几何结构不合适时，需要重参数化、改进提议、分块更新或使用其他采样算法。

Gibbs 采样#

Gibbs 采样把多元目标分布拆成一组完整条件分布。设

$$x=(x_1,x_2,\ldots,x_D),$$

一次系统扫描依次生成

$$x_1^{(t+1)} \sim\pi(x_1\mid x_2^{(t)},\ldots,x_D^{(t)}),$$ $$x_2^{(t+1)} \sim\pi(x_2\mid x_1^{(t+1)},x_3^{(t)},\ldots,x_D^{(t)}),$$

直到

$$x_D^{(t+1)} \sim\pi(x_D\mid x_1^{(t+1)},\ldots,x_{D-1}^{(t+1)}).$$

每个坐标更新都从目标的精确条件分布抽样，因此可以看作接受概率恒为 1 的 MH 更新。单个更新保持 $\pi$ 不变，多个更新的组合也保持 $\pi$ 不变。系统扫描形成的整轮转移未必可逆，但平稳性不依赖整轮详细平衡。

以相关系数为 $\rho$ 的二维标准高斯为例，完整条件分布为

$$X\mid Y=y \sim\mathcal{N}(\rho y,1-\rho^2),$$ $$Y\mid X=x \sim\mathcal{N}(\rho x,1-\rho^2).$$

这两个一维高斯都容易直接抽样。$|\rho|$ 接近 1 时，条件方差 $1-\rho^2$ 很小，坐标更新会沿狭长联合分布缓慢前进，样本自相关随之升高。接受率为 1 并不等于混合良好。

块 Gibbs 采样#

若一组高度相关的变量可以从联合条件分布中一起抽样，可以把它们合成一个块。块更新减少逐坐标移动造成的慢混合，但条件分布的维数和抽样成本会增加。层次贝叶斯模型中，中心化与非中心化参数化也会改变条件相关性，因而显著影响 Gibbs 链的效率。

Gibbs 采样要求完整条件分布能够准确抽样。只知道条件密度的未归一化形式时，可以在某个坐标更新内部再使用 MH，形成 Metropolis-within-Gibbs。此时需要同时检查每个内部更新的遍历性与整条链的混合情况。

祖先采样#

祖先采样（ancestral sampling）利用有向无环图的生成顺序。若贝叶斯网络按拓扑顺序写成

$$p(x_1,\ldots,x_D) =\prod_{i=1}^{D}p(x_i\mid x_{\operatorname{pa}(i)}),$$

便可以从没有父节点的变量开始，再依次从每个局部条件分布抽样：

$$x_i\sim p(x_i\mid x_{\operatorname{pa}(i)}).$$

当处理 $X_i$ 时，它的父节点已经取得数值。把所有节点走完一次，就得到联合分布的一个精确样本；重新使用独立随机数重复整套过程，可以得到相互独立的联合样本。祖先采样不需要马尔可夫链、接受概率或烧入，它依赖的是 DAG 分解中的每个局部条件分布都已归一化且能够抽样。

隐马尔可夫模型的前向生成就是祖先采样。先生成

$$Z_1\sim p(z_1), \qquad X_1\sim p(x_1\mid z_1),$$

之后对 $t=2,\ldots,T$ 依次生成

$$Z_t\sim p(z_t\mid z_{t-1}), \qquad X_t\sim p(x_t\mid z_t).$$

这些样本适合检查模型能够生成怎样的序列，也可用于估计联合分布下的期望。若已经观察到 $X_{1:T}=x_{1:T}$，目标会变成后验 $p(z_{1:T}\mid x_{1:T})$。把观测节点固定后继续按原生成方向抽样，通常不会得到这个后验；观测会把信息传向祖先和同层变量，需要使用前向后向、条件采样、重要性方法或 MCMC 等推断算法。

采样方法#

这些方法的差别可以归结为目标分布提供了哪些可用操作，以及所得样本是否独立：

若目标分布可以直接抽样，独立精确样本通常最容易分析。只能计算未归一化密度时，低维且有紧包络的目标可用拒绝采样；存在合适全局提议时，可用重要性采样估计多个期望或归一化常数；难以构造全局提议时，MCMC 用局部转移换取渐近正确性。概率图模型还提供两种结构化选择：完整条件分布易抽样时使用 Gibbs，按生成方向模拟联合分布时使用祖先采样。

采样程序输出大量数值，不等于积分已经准确。独立样本需要报告样本方差和 MCSE，重要性采样需要检查最大权重与权重有效样本量，MCMC 还要检查多链一致性、自相关、有效样本量和目标区域覆盖情况。只有误差诊断与目标期望对应，随机样本才真正成为可解释的数值估计。

参考文献

1. N. Metropolis and S. Ulam. The Monte Carlo Method. Journal of the American Statistical Association, 1949. DOI: 10.1080/01621459.1949.10483310
2. H. Kahn and A. W. Marshall. Methods of Reducing Sample Size in Monte Carlo Computations. Journal of the Operations Research Society of America, 1953. DOI: 10.1287/opre.1.5.263
3. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller. Equation of State Calculations by Fast Computing Machines. The Journal of Chemical Physics, 1953. DOI: 10.1063/1.1699114
4. W. K. Hastings. Monte Carlo Sampling Methods Using Markov Chains and Their Applications. Biometrika, 1970. DOI: 10.1093/biomet/57.1.97
5. S. Geman and D. Geman. Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1984. DOI: 10.1109/TPAMI.1984.4767596
6. L. Tierney. Markov Chains for Exploring Posterior Distributions. The Annals of Statistics, 1994. DOI: 10.1214/aos/1176325750
7. C. P. Robert and G. Casella. Monte Carlo Statistical Methods. Springer Texts in Statistics, 2004. DOI: 10.1007/978-1-4757-4145-2
8. A. Vehtari, A. Gelman, D. Simpson, B. Carpenter, and P.-C. Bürkner. Rank-Normalization, Folding, and Localization: An Improved R-hat for Assessing Convergence of MCMC. Bayesian Analysis, 2021. DOI: 10.1214/20-BA1221