alterlab-pymc

PyMC 贝叶斯建模

概述

PyMC 是一个用于贝叶斯建模和概率编程的 Python 库。使用当前 API（PyMC ≥ 6，ArviZ ≥ 1.1）构建、拟合、验证和比较贝叶斯模型，包括层级模型、MCMC 抽样（NUTS）、变分推断和 PSIS-LOO 模型比较。

uv pip install "pymc[nutpie]" "arviz[matplotlib,h5netcdf]"   # nutpie sampler, plotting backend, NetCDF I/O

PyMC 6 更改了许多默认设置和名称，而大多数教程仍在使用旧版内容。在复用 PyMC 5 / ArviZ 0.x 代码前，请阅读文末的 版本说明。

何时使用此技能

以下情况应使用此技能：

• 构建贝叶斯模型（线性 / 逻辑回归、层级模型、时间序列等）
• 执行 MCMC 抽样或变分推断
• 进行先验 / 后验预测检验
• 诊断抽样问题（发散、收敛、ESS）
• 使用 PSIS-LOO 交叉验证比较多个模型
• 通过贝叶斯方法实现不确定性量化
• 处理层级 / 多层数据结构
• 以合理严谨的方式处理缺失数据或测量误差

不适用场景

场景	改用
使用 p 值、标准误和残差诊断的频率学派回归	alterlab-statsmodels
针对嵌套数据、以报告为重点的混合效应模型（lme4/glmmTMB/brms、statsmodels MixedLM 或 bambi 公式），并报告 ICC 和方差分量	alterlab-multilevel-models
汇总多项研究中的效应量（固定 / 随机效应、异质性、漏斗图）	alterlab-meta-analysis
点预测机器学习（随机森林、梯度提升、交叉验证调参）	alterlab-scikit-learn

标准贝叶斯工作流

构建和验证贝叶斯模型时，请遵循以下 8 步工作流：

1. 数据准备 — 标准化预测变量、处理缺失数据、设置 coords
2. 构建模型 — 使用弱信息先验、命名的 dims，并为可预测变量使用 pm.Data()
3. 先验预测检验 — 使用 pm.sample_prior_predictive；拟合前验证先验
4. 拟合 — 使用 pm.sample(draws=2000, tune=1000, chains=4, target_accept=0.9)，如果要比较模型，接着运行 pm.compute_log_likelihood(idata)
5. 诊断 — R-hat < 1.01，ESS > 400，无发散，迹线混合良好
6. 后验预测检验 — 使用 pm.sample_posterior_predictive；检查拟合结果与观测数据是否相符
7. 分析 — 使用 az.summary、az.plot_dist（曾为 plot_posterior）、az.plot_forest
8. 预测 — 使用 pm.set_data，然后运行 pm.sample_posterior_predictive；提取 HDI 区间

完整的分步代码：references/workflow_examples.md。

常见模型模式

PyMC 支持线性 / 逻辑 / 泊松回归、层级（多层）模型和时间序列（AR）。 可直接调整使用的代码见 references/model_patterns.md。

当组内观测数量较少时，层级模型应使用非中心化参数化：中心化形式会产生 NUTS 难以探索的漏斗几何结构，从而导致发散。 模板：assets/linear_regression_template.py、assets/hierarchical_model_template.py。

分布选择

选择先验和似然是建模中影响最大的决策。快速选择指南 （尺度参数、无界、正值、概率、相关矩阵；连续 / 计数 / 二元 / 分类似然）见 references/distribution_selection.md。完整目录 见 references/distributions.md。

模型比较与采样

• 比较模型时使用 PSIS-LOO（scripts/model_comparison.py）；结合 dse 解读 elpd_diff，并检查 Pareto-k。ArviZ 1.x 不支持 WAIC。
• 默认使用 NUTS 进行采样（安装了 nutpie 时则使用 nutpie）；出现发散时提高 target_accept；快速近似可使用 ADVI。
• 使用 scripts/model_diagnostics.py 进行诊断（check_diagnostics、create_diagnostic_report）。
• 排查发散、ESS 偏低、R-hat 偏高和采样缓慢等问题。

代码和决策规则：references/model_comparison.md。采样算法详细指南： references/sampling_inference.md。

最佳实践

构建模型

1. 始终对预测变量进行标准化，以改善采样效果
2. 使用弱信息先验（而非平坦先验）
3. 使用命名维度（dims），以提高清晰度
4. 对分层模型使用非中心化参数化
5. 拟合前检查先验预测

采样

1. 运行多个链（至少 4 条），以确保收敛
2. 将 target_accept=0.9 作为基准（必要时调高）
3. 比较模型前调用 pm.compute_log_likelihood(idata)
4. 设置随机种子，以确保可复现

验证

1. 解读结果前先检查诊断指标（R-hat、ESS、发散情况）
2. 使用后验预测检查验证模型
3. 适当时比较多个模型
4. 报告不确定性（HDI 区间，而不只是点估计）

从简单模型开始，逐步增加复杂度，并根据每次预测检查的结果迭代调整模型（参见上文的 8 步工作流）。

资源

此技能包含：

参考资料（references/）

• workflow_examples.md：8 步贝叶斯工作流（数据准备 → 预测）的完整分步代码。
• model_patterns.md：可直接调整使用的线性／逻辑／Poisson 回归、分层模型和 AR 时间序列模型代码。
• distribution_selection.md：按参数／结果类型快速选择先验和似然函数。
• model_comparison.md：PSIS-LOO 比较、诊断脚本和故障排查（发散、ESS、R-hat、采样缓慢）。
• distributions.md：按类别整理的 PyMC 分布全面目录（连续、离散、多变量、混合、时间序列）。选择先验或似然函数时可参考。
• sampling_inference.md：采样算法（NUTS、Metropolis、SMC）、变分推断（ADVI、SVGD）和采样问题处理的详细指南。遇到收敛问题或选择推断方法时可参考。
• workflows.md：单个可运行的端到端脚本（数据准备 → 保存），以及 workflow_examples.md 中未包含的额外实用示例，包括缺失数据插补、QR 重新参数化、混合模型和模型平均辅助函数。

脚本（scripts/）

• model_diagnostics.py：自动检查诊断指标并生成报告（同时支持 PyMC-NUTS 和 nutpie 采样统计数据）。函数：check_diagnostics() 用于快速检查，create_diagnostic_report() 用于生成包含图表的全面分析报告。

• model_comparison.py：PSIS-LOO 模型比较工具。函数：compare_models()、check_loo_reliability()、model_averaging()。

模板（assets/）

• linear_regression_template.py：贝叶斯线性回归完整模板，包含完整工作流（数据准备、先验检查、拟合、诊断、预测）。

• hierarchical_model_template.py：分层／多层模型完整模板，使用非中心化参数化并进行组级分析。

快速参考

构建模型

with pm.Model(coords={'var': names}) as model:
    # Priors
    param = pm.Normal('param', mu=0, sigma=1, dims='var')
    # Likelihood
    y = pm.Normal('y', mu=..., sigma=..., observed=data)

采样

idata = pm.sample(draws=2000, tune=1000, chains=4, target_accept=0.9)  # returns xarray.DataTree
pm.compute_log_likelihood(idata)  # needed for LOO / az.compare

诊断

from scripts.model_diagnostics import check_diagnostics
check_diagnostics(idata)

模型比较

from scripts.model_comparison import compare_models
compare_models({'m1': idata1, 'm2': idata2})  # PSIS-LOO; higher elpd is better

预测

# X must have been wrapped at build time: pm.Data('X', X, dims=('obs', 'predictors'))
with model:
    pm.set_data({'X': X_new}, coords={'obs': range(len(X_new))})
    pm.sample_posterior_predictive(idata, predictions=True, extend_inferencedata=True)
# predictions land in idata.predictions

附加说明

• PyMC 与 ArviZ 集成，用于可视化和诊断。
• 使用 pm.model_to_graphviz(model) 可视化模型结构（需要 graphviz 包）。
• 使用 idata.to_netcdf('results.nc') 保存结果（需要 h5netcdf 或 netCDF4，例如 arviz[h5netcdf]）；使用 az.from_netcdf('results.nc') 加载。
• 对于非常大的模型，可考虑使用小批量 ADVI 或对数据进行子采样。

版本说明（PyMC 6 / ArviZ 1.x）

PyMC 6.0（2026 年 5 月）迁移到了 ArviZ 1.x 和 PyTensor 3。为 PyMC 5 / ArviZ 0.x 编写的代码会在以下方面失效。

• 结果是 xarray.DataTree，而不是 az.InferenceData。仍可通过组访问结果（idata.posterior、idata["posterior"]）；idata.extend(other) 已改为 idata.update(other)；使用 list(idata.children) 列出组。
• 对数似然：采样后调用 pm.compute_log_likelihood(idata)。pm.sample(idata_kwargs={"log_likelihood": True}) 已弃用。
• 采样器：安装 nutpie 后，它会成为默认的 NUTS 采样器（默认进行 400 次调参；PyMC 自带的 NUTS 仍进行 1000 次）。使用 nuts_sampler="pymc" | "nutpie" | "numpyro" | "blackjax" 显式选择。采样统计名称有所不同（tree_depth 与 depth），因此诊断代码应检查实际存在的名称。初始值应传入 initvals=（使用 start= 会失败）。
• 后端：PyTensor 3 默认使用 Numba 编译，因此首次调用会包含 JIT 编译时间。
• 汇总统计：az.summary 报告 89% 等尾区间（eti89_lb、eti89_ub），而不是 94% HDI（hdi_3%、hdi_97%）。使用 az.summary(idata, ci_kind="hdi", ci_prob=0.94)，或设置 az.rcParams["stats.ci_kind"] / ["stats.ci_prob"]。az.hdi(..., prob=0.95) 取代了 hdi_prob=。
• 绘图返回 PlotCollection（使用 pc.savefig("file.png") 保存，使用 pc.show() 显示）；ax= / axes= 参数已移除。名称变更：plot_posterior → plot_dist、plot_ppc → plot_ppc_dist（先验检查：group="prior_predictive"，使用 num_samples= 代替 num_pp_samples=）、plot_trace → plot_trace_dist（轨迹加密度）、plot_dist_comparison → plot_prior_posterior。使用 arviz[matplotlib] 安装绘图后端。
• 模型比较：az.waic 已移除，az.compare(dict) 不再接受 ic= / scale=；它会按 PSIS-LOO elpd 排名（越高越好），使用 stacking 权重，并报告 elpd_diff、dse、p_worse 和 diag_* 列。
• 若要对新数据进行预测，构建模型时必须将预测变量封装在 pm.Data('X', X, dims=...) 中——只有这样才能通过 pm.set_data({'X': X_new}, coords={...}) 替换它们。直接嵌入计算图的普通 NumPy 数组无法替换。
• pm.sample_prior_predictive 使用 draws=（旧的 samples= 关键字已移除）。
• 对样本外数据进行预测时，调用 pm.sample_posterior_predictive(idata, predictions=True, extend_inferencedata=True, ...)；结果随后保存在 idata.predictions 中，而不是 idata.posterior_predictive 中。