简介本资源是一份面向机器学习进阶学习者与贝叶斯统计实践者的Python开源实现聚焦马尔可夫链蒙特卡洛MCMC算法在贝叶斯推断与后验采样中的核心应用。涵盖Metropolis-Hastings、Gibbs与哈密顿蒙特卡洛HMC三大主流采样器集成贝叶斯线性回归、分布参数估计、多维度收敛诊断如Gelman-Rubin、Geweke、有效样本量及轨迹图、角点图、自相关分析等可视化模块配套6类交互式演示场景支持从原理理解到工程复现的完整学习闭环。压缩包共21个文件以11个核心Python源码含采样器、推断、可视化、示例模块、2个CSV测试数据集、1个PNG结果图及README、requirements等辅助文件为主结构清晰、模块解耦便于按需调用或二次开发整体体积仅644KB轻量易部署。目前已有187人下载学习适合具备Python基础与概率统计知识的学习者系统掌握MCMC实战方法。1. MCMC 不是黑匣子用 Python 手撕 Metropolis-Hastings把后验分布“采”出来而不是“猜”出来你手头有一组观测数据想估计某个参数的真实分布——比如某药物疗效的均值、某设备故障率的置信区间、或者一个复杂模型里十几个超参的联合后验。这时候贝叶斯推断告诉你答案藏在后验分布 $p(\theta \mid D)$ 里。但问题来了这个分布往往没有解析解积分算不出来梯度求不了连画个图都得先“猜个形状”。MCMC 就是那个不靠猜、不靠近似、而是靠“走动”来逼近真相的硬核方法。它不直接生成样本而是构造一条马尔可夫链让链的状态逐步“游荡”到高概率区域并在稳定后留下足够多的落点——这些落点就是后验分布的忠实代表。本文不是讲理论证明而是给你一份可运行、可调试、可改参数、可加模型的 Python 实现从零手写 Metropolis-Hastings最基础也最透明的 MCMC 算法配套真实数据模拟、收敛诊断、轨迹可视化以及最关键的——为什么你的链跑飞了、卡住了、烧 CPU 却不出结果适合正在学贝叶斯统计、做参数估计、或需要在科研/工程中落地不确定量化UQ的 Python 工程师和数据分析师。别再调pymc或stan的黑盒接口了先亲手走一遍采样器的每一步你才能真正看懂 traceplot 里那条线到底在说什么。2. 从原理到代码Metropolis-Hastings 算法的四步拆解与 Python 实现MCMC 的核心不是数学魔术而是一套可编程的决策逻辑当前在哪下一步往哪走走过去值不值得要不要接受这四个问题决定了整个采样链的质量。Metropolis-HastingsMH是最易理解、最易调试的起点。它不要求你掌握复杂的跳跃分布设计只要能写出目标分布即未归一化的后验密度和一个对称建议分布比如高斯就能跑起来。下面我们就把 MH 拆成四步每一步对应一段可执行、可打断、可打印中间状态的 Python 代码。2.1 目标分布后验密度的“能量地图”只写 unnormalized 版本贝叶斯后验 $p(\theta \mid D) \propto p(D \mid \theta) p(\theta)$其中归一化常数 $\int p(D \mid \theta) p(\theta) d\theta$ 通常不可算。幸运的是MH 只需计算两个状态之间的比值常数自动约掉。因此我们只需实现unnormalized log-posterior——用对数形式避免下溢且更稳定。import numpy as np from scipy.stats import norm, uniform def log_posterior(theta, data, prior_mean0.0, prior_std10.0): 未归一化的对数后验密度正态似然 正态先验 theta: 标量参数如均值 data: 观测数据数组一维 返回: floatlog(p(data|theta) * p(theta)) # 似然假设数据 ~ N(theta, 1.0)固定方差为1 log_likelihood np.sum(norm.logpdf(data, loctheta, scale1.0)) # 先验theta ~ N(prior_mean, prior_std) log_prior norm.logpdf(theta, locprior_mean, scaleprior_std) return log_likelihood log_prior注意这里scale1.0是似然的标准差prior_std10.0是先验的宽泛程度。这两个值必须与你的实际问题匹配。如果数据量很小比如 n5先验太强会压倒数据如果数据量大n1000先验影响几乎为零。这不是“调参”而是建模选择——你得清楚自己在假设什么。2.2 建议分布小步试探的“脚”选高斯最稳妥MH 要求一个建议proposal分布 $q(\theta \mid \theta)$用于从当前状态 $\theta$ 生成候选 $\theta$。最常用、最安全的选择是以当前点为中心的高斯分布$q(\theta \mid \theta) \mathcal{N}(\theta, \sigma^2)$。它的对称性$q(\theta \mid \theta) q(\theta \mid \theta)$让接受率公式大幅简化。def propose(theta_current, proposal_std0.5): 从当前theta生成候选theta proposal_std: 建议分布标准差控制步长大小 返回: float候选参数值 return np.random.normal(loctheta_current, scaleproposal_std)关键参数说明proposal_std是 MH 最关键的调优参数。它太小如 0.01链移动缓慢收敛极慢太大如 5.0大部分提议被拒绝链原地踏步。经验法则是让接受率落在 0.2–0.5 区间一维问题可放宽到 0.15–0.6。我们后面会监控它。2.3 接受率计算决定“留”还是“走”的唯一判据MH 的灵魂在于接受率 $\alpha$ $$ \alpha \min\left(1,\ \frac{p(\theta)}{p(\theta)} \cdot \frac{q(\theta \mid \theta)}{q(\theta \mid \theta)}\right) $$ 由于我们用了对称建议分布第二项恒为 1所以只需比较目标分布的比值。用对数形式更稳定def metropolis_hastings_step(theta_current, data, proposal_std0.5): 执行单步 MH生成候选、计算接受率、决定是否接受 返回: (theta_next, accepted: bool, log_ratio: float) # 1. 提议新位置 theta_proposed propose(theta_current, proposal_std) # 2. 计算未归一化后验的对数比值 log_post_current log_posterior(theta_current, data) log_post_proposed log_posterior(theta_proposed, data) log_ratio log_post_proposed - log_post_current # 3. 计算接受概率避免 exp overflow if log_ratio 0.0: alpha 1.0 else: alpha np.exp(log_ratio) # log_ratio 0 exp 1 # 4. 随机决定是否接受 u np.random.uniform() accepted u alpha theta_next theta_proposed if accepted else theta_current return theta_next, accepted, log_ratio逻辑说明这段代码返回三个值不只是theta_next。accepted用于统计接受率log_ratio是诊断链行为的核心——如果它长期远小于 0说明提议总在低概率区得调小proposal_std如果它常接近 0说明提议恰在高概率区边缘接受率可能理想。2.4 主循环烧录 采样 存储三阶段不可跳过一个完整的 MH 运行包含三个阶段burn-in烧录期、thinning薄化、sampling正式采样。新手常犯的错误是直接从第 1 步开始存样或把所有迭代都当有效样本。正确做法如下def run_mh_sampler(data, n_iter10000, burn_in2000, thinning2, init_theta0.0, proposal_std0.5): 运行完整 MH 采样器 n_iter: 总迭代次数含 burn-in burn_in: 烧录期长度丢弃前 burn_in 个样本 thinning: 薄化间隔每 thinning 步存 1 个防自相关 返回: numpy array, shape (n_samples,) # 初始化 theta init_theta samples [] acceptance_count 0 # 主循环 for i in range(n_iter): theta, accepted, _ metropolis_hastings_step( theta, data, proposal_std ) if accepted: acceptance_count 1 # 烧录期不存样 if i burn_in: # 薄化只在 i-burn_in 为 thinning 倍数时存 if (i - burn_in) % thinning 0: samples.append(theta) acceptance_rate acceptance_count / n_iter print(fMH 完成总迭代 {n_iter}接受率 {acceptance_rate:.3f}) return np.array(samples) # 示例模拟数据并运行 np.random.seed(42) true_theta 2.5 data np.random.normal(loctrue_theta, scale1.0, size50) # 50 个观测 samples run_mh_sampler(data, n_iter20000, burn_in5000, thinning5)参数说明burn_in5000前 5000 步用于让链“忘记”初始值进入平稳分布。不能凭感觉设需用收敛诊断见第 4 章验证。thinning5每 5 步存 1 个降低样本间自相关。若链混合好autocorrelation decay 快可设为 1若慢则需加大。init_theta0.0初始值可以任意但离真值太远会延长烧录时间。实践中可用最大似然估计MLE或先验均值初始化。3. 多参数扩展从单变量到多元正态后验向量版 MH 实现现实问题极少只有一个参数。比如线性回归$y_i \beta_0 \beta_1 x_i \varepsilon_i$你要同时估计 $\beta_0, \beta_1, \sigma^2$。此时参数 $\theta$ 是向量目标分布是联合后验 $p(\beta_0, \beta_1, \sigma^2 \mid D)$。MH 的框架不变但实现细节必须升级建议分布要支持向量、目标函数要处理多维输入、接受率计算要保持向量化。3.1 向量目标函数用 numpy 广播避免 for 循环我们以二元正态均值估计为例数据 $D {x_i}_{i1}^n$其中 $x_i \in \mathbb{R}^2$假设 $x_i \sim \mathcal{N}(\mu, \Sigma)$$\Sigma$ 已知如单位阵先验 $\mu \sim \mathcal{N}(\mu_0, \Lambda)$。目标是估计二维均值向量 $\mu [\mu_1, \mu_2]^T$。def log_posterior_2d(mu, data, Sigma_invnp.eye(2), prior_mu0np.array([0.0, 0.0]), prior_Lambda_invnp.eye(2)): 二维均值的未归一化对数后验 mu: shape (2,)当前均值向量 data: shape (n, 2)观测数据矩阵 Sigma_inv: 2x2已知协方差逆矩阵此处为单位阵 prior_Lambda_inv: 2x2先验协方差逆矩阵 返回: float n len(data) # 似然sum_i log N(x_i | mu, Sigma) # -0.5 * sum_i (x_i - mu).T Sigma_inv (x_i - mu) const diff data - mu # (n, 2) # 使用 einsum 避免显式循环对每个样本计算二次型 log_likelihood -0.5 * np.einsum(ij,jk,ik-, diff, Sigma_inv, diff) # 先验log N(mu | prior_mu0, Lambda) diff_prior mu - prior_mu0 log_prior -0.5 * diff_prior prior_Lambda_inv diff_prior return log_likelihood log_prior技巧说明np.einsum(ij,jk,ik-, ...)是计算 $\sum_i (x_i - \mu)^T \Sigma^{-1} (x_i - \mu)$ 的高效方式比for循环快 10 倍以上且完全向量化。Sigma_inv和prior_Lambda_inv用逆矩阵而非原矩阵避免求逆开销。3.2 多元建议分布各维度独立高斯或带相关性的多元高斯最简单做法是各维度独立提议对 $\mu_1$ 和 $\mu_2$ 分别加独立高斯噪声。这等价于用对角协方差矩阵的多元高斯。def propose_2d(mu_current, proposal_stdsnp.array([0.3, 0.3])): 二维独立高斯提议 proposal_stds: shape (2,)每个维度的标准差 返回: shape (2,)候选均值向量 return mu_current np.random.normal(scaleproposal_stds)更高级的做法是使用相关多元高斯其协方差矩阵可学习如用前期样本估计但初学者务必从独立开始——相关结构会显著增加调参难度。3.3 向量版 MH 主循环状态更新与存储逻辑不变但需检查维度def run_mh_2d_sampler(data, n_iter20000, burn_in5000, thinning10, init_munp.array([0.0, 0.0]), proposal_stdsnp.array([0.3, 0.3])): 二维 MH 采样器主函数 返回: samples, shape (n_samples, 2) mu init_mu.copy() samples [] acceptance_count 0 for i in range(n_iter): # 提议 mu_proposed propose_2d(mu, proposal_stds) # 计算对数后验比值 log_post_current log_posterior_2d(mu, data) log_post_proposed log_posterior_2d(mu_proposed, data) log_ratio log_post_proposed - log_post_current # 接受/拒绝 alpha 1.0 if log_ratio 0 else np.exp(log_ratio) if np.random.uniform() alpha: mu mu_proposed acceptance_count 1 # 烧录与薄化 if i burn_in and (i - burn_in) % thinning 0: samples.append(mu.copy()) # 注意 copy否则引用同一数组 acceptance_rate acceptance_count / n_iter print(f2D MH 完成接受率 {acceptance_rate:.3f}) return np.array(samples) # 生成二维数据并运行 np.random.seed(42) true_mu np.array([1.8, -0.7]) data_2d np.random.multivariate_normal(meantrue_mu, covnp.eye(2), size100) samples_2d run_mh_2d_sampler(data_2d, proposal_stdsnp.array([0.4, 0.4]))血泪经验samples.append(mu.copy())中的.copy()绝对不能省否则samples列表里存的全是同一个内存地址的引用最后所有行都等于最后一次迭代的mu。这是 Python 初学者在向量化采样中最常见的翻车点之一。3.4 可视化用 scatter kde 展示联合后验比直方图更有信息量单变量用直方图双变量必须用二维 KDE 或等高线图才能看出参数间的相关性import matplotlib.pyplot as plt import seaborn as sns plt.figure(figsize(10, 4)) # 左图散点图 核密度估计 plt.subplot(1, 2, 1) sns.kdeplot(xsamples_2d[:, 0], ysamples_2d[:, 1], fillTrue, cmapBlues, levels6, alpha0.8) plt.scatter(true_mu[0], true_mu[1], cred, s100, markerx, linewidths3, labelfTrue μ {true_mu}) plt.xlabel(r$\mu_1$) plt.ylabel(r$\mu_2$) plt.title(Joint Posterior KDE) plt.legend() # 右图边际分布 plt.subplot(1, 2, 2) plt.hist(samples_2d[:, 0], bins50, alpha0.6, labelr$\mu_1$, densityTrue) plt.hist(samples_2d[:, 1], bins50, alpha0.6, labelr$\mu_2$, densityTrue) plt.xlabel(Parameter Value) plt.ylabel(Density) plt.title(Marginal Posterior Densities) plt.legend() plt.tight_layout() plt.show()为什么重要左图的椭圆形状直接反映了 $\mu_1$ 和 $\mu_2$ 的后验相关性。如果数据中 $x_1$ 和 $x_2$ 高度相关KDE 椭圆会倾斜如果独立则接近圆形。这是任何单变量摘要如均值、标准差都无法提供的信息。4. 避坑指南MH 采样链的五大典型翻车现场与排查路径MCMC 不是“跑起来就完事”。90% 的失败不在代码语法而在采样链的行为异常。以下是我在线性回归、分层模型、生存分析等多个项目中踩过的坑按现象→原因→解决的结构整理每一条都配可复现的诊断代码。4.1 现象链在某个值附近“死锁”几十步都不动原因proposal_std过大导致提议总落在后验密度极低的区域alpha ≈ 0永远被拒绝。诊断打印每步的log_ratio发现它持续 -10即alpha 4.5e-5。解决立即减小proposal_std例如从 2.0 降到 0.2监控接受率若 0.05大概率是步子太大临时加日志if i % 1000 0: print(fStep {i}: log_ratio {log_ratio:.2f})。4.2 现象链缓慢爬行像蜗牛一样一点点挪动原因proposal_std过小每次只挪一点点自相关极高需要上万步才能获得一个有效样本。诊断计算样本自相关函数ACF——用statsmodels.tsa.stattools.acf发现 lag10 时 ACF 仍 0.8。解决增大proposal_std例如从 0.05 到 0.5目标接受率调至 0.3–0.4一维或 0.2–0.3高维进阶用 adaptive MH在 burn-in 期动态调整proposal_std使接受率趋近目标见第 6 章。4.3 现象链发散theta值爆炸式增长如1e10,-1e12原因目标函数log_posterior在某些区域未定义或返回-inf/nan但 MH 仍尝试提议。常见于先验或似然中用了log函数输入 ≤ 0方差参数 $\sigma^2$ 未加约束应 0却提议负值数值下溢/上溢如exp(1000)。诊断在metropolis_hastings_step中加断言assert not np.isnan(log_post_current) and not np.isinf(log_post_current), \ fNaN/Inf in log_post at theta{theta_current}解决对无约束参数如方差做变换用 $\log(\sigma^2)$ 作为采样变量再映射回 $\sigma^2 \exp(\text{log_sigma2})$在log_posterior开头加保护if np.any(np.isnan(data)) or np.any(np.isinf(data)): return -np.inf用np.clip限制提议范围谨慎使用会引入边界效应。4.4 现象烧录期后链看起来“稳定”但边际分布严重偏斜与真值偏差大原因链未真正收敛只是暂时停留在局部高概率峰。尤其在多峰后验中MH 容易困在次优峰。诊断运行多条独立链不同初值用 Gelman-Rubin $\hat{R}$ 统计量arviz.rhat若 $\hat{R} 1.05$说明链间差异显著未收敛画 traceplotarviz.plot_trace观察多条链是否重叠。解决延长burn_in至 3× 原设定或改用更鲁棒的收敛诊断如 Geweke换用并行 tempering或Hamiltonian Monte CarloHMC需pymc或numpyro最实用用多个初值如init_theta [0, 1, -1, 5]跑多条链取所有链合并后的样本。4.5 现象接受率正常~0.3但有效样本量ESS极低原因高自相关 薄化不足。即使接受率 OK若提议方向与后验主轴不一致链仍会 zigzag 前进。诊断计算 ESSarviz.ess(samples)若ESS / n_samples 0.01说明效率低下。解决立即行动增大thinning如从 1 到 20根本解决用协方差自适应提议——用 burn-in 期样本估计后验协方差 $\hat{\Sigma}$设proposal_cov 2.38^2 / d * \hat{\Sigma}d 为维度替代方案换用pymc的NUTS采样器它自动学习后验几何结构。提示所有诊断都应在burn_in之后进行。切勿用前 1000 步算 ACF 或 ESS——那只是热身不是稳态。5. 收敛诊断实战用 Gelman-Rubin、Geweke 和 Effective Sample Size 三重验证“跑够 10000 步”不是收敛的保证。真正的收敛意味着链已进入平稳分布且样本能无偏估计后验期望。这需要量化验证。本章不讲公式推导只给可粘贴、可运行、可解读的诊断代码覆盖最常用的三种方法。5.1 Gelman-Rubin $\hat{R}$多链一致性检验推荐首选原理比较链内方差 vs 链间方差。若 $\hat{R} \approx 1$说明所有链都探索了同一分布。import arviz as az import numpy as np # 运行 4 条独立链不同初值 chains [] inits [-2.0, 0.0, 2.0, 4.0] for init in inits: chain run_mh_sampler(data, n_iter15000, burn_in5000, thinning5, init_thetainit, proposal_std0.4) chains.append(chain) # 转为 arviz InferenceData 格式必需 idata az.from_numpyro( posterior{theta: np.array(chains).T[:, :, None]}, # shape (n_draws, n_chains, n_vars) observed_data{data: data} ) # 计算 R-hat rhat az.rhat(idata) print(fGelman-Rubin R-hat {rhat[theta].item():.3f}) # 可视化 traceplot az.plot_trace(idata, var_names[theta]) plt.show()解读规则$\hat{R} 1.01$优秀$1.01 \le \hat{R} 1.05$勉强可接受需谨慎$\hat{R} \ge 1.05$未收敛必须重跑加 burn-in、换初值、调 proposal。注意az.from_numpyro是占位名实际需用az.convert_to_inference_data或手动构建字典。关键是chains必须是 list of arrays每个 array 是一条链的样本。5.2 Geweke分段均值差异检验检测早期漂移原理将链切成前 10% 和后 50%检验两段均值是否显著不同。适用于单链快速筛查。from statsmodels.tsa.stattools import acf def geweke_diagnostic(samples, first0.1, last0.5, intervals20): Geweke 检验z-score 形式 first: 前段比例默认 10% last: 后段比例默认 50% 返回: z-score|z| 2 表示未收敛 n len(samples) n_first int(first * n) n_last int(last * n) first_mean np.mean(samples[:n_first]) last_mean np.mean(samples[-n_last:]) # 估算方差考虑自相关 def autocovariance(x, max_lag100): acf_vals acf(x, nlagsmax_lag, fftTrue) return np.sum(acf_vals) * (x.var() / len(x)) var_first autocovariance(samples[:n_first]) var_last autocovariance(samples[-n_last:]) z_score (first_mean - last_mean) / np.sqrt(var_first var_last) return z_score z geweke_diagnostic(samples) print(fGeweke z-score {z:.3f} (|z| 2 indicates non-convergence))实操建议对每条链单独运行 Geweke。若某条链|z| 2说明它没走出烧录期应丢弃该链或延长其burn_in。5.3 Effective Sample Size (ESS)衡量“真正独立样本”数量原理ESS $N / (1 2\sum_{k1}^\infty \rho_k)$其中 $\rho_k$ 是 lag-k 自相关。ESS 越接近总样本数链越高效。def effective_sample_size(samples, max_lag100): 手动计算 ESS简化版忽略 FFT 加速 n len(samples) # 计算自相关函数 acf_vals acf(samples, nlagsmax_lag, fftTrue) # 截断首个负值处停止求和Bartletts formula tau 1 2 * np.sum(acf_vals[1:]) ess n / tau return int(ess) ess effective_sample_size(samples) print(fTotal samples: {len(samples)}, ESS: {ess}, ESS ratio: {ess/len(samples):.3f})阈值参考ESS ≥ 100可做粗略推断ESS ≥ 1000可做可信区间估计ESS / n_samples 0.05链效率极低必须优化提议或换算法。5.4 一张表看清诊断结果与行动指南诊断方法正常值范围异常表现应对动作Gelman-Rubin$\hat{R} 1.01$$\hat{R} 1.05$增加 burn-in换初值重跑多链Geweke z-score$|z| 2$$|z| 2.5$单独延长该链 burn-inESS ratio$ 0.1$一维$ 0.02$增大 thinning调 proposal_std用自适应 MH接受率0.2–0.5一维 0.05 或 0.9减小/增大 proposal_stdTraceplot 形态平稳波动无趋势持续上升/下降或长平台检查目标函数数值稳定性注意这五项不是“任选其一”而是必须全部通过。就像体检不能只查血压——R-hat 合格但 ESS 极低说明链虽收敛但效率为零ESS 高但 Geweke 异常说明链后期才稳定前期样本全废。6. 进阶技巧自适应 Metropolis-Hastings 与后验预测检查PPC当你已经能稳定跑出收敛链下一步是让采样器更智能、结果更可信。本章聚焦两个实战价值最高的进阶点自适应提议让 MH 自己学会怎么走和后验预测检查用数据本身验证模型是否合理。它们不增加理论负担但极大提升结果可信度。6.1 自适应 MHburn-in 期动态学习提议协方差标准 MH 的proposal_std是固定超参需人工试错。自适应 MHAM在 burn-in 期边跑边学用已生成的样本估计后验协方差 $\hat{\Sigma}$并据此缩放提议尺度。Roberts et al. (1997) 证明最优尺度为 $2.38^2 / d$d 为维度。def run_adaptive_mh(data, n_iter20000, burn_in5000, thinning5, init_theta0.0, adapt_start1000, adapt_interval100): 自适应 MHburn-in 期每 adapt_interval 步更新一次提议尺度 adapt_start: 从第几步开始适应避开初始不稳定期 theta init_theta samples [] acceptance_count 0 # 存储 burn-in 期样本用于协方差估计 burn_in_samples [] # 初始提议标准差保守起见设小些 proposal_std 0.1 for i in range(n_iter): # 提议 接受 theta_proposed np.random.normal(loctheta, scaleproposal_std) log_post_current log_posterior(theta, data) log_post_proposed log_posterior(theta_proposed, data) log_ratio log_post_proposed - log_post_current alpha 1.0 if log_ratio 0 else np.exp(log_ratio) if np.random.uniform() alpha: theta theta_proposed acceptance_count 1 # 烧录期存储 if i burn_in: burn_in_samples.append(theta) else: if (i - burn_in) % thinning 0: samples.append(theta) # 自适应仅在 burn-in 期内且达到 adapt_start 后 if i burn_in and i adapt_start and i % adapt_interval 0: # 用 burn_in_samples 估计标准差一维 std_est np.std(burn_in_samples[-adapt_interval:]) 1e-6 # Roberts 公式最优尺度 2.38 / sqrt(d)d1 2.38 proposal_std 2.38 * std_est # print(fStep {i}: updated proposal_std {proposal_std:.3f}) acceptance_rate acceptance_count / n_iter print(fAdaptive MH 完成接受率 {acceptance_rate:.3f}, final proposal_std {proposal_std:.3f}) return np.array(samples) # 运行对比 samples_adapt run_adaptive_mh(data, adapt_start2000)优势无需反复试proposal_std对多峰、非球形后验更鲁棒接受率自动趋近理论最优值。注意自适应只能在 burn-in 期进行正式采样期必须冻结提议分布否则破坏马尔可夫性。6.2 后验预测检查PPC本文还有配套的精品资源点击获取
手撕 Metropolis-Hastings：从原理到可调试 Python 实现
简介本资源是一份面向机器学习进阶学习者与贝叶斯统计实践者的Python开源实现聚焦马尔可夫链蒙特卡洛MCMC算法在贝叶斯推断与后验采样中的核心应用。涵盖Metropolis-Hastings、Gibbs与哈密顿蒙特卡洛HMC三大主流采样器集成贝叶斯线性回归、分布参数估计、多维度收敛诊断如Gelman-Rubin、Geweke、有效样本量及轨迹图、角点图、自相关分析等可视化模块配套6类交互式演示场景支持从原理理解到工程复现的完整学习闭环。压缩包共21个文件以11个核心Python源码含采样器、推断、可视化、示例模块、2个CSV测试数据集、1个PNG结果图及README、requirements等辅助文件为主结构清晰、模块解耦便于按需调用或二次开发整体体积仅644KB轻量易部署。目前已有187人下载学习适合具备Python基础与概率统计知识的学习者系统掌握MCMC实战方法。1. MCMC 不是黑匣子用 Python 手撕 Metropolis-Hastings把后验分布“采”出来而不是“猜”出来你手头有一组观测数据想估计某个参数的真实分布——比如某药物疗效的均值、某设备故障率的置信区间、或者一个复杂模型里十几个超参的联合后验。这时候贝叶斯推断告诉你答案藏在后验分布 $p(\theta \mid D)$ 里。但问题来了这个分布往往没有解析解积分算不出来梯度求不了连画个图都得先“猜个形状”。MCMC 就是那个不靠猜、不靠近似、而是靠“走动”来逼近真相的硬核方法。它不直接生成样本而是构造一条马尔可夫链让链的状态逐步“游荡”到高概率区域并在稳定后留下足够多的落点——这些落点就是后验分布的忠实代表。本文不是讲理论证明而是给你一份可运行、可调试、可改参数、可加模型的 Python 实现从零手写 Metropolis-Hastings最基础也最透明的 MCMC 算法配套真实数据模拟、收敛诊断、轨迹可视化以及最关键的——为什么你的链跑飞了、卡住了、烧 CPU 却不出结果适合正在学贝叶斯统计、做参数估计、或需要在科研/工程中落地不确定量化UQ的 Python 工程师和数据分析师。别再调pymc或stan的黑盒接口了先亲手走一遍采样器的每一步你才能真正看懂 traceplot 里那条线到底在说什么。2. 从原理到代码Metropolis-Hastings 算法的四步拆解与 Python 实现MCMC 的核心不是数学魔术而是一套可编程的决策逻辑当前在哪下一步往哪走走过去值不值得要不要接受这四个问题决定了整个采样链的质量。Metropolis-HastingsMH是最易理解、最易调试的起点。它不要求你掌握复杂的跳跃分布设计只要能写出目标分布即未归一化的后验密度和一个对称建议分布比如高斯就能跑起来。下面我们就把 MH 拆成四步每一步对应一段可执行、可打断、可打印中间状态的 Python 代码。2.1 目标分布后验密度的“能量地图”只写 unnormalized 版本贝叶斯后验 $p(\theta \mid D) \propto p(D \mid \theta) p(\theta)$其中归一化常数 $\int p(D \mid \theta) p(\theta) d\theta$ 通常不可算。幸运的是MH 只需计算两个状态之间的比值常数自动约掉。因此我们只需实现unnormalized log-posterior——用对数形式避免下溢且更稳定。import numpy as np from scipy.stats import norm, uniform def log_posterior(theta, data, prior_mean0.0, prior_std10.0): 未归一化的对数后验密度正态似然 正态先验 theta: 标量参数如均值 data: 观测数据数组一维 返回: floatlog(p(data|theta) * p(theta)) # 似然假设数据 ~ N(theta, 1.0)固定方差为1 log_likelihood np.sum(norm.logpdf(data, loctheta, scale1.0)) # 先验theta ~ N(prior_mean, prior_std) log_prior norm.logpdf(theta, locprior_mean, scaleprior_std) return log_likelihood log_prior注意这里scale1.0是似然的标准差prior_std10.0是先验的宽泛程度。这两个值必须与你的实际问题匹配。如果数据量很小比如 n5先验太强会压倒数据如果数据量大n1000先验影响几乎为零。这不是“调参”而是建模选择——你得清楚自己在假设什么。2.2 建议分布小步试探的“脚”选高斯最稳妥MH 要求一个建议proposal分布 $q(\theta \mid \theta)$用于从当前状态 $\theta$ 生成候选 $\theta$。最常用、最安全的选择是以当前点为中心的高斯分布$q(\theta \mid \theta) \mathcal{N}(\theta, \sigma^2)$。它的对称性$q(\theta \mid \theta) q(\theta \mid \theta)$让接受率公式大幅简化。def propose(theta_current, proposal_std0.5): 从当前theta生成候选theta proposal_std: 建议分布标准差控制步长大小 返回: float候选参数值 return np.random.normal(loctheta_current, scaleproposal_std)关键参数说明proposal_std是 MH 最关键的调优参数。它太小如 0.01链移动缓慢收敛极慢太大如 5.0大部分提议被拒绝链原地踏步。经验法则是让接受率落在 0.2–0.5 区间一维问题可放宽到 0.15–0.6。我们后面会监控它。2.3 接受率计算决定“留”还是“走”的唯一判据MH 的灵魂在于接受率 $\alpha$ $$ \alpha \min\left(1,\ \frac{p(\theta)}{p(\theta)} \cdot \frac{q(\theta \mid \theta)}{q(\theta \mid \theta)}\right) $$ 由于我们用了对称建议分布第二项恒为 1所以只需比较目标分布的比值。用对数形式更稳定def metropolis_hastings_step(theta_current, data, proposal_std0.5): 执行单步 MH生成候选、计算接受率、决定是否接受 返回: (theta_next, accepted: bool, log_ratio: float) # 1. 提议新位置 theta_proposed propose(theta_current, proposal_std) # 2. 计算未归一化后验的对数比值 log_post_current log_posterior(theta_current, data) log_post_proposed log_posterior(theta_proposed, data) log_ratio log_post_proposed - log_post_current # 3. 计算接受概率避免 exp overflow if log_ratio 0.0: alpha 1.0 else: alpha np.exp(log_ratio) # log_ratio 0 exp 1 # 4. 随机决定是否接受 u np.random.uniform() accepted u alpha theta_next theta_proposed if accepted else theta_current return theta_next, accepted, log_ratio逻辑说明这段代码返回三个值不只是theta_next。accepted用于统计接受率log_ratio是诊断链行为的核心——如果它长期远小于 0说明提议总在低概率区得调小proposal_std如果它常接近 0说明提议恰在高概率区边缘接受率可能理想。2.4 主循环烧录 采样 存储三阶段不可跳过一个完整的 MH 运行包含三个阶段burn-in烧录期、thinning薄化、sampling正式采样。新手常犯的错误是直接从第 1 步开始存样或把所有迭代都当有效样本。正确做法如下def run_mh_sampler(data, n_iter10000, burn_in2000, thinning2, init_theta0.0, proposal_std0.5): 运行完整 MH 采样器 n_iter: 总迭代次数含 burn-in burn_in: 烧录期长度丢弃前 burn_in 个样本 thinning: 薄化间隔每 thinning 步存 1 个防自相关 返回: numpy array, shape (n_samples,) # 初始化 theta init_theta samples [] acceptance_count 0 # 主循环 for i in range(n_iter): theta, accepted, _ metropolis_hastings_step( theta, data, proposal_std ) if accepted: acceptance_count 1 # 烧录期不存样 if i burn_in: # 薄化只在 i-burn_in 为 thinning 倍数时存 if (i - burn_in) % thinning 0: samples.append(theta) acceptance_rate acceptance_count / n_iter print(fMH 完成总迭代 {n_iter}接受率 {acceptance_rate:.3f}) return np.array(samples) # 示例模拟数据并运行 np.random.seed(42) true_theta 2.5 data np.random.normal(loctrue_theta, scale1.0, size50) # 50 个观测 samples run_mh_sampler(data, n_iter20000, burn_in5000, thinning5)参数说明burn_in5000前 5000 步用于让链“忘记”初始值进入平稳分布。不能凭感觉设需用收敛诊断见第 4 章验证。thinning5每 5 步存 1 个降低样本间自相关。若链混合好autocorrelation decay 快可设为 1若慢则需加大。init_theta0.0初始值可以任意但离真值太远会延长烧录时间。实践中可用最大似然估计MLE或先验均值初始化。3. 多参数扩展从单变量到多元正态后验向量版 MH 实现现实问题极少只有一个参数。比如线性回归$y_i \beta_0 \beta_1 x_i \varepsilon_i$你要同时估计 $\beta_0, \beta_1, \sigma^2$。此时参数 $\theta$ 是向量目标分布是联合后验 $p(\beta_0, \beta_1, \sigma^2 \mid D)$。MH 的框架不变但实现细节必须升级建议分布要支持向量、目标函数要处理多维输入、接受率计算要保持向量化。3.1 向量目标函数用 numpy 广播避免 for 循环我们以二元正态均值估计为例数据 $D {x_i}_{i1}^n$其中 $x_i \in \mathbb{R}^2$假设 $x_i \sim \mathcal{N}(\mu, \Sigma)$$\Sigma$ 已知如单位阵先验 $\mu \sim \mathcal{N}(\mu_0, \Lambda)$。目标是估计二维均值向量 $\mu [\mu_1, \mu_2]^T$。def log_posterior_2d(mu, data, Sigma_invnp.eye(2), prior_mu0np.array([0.0, 0.0]), prior_Lambda_invnp.eye(2)): 二维均值的未归一化对数后验 mu: shape (2,)当前均值向量 data: shape (n, 2)观测数据矩阵 Sigma_inv: 2x2已知协方差逆矩阵此处为单位阵 prior_Lambda_inv: 2x2先验协方差逆矩阵 返回: float n len(data) # 似然sum_i log N(x_i | mu, Sigma) # -0.5 * sum_i (x_i - mu).T Sigma_inv (x_i - mu) const diff data - mu # (n, 2) # 使用 einsum 避免显式循环对每个样本计算二次型 log_likelihood -0.5 * np.einsum(ij,jk,ik-, diff, Sigma_inv, diff) # 先验log N(mu | prior_mu0, Lambda) diff_prior mu - prior_mu0 log_prior -0.5 * diff_prior prior_Lambda_inv diff_prior return log_likelihood log_prior技巧说明np.einsum(ij,jk,ik-, ...)是计算 $\sum_i (x_i - \mu)^T \Sigma^{-1} (x_i - \mu)$ 的高效方式比for循环快 10 倍以上且完全向量化。Sigma_inv和prior_Lambda_inv用逆矩阵而非原矩阵避免求逆开销。3.2 多元建议分布各维度独立高斯或带相关性的多元高斯最简单做法是各维度独立提议对 $\mu_1$ 和 $\mu_2$ 分别加独立高斯噪声。这等价于用对角协方差矩阵的多元高斯。def propose_2d(mu_current, proposal_stdsnp.array([0.3, 0.3])): 二维独立高斯提议 proposal_stds: shape (2,)每个维度的标准差 返回: shape (2,)候选均值向量 return mu_current np.random.normal(scaleproposal_stds)更高级的做法是使用相关多元高斯其协方差矩阵可学习如用前期样本估计但初学者务必从独立开始——相关结构会显著增加调参难度。3.3 向量版 MH 主循环状态更新与存储逻辑不变但需检查维度def run_mh_2d_sampler(data, n_iter20000, burn_in5000, thinning10, init_munp.array([0.0, 0.0]), proposal_stdsnp.array([0.3, 0.3])): 二维 MH 采样器主函数 返回: samples, shape (n_samples, 2) mu init_mu.copy() samples [] acceptance_count 0 for i in range(n_iter): # 提议 mu_proposed propose_2d(mu, proposal_stds) # 计算对数后验比值 log_post_current log_posterior_2d(mu, data) log_post_proposed log_posterior_2d(mu_proposed, data) log_ratio log_post_proposed - log_post_current # 接受/拒绝 alpha 1.0 if log_ratio 0 else np.exp(log_ratio) if np.random.uniform() alpha: mu mu_proposed acceptance_count 1 # 烧录与薄化 if i burn_in and (i - burn_in) % thinning 0: samples.append(mu.copy()) # 注意 copy否则引用同一数组 acceptance_rate acceptance_count / n_iter print(f2D MH 完成接受率 {acceptance_rate:.3f}) return np.array(samples) # 生成二维数据并运行 np.random.seed(42) true_mu np.array([1.8, -0.7]) data_2d np.random.multivariate_normal(meantrue_mu, covnp.eye(2), size100) samples_2d run_mh_2d_sampler(data_2d, proposal_stdsnp.array([0.4, 0.4]))血泪经验samples.append(mu.copy())中的.copy()绝对不能省否则samples列表里存的全是同一个内存地址的引用最后所有行都等于最后一次迭代的mu。这是 Python 初学者在向量化采样中最常见的翻车点之一。3.4 可视化用 scatter kde 展示联合后验比直方图更有信息量单变量用直方图双变量必须用二维 KDE 或等高线图才能看出参数间的相关性import matplotlib.pyplot as plt import seaborn as sns plt.figure(figsize(10, 4)) # 左图散点图 核密度估计 plt.subplot(1, 2, 1) sns.kdeplot(xsamples_2d[:, 0], ysamples_2d[:, 1], fillTrue, cmapBlues, levels6, alpha0.8) plt.scatter(true_mu[0], true_mu[1], cred, s100, markerx, linewidths3, labelfTrue μ {true_mu}) plt.xlabel(r$\mu_1$) plt.ylabel(r$\mu_2$) plt.title(Joint Posterior KDE) plt.legend() # 右图边际分布 plt.subplot(1, 2, 2) plt.hist(samples_2d[:, 0], bins50, alpha0.6, labelr$\mu_1$, densityTrue) plt.hist(samples_2d[:, 1], bins50, alpha0.6, labelr$\mu_2$, densityTrue) plt.xlabel(Parameter Value) plt.ylabel(Density) plt.title(Marginal Posterior Densities) plt.legend() plt.tight_layout() plt.show()为什么重要左图的椭圆形状直接反映了 $\mu_1$ 和 $\mu_2$ 的后验相关性。如果数据中 $x_1$ 和 $x_2$ 高度相关KDE 椭圆会倾斜如果独立则接近圆形。这是任何单变量摘要如均值、标准差都无法提供的信息。4. 避坑指南MH 采样链的五大典型翻车现场与排查路径MCMC 不是“跑起来就完事”。90% 的失败不在代码语法而在采样链的行为异常。以下是我在线性回归、分层模型、生存分析等多个项目中踩过的坑按现象→原因→解决的结构整理每一条都配可复现的诊断代码。4.1 现象链在某个值附近“死锁”几十步都不动原因proposal_std过大导致提议总落在后验密度极低的区域alpha ≈ 0永远被拒绝。诊断打印每步的log_ratio发现它持续 -10即alpha 4.5e-5。解决立即减小proposal_std例如从 2.0 降到 0.2监控接受率若 0.05大概率是步子太大临时加日志if i % 1000 0: print(fStep {i}: log_ratio {log_ratio:.2f})。4.2 现象链缓慢爬行像蜗牛一样一点点挪动原因proposal_std过小每次只挪一点点自相关极高需要上万步才能获得一个有效样本。诊断计算样本自相关函数ACF——用statsmodels.tsa.stattools.acf发现 lag10 时 ACF 仍 0.8。解决增大proposal_std例如从 0.05 到 0.5目标接受率调至 0.3–0.4一维或 0.2–0.3高维进阶用 adaptive MH在 burn-in 期动态调整proposal_std使接受率趋近目标见第 6 章。4.3 现象链发散theta值爆炸式增长如1e10,-1e12原因目标函数log_posterior在某些区域未定义或返回-inf/nan但 MH 仍尝试提议。常见于先验或似然中用了log函数输入 ≤ 0方差参数 $\sigma^2$ 未加约束应 0却提议负值数值下溢/上溢如exp(1000)。诊断在metropolis_hastings_step中加断言assert not np.isnan(log_post_current) and not np.isinf(log_post_current), \ fNaN/Inf in log_post at theta{theta_current}解决对无约束参数如方差做变换用 $\log(\sigma^2)$ 作为采样变量再映射回 $\sigma^2 \exp(\text{log_sigma2})$在log_posterior开头加保护if np.any(np.isnan(data)) or np.any(np.isinf(data)): return -np.inf用np.clip限制提议范围谨慎使用会引入边界效应。4.4 现象烧录期后链看起来“稳定”但边际分布严重偏斜与真值偏差大原因链未真正收敛只是暂时停留在局部高概率峰。尤其在多峰后验中MH 容易困在次优峰。诊断运行多条独立链不同初值用 Gelman-Rubin $\hat{R}$ 统计量arviz.rhat若 $\hat{R} 1.05$说明链间差异显著未收敛画 traceplotarviz.plot_trace观察多条链是否重叠。解决延长burn_in至 3× 原设定或改用更鲁棒的收敛诊断如 Geweke换用并行 tempering或Hamiltonian Monte CarloHMC需pymc或numpyro最实用用多个初值如init_theta [0, 1, -1, 5]跑多条链取所有链合并后的样本。4.5 现象接受率正常~0.3但有效样本量ESS极低原因高自相关 薄化不足。即使接受率 OK若提议方向与后验主轴不一致链仍会 zigzag 前进。诊断计算 ESSarviz.ess(samples)若ESS / n_samples 0.01说明效率低下。解决立即行动增大thinning如从 1 到 20根本解决用协方差自适应提议——用 burn-in 期样本估计后验协方差 $\hat{\Sigma}$设proposal_cov 2.38^2 / d * \hat{\Sigma}d 为维度替代方案换用pymc的NUTS采样器它自动学习后验几何结构。提示所有诊断都应在burn_in之后进行。切勿用前 1000 步算 ACF 或 ESS——那只是热身不是稳态。5. 收敛诊断实战用 Gelman-Rubin、Geweke 和 Effective Sample Size 三重验证“跑够 10000 步”不是收敛的保证。真正的收敛意味着链已进入平稳分布且样本能无偏估计后验期望。这需要量化验证。本章不讲公式推导只给可粘贴、可运行、可解读的诊断代码覆盖最常用的三种方法。5.1 Gelman-Rubin $\hat{R}$多链一致性检验推荐首选原理比较链内方差 vs 链间方差。若 $\hat{R} \approx 1$说明所有链都探索了同一分布。import arviz as az import numpy as np # 运行 4 条独立链不同初值 chains [] inits [-2.0, 0.0, 2.0, 4.0] for init in inits: chain run_mh_sampler(data, n_iter15000, burn_in5000, thinning5, init_thetainit, proposal_std0.4) chains.append(chain) # 转为 arviz InferenceData 格式必需 idata az.from_numpyro( posterior{theta: np.array(chains).T[:, :, None]}, # shape (n_draws, n_chains, n_vars) observed_data{data: data} ) # 计算 R-hat rhat az.rhat(idata) print(fGelman-Rubin R-hat {rhat[theta].item():.3f}) # 可视化 traceplot az.plot_trace(idata, var_names[theta]) plt.show()解读规则$\hat{R} 1.01$优秀$1.01 \le \hat{R} 1.05$勉强可接受需谨慎$\hat{R} \ge 1.05$未收敛必须重跑加 burn-in、换初值、调 proposal。注意az.from_numpyro是占位名实际需用az.convert_to_inference_data或手动构建字典。关键是chains必须是 list of arrays每个 array 是一条链的样本。5.2 Geweke分段均值差异检验检测早期漂移原理将链切成前 10% 和后 50%检验两段均值是否显著不同。适用于单链快速筛查。from statsmodels.tsa.stattools import acf def geweke_diagnostic(samples, first0.1, last0.5, intervals20): Geweke 检验z-score 形式 first: 前段比例默认 10% last: 后段比例默认 50% 返回: z-score|z| 2 表示未收敛 n len(samples) n_first int(first * n) n_last int(last * n) first_mean np.mean(samples[:n_first]) last_mean np.mean(samples[-n_last:]) # 估算方差考虑自相关 def autocovariance(x, max_lag100): acf_vals acf(x, nlagsmax_lag, fftTrue) return np.sum(acf_vals) * (x.var() / len(x)) var_first autocovariance(samples[:n_first]) var_last autocovariance(samples[-n_last:]) z_score (first_mean - last_mean) / np.sqrt(var_first var_last) return z_score z geweke_diagnostic(samples) print(fGeweke z-score {z:.3f} (|z| 2 indicates non-convergence))实操建议对每条链单独运行 Geweke。若某条链|z| 2说明它没走出烧录期应丢弃该链或延长其burn_in。5.3 Effective Sample Size (ESS)衡量“真正独立样本”数量原理ESS $N / (1 2\sum_{k1}^\infty \rho_k)$其中 $\rho_k$ 是 lag-k 自相关。ESS 越接近总样本数链越高效。def effective_sample_size(samples, max_lag100): 手动计算 ESS简化版忽略 FFT 加速 n len(samples) # 计算自相关函数 acf_vals acf(samples, nlagsmax_lag, fftTrue) # 截断首个负值处停止求和Bartletts formula tau 1 2 * np.sum(acf_vals[1:]) ess n / tau return int(ess) ess effective_sample_size(samples) print(fTotal samples: {len(samples)}, ESS: {ess}, ESS ratio: {ess/len(samples):.3f})阈值参考ESS ≥ 100可做粗略推断ESS ≥ 1000可做可信区间估计ESS / n_samples 0.05链效率极低必须优化提议或换算法。5.4 一张表看清诊断结果与行动指南诊断方法正常值范围异常表现应对动作Gelman-Rubin$\hat{R} 1.01$$\hat{R} 1.05$增加 burn-in换初值重跑多链Geweke z-score$|z| 2$$|z| 2.5$单独延长该链 burn-inESS ratio$ 0.1$一维$ 0.02$增大 thinning调 proposal_std用自适应 MH接受率0.2–0.5一维 0.05 或 0.9减小/增大 proposal_stdTraceplot 形态平稳波动无趋势持续上升/下降或长平台检查目标函数数值稳定性注意这五项不是“任选其一”而是必须全部通过。就像体检不能只查血压——R-hat 合格但 ESS 极低说明链虽收敛但效率为零ESS 高但 Geweke 异常说明链后期才稳定前期样本全废。6. 进阶技巧自适应 Metropolis-Hastings 与后验预测检查PPC当你已经能稳定跑出收敛链下一步是让采样器更智能、结果更可信。本章聚焦两个实战价值最高的进阶点自适应提议让 MH 自己学会怎么走和后验预测检查用数据本身验证模型是否合理。它们不增加理论负担但极大提升结果可信度。6.1 自适应 MHburn-in 期动态学习提议协方差标准 MH 的proposal_std是固定超参需人工试错。自适应 MHAM在 burn-in 期边跑边学用已生成的样本估计后验协方差 $\hat{\Sigma}$并据此缩放提议尺度。Roberts et al. (1997) 证明最优尺度为 $2.38^2 / d$d 为维度。def run_adaptive_mh(data, n_iter20000, burn_in5000, thinning5, init_theta0.0, adapt_start1000, adapt_interval100): 自适应 MHburn-in 期每 adapt_interval 步更新一次提议尺度 adapt_start: 从第几步开始适应避开初始不稳定期 theta init_theta samples [] acceptance_count 0 # 存储 burn-in 期样本用于协方差估计 burn_in_samples [] # 初始提议标准差保守起见设小些 proposal_std 0.1 for i in range(n_iter): # 提议 接受 theta_proposed np.random.normal(loctheta, scaleproposal_std) log_post_current log_posterior(theta, data) log_post_proposed log_posterior(theta_proposed, data) log_ratio log_post_proposed - log_post_current alpha 1.0 if log_ratio 0 else np.exp(log_ratio) if np.random.uniform() alpha: theta theta_proposed acceptance_count 1 # 烧录期存储 if i burn_in: burn_in_samples.append(theta) else: if (i - burn_in) % thinning 0: samples.append(theta) # 自适应仅在 burn-in 期内且达到 adapt_start 后 if i burn_in and i adapt_start and i % adapt_interval 0: # 用 burn_in_samples 估计标准差一维 std_est np.std(burn_in_samples[-adapt_interval:]) 1e-6 # Roberts 公式最优尺度 2.38 / sqrt(d)d1 2.38 proposal_std 2.38 * std_est # print(fStep {i}: updated proposal_std {proposal_std:.3f}) acceptance_rate acceptance_count / n_iter print(fAdaptive MH 完成接受率 {acceptance_rate:.3f}, final proposal_std {proposal_std:.3f}) return np.array(samples) # 运行对比 samples_adapt run_adaptive_mh(data, adapt_start2000)优势无需反复试proposal_std对多峰、非球形后验更鲁棒接受率自动趋近理论最优值。注意自适应只能在 burn-in 期进行正式采样期必须冻结提议分布否则破坏马尔可夫性。6.2 后验预测检查PPC本文还有配套的精品资源点击获取
