RJMCMCReversible Jump Markov Chain Monte Carlo在贝叶斯计算圈子里一直是个听着简单做起来想摔键盘的算法。我最早接触它是在做高斯混合模型分量数估计的时候跑了三天三夜后验分布还是跟心电图似的乱跳后来才发现是雅可比行列式算错了。这算法解决的问题很直接当你不知道模型长什么样——不知道有几个分量、几个变量、几个变化点——但又想同时做模型选择和参数估计时传统MCMC只能在一棵树上吊死RJMCMC却能从一棵树跳到另一棵树。这篇文章我会把RJMCMC从原理到代码一步步拆开带着你写出一个能跑的高斯混合模型RJMCMC实现也把那些文档里不会写的坑一并交代清楚。适合正在做贝叶斯模型选择、变化点检测、变量选择或者刚被可逆跳跃折磨得睡不着的朋友参考。1. 从MCMC到RJMCMC为什么要做可变维采样1.1 模型选择问题的本质先想清楚一件事常规MCMC到底在干嘛。假设你有一个固定模型参数是θ似然是p(y|θ)先验是p(θ)你想估计θ的后验分布p(θ|y)。Metropolis-Hastings算法做的事情是在参数空间里随机游走接受或拒绝新提议最终让采样点逼近后验分布。这个过程能工作的前提是θ的维度是固定的。你从一个二维高斯里提一个样本跳到三维空间去MH算法根本没有这个机制。但实际统计问题里模型本身往往是未知的。拿高斯混合模型来说数据到底有几个分量这个k本身就是我们要推断的量。再举个例子时间序列的变化点检测——到底有几个变化点每个变化点的位置和跳跃幅度是多少如果你把k固定成3那拟合出来效果可能还行但你怎么知道k4不会更好这里就涉及一个经典问题模型选择。传统做模型选择有两条路。第一算贝叶斯因子对每个候选模型分别跑MCMC算边际似然然后比大小。问题在于边际似然本身很难算而且当候选模型数量很多比如k从1到100时计算量爆炸。第二用AIC/BIC之类的信息准则但这些都是点估计忽略了参数的不确定性而且惩罚项的选择有相当的主观性。RJMCMC的思路完全不一样它把模型索引也当成一个随机变量扔进采样过程。也就是说状态空间不再是某个单独的参数空间而是所有模型参数的并集。每迭代一步你不仅更新参数还可能从当前模型跳到一个维度更高或更低的模型。这样一来你最终得到的一整条采样链天然就包含了各个模型的后验概率——模型选择问题和参数估计问题一次性解决。1.2 跨越不同维度的状态空间但要实现在不同模型之间跳跃第一个拦路虎就来了维度不匹配。假设当前模型M₁的参数是θ₁维度是3你想跳到模型M₂参数θ₂维度是5。从3维空间跳到一个5维空间套用MH的接受率公式你必须计算提议密度之比但两个不同维度空间里的概率密度根本没法直接比——数学上这俩是可测空间但密度定义在完全不同的测度上。Green在1995年那篇论文里的核心贡献就是把这个问题从数学层面破解了。他的思路非常巧妙与其在两个不同维度的空间之间硬跳不如把跳跃过程放到一个维度匹配的扩展空间里来做。简单说你把低维的参数加上一个辅助变量u让(θ₁,u)和θ₂的维度一样然后构造一个可逆的双射函数θ₂ g(θ₁, u)。这样低维到高维的跳跃就变成了同维度空间里的一个确定性映射加一个随机抽样维度的问题被绕过去了。这就是可逆跳跃这个名字的来历。不是跳跃本身可逆而是构建了一个可逆的映射结构让正向跳跃和反向跳跃在扩展空间里形成一一对应的关系。理解这一点后面接受率的推导就顺理成章。2. 算法核心原理可逆跳跃的设计逻辑2.1 可逆跳跃与维数匹配现在我用一个具体的例子把维数匹配讲透。假设你要从模型M₁参数θ₁维度d₁跳到模型M₂参数θ₂维度d₂且d₁ d₂。维数匹配的标准做法是从某个辅助分布q(u)里抽一个随机数uu的维度等于d₂ - d₁。然后构造函数g(θ₁, u) → θ₂。这个函数必须是一个微分同胚也就是可逆且光滑。反过来从M₂跳回M₁的时候你必须能从θ₂里恢复出θ₁和u那就是g的逆映射。注意一个细节反向跳跃的时候u不是随便抽出来的。它是由θ₂通过逆映射唯一确定的所以反向提议的密度里u对应的那部分实际上是由雅可比行列式贡献的。这就是为什么RJMCMC的接受率里会莫名其妙多出一个雅可比项。很多人第一次写代码时把这一项漏掉结果高维模型的采样结果完全偏离真实后验但又不完全错到处找bug都找不出来。我自己常用的一个验证方法写一个维度差为1的简单跳跃比如从常数函数跳到一个分段常数函数。这样辅助变量只是一个均匀随机数雅可比的计算也简单清晰可以手动验证接受率是否正确然后再上复杂模型。2.2 接受率推导与Hastings比RJMCMC的接受率完整形式长这样α min(1, 后验比 × 提议比 × 雅可比)其中后验比 [p(y|θ₂, M₂) p(θ₂|M₂) p(M₂)] / [p(y|θ₁, M₁) p(θ₁|M₁) p(M₁)]提议比 [q_M₁₂(M₂→M₁的跳转概率) × q(u的反向密度)] / [q_M₁→M₂(跳转概率) × q(u的提议密度)]雅可比 |∂g(θ₁,u) / ∂(θ₁,u)|每一项的意义都不难理解。后验比是常规的似然和先验之比对应MH里的后验项。提议比是跳转概率之比也就是你以多大的概率选择从M₁跳到M₂又选择了什么样的辅助变量分布。最关键、也最容易被忽略的是雅可比项它代表了从(θ₁,u)到θ₂的坐标变换对体积元素的影响。举个例子。你把一个权重为w的分量分裂成两个权重为w₁和w₂的分量并令w₁ w·u, w₂ w·(1-u)u是[0,1]上抽的。这个映射的雅可比行列式计算出来就是w。如果忘了乘这个w接受率就偏了分裂移动的频率就会失真。还有一点容易忽视在做分裂/合并这类移动时提议分布的设计不是随便选的。最常见的是用矩匹配的方法保证分裂前后的前几阶矩一致。比如均值、方差、权重都要有对应的映射关系。如果设计得不好接受率会低得离谱链一直在原地打转。3. 实操设计以高斯混合模型为例3.1 模型设定与移动类型划分理论说再多不如动手写个完整案例。高斯混合模型是RJMCMC最经典的练手场景因为它有非常自然的生灭和分裂合并操作。模型设定如下假设观测数据y (y₁, ..., yₙ)服从一个k分量高斯混合分布其中k未知。每个分量的参数是(μⱼ, σⱼ², wⱼ)wⱼ是混合权重满足Σwⱼ1。模型M_k的参数空间维度是3k-1因为权重有和约束所以少一个自由度。我们需要设计的移动类型有四类参数更新Param Update在给定k的情况下用MH或Gibbs更新每个分量的参数。生灭移动Birth Death增加或删除一个分量。分裂合并移动Split Merge把一个分量劈成两个或把两个分量合成一个。权重调整Weight Move有时候单独放在生灭里一起做有时设计成独立步骤让权重空间探索更充分。四类移动缺一不可。生灭移动负责探索整个模型空间的大尺度结构分裂合并则在局部细节上更高效。两者配合才能保证可逆性。比如你做了一个Birth理论上对应的是Death但Birth从一个空位置生成一个新分量直接删掉它几乎不可能被接受。这时候就需要Split/Merge来互补——Split是在现有分量的基础上做精细调整Merge是它的逆操作两者形成高质量的可逆对。3.2 分裂/合并与生灭移动的实现细节先看Birth怎么做。我说说最常见的实现方式在数据范围内生成一个新分量。具体来说从某个先验分布里抽取新分量的权重w_new再在数据范围内均匀抽取均值μ_new方差σ_new²从某个逆Gamma先验里抽。但这里有个容易掉的坑直接生成w_new会遇到约束问题。混合权重有和等于1的约束插进来一个新权重已有的权重必须按比例缩放。这就要用到可逆映射的设计思路。我在代码里常用的做法是这样的。假设当前有k个分量权重为w₁, ..., wₖ。Birth时抽一个u ~ Beta(1, k)新分量的权重设为w_new u现有权重缩放为wⱼ wⱼ·(1-u)。这样就保证了Σwⱼ w_{new} 1。在计算雅可比时要考虑到这个从(w₁,...,wₖ, u)到(w₁,...,wₖ, w_{new})的变换。具体算出来会发现雅可比里含有一个(1-u)^{k-1}项这个项在双向跳跃中都需要仔细处理。Split移动则更精细。假设你要把第j个分量分裂成两个分量可逆映射的设计要保证分裂前后的一些重要特征尽可能保持一致。常用的方案是矩匹配令分裂前后总体的均值、方差、权重之和相等。具体来说假设原分量参数为(w, μ, σ²)分裂成两个分量(w₁, μ₁, σ₁²)和(w₂, μ₂, σ₂²)。先抽两个辅助变量u₁, u₂ ~ Beta(2,2)然后做映射w₁ w·u₁, w₂ w·(1-u₁) μ₁ μ - u₂·σ√(w₂/w₁), μ₂ μ u₂·σ√(w₁/w₂) σ₁² u₁(1-u₂²)·σ²·(w/w₁), σ₂² (1-u₁)(1-u₂²)·σ²·(w/w₂)这样设计的妙处在于它保证了分裂前后一阶矩总均值和二阶中心矩总方差不变。你可以在代码里用数值验证一下分裂前后的期望和方差确实保持不变。这在很大程度上提高了Split移动的接受率因为数据下的似然不会因为分裂而产生剧烈变化。4. 代码级实现核心函数与关键参数4.1 提议分布与辅助变量采样整个RJMCMC最核心的部分就是设计提议分布。提议分布直接决定了链的混合效率和接受率。常见的提议分布选择如下辅助变量u一般选Beta分布。Beta参数的选择会显著影响分裂/合并的接受率。我试过Beta(2,2)和Beta(1,1)即均匀分布前者偏重中间值更保守接受率略高但移动步长小后者更容易产生极端的权重比活跃但接受率低。实际中要根据数据量调。新分量的均值在观测数据范围内均匀抽样或者用经验分布抽样。如果数据有明显的聚簇特征可以考虑用数据的某个随机子集的均值作为提议均值这能显著提高Birth的接受率。新分量的方差从逆Gamma(a, b)先验里抽a和b根据自己的先验知识设定。如果完全没有信息也可以直接用数据整体方差的某个比例。再提醒一个关键点辅助变量的维度必须严格等于d₂ - d₁。Split时d₂ - d₁ 3因为每个高斯分量贡献3个参数权重、均值、方差所以辅助变量需要是3维。但前面矩匹配方案里我用的是(μ₁, μ₂, σ₁², σ₂²)共4个参数对应原分量3个参数此时辅助变量维度是1还是别的这里要仔细建模。其实在这个特定的Split映射里我用的是(u₁, u₂)两个辅助变量。但注意分裂后参数是4个原分量是3个差值是1。维度严格来说差值是1需要用1个辅助变量。但矩匹配的映射并不总是使用最少数量的辅助变量。这就是RJMCMC实现容易出错的地方——在做映射设计时必须核对一下输入(θ, u)的维度与输出θ的维度是否相等。我记得第一次写这个映射的时候心里想的是两个辅助变量看起来更灵活但后来推导雅可比时发现维度不匹配整个马尔可夫链的细致平衡被破坏结果采出来的k分布忽高忽低怀疑人生好几天。4.2 接受率计算与完整迭代逻辑光说不练假把式直接上一段核心代码。这段代码展示的是单步迭代的逻辑包括参数更新、分裂、合并、生灭四种移动以及雅可比行列式的计算。import numpy as np from scipy.stats import norm, invgamma, beta def rjmcmc_step(y, state, move_probs): 执行一步RJMCMC迭代。 参数: y: 观测数据shape (n,) state: 包含当前模型状态的dict键包括: - k: 分量个数 - mu: 均值列表长度k - sigma2: 方差列表长度k - w: 权重列表长度k move_probs: 四种移动的概率[param_update, birth, death, split_merge] 返回: 更新后的state n len(y) move_type np.random.choice([param, birth, death, split_merge], pmove_probs) if move_type param: # 固定模型用MH分别更新每个分量的 (mu, sigma2) for j in range(state[k]): mu_new state[mu][j] 0.1 * state[sigma2][j]**0.5 * np.random.randn() log_ratio compute_log_likelihood(y, state, j, mu_new) - compute_log_likelihood(y, state, j, state[mu][j]) if np.log(np.random.rand()) log_ratio: state[mu][j] mu_new return state if move_type birth: # 生: 从先验中采样新分量参数 u np.random.beta(1, state[k]) w_new u mu_new np.random.uniform(y.min(), y.max()) sigma2_new invgamma.rvs(a2, scalenp.var(y) / 2) # 新状态的构建 state_new state.copy() state_new[k] 1 state_new[w] [w * (1 - u) for w in state[w]] [w_new] state_new[mu] state[mu] [mu_new] state_new[sigma2] state[sigma2] [sigma2_new] # 计算接受率(对数形式) # 后验比(log likelihood ratio log prior ratio) log_lik_ratio compute_log_likelihood(y, state_new) - compute_log_likelihood(y, state) log_prior_ratio compute_log_prior(state_new) - compute_log_prior(state) # 提议比 log_prop_ratio np.log(state[k]) - np.log(state_new[k]) # 雅可比: 从(w_old, u) 到 (w_new, w_old_scaled) 的变换 log_jacobian np.log(1 - u) ** (state[k] - 1) # 简化实际需要考虑完整矩阵 log_alpha log_lik_ratio log_prior_ratio log_prop_ratio log_jacobian if np.log(np.random.rand()) log_alpha: return state_new return state # 其他移动类型类似death、split_merge # ...因为这里的代码是简化版重点展示结构。实际工程实现时我强烈建议把所有对数计算都封装成函数并且每一步都单独记录接受率、雅可比等中间值便于调试。完整代码里log_prior和compute_log_likelihood都要仔细处理。先验部分对k可以设一个Poisson先验截断到某个最大值对每个分量的参数设共轭先验。似然部分要用log-sum-exp技巧计算防止数值下溢。工程上还有几个细节值得注意。步长/提议方差的设定参数更新移动里的提议方差上面代码里我用的0.1 × σ需要自适应调整。目标是让接受率维持在0.2~0.5之间。我通常跑一个几百步的预热阶段用简单的自适应算法在线调整方差。移动概率分配birth、death、split、merge、parameter update这五类移动的概率要合理分配。我一般设置成param update占50%其余四类合计50%。如果k的探索过慢就加大birth/death的比例如果混合度差就加大split/merge的比例。随机数种子RJMCMC的随机性大跑一次可能结果不稳定。复制多个种子运行看k的边缘后验分布是否一致是判断实现正确性的一个有效方法。5. 常见问题与排查技巧实录5.1 接受率过低或过高怎么办RJMCMC最常见的毛病就是接受率要么低得让人绝望要么高到不太正常。这两种情况的处理策略完全不同。接受率过低通常是提议尺度不合适。比如Birth提议的新方差异常大导致新分量覆盖范围太大似然不升反降。解决办法就是让提议分布尽量贴近数据范围或者改用经验分布抽样。还有可能是Split的矩匹配做得不够好导致分裂后整体矩大变接受率骤降。可以先在模拟数据上测试看看split的成功率低于5%就该重新设计映射了。我一般期望的参数更新接受率是20%~50%Split/Merge是5%~20%Birth/Death是1%~10%这个区间仅供参考。接受率过高反而可能是链在局部模型空间里打转。比如你只在某个固定的k附近做小幅度的移动样本全挤在一起接受率自然高但模型空间的探索严重不足。这时候要增加跨空间移动的比例尤其是死亡和合并的频率。另外一个原因是辅助变量的分布选得太集中提议太保守接受率虚高但每一步的影响太小链的有效样本量其实很低。还要检查雅可比行列式是否正确。我踩过的坑是分裂映射的雅可比算成了绝对值忘了它是负值——雅可比行列式可正可负但受率里的变换因子要取绝对值。这个错误会让反向跳跃时的接受率算错导致k的分布产生系统性偏差。检验方法很简单构造两组参数一组正向一组逆向做数值验证看看两个方向的Hastings比乘积是否等于1。5.2 后验分布的多峰性与标签切换高斯混合模型里有个臭名昭著的标签切换问题。如果不对分量做约束采样过程中分量的标签会随机交换后验分布会出现k个对称的模式最终导致均值和方差的边际估计失真。RJMCMC在处理这个问题时格外头疼——因为你还要在变维空间里做标签对齐。我推荐的做法是对参数施加必要的排序约束。最常用的是按均值排序μ₁ μ₂ ... μₖ或者在更新后做排序修正。但这个做法在RJMCMC里要小心排序本身会改变雅可比项吗如果排序函数是分段线性且几乎处处可微雅可比贡献是1所以可以忽略但要确保在边界点不会出问题。另一种做法是做标签的relabeling后处理采样结束后把每次迭代的标签对齐到某个参考模型。这个方案更稳健但实现稍复杂。除了标签问题多峰性还有另一个挑战RJMCMC在模型空间跳跃时可能从一个高概率的盆地跳不到另一个盆地导致后验分布在不同链上收敛到不同模式。我的经验是用并行链策略开多条链从不同初始值出发然后用Gelman-Rubin诊断配合k的边缘分布比较来看是否收敛。必要时使用模拟回火simulated tempering或并行回火parallel tempering来穿越能量壁垒。再分享一个调试技巧在实现完整RJMCMC之前先固定k的取值跑一遍普通MCMC把每个k下的参数后验分布都采出来。一方面可以作为正确性的基准另一方面可以为RJMCMC里的提议分布提供参考。如果RJMCMC在k3附近的参数估计和单独跑3分量的标准MCMC结果差距太大说明你的实现里大概率有bug——这是我验证过无数次的通用检查套路。最后说一个容易忽略的细节结果汇报时k的后验分布往往会受先验影响。如果你给k设了一个比较弱的先验而数据恰好在k2和k3之间模糊不清后验概率可能很接近。这时候不要只看最多的一个k要把整个后验分布画出来。RJMCMC的优势本来就在这里给你的是完整的模型后验分布而不是一个硬生生的模型编号。我在实际项目中碰到过不少情况k2和k3都有接近40%的后验概率这时候只报一个最优模型其实是在抛信息。你应该报告两边的概率和对应的参数后验让用数据的人自己判断。RJMCMC的实现难度确实不低但一旦攻下来了你在贝叶斯计算这块的功力会明显提升。我个人的体会是先把维度差1的简单问题彻底搞清楚再去碰高维的复杂场景每次改动映射函数都要跑模拟数据验证一遍接受率和边际分布。这套流程走顺了以后什么变化点检测、变量选择、稀疏回归都是一个套路换个模型设计提议分布和雅可比而已。