1. 从算不动说起为什么需要MCMC做概率图模型的人迟早会撞上一堵墙贝叶斯网络的后验分布往往没有解析解。你搭好了一个结构清晰的网络条件概率表也填好了观测数据也进来了结果一到要算后验概率或者做参数估计的时候发现积分根本积不出来。这不是你数学不好而是高维空间里的归一化常数本身就是个噩梦。我最早接触这个问题是在做一个故障诊断的小项目网络里也就十来个节点但变量之间的耦合关系一多联合分布的维度就上去了。当时天真地想着用网格法近似每个维度切100个点10个维度就是100的10次方这个数字大到连想都不用想。后来才明白这就是所谓的维数灾难也是MCMC方法存在的根本理由。MCMC全称是马尔可夫链蒙特卡洛Markov Chain Monte Carlo。拆开看就是两个东西的组合用马尔可夫链去构造一个能收敛到目标分布的随机过程再用蒙特卡洛方法从这个链上采样做估计。它的核心价值在于你不需要知道目标分布的归一化常数只需要知道它正比于某个容易计算的函数就行。这一点对贝叶斯推断来说简直是救命稻草因为后验分布正好就是似然乘以先验再除以一个我们算不出来的归一化项。这篇文章我打算把MCMC的抽样原理从头到尾捋一遍重点讲清楚三件事为什么直接采样不行、马尔可夫链凭什么能收敛到目标分布、以及最常用的两种算法Metropolis-Hastings和Gibbs采样到底是怎么工作的。适合已经了解贝叶斯网络基本概念、但一看到MCMC就头大的读者。我会尽量用生活化的类比和具体的计算例子把那些看起来吓人的公式拆开揉碎。2. 直接采样为什么行不通归一化常数的诅咒2.1 贝叶斯推断里的那个除不掉的项贝叶斯定理的形式大家都熟后验正比于似然乘以先验。写成公式就是P(θ|D) P(D|θ)P(θ) / P(D)。问题就出在分母那个P(D)上它等于对所有的θ积分P(D|θ)P(θ)dθ。在低维情况下这个积分还能算比如θ是一维或者二维的时候数值积分完全够用。但贝叶斯网络里θ通常是一组参数维度可能几十上百这个积分就变成了高维空间上的多重积分。我见过很多初学者在这里卡住觉得不就是个积分吗数值方法那么多。但高维积分和低维积分完全是两码事。一维的时候你用梯形法切1000个点就能得到很高精度十维的时候每个维度切10个点就是10的10次方个点计算量直接爆炸。更麻烦的是高维空间里大部分体积都集中在角落附近均匀网格根本抓不住概率质量集中的区域。MCMC的巧妙之处就在于它绕开了这个积分。它不需要你算出P(D)只需要你能计算未归一化的后验P(D|θ)P(θ)。因为马尔可夫链的构造只依赖于目标分布之间的比值归一化常数在比值里自然消掉了。这就像是你不需要知道一座山的具体海拔高度只需要知道哪里比哪里高就能顺着坡度走到山顶。2.2 蒙特卡洛估计的基本逻辑蒙特卡洛方法的核心思想其实很朴素如果你能从某个分布里采样那就可以用样本均值去估计这个分布下的期望。比如你想算E[f(θ)]其中θ服从后验分布P(θ|D)那你就从后验里抽N个样本θ1到θN然后用(1/N)Σf(θi)去估计。根据大数定律N足够大的时候这个估计会收敛到真实值。这个逻辑本身没问题问题在于如果你能从某个分布里采样这个前提。对于简单的分布比如正态分布、均匀分布我们有现成的采样方法。但对于贝叶斯网络的后验分布你根本不知道它长什么样更别说直接采样了。MCMC解决的就是这个怎么从复杂分布里采样的问题。它构造一条马尔可夫链让这条链的平稳分布恰好就是你要的目标分布然后你沿着链走走够了之后收集到的样本就近似服从目标分布。这里有个关键点需要强调MCMC产生的样本不是独立的。相邻的样本之间有相关性因为马尔可夫链的下一步依赖于当前状态。这跟独立同分布的理想情况有差距但在实际应用中只要链混合得足够好样本之间的相关性衰减得够快估计结果依然可用。只不过有效样本量会比实际样本数小这个后面会详细说。2.3 高维空间带来的直觉失效在低维空间里培养的直觉到了高维空间往往会骗你。举个经典的例子一个高维球体大部分体积都集中在表面附近的一层薄壳里。如果你在球内均匀采样绝大多数点都会靠近表面。类似地高维概率分布的质量往往集中在一些低维的流形附近而不是均匀散布在整个空间。这意味着基于均匀网格或者简单拒绝采样的方法在高维下效率极低。拒绝采样需要你找一个覆盖目标分布的提议分布高维下这个提议分布很难找而且接受率会随着维度增加指数下降。MCMC通过让样本在状态空间里游走自动地探索概率质量高的区域不需要你事先知道这些区域在哪里。这是它相比其他采样方法的根本优势。3. 马尔可夫链的平稳分布MCMC的理论基石3.1 马尔可夫链的直观理解马尔可夫链就是一个状态序列下一个状态只依赖于当前状态跟更早的历史无关。这个性质叫无记忆性。你可以把它想象成一个喝醉的人在走路他下一步往哪走只取决于他现在站在哪不取决于他是怎么走到这里的。在MCMC里这个喝醉的人就是在参数空间里游走。我们设计一套转移规则使得这个游走过程最终会稳定下来稳定之后他在各个区域停留的时间比例正好等于目标分布的概率。换句话说链的平稳分布就是我们的目标后验分布。这里需要区分几个概念。转移核T(x|x)表示从状态x转移到x的概率密度。如果存在一个分布π使得π(x) ∫π(x)T(x|x)dx那π就是这个链的平稳分布。MCMC的目标就是构造一个转移核让它的平稳分布恰好是我们想要的后验分布。3.2 细致平衡条件构造平稳分布的捷径直接验证一个分布是不是平稳分布需要做积分这往往很难。但有一个充分条件叫细致平衡条件验证起来容易得多。细致平衡说的是π(x)T(x|x) π(x)T(x|x)也就是从x转移到x的概率流等于从x转移到x的概率流。如果细致平衡成立对两边同时关于x积分左边得到∫π(x)T(x|x)dx右边得到π(x)∫T(x|x)dx π(x)。所以细致平衡确实能推出平稳性。Metropolis-Hastings算法的设计就是围绕如何构造满足细致平衡的转移核来展开的。我个人的理解是细致平衡是一种可逆性的体现。链在平稳分布下任意两个状态之间的双向流量相等整个系统处于动态平衡。这个条件比平稳性更强但更容易构造和验证所以成了MCMC算法设计的主要工具。3.3 收敛性链什么时候才算热透理论上只要链是不可约的从任何状态都能到达任何正概率状态和非周期的不会陷入固定周期的循环它就会收敛到平稳分布。但理论保证的是足够长时间后会收敛实际中你永远不知道这个足够长到底是多长。这就引出了MCMC实践中最头疼的问题收敛诊断。你跑了一条链怎么知道它已经收敛到目标分布了常用的方法包括看迹图trace plot是否稳定、计算Gelman-Rubin统计量多条链之间的方差与链内方差的比值、检查自相关函数衰减速度等。但这些方法都只是启发式的没有哪个能给出绝对的保证。我的经验是宁可多跑一段burn-in预热期把前面的样本丢掉。因为链从初始状态出发需要一段时间才能忘记初始值的影响走到目标分布的高概率区域。burn-in的长度没有固定标准取决于初始值和目标分布的距离以及链的混合速度。一般至少丢掉前10%到20%的样本复杂问题可能要丢掉一半。4. Metropolis-Hastings算法拒绝采样与马尔可夫链的结合4.1 提议分布与接受率的博弈Metropolis-Hastings简称MH算法的思路很巧妙既然不能直接从目标分布π里采样那就从一个容易采样的提议分布q里生成候选样本然后以一定的概率决定接受还是拒绝这个候选。接受的概率设计得恰好使得最终的转移核满足细致平衡。具体来说假设当前状态是x从q(x|x)里抽一个候选x。然后计算接受率α min(1, [π(x)q(x|x)] / [π(x)q(x|x)])。以概率α接受x否则留在x。注意这里π的归一化常数在比值里消掉了所以只需要计算未归一化的π就行。这个接受率的设计不是随便来的。它恰好使得π(x)q(x|x)α(x,x) π(x)q(x|x)α(x,x)满足细致平衡。直观理解就是如果候选状态x的概率密度比当前状态x高那就一定接受如果比x低那就以概率比值接受避免链总是往低概率区域跑。4.2 对称提议分布下的简化Metropolis算法如果提议分布是对称的即q(x|x) q(x|x)那接受率就简化为α min(1, π(x)/π(x))。这就是原始的Metropolis算法比MH更简单。对称提议分布最常见的例子是以当前状态为中心的正态分布也就是x x ε其中ε服从均值为零的正态分布。这种随机游走式的提议分布实现起来很简单但有一个明显的缺点步长不好选。步长太小链走得慢自相关高需要很多样本才能有效探索整个空间步长太大接受率低链经常原地不动同样效率低下。实践中一般调整步长使得接受率在20%到50%之间这是一个经验性的平衡点。我刚开始用MH的时候经常犯的一个错误是步长设得太小结果迹图看起来平滑得很以为收敛得很好实际上链根本没怎么动有效样本量少得可怜。后来学会了看自相关函数才发现相邻样本的相关性高得离谱。所以光看迹图不够一定要结合自相关和有效样本量来判断。4.3 一个手算例子从混合高斯后验中采样为了把MH的流程说清楚我构造一个简单的例子。假设目标分布是一个一维的混合高斯0.3×N(0,1) 0.7×N(5,1)。这个分布的归一化常数是1但我们可以假装不知道只用未归一化的形式π(x) 0.3×exp(-x²/2) 0.7×exp(-(x-5)²/2)。提议分布用对称的正态分布q(x|x) N(x, 2²)步长2。从x0开始第一步从N(0,4)里抽一个候选假设抽到x1.5。计算接受率π(1.5) 0.3×exp(-1.125) 0.7×exp(-(1.5-5)²/2) 0.3×0.325 0.7×exp(-6.125) ≈ 0.0975 0.0015 0.099。π(0) 0.3×1 0.7×exp(-12.5) ≈ 0.3。接受率α min(1, 0.099/0.3) 0.33。然后生成一个0到1之间的均匀随机数如果小于0.33就接受否则拒绝。这个例子里从x0跳到x1.5是往低概率区域走所以接受率小于1。如果候选是x4.8π(4.8) ≈ 0.3×exp(-11.52) 0.7×exp(-0.02) ≈ 0.686接受率就是min(1, 0.686/0.3) 1一定接受。这样链就会倾向于在第二个高斯峰附近停留最终两个峰之间的样本比例会接近0.3:0.7。实际跑的时候你会发现链偶尔会从一个峰跳到另一个峰但跳跃的频率取决于两个峰之间的低谷有多深。如果两个峰离得远、方差小链可能在一个峰里待很久才跳过去一次。这就是所谓的混合不良是多峰分布采样的经典难题。5. Gibbs采样条件分布已知时的优雅方案5.1 逐坐标更新的核心思想Gibbs采样是MH的一个特例但它的形式更简洁不需要设计提议分布也不需要计算接受率。它的做法是每次只更新一个变量或一组变量从该变量在给定其他变量当前值下的条件分布中采样。因为条件分布是精确采样的所以接受率恒为1每次采样都被接受。具体来说假设参数向量θ (θ1, θ2, ..., θd)Gibbs采样的一轮迭代是从P(θ1|θ2,...,θd, D)中抽θ1从P(θ2|θ1,θ3,...,θd, D)中抽θ2依此类推直到更新完所有变量。然后重复这个过程。这个方法的优雅之处在于它把高维联合分布的采样问题分解成了一维条件分布的采样问题。而条件分布往往比联合分布简单得多很多时候有解析形式可以直接采样。在贝叶斯网络里一个节点的条件分布只依赖于它的父节点、子节点和子节点的其他父节点这个局部性使得Gibbs采样特别自然。5.2 贝叶斯网络中的Gibbs采样实操在贝叶斯网络里做Gibbs采样关键是要能写出每个节点的全条件分布。对于一个节点Xi它的全条件分布正比于它自己的条件概率P(Xi|Pa(Xi))乘以所有以Xi为父节点的子节点的条件概率P(Xj|Pa(Xj))。这个乘积通常只涉及少数几个节点计算量很小。我拿一个简单的例子说明。假设有三个二值变量A、B、C结构是A→B→C。联合分布是P(A)P(B|A)P(C|B)。Gibbs采样的一步是固定B和C的当前值从P(A|B,C)中抽A。根据贝叶斯定理P(A|B,C) ∝ P(A)P(B|A)P(C|B)。注意P(C|B)跟A无关所以在给定B的情况下它是个常数可以忽略。所以P(A|B,C) ∝ P(A)P(B|A)这就是A的全条件分布。实际实现的时候你只需要遍历每个节点计算它的全条件分布通常就是几个概率值的乘积归一化后采样。这个过程不需要计算任何边缘概率或者归一化常数非常高效。但有一个前提全条件分布必须能采样。对于离散变量直接按概率分布抽就行对于连续变量如果条件分布是标准形式比如正态、Gamma也可以直接抽如果条件分布形式复杂可能还需要嵌套一个MH步骤这就变成了Metropolis-within-Gibbs。5.3 Gibbs采样的软肋强相关变量Gibbs采样有一个致命的弱点当变量之间高度相关时逐坐标更新会非常慢。因为每次只动一个坐标如果两个变量几乎同步变化链就会沿着一条狭窄的对角线缓慢移动需要很多步才能探索整个后验分布。我遇到过一个典型的例子两个参数θ1和θ2的后验分布是一个狭长的椭圆长轴沿着θ1θ2的方向。Gibbs采样每次只更新θ1或θ2链就会在椭圆上之字形前进效率极低。这种情况下MH用各向异性的提议分布反而可能更好或者需要对变量做重参数化让它们在新的坐标系下近似独立。解决这个问题的常用方法包括对变量做分组更新把强相关的变量放在一组一起采样、使用辅助变量方法、或者直接换用Hamiltonian Monte CarloHMC这类利用梯度信息的算法。HMC不在本文的讨论范围内但它是现代贝叶斯计算的主流工具值得后续单独展开。6. 收敛诊断与样本处理那些文档里不写的经验6.1 burn-in到底该丢多少Burn-in预热期是指链从初始状态走到目标分布高概率区域所需的时间。这段时间的样本不服从目标分布必须丢掉。但丢多少没有理论上的确定答案取决于初始值和目标分布的距离以及链的混合速度。我的做法是先跑一条短链比如几千步看迹图什么时候稳定下来把稳定点之前的样本丢掉。然后正式跑的时候burn-in长度设为这个稳定点的两倍保险一些。如果初始值离高概率区域很远burn-in可能需要几万甚至几十万步。一个实用的技巧是初始化在多个不同的起点如果不同起点的链最终都收敛到同一个区域说明burn-in够了。还有一个常见的误区是认为burn-in之后样本就独立了。不是的burn-in只是去掉了初始值的 transient 影响样本之间的自相关依然存在。所以burn-in和 thinning稀释是两回事。Thinning是每隔k步取一个样本降低自相关但代价是丢弃了大量样本。现代观点认为如果计算资源允许与其thinning不如多跑几步保留所有样本因为thinning虽然降低了自相关但总信息量可能反而减少。6.2 有效样本量比样本总数更重要的指标有效样本量ESS衡量的是你的MCMC样本相当于多少个独立样本。计算公式是ESS N / (1 2Σρk)其中ρk是滞后k阶的自相关。如果样本完全独立ESS等于N如果高度相关ESS可能只有N的百分之几。我见过很多人跑了一百万步的链以为样本量很大结果一算ESS只有几千估计的蒙特卡洛误差大得吓人。所以看ESS比看总步数靠谱得多。一般来说对于后验均值这样的简单统计量ESS达到几千就够用了对于尾部分位数或者极端事件概率可能需要ESS上万。提高ESS的方法包括优化提议分布的步长、对参数做重参数化、使用更高效的采样算法如HMC、或者并行跑多条链。其中重参数化往往是最容易被忽视但效果最显著的手段。比如对于层级模型把中心化参数化改成非中心化参数化ESS可能提高一个数量级。6.3 多链诊断Gelman-Rubin统计量单条链的迹图看起来稳定不代表收敛了。可能链只是在一个局部区域里慢慢游走还没探索到其他区域。多链诊断是更可靠的方法从不同的初始值跑多条链比较链间方差和链内方差。Gelman-Rubin统计量也叫R-hat的定义是R-hat sqrt( (N-1)/N × W B/N ) / W其中W是链内方差B是链间方差N是每条链的长度。如果所有链都收敛到同一个分布R-hat应该接近1。一般要求R-hat小于1.1严格一点要求小于1.05。我自己的习惯是至少跑4条链初始值分散在参数空间的不同区域。如果R-hat大于1.1说明链之间还没达成一致需要继续跑或者检查模型是否有问题。R-hat的一个局限是它只检测均值层面的收敛对于方差或者尾部的收敛不敏感。所以最好结合ESS和迹图一起看。7. 从理论到代码一个完整的MH采样实现7.1 目标分布与提议分布的定义我用Python写一个完整的MH采样例子目标分布还是那个混合高斯π(x) ∝ 0.3×exp(-x²/2) 0.7×exp(-(x-5)²/2)。提议分布用对称正态N(x, 2²)。代码里我会把每一步的计算都写清楚方便对照前面的理论。import numpy as np import matplotlib.pyplot as plt def target_unnormalized(x): return 0.3 * np.exp(-x**2 / 2) 0.7 * np.exp(-(x-5)**2 / 2) def metropolis_hastings(n_samples, step_size, x_init): samples np.zeros(n_samples) x_current x_init accept_count 0 for i in range(n_samples): x_proposal np.random.normal(x_current, step_size) ratio target_unnormalized(x_proposal) / target_unnormalized(x_current) accept_prob min(1.0, ratio) if np.random.uniform() accept_prob: x_current x_proposal accept_count 1 samples[i] x_current accept_rate accept_count / n_samples return samples, accept_rate这段代码里target_unnormalized就是未归一化的目标分布注意没有除以任何常数。metropolis_hastings函数里每次从当前状态出发抽一个候选计算接受概率然后决定是否接受。接受率会被记录下来用来判断步长是否合适。7.2 步长选择与接受率的调参过程步长是MH算法里最重要的调参。步长太小接受率高但链移动慢步长太大链移动快但接受率低。我一般先跑一个短链试几个步长看哪个的接受率在30%到40%之间。for step in [0.5, 1.0, 2.0, 4.0, 8.0]: _, acc metropolis_hastings(5000, step, 0.0) print(fstep_size{step}, accept_rate{acc:.3f})跑出来大概是步长0.5时接受率约0.95步长1.0时约0.85步长2.0时约0.55步长4.0时约0.25步长8.0时约0.08。根据经验步长2.0到4.0之间比较合适。我选步长3.0接受率大概在0.35左右。这里有个细节接受率不是越高越好。接受率太高说明步长太小链在原地打转接受率太低说明步长太大链经常被拒绝。最优接受率取决于目标分布的形状和维度一维情况下理论最优大约是0.44高维情况下会降到0.23左右。但这些都是渐近结果实际中不用太纠结30%到50%都是可接受的。7.3 结果验证迹图、直方图与真实分布对比跑完采样后我习惯做三个检查迹图看链是否稳定、直方图看样本分布是否接近目标、自相关函数看样本相关性。samples, acc metropolis_hastings(50000, 3.0, 0.0) burn_in 5000 samples_kept samples[burn_in:] fig, axes plt.subplots(1, 3, figsize(15, 4)) axes[0].plot(samples) axes[0].set_title(Trace Plot) axes[1].hist(samples_kept, bins50, densityTrue, alpha0.6) x_grid np.linspace(-5, 12, 500) axes[1].plot(x_grid, target_unnormalized(x_grid) / 3.0, r-) axes[1].set_title(Histogram vs Target) axes[2].acorr(samples_kept, maxlags100) axes[2].set_title(Autocorrelation) plt.show()迹图应该看起来像一条稳定的噪声没有明显的趋势。直方图应该跟红色的目标曲线吻合。自相关函数应该快速衰减到零附近如果衰减很慢说明混合不好。这个例子里因为两个峰离得比较远链偶尔会在两个峰之间跳跃自相关可能在滞后几十步的地方还有小幅波动这是正常的。8. 写在最后一些踩过的坑和实用建议MCMC的理论看起来优雅但实际用起来坑不少。我最早做贝叶斯推断的时候最大的错误是只看迹图就下结论。迹图稳定不代表收敛可能只是链在一个局部区域里慢慢晃。后来养成了多链诊断的习惯才避免了几个严重的错误。另一个常见的坑是忽略参数的相关性。Gibbs采样在强相关变量上效率极低MH在高维空间里也很难设计好的提议分布。这时候要么做重参数化要么换用HMC或者变分推断。没有哪种方法是万能的关键是根据问题的结构选择合适的工具。还有一个经验是MCMC的计算成本往往被低估。跑一条链可能需要几分钟到几小时调参和诊断又要花好几倍的时间。所以在实际项目里如果后验分布可以用变分推断或者拉普拉斯近似快速得到先用这些方法做快速迭代等模型结构稳定了再用MCMC做精细推断。这样能省下大量时间。最后说一个诊断的小技巧如果你不确定burn-in够不够可以跑两条链一条从合理的初始值出发一条从极端的初始值出发比如参数空间的角落。如果两条链的迹图在某个点之后重合了说明burn-in够了。这个方法虽然笨但很可靠。