14.4.1 不变分布、时间平均与算法目标
如果目标分布 \(\pi\) 容易计算比例,却难以直接生成独立样本,可以构造以 \(\pi\) 为平稳分布的马尔可夫链 \(X_0,X_1,\ldots\),再用时间平均估计 \(\pi(h)=\sum_x\pi(x)h(x)\)。本节先在有限状态空间上给出完整推导,连续情形的更新公式随后说明。第12.1 节至12.4 节提供转移矩阵、常返性、可逆性与耦合的基础。
不变性是 \(\pi P=\pi\),它说明从 \(\pi\) 出发时每一时刻的边缘仍为 \(\pi\)。分布收敛是从其他初始分布出发时 \(\mathcal L(X_n)\to\pi\)。时间平均收敛则是 \(n^{-1}\sum_{k=0}^{n-1}h(X_k)\to\pi(h)\)。不变性本身不足以推出后两项;例如恒等转移 \(P=I\) 保持任何分布,却从不离开初始状态。
平稳也不等于独立。即便 \(X_0\sim\pi\),相邻输出仍可能高度相关。马尔可夫链蒙特卡洛(Markov chain Monte Carlo,MCMC)的误差分析必须处理这种依赖,不能把样本数量直接代入上一节的 iid 标准误公式。
14.4.2 Metropolis–Hastings 接受率从何而来
令 \(S_+=\{x:\pi(x)>0\}\),从其中一个状态开始。提议转移为 \(q(x,y)\)。当前位于 \(x\) 时先抽 \(Y\sim q(x,\cdot)\),再取一个独立均匀数决定是否接受;拒绝时仍记录当前状态,不能只保留被接受的状态。
对实际可能提出的 \(y\),即 \(q(x,y)>0\),取
当 \(\pi(y)=0\) 或 \(q(y,x)=0\) 时拒绝这个提议。\(q(x,y)=0\) 的点不会被提出,可把接受率定义为零。这样避免了用没有定义的 \(0/0\) 作判断。若 \(\pi=\widetilde\pi/Z\),未知常数 \(Z\) 在比值中约去。
对 \(x\ne y\),转移概率为 \(P(x,y)=q(x,y)\alpha(x,y)\),其余概率放在 \(P(x,x)=1-\sum_{y\ne x}P(x,y)\)。因为 \(0\le\alpha\le1\),这是合法转移矩阵。并且
对角项自动满足同一等式。固定 \(y\) 并对 \(x\) 求和,得 \(\sum_x\pi(x)P(x,y)=\pi(y)\sum_xP(y,x)=\pi(y)\),所以 \(\pi\) 不变,且链可逆。上述最小值形式说明接受率调节的是两个方向的平稳概率流。
若提议对称,接受率简化为 \(\min\{1,\pi(y)/\pi(x)\}\);非对称提议不能省掉 \(q(y,x)/q(x,y)\)。有限正支撑上,只要双向提议都为正的边组成连通图,接受后的链就不可约。非周期性可另外检查自环;也可将转移替换为 \((I+P)/2\),但增加停留会改变统计效率。
连续空间中,以目标密度和提议密度 \(q(y\mid x)\) 替代质量函数,仍使用同样的比值。接受部分满足对称密度恒等式,拒绝部分则是当前点的点质量。目标密度为零的提议仍拒绝;在目标的零测例外处选择版本不影响不变性。一般状态空间的收敛还需不可约性、Harris 常返等条件,单靠这个代数恒等式不能保证从任意起点收敛。
14.4.3 一个可手算的两状态算法
在状态 \(0,1\) 上取 \(\pi=(3/4,1/4)\),每次提议翻转状态。由式 (14.4-1),从 \(0\) 到 \(1\) 的接受率为 \(1/3\),反方向为 \(1\),所以
两方向的平稳流都为 \(1/4\)。设 \(a_n=P(X_n=1)\),则 \(a_{n+1}=(1-a_n)/3\),解得 \(a_n-1/4=(-1/3)^n(a_0-1/4)\)。例如从 \(X_0=0\) 出发,\(a_1=1/3\)、\(a_2=2/9\)、\(a_3=7/27\),它们围绕 \(1/4\) 交替趋近。平稳分布已经正确,并不意味着第一步就得到目标分布。
若只保留接受的翻转,记录就成为 \(0,1,0,1,\ldots\),长期比例变成 \(1/2,1/2\)。这不是原链的目标分布。保留拒绝时的重复值,才能让停留时间正确参与平均。
14.4.4 Gibbs 条件更新为何保持联合分布
设目标是两变量的联合分布 \(\pi(x,y)\),并能从完整条件分布 \(\pi(x\mid y)\)、\(\pi(y\mid x)\) 抽样。Gibbs 更新保持一个坐标不动,按另一个坐标的完整条件分布重新生成。零边缘质量的条件值可选任意版本,但它们不影响平稳测度下的不变性。
若先更新 \(x\),从平稳联合分布出发得到新状态 \((x',y)\) 的概率为
因此单个更新核 \(K_x\) 保持 \(\pi\),\(K_y\) 也一样。顺序更新的核 \(K_xK_y\) 满足 \(\pi K_xK_y=(\pi K_x)K_y=\pi\);以固定概率随机选择坐标的混合核也不变。连续情形把求和改为对条件分布核积分,证明相同。单坐标核可逆,固定顺序的复合一般未必可逆。
令 \(\pi(0,0)=\pi(1,1)=0.4\),\(\pi(0,1)=\pi(1,0)=0.1\)。两个边缘都公平,但给定另一个坐标时,取相同值的条件概率为 \(0.4/(0.4+0.1)=0.8\)。从 \((0,0)\) 先更新 \(x\)、再更新 \(y\),到 \((1,1)\) 的概率为 \(0.2\times0.8=0.16\)。这一整轮的四个终点概率依次为 \(0.64,0.16,0.04,0.16\),和为 \(1\),却不等于目标概率表。仍需运行多步。
若目标仅在 \((0,0)\) 和 \((1,1)\) 上各放 \(1/2\) 的质量,单坐标 Gibbs 更新会永远固定在初始点。完整条件分布完全正确,整个链仍不可约性失败。这也是检查目标支撑结构的必要性所在。
14.4.5 有限链的平均收敛与相关误差
有限不可约链有唯一平稳分布 \(\pi\)。对任何实函数 \(h\) 和任意初始分布,\(n^{-1}\sum_{k=0}^{n-1}h(X_k)\to\pi(h)\) 几乎必然。这个结论不要求非周期;若还要 \(P^n(x,\cdot)\to\pi\),则需要非周期性。
记 \(f=h-\pi(h)\)。在满足 \(\pi(g)=0\) 的有限维空间上,\(I-P\) 可逆:若 \(Pg=g\),在 \(g\) 的最大值点,相邻可达点必须也取最大值;不可约性使 \(g\) 恒定,再由均值为零知 \(g=0\)。因此存在有界的 \(g\) 使 \(g-Pg=f\)。令
各增量对过去的条件均值为零,且一致有界。于是 \(\sum_{k=1}^n(M_k-M_{k-1})/k\) 是 \(L^2\) 有界鞅:不同增量正交,二阶矩之和由常数倍 \(\sum_k k^{-2}\) 控制。由鞅收敛定理该级数几乎必然收敛,再由Kronecker 引理得 \(M_n/n\to0\)。有界端点差除以 \(n\) 也趋零,遂得所需平均收敛。
为说明误差,先假设链已平稳,记 \(H_k=h(X_k)\)、\(\gamma_j=\operatorname{Cov}_\pi(H_0,H_j)\)。直接展开有限和的方差,得到精确等式
如果 \(\sum_{j\ge1}|\gamma_j|<\infty\),则括号趋于 \(\sigma_{\rm as}^2=\gamma_0+2\sum_{j\ge1}\gamma_j\ge0\)。协方差可求和给出这个方差极限,本身并不在任意平稳过程中自动证明 CLT。对有限不可约非周期链,相关性几何衰减,且马尔可夫链 CLT 成立;这里引用该CLT,相关条件与误差估计见文献 [3];若渐近方差为零则对应退化极限。
当 \(\gamma_0>0\) 且 \(\sigma_{\rm as}^2>0\) 时,定义积分自相关时间 \(\tau=\sigma_{\rm as}^2/\gamma_0=1+2\sum_{j\ge1}\rho_j\),以及对这个观测函数的有效样本量 \(n_{\rm eff}\approx n/\tau\)。它来自渐近方差匹配,数值随观测函数而变化;原始输出仍是相关样本。
例 1 对 \(h(x)=\mathbf1_{\{x=1\}}\),有 \(\gamma_0=3/16\)。条件均值递推给出 \(E[h(X_k)-1/4\mid X_0]=(-1/3)^k(h(X_0)-1/4)\),所以 \(\rho_k=(-1/3)^k\)、\(\tau=1+2\sum_{k\ge1}(-1/3)^k=1/2\),\(\sigma_{\rm as}^2=3/32\)。负相关使此函数的有效样本量约为 \(2n\),这不违背样本依赖的事实。
14.4.6 预热、诊断与习题
丢弃最初若干步称为预热(burn-in)。丢弃一个固定长度的前缀不改变上述时间平均极限,但没有适用于所有目标分布的通用预热步数。有限长度的轨迹可能长期困在一个区域,表面稳定也无法证明已经探索全部支撑。多起点比较、轨迹图、自相关和分批均值可发现问题;不能仅凭其中一项“看起来正常”就断言收敛。
估计误差时,可以在适用的链 CLT 和一致方差估计条件下采用分批均值或谱方差方法。直接删去大部分样本并不自动提高单位计算量的效率;相关性较低和总信息更多是两个不同的比较目标。关于一般状态空间条件与实践误差估计,见文献 [3]。
习题 14.4-1 · 非对称提议不能遗漏修正比
目标为两个状态上的均匀分布,提议满足 \(q(0,1)=0.9\)、\(q(1,0)=0.3\)。求接受率及真实转移;若错误地总是接受,会得到什么平稳分布?
解答。接受率依次为 \(1/3\) 与 \(1\),所以实际转移两方向都为 \(0.3\),自环都为 \(0.7\),均匀分布不变。若总是接受,则转移就是提议核,平衡方程为 \(0.9\pi_0=0.3\pi_1\),故平稳分布是 \((1/4,3/4)\),与目标不同。
习题 14.4-2 · 顺序 Gibbs 可以不满足详细平衡
在例 2 中令 \(K\) 表示先更新 \(x\) 再更新 \(y\) 的一整轮,计算 \(\pi(0,0)K((0,0),(1,0))\) 及反方向的平稳流。
解答。正方向转移为 \(0.2\times0.2=0.04\),流量为 \(0.4\times0.04=0.016\)。反方向先把 \(x\) 更新为 \(0\) 的概率为 \(0.8\),再把 \(y\) 更新为 \(0\) 的概率也为 \(0.8\),所以流量为 \(0.1\times0.64=0.064\)。两者不等,故整轮核不满足详细平衡。然而每个单坐标更新都保持目标,复合仍然不变;可逆性是比不变性更强的性质。
习题 14.4-3 · 相同目标也可以混合得很慢
在例 1 的每一步,以概率 \(\varepsilon\) 执行原转移,否则不动。求新的滞后相关系数与积分自相关时间,计算 \(\varepsilon=0.03\) 的值。
解答。新核为 \(P_\varepsilon=(1-\varepsilon)I+\varepsilon P\),故仍保持 \(\pi\)。中心化观测原来的特征值为 \(-1/3\),现在变成 \(r=1-4\varepsilon/3\)。当 \(0<\varepsilon\le1\) 时 \(|r|<1\),于是 \(\rho_k=r^k\)、\(\tau=(1+r)/(1-r)\)。取 \(0.03\),得 \(r=0.96\)、\(\tau=49\),有效样本量约 \(n/49\)。目标分布没有改变,但大量停留显著降低这个观测均值的效率。
习题 14.4-4 · 周期性与时间平均
链在 \(0,1\) 间确定性交替。说明从状态 \(0\) 出发的分布是否收敛,以及状态 \(1\) 的时间比例是否收敛。能否使用本节绝对可求和的相关公式?
解答。转移矩阵为 \(\begin{pmatrix}0&1\\1&0\end{pmatrix}\),不可约但周期为 \(2\)。分布在两点质量之间交替,不收敛;前 \(n\) 步中状态 \(1\) 的计数与 \(n/2\) 相差至多 \(1\),故时间比例趋于 \(1/2\)。若改从平稳均匀分布出发,中心化指标的相关系数为 \((-1)^k\),绝对值之和发散,所以不能套用要求绝对可求和的无限级数公式。有限和的方差等式仍成立。
关键术语
- MCMC:用具有目标平稳分布的链估计期望。
- Metropolis–Hastings:通过接受拒绝修正提议的方向概率流。
- 详细平衡:每对状态之间的平稳流相等。
- Gibbs 更新:从完整条件分布重新生成一个坐标。
- 泊松方程:把中心化观测表示为 \(g-Pg\)。
- 渐近方差:长时间平均的方差乘以样本数后的极限。
- 积分自相关时间:把全部相关滞后累积后的方差比例。
- 有效样本量:针对指定观测函数按方差匹配定义的 iid 等价数量。
- 预热:丢弃初始前缀以缓和初始分布影响,需单独评估。
参考文献
- Aldous, D. & Fill, J. A. Reversible Markov Chains and Random Walks on Graphs, §11.2. 作者在线章节:Metropolis–Hastings 构造、连通条件与条件更新。
- Durrett, R. (2019). Probability: Theory and Examples, 5th ed., §5.5–5.8。作者版教材:平稳分布、渐近行为及一般状态空间的收敛条件。
- Geyer, C. J. Stat 5421 Lecture Notes: Bayesian Inference via Markov Chain Monte Carlo. 明尼苏达大学官方讲义:MCMC、相关输出的误差估计和诊断边界。