pymc

PyMC 贝叶斯建模

概述

PyMC 是一个用于贝叶斯建模和概率编程的 Python 库。使用 PyMC 的现代 API（6.x 及以上版本）构建、拟合、验证和比较贝叶斯模型，包括分层模型、MCMC 采样（NUTS）、变分推断、后验预测检验和模型比较（LOO、WAIC）。

当前版本和设置

截至 2026 年 6 月，PyMC 6.0.1 是当前稳定版本。它要求 Python 3.12+，使用 PyTensor 3 作为计算图后端，并默认使用 Numba 等编译后端。为确保本地环境可复现，请固定版本：

uv pip install "pymc[nutpie]==6.0.1"

nutpie 扩展可启用更快的 Rust/Numba NUTS 实现。如果使用 NumPyro 或 BlackJAX，请在同一环境中安装这些可选采样器依赖，并在项目锁定文件中固定版本。

何时使用此技能

在以下情况下应使用此技能：

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

标准贝叶斯工作流

切勿先采样、后检查。包含代码的八步工作流记录于 references/standard_workflow.md：

1. 数据准备——包括对预测变量进行标准化，以便先验更易解释。
2. 模型构建——在 pm.Model 上下文中设置先验和似然。
3. 先验预测检验——在拟合前确认先验所隐含的数据是合理的。
4. 拟合模型——使用显式随机种子调用 pm.sample()。
5. 检查诊断结果——R-hat、ESS、发散情况。发散会使拟合结果无效；应修正模型或重新参数化，而不是调高 target_accept 后寄希望于问题消失。
6. 后验预测检验——拟合后的模型能否重现观测数据？
7. 分析结果——根据后验分布生成摘要和区间。
8. 进行预测——通过 pm.set_data 和后验预测采样对新数据进行预测。

可复用的模型结构和模型比较见 references/model_patterns.md。

分布选择指南

先验分布

尺度参数（σ、τ）：

• pm.HalfNormal('sigma', sigma=1)——默认选择
• pm.Exponential('sigma', lam=1)——替代方案
• pm.Gamma('sigma', alpha=2, beta=1)——信息量更大的选择

无界参数：

• pm.Normal('theta', mu=0, sigma=1)——适用于标准化数据
• pm.StudentT('theta', nu=3, mu=0, sigma=1)——对异常值具有稳健性

正参数：

• pm.LogNormal('theta', mu=0, sigma=1)
• pm.Gamma('theta', alpha=2, beta=1)

概率：

• pm.Beta('p', alpha=2, beta=2)——弱信息先验
• pm.Uniform('p', lower=0, upper=1)——非信息先验（谨慎使用）

相关矩阵：

• pm.LKJCholeskyCov('chol', n=n_vars, eta=2, sd_dist=pm.HalfNormal.dist(1))——首选的协方差先验
• pm.LKJCorr('corr', n=n_vars, eta=2)——仅用于相关性的先验；eta=1 时为均匀分布，eta>1 时更偏好单位矩阵

似然函数

连续结果：

• pm.Normal('y', mu=mu, sigma=sigma) - 连续数据的默认选择
• pm.StudentT('y', nu=nu, mu=mu, sigma=sigma) - 对异常值具有鲁棒性

计数数据：

• pm.Poisson('y', mu=lambda) - 等离散计数
• pm.NegativeBinomial('y', mu=mu, alpha=alpha) - 过度离散计数
• pm.ZeroInflatedPoisson('y', psi=psi, mu=mu) - 零值过多
• pm.HurdleNegativeBinomial('y', psi=psi, mu=mu, alpha=alpha) - 零值过多且过度离散

二元结果：

• pm.Bernoulli('y', p=p) 或 pm.Bernoulli('y', logit_p=logit_p)

分类结果：

• pm.Categorical('y', p=probs)

参见： references/distributions.md，了解完整的分布参考

采样与推断

使用 NUTS 进行 MCMC

大多数模型的默认推荐方法：

idata = pm.sample(
    draws=2000,
    tune=1000,
    chains=4,
    target_accept=0.9,
    random_seed=42
)

必要时进行调整：

• 发散 → target_accept=0.95 或更高
• 采样缓慢 → 使用 ADVI 进行初始化
• 离散参数 → 对离散变量使用 pm.Metropolis()

变分推断

用于探索或初始化的快速近似方法：

with model:
    approx = pm.fit(n=20000, method='advi')

    # Use for initialization
    initvals = approx.sample(return_inferencedata=False)[0]
    idata = pm.sample(initvals=initvals)

权衡：

• 比 MCMC 快得多
• 属于近似方法（可能低估不确定性）
• 适用于大型模型或快速探索

参见： references/sampling_inference.md，了解详细的采样指南

诊断脚本

全面诊断

from scripts.model_diagnostics import create_diagnostic_report

create_diagnostic_report(
    idata,
    var_names=['alpha', 'beta', 'sigma'],
    output_dir='diagnostics/'
)

生成：

• 轨迹图
• 秩图（检查混合情况）
• 自相关图
• 能量图
• 局部 ESS 图
• 汇总统计 CSV

快速诊断检查

from scripts.model_diagnostics import check_diagnostics

results = check_diagnostics(idata)

检查 R-hat、ESS、发散情况和树深度。

常见问题与解决方案

发散

症状： idata.sample_stats.diverging.sum() > 0

解决方案：

1. 将 target_accept=0.95 或 0.99 调高
2. 使用非中心化参数化（层级模型）
3. 增加强先验以约束参数
4. 检查模型设定是否有误

有效样本量偏低

症状： ESS < 400

解决方案：

1. 增加采样次数：draws=5000
2. 重新参数化以降低后验相关性
3. 对预测变量相关的回归使用 QR 分解

R-hat 偏高

症状： R-hat > 1.01

解决方案：

1. 延长链的运行时间：tune=2000, draws=5000
2. 检查是否存在多峰性
3. 使用 ADVI 改进初始化

采样缓慢

解决方案：

1. 使用 ADVI 初始化
2. 降低模型复杂度
3. 增加并行化程度：cores=8, chains=8
4. 适用时使用变分推断

最佳实践

模型构建

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

采样

1. 运行多条链（至少 4 条），以检查收敛情况
2. 将 target_accept=0.9 作为基准值（必要时调高）
3. 设置 log_likelihood=True，以便比较模型
4. 设置随机种子，以确保结果可复现

验证

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

工作流

1. 从简单模型开始，逐步增加复杂度
2. 先验预测检查 → 拟合 → 诊断 → 后验预测检查
3. 根据检查结果迭代调整模型设定
4. 记录假设和先验选择

资源

此技能包含：

参考资料（references/）

• distributions.md：按类别（连续、离散、多变量、混合、时间序列）整理的 PyMC 分布完整目录。选择先验或似然时使用。

• sampling_inference.md：关于采样算法（NUTS、Metropolis、SMC）、变分推断（ADVI、SVGD）以及处理采样问题的详细指南。遇到收敛问题或选择推断方法时使用。

• workflows.md：常见模型类型、数据准备、先验选择和模型验证的完整工作流示例与代码模式。可作为标准贝叶斯分析的操作手册。

脚本（scripts/）

• model_diagnostics.py：自动诊断检查与报告生成。函数：check_diagnostics() 用于快速检查，create_diagnostic_report() 用于生成包含图表的综合分析报告。

• model_comparison.py：基于 PSIS-LOO ELPD 的模型比较工具；这是 ArviZ 1.x compare() 唯一用于排序的准则。函数：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)

诊断

from scripts.model_diagnostics import check_diagnostics
check_diagnostics(idata)

模型比较

from scripts.model_comparison import compare_models
compare_models({'m1': idata1, 'm2': idata2}, ic='loo')

预测

with model:
    pm.set_data({'X_data': X_new})
    pred = pm.sample_posterior_predictive(idata, predictions=True)

补充说明

• PyMC 与 ArviZ 集成，用于可视化和诊断；PyMC 6 / ArviZ 1 使用 xarray DataTree，同时保留 .posterior 和 .posterior_predictive 等熟悉的组
• 使用 pm.model_to_graphviz(model) 可视化模型结构
• 使用 idata.to_netcdf('results.nc') 保存结果
• 使用 az.from_netcdf('results.nc') 加载
• 对于非常大的模型，可考虑使用 minibatch ADVI 或数据子采样

Agent 操作流程

1. 检查环境。 确认 Python 环境和库版本（python -c "import pkg; print(pkg.__version__)"），并检查数据的形状、类型和缺失值。
2. 明确输入。 确认数据格式、标识符和参数。对于会影响结果的值，应询问而非猜测。
3. 先运行小规模版本。 先在样本或单个折上运行，并检查运行时间和内存。
4. 按照上述说明和参考资料执行完整任务。
5. 验证结果。 使用留出数据、固定随机种子和适当的指标；检查是否存在数据泄漏；报告不确定性（置信区间、不同随机种子的标准差）。
6. 报告。 说明运行了什么（版本、命令、参数）、检查了什么，以及仍有哪些不确定之处。

如果发生这种情况	请这样做
内存不足或运行非常缓慢	进行子采样，使用分块或惰性计算，或缩小模型规模，并告知用户做了哪些更改。
这些说明中提到的函数、标志或端点在已安装版本中不存在	查阅已安装版本自带的文档（help()、--help、官方文档），进行调整，并告知用户。切勿编造 API。
必需的输入、标识符或参数不明确	询问用户，或在运行前明确说明假设。

诚信规则

• 切勿编造结果、参数、标识符、引文或统计数据。如果某项内容无法运行或验证，应如实说明。
• 切勿报告本次会话中未计算的指标；应展示生成每个数字的代码路径。
• 请将此处与版本相关的细节视为可能已过时，并根据已安装版本的官方文档进行确认。
• 对于需要付费、占用共享 GPU 或云配额、涉及个人或患者数据，或无法撤销的操作，务必事先询问。

相关技能

• statistical-analysis：针对研究数据的引导式统计分析——检验选择、假设检查、效应量、功效分析、贝叶斯替代方案，等等……