简介在机器学习与统计建模中贝叶斯推断提供了一套将先验知识与观测数据融合的完备框架然而真实模型的后验分布往往没有解析解归一化常数也难以计算。蒙特卡洛采样与马尔可夫链的结合——MCMC方法成为解决该问题的核心工具。通过构造一条以目标后验为平稳分布的随机游走链MCMC能够在高维空间中高效采集样本从而近似后验的均值、方差与分位数。从基础的Metropolis-Hastings算法到现代概率编程框架PyMC中内置的NUTS采样器MCMC已广泛用于回归分析、分层模型、深度学习不确定性量化等场景。理解其原理与诊断指标如Rhat、ESS是实践的关键。本文从手写采样器出发系统地讲解MCMC算法家族、Python实现及实战调参经验帮助读者真正掌握贝叶斯推断的计算逻辑。 做贝叶斯推断的人迟早都会撞上MCMC。我最早接触这个概念是在看一篇分层模型的论文时模型本身很简单但后验没有解析形式作者直接用了MCMC采样。当时我天真地想不就是采样嘛直接从分布里抽不就行了真动起手来才发现难点在于我们根本不知道后验分布长什么样连归一化常数都算不出来。换句话说你知道目标密度的大小关系却不知道它确切的值。MCMC就是在这种情况下硬生生开辟出一条路构造一条马尔可夫链让它的平稳分布等于目标后验然后从链上取样本用样本近似这个分布。本文会把MCMC的原理、选型、Python手写实现以及实际跑PyMC时踩过的坑按我自己的操作路径完整过一遍。适合刚入门贝叶斯推断或者已经会用相关库但想搞清楚内部逻辑的人。1. 为什么要用MCMC做贝叶斯推断1.1 贝叶斯推断的计算瓶颈贝叶斯推断的核心公式其实非常简短后验分布 \(p(\theta | y) \propto p(y | \theta) \, p(\theta)\)右边的分子是似然乘以先验计算成本不高。真正的麻烦在分母就是那个边际似然 \(p(y) \int p(y | \theta) p(\theta) d\theta\)。很多入门教程会轻描淡写地说“分母是归一化常数可以忽略”但如果你要的是概率值而不是一个相对大小这个积分就躲不掉。问题在于这个积分大多数情况下根本没有解析解。简单模型比如共轭先验的线性回归可以直接手推一旦换成逻辑回归、混合模型、随机效应模型后验立刻变成一个高维、非标准、形状复杂的分布。积分在高维空间里几乎不可能用数值方法精确求解。有人说贝叶斯统计的“贝叶斯”部分其实不难难的是“算”这句话我深有体会。打个比方你想知道某个小区里穿红色衣服的人占多少比例但你手里没有小区住户名单。你只能在小区里来回逛遇到一个人就记一下他的衣服颜色逛得足够久统计出来的比例就会接近真实值。MCMC干的就是这件事——它不直接计算分布而是在分布里“逛街”用逛出来的样本频率近似真实概率。1.2 蒙特卡洛与马尔可夫链的碰撞既然直接算积分不行那就采样。蒙特卡洛方法的核心思想非常简单用样本均值近似期望。如果我能从后验分布里抽出 N 个样本那么后验的均值、方差、分位数都能用样本统计量来估计。问题是怎么从“只知道相对密度、不知道归一化常数”的分布里抽样直接采样需要分布有解析的累积分布函数拒绝采样在高维空间里效率会指数级下降。这时候马尔可夫链就派上用场了。马尔可夫链是一个随机过程下一个状态只依赖当前状态和更早的历史无关。我们构造一个转移核让状态按照某种规则不断移动。关键的设计目标是这条链长期运行后状态出现的频率恰好等于目标分布。只要满足这个条件我就不需要知道归一化常数也能得到目标分布的样本。这个想法第一次接触会觉得有点绕但一旦想通会觉得很漂亮。链条在状态空间里游走走到高概率区域的时间长在低概率区域一晃而过。最终记录下来的状态轨迹就像是目标分布的一个“指纹”它天然带有了目标分布的形态信息。1.3 三个基础概念先理清开始上手前有几个高频概念必须搞清楚不然后面看轨迹图、看诊断指标会一头雾水。第一个是平稳分布。马尔可夫链在转移核作用下状态分布会逐渐稳定下来达到一个不再变化的分布这就是平稳分布。MCMC的目标就是构造转移核让目标后验成为这条链的平稳分布。第二个是 burn-in中文常叫预烧期或老化期。链条起点是任意选的初始状态大概率不在高概率区域前面的样本不能反映目标分布需要扔掉。扔掉的这部分就叫 burn-in。我习惯把 warmup 和 burn-in 混着叫但在PyMC里tune 阶段就是专门用来让链条收敛和自适应步长的这个阶段的样本不会保留。第三个是自相关与有效样本量 ESSEffective Sample Size。MCMC的样本不是独立的相邻样本之间有相关性。自相关性越高样本携带的信息越少有效样本量就越低。你可以有10000个样本但 ESS 只有500那实际效果相当于只有500个独立样本。这直接关系到后面说到的收敛性诊断。采样方法原理主要限制直接采样从目标分布直接生成独立样本需要知道分布且能解析采样接受-拒绝采样用提议分布包住目标分布高维时接受率指数下降重要性采样用加权平均近似期望权重容易发散局部失控MCMC构造马尔可夫链逼近目标分布样本自相关需要诊断收敛性2. MCMC算法家族与选型思路2.1 Metropolis-Hastings理解MCMC的入门钥匙Metropolis-HastingsMH是几乎所有MCMC入门教程的起点也是我自己手写过的第一个采样器。它的逻辑很直观从当前状态出发按一个提议分布随机跳到一个新状态然后以一定概率接受这个新状态。接受概率写成\[\alpha \min\left(1, \frac{p(\theta_{\text{new}})}{p(\theta_{\text{old}})} \times \frac{q(\theta_{\text{old}} | \theta_{\text{new}})}{q(\theta_{\text{new}} | \theta_{\text{old}})}\right)\]这里的 q 是提议分布。如果提议分布是对称的比如用正态分布 \(N(0, \sigma^2)\) 作为增量那后面那个比值正好抵消接受率就简化为目标密度的比值。注意这句话我们只需要知道目标密度的比值所以分母归一化常数根本不用算。这就是MH能工作的核心原因。实际写代码时我强烈建议在 log 域里计算。因为直接算概率密度高维后验很容易产生极小值下溢到0然后 log 0 直接报错。用对数概率相减数值稳定得多。这个细节我第一版代码就踩过坑后面会专门说。MH采样器的效果高度依赖提议分布的步长。步长太小接受率高但链条像蜗牛一样爬自相关巨大要很久才能探索完整分布。步长太大提议点总跳到概率极低的区域接受率低链条经常原地踏步。经验上对称随机游走MH的接受率在20%到50%之间比较合理。当然如果用了自适应步长或者更先进的算法这个经验值就不一定适用了。2.2 Gibbs采样当条件分布好抽的时候Gibbs采样是MH的一个特例它的核心思想是与其在高维空间里直接跳不如一次只更新一个变量其他变量固定不变。每次更新时从当前变量的条件后验分布中抽取新值。Gibbs采样的优势在于它不需要手动调步长也没有接受/拒绝环节。你只需要能写出每个变量的条件分布。比如在二元正态分布中如果 \((X, Y)\) 服从均值为0、相关系数为 \(\rho\) 的标准二元正态那么条件分布是\[X | Y y \sim N(\rho y, \, 1 - \rho^2)\]这个条件分布是解析的直接抽就行。经典的蒙特卡洛教程里都会用这个例子演示Gibbs采样。我自己跑过一次之后最大的感受是Gibbs每次只朝一个方向移动虽然探索效率不算高但每步都落在高概率区域不会浪费样本。它是理解“分而治之”思想的最佳教材。但现实世界的模型没那么温柔。很多参数的条件分布根本抽不了这时要么用MH对单个参数更新要么直接用下一节说的HMC。2.3 HMC与NUTS现代贝叶斯推断的主力如果你只用过PyMC或者Stan那你实际用的多半不是MH而是HMCHamiltonian Monte Carlo或者它的自动化版本NUTSNo-U-Turn Sampler。HMC和MH最大的区别是它利用了目标分布的梯度信息。想象一个无摩擦的小球在势能场里滚动势能就是负的log后验概率密度。小球从当前位置出发沿梯度方向积累动量滑向能量低、概率高的区域。因为这个过程带有“惯性”它可以沿同一方向连续移动好几步不会像随机游走那样来回震荡。这正是普通MH在高维空间里表现不佳的根本原因——随机游走需要 \(O(n^2)\) 甚至更多步才能探索完空间而HMC可以一路滑过去。NUTS是HMC的改进它自动选择轨迹长度不需要手动调。PyMC和Stan把NUTS默认配置得已经很好了但对复杂模型仍然需要通过 target_accept 参数来控制步长。我习惯把 target_accept 提到0.9以上尤其在数据量不大但模型后验比较尖锐的时候。2.4 算法选型建议算法适合场景优点缺点Metropolis-Hastings低维、分布形状简单实现简单通用性强高维效率低需调步长Gibbs采样条件分布可解析无接受率参数少不适用所有模型HMC/NUTS中等维度、梯度可计算效率高适合复杂后验需要自动微分支持3. Python代码实现从手写MH到PyMC实战3.1 环境准备我使用的环境是 Python 3.9 及以上版本。基础包是 NumPy 和 SciPy画图用 Matplotlib做贝叶斯模型用 PyMC 和 ArviZ。pip install numpy scipy matplotlib arviz pymc如果你还没装 Python建议直接用 Anaconda 或者官方 Python 安装包都行。安装的时候一定注意 Python 版本PyMC 对版本有要求太老的 3.6 版本装不上新版 PyMC。装完可以用python -c import pymc; print(pymc.__version__)验证一下。3.2 手写 Metropolis-Hastings从一个混合高斯分布采样我第一步练习用的目标分布是一个混合高斯分布密度函数长这样import numpy as np from scipy.stats import norm np.random.seed(42) def target_log_prob(x): # 双峰混合高斯0.4 * N(-3, 0.8) 0.6 * N(2, 1.2) log_density np.log( 0.4 * norm.pdf(x, loc-3, scale0.8) 0.6 * norm.pdf(x, loc2, scale1.2) ) return log_density然后实现一个最简单的随机游走MHdef metropolis_hastings(n_samples10000, warmup1000, init0.0, proposal_scale1.0): samples [] x_current init accepted 0 for _ in range(warmup n_samples): # 提议分布以当前点为中心的正态分布 x_proposal x_current np.random.normal(0, proposal_scale) # 计算接受概率对数域 log_alpha target_log_prob(x_proposal) - target_log_prob(x_current) # 接受或拒绝 if np.log(np.random.rand()) log_alpha: x_current x_proposal accepted 1 if _ warmup: samples.append(x_current) return np.array(samples), accepted / (warmup n_samples) samples, accept_rate metropolis_hastings(proposal_scale1.2) print(fAcceptance rate: {accept_rate:.2f})运行完你会发现即便这么简单的模型如果 proposal_scale 设得过小比如0.1链条要跑很久才能在两个峰之间跳转一次设得过大比如5大部分提议点都落在密度极低的区域接受率会掉到10%以下。我调了几次1.2左右在这个例子里表现不错接受率大约在40%上下。把样本画成直方图能看到两个峰的位置和高度比例都和真实分布吻合。这个练习强烈建议亲手做一遍它能让你直观理解MH内部发生了什么。3.3 手写MH实现贝叶斯线性回归等你对单变量MH熟了就可以直接上贝叶斯线性回归。生成模拟数据np.random.seed(123) true_intercept, true_slope, true_sigma 1.5, 2.0, 1.0 x np.random.uniform(-3, 3, size100) y true_intercept true_slope * x np.random.normal(0, true_sigma, sizelen(x))我们要推断的参数是截距 intercept、斜率 slope 和误差标准差 sigma。为了方便采样我让 sigma 在对数空间里计算这样它恒为正。def log_posterior(params, x, y): intercept, slope, log_sigma params sigma np.exp(log_sigma) resid y - (intercept slope * x) log_likelihood -len(x) * np.log(sigma) - 0.5 * np.sum((resid / sigma) ** 2) log_prior -0.5 * (intercept ** 2 slope ** 2) / 100 # 弱先验 return log_likelihood log_prior采样时直接用多元正态提议分布def bayesian_linear_regression_mh(n_samples10000, warmup1000, step_size0.1): params np.array([0.0, 0.0, 0.0]) # intercept, slope, log_sigma samples [] accepted 0 for i in range(warmup n_samples): proposal params np.random.normal(0, step_size, size3) log_alpha (log_posterior(proposal, x, y) - log_posterior(params, x, y)) if np.log(np.random.rand()) log_alpha: params proposal accepted 1 if i warmup: samples.append(params) return np.array(samples), accepted / (warmup n_samples) samples, accept_rate bayesian_linear_regression_mh(step_size0.1) print(fAcceptance rate: {accept_rate:.2f})这里有几个细节值得注意。第一提议分布的对角协方差 0.1 是手工调的如果你换数据这个值很可能要重调。第二我一开始没对 y 和 x 做标准化导致 slope 的后验和 intercept 的后验尺度差异很大同一个 step_size 对两个参数的感受完全不同采样效率很低。后来我先对 x 做了标准化再把结果反算回去情况改善很多。这就是为什么 PyMC 和 Stan 内部会对参数空间做各种变换——它们在做和人类一样的调整。3.4 用PyMC实现同一个模型有了手写MH的基础再用PyMC就会觉得亲切很多。PyMC内部用NUTS采样器代码量少一个量级import pymc as pm import arviz as az with pm.Model() as model: intercept pm.Normal(intercept, mu0, sigma10) slope pm.Normal(slope, mu0, sigma10) sigma pm.HalfNormal(sigma, sigma5) mu intercept slope * x obs pm.Normal(obs, mumu, sigmasigma, observedy) trace pm.sample( draws1000, tune1000, chains4, target_accept0.9, random_seed42 )这里的参数解释一下draws 是正式采样的样本数tune 是预热样本数chains 是并行跑的独立链数target_accept 是目标接受率。跑完以后用 ArviZ 看结果az.summary(trace)你会看到每个参数的均值、标准差、HDI区间和 Rhat。Rhat 的一项核心判断标准是接近1如果明显大于1.01就说明链条还没收敛。3.5 后验预测与不确定度贝叶斯推断的一个独特优势是能给出预测区间而不只是点估计。采样完成后可以从后验预测分布中抽取新数据with model: posterior_predictive pm.sample_posterior_predictive( trace, random_seed42 )这里的预测样本会体现参数不确定性和观测噪声两部分。画出来以后你能看到一条回归线周围一个完整的带状区间。我第一次跑通这个流程时最大的震撼就是频率派的置信区间和贝叶斯的可信区间在解释上完全不同但数据结构好的时候数值往往很接近。理解这个差异你才能真正理解贝叶斯推断的意义。4. 收敛性诊断与采样细节4.1 轨迹图怎么看才算正常轨迹图trace plot是判断采样是否正常的第一个依据。把采样链上的值按迭代次数画出来你会看到几种典型形态。正常的轨迹图像一条挤满小毛刺的毛毛虫波动平稳没有明显的趋势也没有长时间停留在某个区域。如果轨迹图出现明显的分段跳跃或者一条链稳定呆在一个区域、另一条链呆在另一个区域那大概率是收敛出了问题。我还有一个习惯第一次跑模型时会故意把 draws 设小一点先跑几百样本快速看一眼轨迹图。如果图里出现明显的单峰抖动或者断裂我再决定是调整模型还是调高 target_accept。这比一上来就跑几千样本节省很多时间。burn-in 的长度并没有固定的公式我通常的做法是让 tune 和 draws 等长比如各1000。如果模型复杂后验相关性很高我会把 tune 提到 draws 的两倍。4.2 Rhat 和有效样本量 ESSRhat 也叫 Gelman-Rubin 统计量核心逻辑是比较多重独立链之间的方差和单条链内部的方差。如果多条链都收敛到了同一个分布那么链间方差应该很小Rhat 接近 1。如果有的链跑到了不同的区域链间方差就会特别大Rhat 也会明显变大。我在实际项目里的阈值是Rhat 小于 1.01 才放心。有些教程说小于 1.1 就行但对比较重要的结论我倾向于更严格的标准。ESS 这个指标同样重要。ArviZ 输出的 summarize 表里ess_bulk 和 ess_tail 分别代表后验主体和尾部的有效样本量。如果 ESS 太低后验的分位数估计就会不稳定。增加 draws 当然能提高 ESS但更聪明的做法是先检查自相关性。如果自相关系数衰减得很慢那问题可能出在模型参数化上而不是样本量不够。4.3 target_accept 与自适应步长MH 需要手动调步长HMC 也不是完全不需要管步长。NUTS 内部会自适应调整步长调整的依据就是拒绝率和发散数量。PyMC 默认的 target_accept 是 0.8但对后验比较尖锐或者参数尺度差异大的模型默认值很容易产生发散。我跑这类模型时会把 target_accept 提到 0.9 甚至 0.95代价是采样时间明显增加但换来的是更少的发散和更稳定的结果。一个需要特别注意的坑当你把 target_accept 调到很高采样速度会慢很多但 Rhat 依然不满意。这时候不要再无脑加 target_accept 了应该回去检查模型设定比如是不是先验太强、数据是否有异常值、参数化是不是有问题。4.4 采样实战心得这几条经验是多次跑崩以后总结出来的写在这里供参考。第一对数据进行标准化。把 x 和 y 的尺度大致统一到一个量级后验的曲率更均匀采样器更容易探索。用完模型以后再把斜率结果反变换回原始尺度这个操作很简单但收益很大。第二参数化要尽量“去相关”。参数之间相关性越强采样器越难走。比如随机效应模型直接抽随机截距和随机斜率相关性经常很高改用非中心化参数化后采样效率能提升好几倍。这个概念需要单独写一篇但记住一个方向把模型里的强相关参数想办法解耦。第三跑多链、留种子。我至少跑4条链并且固定 random_seed这样别人可以复现我的结果我自己排查问题也方便。第四先小规模试跑再大规模采样。我见过很多人在第一次跑模型时就设 draws20000结果等了十分钟发现 Rhat 爆炸。先跑 200 样本快速定位问题比硬跑长链高效得多。5. 常见问题与排查技巧实录5.1 Rhat 不收敛怎么办这个是我被问过最多的问题。Rhat 大于 1.01 时先别急着加样本量。按照这个顺序排查先看轨迹图如果多条链各自漂移在不同区域说明后验可能是多峰的或者链没有跳到高概率区域。此时可以增加 tune 阶段、提高 target_accept或者改用更宽的先验。如果轨迹图看起来正常但 Rhat 依然偏高那很可能是参数化导致的采样效率差。尝试把相关参数转换成弱相关的表示形式。我记得有一次跑一个带随机斜率的模型Rhat 死活到不了 1.01后来改用非中心化参数化一次就收敛了。还有一种情况是后验本身就存在多个不连通的模态。这是MCMC的天然难点任何诊断方法都很难解决。这种情况只能改用其他方法比如做退火采样或者直接换模型。5.2 接受率过高或过低怎么调如果你是手写MH接受率就是直接衡量步长大小的指标。接受率过高说明提议分布步长太小链条自相关太高接受率过低说明提议步长太大链条在拒绝中空转。经验上随机游走MH的接受率保持在 20%-50% 附近比较理想。调法很简单接受率偏高就增大提议分布的标准差偏低就减小。如果一开始对尺度没有概念可以先从 0.1 开始观察几十步再逐步调整。你还可以在迭代过程中采样自适应调节步长但要注意不要破坏马尔可夫性的推导基础——直接套现成的自适应库会减少很多麻烦。5.3 报错Initial evaluation of model at starting point failed这个报错在跑PyMC时很常见字面意思是在初始点评估模型时失败。通常原因在于参数初始值设置了不支持的取值或者方差参数初始化为负数之类。解决方法是检查模型里的参数约束。例如sigma pm.HalfNormal(...)就隐含了非负约束不需要再手动限制。还有一次我遇到这个报错原因是数据里有 NaN。检查数据、删除异常样本问题就解决了。这个报错虽然吓人但一般都不是模型设计的大问题而是数据或初始化的小问题。5.4 常用问题速查表问题现象可能原因常用处理手段Rhat 1.01收敛不充分链未混合增加 tune提升 target_accept改参数化轨迹图有明显分段跳跃多峰后验或步长过大换采样器、检查后验模态采样速度极慢后验相关性强、数据量太大标准化、参数去相关、减少样本量ESS 很低自相关强有效样本量不足增加 draws、降低自相关、调整参数化出现发散divergence步长不适配后验曲率提高 target_accept 到 0.95初始评估失败初始值非法、数据含NaN检查数据、调整初值、检查约束后验结果对初始值敏感链没有充分收敛增大 warmup增加链数我在实际使用中最大的体会是MCMC是一个非常依赖经验的工具理论和代码只是基础真正值钱的技能是“看样本判断模型状态”的能力。手写一次MH亲眼见过什么叫做“链没跑偏”再去看PyMC输出的时候心里会踏实很多。如果你刚开始接触贝叶斯推断建议先不要急着套高级库先用纯NumPy手写一个最简单的采样器跑通混合高斯、跑通线性回归再回到PyMC的舒适环境。这个弯路我走过值得走。本文还有配套的精品资源点击获取 MCMC贝叶斯推断入门：从手写采样器到PyMC实战
简介在机器学习与统计建模中贝叶斯推断提供了一套将先验知识与观测数据融合的完备框架然而真实模型的后验分布往往没有解析解归一化常数也难以计算。蒙特卡洛采样与马尔可夫链的结合——MCMC方法成为解决该问题的核心工具。通过构造一条以目标后验为平稳分布的随机游走链MCMC能够在高维空间中高效采集样本从而近似后验的均值、方差与分位数。从基础的Metropolis-Hastings算法到现代概率编程框架PyMC中内置的NUTS采样器MCMC已广泛用于回归分析、分层模型、深度学习不确定性量化等场景。理解其原理与诊断指标如Rhat、ESS是实践的关键。本文从手写采样器出发系统地讲解MCMC算法家族、Python实现及实战调参经验帮助读者真正掌握贝叶斯推断的计算逻辑。 做贝叶斯推断的人迟早都会撞上MCMC。我最早接触这个概念是在看一篇分层模型的论文时模型本身很简单但后验没有解析形式作者直接用了MCMC采样。当时我天真地想不就是采样嘛直接从分布里抽不就行了真动起手来才发现难点在于我们根本不知道后验分布长什么样连归一化常数都算不出来。换句话说你知道目标密度的大小关系却不知道它确切的值。MCMC就是在这种情况下硬生生开辟出一条路构造一条马尔可夫链让它的平稳分布等于目标后验然后从链上取样本用样本近似这个分布。本文会把MCMC的原理、选型、Python手写实现以及实际跑PyMC时踩过的坑按我自己的操作路径完整过一遍。适合刚入门贝叶斯推断或者已经会用相关库但想搞清楚内部逻辑的人。1. 为什么要用MCMC做贝叶斯推断1.1 贝叶斯推断的计算瓶颈贝叶斯推断的核心公式其实非常简短后验分布 \(p(\theta | y) \propto p(y | \theta) \, p(\theta)\)右边的分子是似然乘以先验计算成本不高。真正的麻烦在分母就是那个边际似然 \(p(y) \int p(y | \theta) p(\theta) d\theta\)。很多入门教程会轻描淡写地说“分母是归一化常数可以忽略”但如果你要的是概率值而不是一个相对大小这个积分就躲不掉。问题在于这个积分大多数情况下根本没有解析解。简单模型比如共轭先验的线性回归可以直接手推一旦换成逻辑回归、混合模型、随机效应模型后验立刻变成一个高维、非标准、形状复杂的分布。积分在高维空间里几乎不可能用数值方法精确求解。有人说贝叶斯统计的“贝叶斯”部分其实不难难的是“算”这句话我深有体会。打个比方你想知道某个小区里穿红色衣服的人占多少比例但你手里没有小区住户名单。你只能在小区里来回逛遇到一个人就记一下他的衣服颜色逛得足够久统计出来的比例就会接近真实值。MCMC干的就是这件事——它不直接计算分布而是在分布里“逛街”用逛出来的样本频率近似真实概率。1.2 蒙特卡洛与马尔可夫链的碰撞既然直接算积分不行那就采样。蒙特卡洛方法的核心思想非常简单用样本均值近似期望。如果我能从后验分布里抽出 N 个样本那么后验的均值、方差、分位数都能用样本统计量来估计。问题是怎么从“只知道相对密度、不知道归一化常数”的分布里抽样直接采样需要分布有解析的累积分布函数拒绝采样在高维空间里效率会指数级下降。这时候马尔可夫链就派上用场了。马尔可夫链是一个随机过程下一个状态只依赖当前状态和更早的历史无关。我们构造一个转移核让状态按照某种规则不断移动。关键的设计目标是这条链长期运行后状态出现的频率恰好等于目标分布。只要满足这个条件我就不需要知道归一化常数也能得到目标分布的样本。这个想法第一次接触会觉得有点绕但一旦想通会觉得很漂亮。链条在状态空间里游走走到高概率区域的时间长在低概率区域一晃而过。最终记录下来的状态轨迹就像是目标分布的一个“指纹”它天然带有了目标分布的形态信息。1.3 三个基础概念先理清开始上手前有几个高频概念必须搞清楚不然后面看轨迹图、看诊断指标会一头雾水。第一个是平稳分布。马尔可夫链在转移核作用下状态分布会逐渐稳定下来达到一个不再变化的分布这就是平稳分布。MCMC的目标就是构造转移核让目标后验成为这条链的平稳分布。第二个是 burn-in中文常叫预烧期或老化期。链条起点是任意选的初始状态大概率不在高概率区域前面的样本不能反映目标分布需要扔掉。扔掉的这部分就叫 burn-in。我习惯把 warmup 和 burn-in 混着叫但在PyMC里tune 阶段就是专门用来让链条收敛和自适应步长的这个阶段的样本不会保留。第三个是自相关与有效样本量 ESSEffective Sample Size。MCMC的样本不是独立的相邻样本之间有相关性。自相关性越高样本携带的信息越少有效样本量就越低。你可以有10000个样本但 ESS 只有500那实际效果相当于只有500个独立样本。这直接关系到后面说到的收敛性诊断。采样方法原理主要限制直接采样从目标分布直接生成独立样本需要知道分布且能解析采样接受-拒绝采样用提议分布包住目标分布高维时接受率指数下降重要性采样用加权平均近似期望权重容易发散局部失控MCMC构造马尔可夫链逼近目标分布样本自相关需要诊断收敛性2. MCMC算法家族与选型思路2.1 Metropolis-Hastings理解MCMC的入门钥匙Metropolis-HastingsMH是几乎所有MCMC入门教程的起点也是我自己手写过的第一个采样器。它的逻辑很直观从当前状态出发按一个提议分布随机跳到一个新状态然后以一定概率接受这个新状态。接受概率写成\[\alpha \min\left(1, \frac{p(\theta_{\text{new}})}{p(\theta_{\text{old}})} \times \frac{q(\theta_{\text{old}} | \theta_{\text{new}})}{q(\theta_{\text{new}} | \theta_{\text{old}})}\right)\]这里的 q 是提议分布。如果提议分布是对称的比如用正态分布 \(N(0, \sigma^2)\) 作为增量那后面那个比值正好抵消接受率就简化为目标密度的比值。注意这句话我们只需要知道目标密度的比值所以分母归一化常数根本不用算。这就是MH能工作的核心原因。实际写代码时我强烈建议在 log 域里计算。因为直接算概率密度高维后验很容易产生极小值下溢到0然后 log 0 直接报错。用对数概率相减数值稳定得多。这个细节我第一版代码就踩过坑后面会专门说。MH采样器的效果高度依赖提议分布的步长。步长太小接受率高但链条像蜗牛一样爬自相关巨大要很久才能探索完整分布。步长太大提议点总跳到概率极低的区域接受率低链条经常原地踏步。经验上对称随机游走MH的接受率在20%到50%之间比较合理。当然如果用了自适应步长或者更先进的算法这个经验值就不一定适用了。2.2 Gibbs采样当条件分布好抽的时候Gibbs采样是MH的一个特例它的核心思想是与其在高维空间里直接跳不如一次只更新一个变量其他变量固定不变。每次更新时从当前变量的条件后验分布中抽取新值。Gibbs采样的优势在于它不需要手动调步长也没有接受/拒绝环节。你只需要能写出每个变量的条件分布。比如在二元正态分布中如果 \((X, Y)\) 服从均值为0、相关系数为 \(\rho\) 的标准二元正态那么条件分布是\[X | Y y \sim N(\rho y, \, 1 - \rho^2)\]这个条件分布是解析的直接抽就行。经典的蒙特卡洛教程里都会用这个例子演示Gibbs采样。我自己跑过一次之后最大的感受是Gibbs每次只朝一个方向移动虽然探索效率不算高但每步都落在高概率区域不会浪费样本。它是理解“分而治之”思想的最佳教材。但现实世界的模型没那么温柔。很多参数的条件分布根本抽不了这时要么用MH对单个参数更新要么直接用下一节说的HMC。2.3 HMC与NUTS现代贝叶斯推断的主力如果你只用过PyMC或者Stan那你实际用的多半不是MH而是HMCHamiltonian Monte Carlo或者它的自动化版本NUTSNo-U-Turn Sampler。HMC和MH最大的区别是它利用了目标分布的梯度信息。想象一个无摩擦的小球在势能场里滚动势能就是负的log后验概率密度。小球从当前位置出发沿梯度方向积累动量滑向能量低、概率高的区域。因为这个过程带有“惯性”它可以沿同一方向连续移动好几步不会像随机游走那样来回震荡。这正是普通MH在高维空间里表现不佳的根本原因——随机游走需要 \(O(n^2)\) 甚至更多步才能探索完空间而HMC可以一路滑过去。NUTS是HMC的改进它自动选择轨迹长度不需要手动调。PyMC和Stan把NUTS默认配置得已经很好了但对复杂模型仍然需要通过 target_accept 参数来控制步长。我习惯把 target_accept 提到0.9以上尤其在数据量不大但模型后验比较尖锐的时候。2.4 算法选型建议算法适合场景优点缺点Metropolis-Hastings低维、分布形状简单实现简单通用性强高维效率低需调步长Gibbs采样条件分布可解析无接受率参数少不适用所有模型HMC/NUTS中等维度、梯度可计算效率高适合复杂后验需要自动微分支持3. Python代码实现从手写MH到PyMC实战3.1 环境准备我使用的环境是 Python 3.9 及以上版本。基础包是 NumPy 和 SciPy画图用 Matplotlib做贝叶斯模型用 PyMC 和 ArviZ。pip install numpy scipy matplotlib arviz pymc如果你还没装 Python建议直接用 Anaconda 或者官方 Python 安装包都行。安装的时候一定注意 Python 版本PyMC 对版本有要求太老的 3.6 版本装不上新版 PyMC。装完可以用python -c import pymc; print(pymc.__version__)验证一下。3.2 手写 Metropolis-Hastings从一个混合高斯分布采样我第一步练习用的目标分布是一个混合高斯分布密度函数长这样import numpy as np from scipy.stats import norm np.random.seed(42) def target_log_prob(x): # 双峰混合高斯0.4 * N(-3, 0.8) 0.6 * N(2, 1.2) log_density np.log( 0.4 * norm.pdf(x, loc-3, scale0.8) 0.6 * norm.pdf(x, loc2, scale1.2) ) return log_density然后实现一个最简单的随机游走MHdef metropolis_hastings(n_samples10000, warmup1000, init0.0, proposal_scale1.0): samples [] x_current init accepted 0 for _ in range(warmup n_samples): # 提议分布以当前点为中心的正态分布 x_proposal x_current np.random.normal(0, proposal_scale) # 计算接受概率对数域 log_alpha target_log_prob(x_proposal) - target_log_prob(x_current) # 接受或拒绝 if np.log(np.random.rand()) log_alpha: x_current x_proposal accepted 1 if _ warmup: samples.append(x_current) return np.array(samples), accepted / (warmup n_samples) samples, accept_rate metropolis_hastings(proposal_scale1.2) print(fAcceptance rate: {accept_rate:.2f})运行完你会发现即便这么简单的模型如果 proposal_scale 设得过小比如0.1链条要跑很久才能在两个峰之间跳转一次设得过大比如5大部分提议点都落在密度极低的区域接受率会掉到10%以下。我调了几次1.2左右在这个例子里表现不错接受率大约在40%上下。把样本画成直方图能看到两个峰的位置和高度比例都和真实分布吻合。这个练习强烈建议亲手做一遍它能让你直观理解MH内部发生了什么。3.3 手写MH实现贝叶斯线性回归等你对单变量MH熟了就可以直接上贝叶斯线性回归。生成模拟数据np.random.seed(123) true_intercept, true_slope, true_sigma 1.5, 2.0, 1.0 x np.random.uniform(-3, 3, size100) y true_intercept true_slope * x np.random.normal(0, true_sigma, sizelen(x))我们要推断的参数是截距 intercept、斜率 slope 和误差标准差 sigma。为了方便采样我让 sigma 在对数空间里计算这样它恒为正。def log_posterior(params, x, y): intercept, slope, log_sigma params sigma np.exp(log_sigma) resid y - (intercept slope * x) log_likelihood -len(x) * np.log(sigma) - 0.5 * np.sum((resid / sigma) ** 2) log_prior -0.5 * (intercept ** 2 slope ** 2) / 100 # 弱先验 return log_likelihood log_prior采样时直接用多元正态提议分布def bayesian_linear_regression_mh(n_samples10000, warmup1000, step_size0.1): params np.array([0.0, 0.0, 0.0]) # intercept, slope, log_sigma samples [] accepted 0 for i in range(warmup n_samples): proposal params np.random.normal(0, step_size, size3) log_alpha (log_posterior(proposal, x, y) - log_posterior(params, x, y)) if np.log(np.random.rand()) log_alpha: params proposal accepted 1 if i warmup: samples.append(params) return np.array(samples), accepted / (warmup n_samples) samples, accept_rate bayesian_linear_regression_mh(step_size0.1) print(fAcceptance rate: {accept_rate:.2f})这里有几个细节值得注意。第一提议分布的对角协方差 0.1 是手工调的如果你换数据这个值很可能要重调。第二我一开始没对 y 和 x 做标准化导致 slope 的后验和 intercept 的后验尺度差异很大同一个 step_size 对两个参数的感受完全不同采样效率很低。后来我先对 x 做了标准化再把结果反算回去情况改善很多。这就是为什么 PyMC 和 Stan 内部会对参数空间做各种变换——它们在做和人类一样的调整。3.4 用PyMC实现同一个模型有了手写MH的基础再用PyMC就会觉得亲切很多。PyMC内部用NUTS采样器代码量少一个量级import pymc as pm import arviz as az with pm.Model() as model: intercept pm.Normal(intercept, mu0, sigma10) slope pm.Normal(slope, mu0, sigma10) sigma pm.HalfNormal(sigma, sigma5) mu intercept slope * x obs pm.Normal(obs, mumu, sigmasigma, observedy) trace pm.sample( draws1000, tune1000, chains4, target_accept0.9, random_seed42 )这里的参数解释一下draws 是正式采样的样本数tune 是预热样本数chains 是并行跑的独立链数target_accept 是目标接受率。跑完以后用 ArviZ 看结果az.summary(trace)你会看到每个参数的均值、标准差、HDI区间和 Rhat。Rhat 的一项核心判断标准是接近1如果明显大于1.01就说明链条还没收敛。3.5 后验预测与不确定度贝叶斯推断的一个独特优势是能给出预测区间而不只是点估计。采样完成后可以从后验预测分布中抽取新数据with model: posterior_predictive pm.sample_posterior_predictive( trace, random_seed42 )这里的预测样本会体现参数不确定性和观测噪声两部分。画出来以后你能看到一条回归线周围一个完整的带状区间。我第一次跑通这个流程时最大的震撼就是频率派的置信区间和贝叶斯的可信区间在解释上完全不同但数据结构好的时候数值往往很接近。理解这个差异你才能真正理解贝叶斯推断的意义。4. 收敛性诊断与采样细节4.1 轨迹图怎么看才算正常轨迹图trace plot是判断采样是否正常的第一个依据。把采样链上的值按迭代次数画出来你会看到几种典型形态。正常的轨迹图像一条挤满小毛刺的毛毛虫波动平稳没有明显的趋势也没有长时间停留在某个区域。如果轨迹图出现明显的分段跳跃或者一条链稳定呆在一个区域、另一条链呆在另一个区域那大概率是收敛出了问题。我还有一个习惯第一次跑模型时会故意把 draws 设小一点先跑几百样本快速看一眼轨迹图。如果图里出现明显的单峰抖动或者断裂我再决定是调整模型还是调高 target_accept。这比一上来就跑几千样本节省很多时间。burn-in 的长度并没有固定的公式我通常的做法是让 tune 和 draws 等长比如各1000。如果模型复杂后验相关性很高我会把 tune 提到 draws 的两倍。4.2 Rhat 和有效样本量 ESSRhat 也叫 Gelman-Rubin 统计量核心逻辑是比较多重独立链之间的方差和单条链内部的方差。如果多条链都收敛到了同一个分布那么链间方差应该很小Rhat 接近 1。如果有的链跑到了不同的区域链间方差就会特别大Rhat 也会明显变大。我在实际项目里的阈值是Rhat 小于 1.01 才放心。有些教程说小于 1.1 就行但对比较重要的结论我倾向于更严格的标准。ESS 这个指标同样重要。ArviZ 输出的 summarize 表里ess_bulk 和 ess_tail 分别代表后验主体和尾部的有效样本量。如果 ESS 太低后验的分位数估计就会不稳定。增加 draws 当然能提高 ESS但更聪明的做法是先检查自相关性。如果自相关系数衰减得很慢那问题可能出在模型参数化上而不是样本量不够。4.3 target_accept 与自适应步长MH 需要手动调步长HMC 也不是完全不需要管步长。NUTS 内部会自适应调整步长调整的依据就是拒绝率和发散数量。PyMC 默认的 target_accept 是 0.8但对后验比较尖锐或者参数尺度差异大的模型默认值很容易产生发散。我跑这类模型时会把 target_accept 提到 0.9 甚至 0.95代价是采样时间明显增加但换来的是更少的发散和更稳定的结果。一个需要特别注意的坑当你把 target_accept 调到很高采样速度会慢很多但 Rhat 依然不满意。这时候不要再无脑加 target_accept 了应该回去检查模型设定比如是不是先验太强、数据是否有异常值、参数化是不是有问题。4.4 采样实战心得这几条经验是多次跑崩以后总结出来的写在这里供参考。第一对数据进行标准化。把 x 和 y 的尺度大致统一到一个量级后验的曲率更均匀采样器更容易探索。用完模型以后再把斜率结果反变换回原始尺度这个操作很简单但收益很大。第二参数化要尽量“去相关”。参数之间相关性越强采样器越难走。比如随机效应模型直接抽随机截距和随机斜率相关性经常很高改用非中心化参数化后采样效率能提升好几倍。这个概念需要单独写一篇但记住一个方向把模型里的强相关参数想办法解耦。第三跑多链、留种子。我至少跑4条链并且固定 random_seed这样别人可以复现我的结果我自己排查问题也方便。第四先小规模试跑再大规模采样。我见过很多人在第一次跑模型时就设 draws20000结果等了十分钟发现 Rhat 爆炸。先跑 200 样本快速定位问题比硬跑长链高效得多。5. 常见问题与排查技巧实录5.1 Rhat 不收敛怎么办这个是我被问过最多的问题。Rhat 大于 1.01 时先别急着加样本量。按照这个顺序排查先看轨迹图如果多条链各自漂移在不同区域说明后验可能是多峰的或者链没有跳到高概率区域。此时可以增加 tune 阶段、提高 target_accept或者改用更宽的先验。如果轨迹图看起来正常但 Rhat 依然偏高那很可能是参数化导致的采样效率差。尝试把相关参数转换成弱相关的表示形式。我记得有一次跑一个带随机斜率的模型Rhat 死活到不了 1.01后来改用非中心化参数化一次就收敛了。还有一种情况是后验本身就存在多个不连通的模态。这是MCMC的天然难点任何诊断方法都很难解决。这种情况只能改用其他方法比如做退火采样或者直接换模型。5.2 接受率过高或过低怎么调如果你是手写MH接受率就是直接衡量步长大小的指标。接受率过高说明提议分布步长太小链条自相关太高接受率过低说明提议步长太大链条在拒绝中空转。经验上随机游走MH的接受率保持在 20%-50% 附近比较理想。调法很简单接受率偏高就增大提议分布的标准差偏低就减小。如果一开始对尺度没有概念可以先从 0.1 开始观察几十步再逐步调整。你还可以在迭代过程中采样自适应调节步长但要注意不要破坏马尔可夫性的推导基础——直接套现成的自适应库会减少很多麻烦。5.3 报错Initial evaluation of model at starting point failed这个报错在跑PyMC时很常见字面意思是在初始点评估模型时失败。通常原因在于参数初始值设置了不支持的取值或者方差参数初始化为负数之类。解决方法是检查模型里的参数约束。例如sigma pm.HalfNormal(...)就隐含了非负约束不需要再手动限制。还有一次我遇到这个报错原因是数据里有 NaN。检查数据、删除异常样本问题就解决了。这个报错虽然吓人但一般都不是模型设计的大问题而是数据或初始化的小问题。5.4 常用问题速查表问题现象可能原因常用处理手段Rhat 1.01收敛不充分链未混合增加 tune提升 target_accept改参数化轨迹图有明显分段跳跃多峰后验或步长过大换采样器、检查后验模态采样速度极慢后验相关性强、数据量太大标准化、参数去相关、减少样本量ESS 很低自相关强有效样本量不足增加 draws、降低自相关、调整参数化出现发散divergence步长不适配后验曲率提高 target_accept 到 0.95初始评估失败初始值非法、数据含NaN检查数据、调整初值、检查约束后验结果对初始值敏感链没有充分收敛增大 warmup增加链数我在实际使用中最大的体会是MCMC是一个非常依赖经验的工具理论和代码只是基础真正值钱的技能是“看样本判断模型状态”的能力。手写一次MH亲眼见过什么叫做“链没跑偏”再去看PyMC输出的时候心里会踏实很多。如果你刚开始接触贝叶斯推断建议先不要急着套高级库先用纯NumPy手写一个最简单的采样器跑通混合高斯、跑通线性回归再回到PyMC的舒适环境。这个弯路我走过值得走。本文还有配套的精品资源点击获取