MCMC采样完全指南：从Metropolis-Hastings到Gibbs实战

从贝叶斯推断的一个具体需求说起。我第一次正儿八经把蒙特卡洛马尔可夫链MCMC捡起来学明白是因为手头有个模型的后验分布根本没有解析解。数据量不大算一个参数的期望却处处碰壁。翻资料的时候满屏的Metropolis-Hastings、Gibbs采样、R-hat、有效样本量每个词单独看都认识连起来就不知道它们在唱哪出戏。这篇小结不是教科书式的定义堆砌而是按我实际学习时“从问题到解法、从原理到踩坑”的顺序来写希望给同样被MCMC绕晕的人一条比较顺的上手路径。1. 从哪开始MCMC到底在解决什么问题学MCMC之前得先搞清楚它到底被造出来干嘛的。很多教程一上来就讲马尔可夫链的数学性质结果初学者连“为什么要采样”都没想明白后面全是空中楼阁。1.1 一个算不动的积分为何拦住所有人在贝叶斯推断里核心任务通常是计算后验期望这种东西E[f(θ) | x] ∫ f(θ) p(θ|x) dθ如果后验分布是正态分布、Beta分布这类有现成公式的情况积分可以手算或查表不需要任何采样方法。但现实模型稍微复杂一点比如参数之间互相嵌套、似然函数长得不友好后验分布的归一化常数就积不出来了。没有归一化常数期望算不了区间估计也算不了整个推断就卡死。更麻烦的是当参数维度到了几十甚至上百即便你知道密度函数的表达式想在整个高维空间做数值积分也几乎不可能。多维求积的格子点数量随维度指数增长这一下就回到了“超级难算”的老路上。MCMC干的事情就一句话不去硬算那个积分而是从目标分布里抽取大量样本用样本均值去逼近期望。只要样本确实是按照目标分布采出来的大数定律就帮你把积分问题变成了统计问题。1.2 直接采样、拒绝采样与重要性采样的天花板既然要采样为什么不直接采答案很现实目标分布的累积分布函数往往没有解析反函数你连怎么把均匀分布随机数变换成目标分布都不知道。拒绝采样大概是大多数人第一次接触的非直接采样方法。逻辑很直观找一个好采样的提议分布把目标分布包住然后在提议分布里随机取点、随机决定要不要接受。理论很漂亮但实际操作里它要求一个足够好的“包络”只要目标分布在某个区域的密度比值很高拒绝率就会暴涨。一维情况下调一调还好到了十维、二十维可行区域在整个状态空间里占的比例小到我们很难提前设计一个覆盖它的提议分布拒绝采样直接就废了。重要性采样则是换了个思路不直接生成目标分布的样本而是给提议分布的样本配上权重用加权求和来逼近期望。问题在于当目标分布和提议分布差别较大的时候权重方差会非常大少数几个样本权重接近1剩下全是0估算结果震荡剧烈完全不可信。这些方法的共同弱点是把“独立性”当成金科玉律总想着一次生成一批互相独立的样本。方向本身没有错但复杂度一高独立采样的成本就让人扛不住了。1.3 MCMC的破局思路不采独立样本采一条链条MCMC愿意牺牲样本之间的独立性换来一个极其重要的好处我们只需要能够计算目标分布的相对密度连归一化常数都不需要就可以设计一个马尔可夫链让它的稳态分布正好等于目标分布。然后沿着这条链一直走走到哪、采到哪最后把链条中段以后的那些样本当作目标分布的近似样本。我刚开始不太理解“稳态分布等于目标分布”这个提法总怀疑这种依赖前一个样本的采样方式会不会有偏向。后来想通了依赖性是现实最终样本集合的分布是理论保证的。MCMC相当于用“走路”代替“投飞镖”只要走路规则设计得当走过足够长时间以后你在各个区域停留的频率就会逼近目标分布。搞清楚这个定位以后再回头看Metropolis-Hastings和Gibbs就没有那么神秘了。2. 马尔可夫链的“稳”为什么链能代表目标分布MCMC里面的马尔可夫链不是随便一条链都行它必须收敛到一个我们想要的稳态分布。这一节把“马尔可夫链凭什么能代表目标分布”这件事讲透。2.1 从简单随机游走出发理解转移核马尔可夫链的核心是一个转移核或者说条件分布T(θ | θ)表示从当前状态θ转移到新状态θ的概率。最简单的情形是随机游走θ θ εε是某个随机扰动。这里有个容易被忽视的点马尔可夫链的下一跳只依赖当前状态不依赖更早的历史。这个性质看起来是很强的限制实际却是MCMC能成立的基石。因为一旦转移规则固定链条的“记忆”就只有一阶经过足够多步以后初始状态的影响会逐渐衰减链条的分布会趋于一个由转移核唯一决定的极限分布。举个例子你从家里的任何一个房间开始可乐朗姆地走只要路线上不出现“只出不进”的陷阱位置走一万步之后在哪个房间的概率分布会稳定下来。这个稳定分布不取决于起点只取决于你“怎么走”的规则。2.2 详细平衡条件一个充分不必要条件MCMC设计转移核最常用的套路是满足详细平衡条件π(θ) T(θ | θ) π(θ) T(θ | θ)这个式子读起来就是从θ流向θ的概率质量等于从θ流回θ的概率质量。如果这个条件对所有的状态对都成立两者一抵消π自然就是这条链的稳态分布。顺带说一句这是充分条件不是必要条件。意思是满足详细平衡就一定收敛到π但有些链不满足详细平衡照样有稳态分布π。只不过工程界99%的MCMC算法都基于详细平衡来设计因为它实在是太方便了。你只要证明或者构造T的时候按这个等式调稳态分布就不用担心。我初学的时候一直纠结一个问题为什么不能直接让T(θ|θ) π(θ)一步就跳到目标分布理论上当然可以但你得能从π里采样才能构造这个转移核——这就是原问题本身。MCMC的巧妙恰恰在于它允许你从一个更容易采样的分布里提出“候选”再用一个接受/拒绝的步骤把详细平衡条件修补起来。这个修补动作就是下一章Metropolis-Hastings的核心。2.3 稳态分布与遍历性初学时最容易混淆的两个概念“稳态分布”说的是链走到无限步以后θ的边际分布稳定在π(θ)不再随时间变化。而“遍历性”说的是时间平均等于空间平均。对MCMC来说我们真正依赖的是后者链上依次访问的样本点长期来看足够“均匀地”跑遍目标分布的支撑区域使得样本均值能逼近期望。这两个概念的区别在应用中的体现是一条链可能稳态分布已经算对了但遍历性很差典型表现就是链在状态空间里磨磨蹭蹭自相关性极高跑十万步实际有效的信息量可能只相当于几百个独立样本。这时候光看链是否“收敛”是不够的还得看它“走得多利索”。所以后面专门有一章讲诊断工具本质都是在衡量遍历性。3. Metropolis-Hastings算法接受-拒绝藏在里面的采样器Metropolis-Hastings简称MH是最基础、最能直观理解MCMC工作原理的算法。它把上一节说的“修补详细平衡条件”落实到具体步骤上。3.1 从提议分布到接受概率的完整推导设想我们已经有一个容易采样的提议分布q(θ | θ)比如以当前θ为中心的正态分布。我们按q提出一个候选位置θ然后决定接受还是拒绝。接受概率a(θ, θ)的设计目标是让最终的转移核满足详细平衡。转移核可以写成T(θ | θ) q(θ | θ) × a(θ, θ)代进详细平衡条件π(θ) q(θ | θ) a(θ, θ) π(θ) q(θ | θ) a(θ, θ)要同时照顾两个方向的等式一个经典取法是令a(θ, θ) min(1, [π(θ) q(θ | θ)] / [π(θ) q(θ | θ)])验证起来非常干净。如果分子比分母大接受率是1从θ到θ的流直接放行此时反方向从θ到θ的流自动压缩成一个比例两边重新平衡。这就是MH里“接受概率”公式的来源不是凭空造出来的而是详细平衡条件下解出来的。一个重要推论公式里π只出现在分子分母的比例中归一化常数因子会上下约掉。这就是MCMC比直接积分高明的地方——哪怕后验分布的归一化常数无法计算只要你能算出未归一化的密度函数值就能跑算法。绝大多数贝叶斯模型恰恰满足这个条件。3.2 一个一维t分布的完整示例含代码拿一个具体例子走一遍流程。假设目标分布是自由度4的学生t分布密度函数大家都会背但累积分布函数没有解析反函数想直接采样并不直观。我们用MH来采。提议分布选正态随机游走θ θ εε ~ N(0, σ²)。因为正态分布对称q(θ|θ) q(θ|θ)接受率简化成a min(1, π(θ) / π(θ))这意味着如果候选点密度比当前点高一定接受如果低按密度比值概率接受。下面是完整代码import numpy as np def target_density(x): # Student-t, df4, 未归一化即可 return (1 x**2 / 4) ** (-(4 1) / 2) def mh_sample(n_samples, init0.0, sigma2.0, warmup1000): samples [] x init for i in range(n_samples warmup): proposal x sigma * np.random.randn() ratio target_density(proposal) / target_density(x) if np.random.rand() ratio: x proposal if i warmup: samples.append(x) return np.array(samples) samples mh_sample(5000) print(均值约为:, samples.mean())跑完以后样本均值应该接近0样本分位数和t分布理论分位数的偏差在采样误差范围内。我第一次跑通这个例子时印象很深因为代码总共不到二十行却解决了教科书上需要查表才能完成的采样问题。3.3 步长与接受率的博弈sigma的取值是整个例子最需要经验的地方。sigma太小候选点总是赖在当前位置附近接受率高但链跑得慢自相关强有效样本量低。sigma太大候选点动不动跳到密度很低的地方接受率暴跌链老是原地踏步。我一开始习惯把sigma调到接受率看上去挺健康的水平比如70%左右后来才发现这个目标要分维度看。对随机游走Metropolis高维情形下的渐近最优接受率大约在0.234一维简单问题则往往可以接受更高一些的接受率。实操时你可以先跑一小段用如下代码估算接受率def acceptance_rate(samples, target_density): # 这里简化示意实际需要记录是否接受 pass经验是一维问题优先把接受率调到40%到70%之间高维问题往20%到35%靠。没有绝对标准但范围可以帮你快速判断步长是否离谱。4. 高维问题随机游走效率崩塌与Gibbs采样MH看起来已经挺通用了但一旦参数维度上来随机游走Metropolis就力不从心。这一节聊聊维数灾难是如何毁掉随机游走的以及Gibbs采样作为补救方案的核心逻辑。4.1 维数灾难如何让Metropolis失效在高维空间里随机游走Metropolis有个著名的效率结局如果你想让链走得足够远就必须把步长拉大但步长一大候选点就会经常落进高密度区域之外接受率暴跌。如果你缩小步长保住接受率链又只能在局部小碎步移动自相关高得吓人。数学上有一个粗糙的标度律要达到同样的移动效率随机游走Metropolis需要的步数会随维度线性甚至更差地增长。直观理解是高维空间中“方向”的数量太多随机选一个方向往前走很大概率是在跟目标分布的等值面横着切而不是往高密度中心走。我在一个30维的简单高斯模型上试过随机游走Metropolis跑了二十万步有效样本量连几千都不到均值估计的方差大得没法用。这个教训让我彻底明白了为什么学术界后来发展出那么多更高级的采样器。4.2 条件分布拆解Gibbs采样的核心逻辑Gibbs采样的想法非常务实既然高维联合分布难采那我就一次只采一个维度其他维度固定住并且利用目标分布的条件分布来采样。它之所以在MCMC里地位特殊是因为当条件分布可以直接采样时根本不涉及接受拒绝的过程。每一轮迭代里给定其他参数当前的取值从每个参数的条件后验分布中抽一个新值然后立刻用它去更新其他参数的条件。从MH的角度看Gibbs其实是接受率恒为1的MCMC特例因为从条件分布采样相当于永远在“往更高概率密度的合理区域移动”详细平衡条件自动满足。这不光让代码简单也避开了步长调节的麻烦。所以但凡模型的条件后验是标准分布比如正态、Gamma、Beta这些Gibbs都是首选。4.3 一个二维例子从条件采样到联合采样用二元正态分布来演示最清晰。标准化以后的二元正态假设相关系数ρ那么条件分布是x1 | x2 ~ N(ρ x2, 1 - ρ²), x2 | x1 ~ N(ρ x1, 1 - ρ²)交替采样就能获得(x1, x2)的联合样本。实现如下import numpy as np def gibbs_2d(n_samples, rho0.8, init(0.0, 0.0), warmup1000): x1, x2 init samples [] for i in range(n_samples warmup): x1 rho * x2 np.sqrt(1 - rho**2) * np.random.randn() x2 rho * x1 np.sqrt(1 - rho**2) * np.random.randn() if i warmup: samples.append((x1, x2)) return np.array(samples) samples gibbs_2d(5000) print(np.cov(samples.T)) # 应接近 [[1, 0.8], [0.8, 1]]我跑这个例子的时候发现rho越大链条走起来越像爬坡每一步的变化都很小。原因是两个变量高度相关时条件分布的方差1-ρ²变得很小每个维度只能在很窄的范围里扰动。这解释了为什么Gibbs面对高度相关的参数时也会出现自相关高的问题——不是算法错了而是模型参数化方式对采样效率影响巨大。5. 链收敛了吗诊断工具与实际使用体验跑到这里你手里已经有一堆样本了。但能不能直接拿去算期望不行前提是链得收敛、样本得足够有效。这一章讲的都是我在实际跑实验时真的会看的诊断指标。5.1 痕迹图、自相关图与有效样本量痕迹图就是把链上每个step的θ画出来。一眼扫过去如果链像一条毛毛虫一样在一个水平带里随机波动说明它比较健康如果长期贴着某个边界往上爬或者往下掉说明要么warm-up不够要么目标分布有多峰结构链还没遍历到另一个峰。自相关图衡量的是间隔k步的两个样本之间的相关性。理想情况下k增大自相关迅速降到0说明每一次跳跃都有信息量。如果k到了50甚至100自相关还在0.5以上就要警惕有效样本量太低。有效样本量ESS是我最看重的单一指标。它告诉你“当前N个样本大约相当于多少个独立样本”。计算公式不复杂但通常依赖统计软件实现。实操经验是ESS至少要有几百个才能保证均值估计的方差在一个可控范围如果你要做分位数估计需要的ESS更高。5.2 R-hat与多链运行R-hat是判断多条链是否收敛到同一分布的经典指标。做法是同时从不同的初始值出发跑多条独立链比较链间方差和链内方差。R-hat接近1说明不同起点没影响链基本收敛。经验阈值通常取1.01或1.05但对高维复杂模型仅靠R-hat还不够最好结合ESS一起看。有意思的是R-hat本身也常常被误用。它诊断的是“几条链是否长得像同分布”而不是“链是否覆盖了整个状态空间”。如果目标分布有一大一小两个峰多链偶尔会集体落在同一个大峰里R-hat照样非常漂亮。这就是为什么多峰分布场景下要额外关注链是否从不同初始值出发都给出了相似的宏观结果。5.3 warm-up和thin到底该不该用warm-up也叫burn-in很好理解初始位置可能远离稳态链条前一段样本地弃之不用通常弃掉前10%到50%。我一般习惯直接跑两段第一段只用来调试步长和看痕迹图确定没问题了再从头正式采样。thin每隔k步保留一个样本的作用是降低样本自相关。以前计算资源紧张时thin几乎是必选项。但现在的共识是thin没有必要盲目用因为它会损失信息量。更好的思路是先把步长和参数化调好让链的自相关本身比较低。除非你有文件大小、后处理耗时方面的限制否则保留全部样本反而更稳妥。6. 实操中踩过的坑与个人心得学MCMC最难的往往不是算法本身而是那些书上不写、但跑起来一定会遇到的坑。最后这部分全是一些具体教训每一条都是真金白银换来的。6.1 接受率不达标的真实代价有一段时间我在一个带约束的参数空间里跑MH发现接受率一直在80%以上还挺高兴结果算了半天均值发现跟预期差很远。后来看了自相关图才知道链条每一步都在原地踏步有效样本量极低十万步的样本还不如别人两千步独立样本可靠。这件事给我的教训是接受率高不等于链好。随机游走Metropolis里接受率过高通常意味着步长太小链在空间里“磨洋工”。正确做法不是单纯看接受率而是结合ESS来评估。如果ESS不高果断增大步长哪怕接受率掉到30%左右只要ESS上来了整体效率反而更高。6.2 高维模型别硬上随机游走三十维高斯模型的教训让我彻底放弃了对高维问题使用朴素随机游走Metropolis的幻想。高维场景优先考虑梯度信息或者模型结构信息。我后来在实际项目中试过利用对数后验梯度构造提议分布的采样器效率比纯随机游走高了一个数量级。如果你能算梯度就不要只用随机游走。如果模型的后验条件分布能算优先试试Gibbs。如果条件分布不是标准分布可以考虑引入辅助变量或者重新参数化让参数之间的相关性降低。一个经验法则链的移动效率很大程度上在模型构建阶段就被决定了一半参数化方式选得好采样器事半功倍。6.3 学习路径建议什么情况下MCMC不是最优解学了MCMC以后很容易走向另一个极端见什么都想采样。实际上MCMC有明确的使用边界。如果后验分布是共轭的直接解析更新公式又快又准如果只是计算某个高维积分但分布比较温顺数值积分或者拉普拉斯近似可能就够了如果模型特别大、样本量特别多变分推断通常比MCMC快好几个数量级虽然牺牲了一些精度但能跑出来的模型就有优势。对于刚入门的人我的路径建议是先搞懂MH拿一维、二维的标准分布跑通再动手实现一遍Gibbs最后才去用现成的概率编程库。因为只有自己手推过接受率公式、亲手调过步长你才能理解那些库里面参数设置的含义。我当年就是跳过推导直接调包结果遇到收敛问题时两眼一抹黑退回来一步步补基础后才真正会用它。回头再看MCMC本身并不算一个“高大上”的算法它更像提供了一套通用的“从不可采样分布中获得样本”的方法论。理解它的关键在于想清楚链为什么能收敛到目标分布以及如何判断链“够不够好”。把这些想通了各种改进版采样器在你眼里都只是同一个框架下的不同工具选择。