SuperAPI AI 工具集
合作ENAI 行业资讯收藏
Skill 广场pymc

pymc

数据与分析

查看项目

更新于 11-04 02:43

技能介绍

PyMC 是用于贝叶斯建模与概率编程的 Python 库。本技能介绍如何使用 PyMC 5.x 及以上版本的现代 API 构建、拟合、验证和比较贝叶斯模型，涉及层次模型、MCMC 采样（NUTS）、变分推断和模型比较（LOO、WAIC）。

适用场景

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

标准工作流

  1. 数据准备：导入 pymc、arviz、numpy；对连续预测变量标准化以改善采样效率；在可能时中心化结果变量；显式处理缺失数据（作为参数对待）；使用 coords 命名维度以保持清晰。
  2. 模型构建：使用 pm.Model(coords=...) 定义模型；先验使用弱信息先验（而非平坦先验）；尺度参数使用 HalfNormal 或 Exponential；尽量使用命名维度（dims）而非 shape；对需要在预测时更新的值使用 pm.Data()。
  3. 先验预测检查：拟合前使用 pm.sample_prior_predictive(samples=1000, random_seed=42) 并用 az.plot_ppc(prior_pred, group='prior') 可视化，检查先验预测是否覆盖合理取值；若先验产生不合理数据则调整并重新检查。
  4. 拟合模型：可选先使用 ADVI 快速探索（pm.fit(n=20000)）；完整 MCMC 推断使用 pm.sample(draws=2000, tune=1000, chains=4, target_accept=0.9, random_seed=42, idata_kwargs={'log_likelihood': True})。参数含义：draws 为每条链样本数，tune 为丢弃的预热样本，chains 为链数（4 条用于收敛检查），target_accept 在困难后验时可提高到 0.95-0.99，log_likelihood=True 用于模型比较。
  5. 检查诊断：使用 scripts/model_diagnostics.py 中的 check_diagnostics(idata, var_names=[...])；检查 R-hat < 1.01（链已收敛）、ESS > 400（有效样本充足）、无发散（NUTS 采样成功）、迹图混合良好。出现问题时的处理：发散 → 提高 target_accept=0.95、使用非中心化参数化；ESS 低 → 增加采样、重参数化降低相关性；R-hat 高 → 运行更久、检查多峰性。
  6. 后验预测检查：pm.sample_posterior_predictive(idata, extend_inferencedata=True, random_seed=42)，并用 az.plot_ppc(idata) 检查后验预测是否捕捉到观测数据模式，是否存在系统性偏差（模型设定错误）；若拟合差则考虑替代模型。
  7. 分析结果：az.summary 输出统计摘要；az.plot_posterior 查看后验分布；az.plot_forest 查看系数估计。
  8. 进行预测：对新预测变量按训练时的均值/标准差标准化；pm.set_data({'X_scaled': X_new_scaled})；pm.sample_posterior_predictive(idata.posterior, var_names=['y_obs'], random_seed=42)；提取预测均值与 az.hdi 预测区间。

常见模型模式

  • 线性回归：连续结果与线性关系；参考 assets/linear_regression_template.py。
  • 逻辑回归：二分类结果，使用 pm.Bernoulli('y', logit_p=logit_p, observed=y_obs)。
  • 层次模型：分组数据使用非中心化参数化（alpha = mu_alpha + sigma_alpha * alpha_offset），避免发散；参考 assets/hierarchical_model_template.py。
  • 泊松回归：计数数据，使用 pm.Poisson；过度分散计数改用 NegativeBinomial。
  • 时间序列：自回归过程可使用 pm.AR('y', rho=rho, sigma=sigma, init_dist=init_dist, observed=y_obs)。

模型比较 使用 scripts/model_comparison.py 中的 compare_models(models, ic='loo') 与 check_loo_reliability(models)。解释：Δloo < 2 模型相似、选更简单模型；2-4 弱证据；4-10 中等证据；>10 强证据。Pareto-k：k < 0.7 时 LOO 可靠，k > 0.7 可考虑 WAIC 或 k 折交叉验证。模型相似时可用 model_averaging(models, var_name='y_obs') 对预测做模型平均。

分布选择

  • 先验：尺度参数默认 pm.HalfNormal('sigma', sigma=1)，也可用 Exponential 或 Gamma；无界参数用 Normal（标准化数据）或 StudentT（抗离群值）；正参数用 LogNormal 或 Gamma；概率用 Beta（弱信息）或 Uniform（少用）；相关矩阵用 LKJCorr。
  • 似然：连续结果默认 Normal，抗离群值用 StudentT；计数用 Poisson（等分散）、NegativeBinomial（过度分散）或 ZeroInflatedPoisson（零过多）；二分类用 Bernoulli；分类用 Categorical。详见 references/distributions.md。

采样与推断

  • MCMC（NUTS）为默认推荐：发散时提高 target_accept；采样慢可用 ADVI 初始化；离散参数可用 pm.Metropolis()。
  • 变分推断（ADVI）：pm.fit(n=20000, method='advi')，远快于 MCMC 但为近似结果（可能低估不确定性），适合大模型快速探索或初始化。详见 references/sampling_inference.md。

诊断脚本

  • create_diagnostic_report(idata, var_names=[...], output_dir='diagnostics/') 生成迹图、秩图、自相关图、能量图、ESS 演化图与统计摘要 CSV。
  • check_diagnostics(idata) 快速检查 R-hat、ESS、发散和树深度。

常见问题与解决方案

  • 发散（diverging.sum() > 0）：提高 target_accept、使用非中心化参数化、加强先验约束、检查模型设定错误。
  • ESS 低（< 400）：增加 draws（如 5000）、重参数化降低后验相关性、对相关预测变量使用 QR 分解。
  • R-hat 高（> 1.01）：延长运行（tune=2000, draws=5000）、检查多峰性、用 ADVI 改善初始化。
  • 采样慢：ADVI 初始化、降低模型复杂度、增加并行（cores=8, chains=8）、在合适时使用变分推断。

最佳实践

  • 建模：始终标准化预测变量；使用弱信息先验；使用命名维度；层次模型使用非中心化参数化；拟合前做先验预测检查。
  • 采样：至少 4 条链；target_accept=0.9 为基线；包含 log_likelihood=True；设置随机种子以保证可复现。
  • 验证：解释前先检查诊断；做后验预测检查；适当时比较多个模型；报告不确定性（HDI 区间而非仅点估计）。
  • 工作流：从简单开始逐步增加复杂度；先验预测检查 → 拟合 → 诊断 → 后验预测检查；根据检查迭代模型设定；记录假设与先验选择。

配套资源

  • references/distributions.md：按类别（连续、离散、多变量、混合、时间序列）组织的 PyMC 分布目录，用于选择先验或似然。
  • references/sampling_inference.md：采样算法（NUTS、Metropolis、SMC）、变分推断（ADVI、SVGD）及采样问题处理指南。
  • references/workflows.md：常见模型类型、数据准备、先验选择和模型验证的完整工作流示例与代码模式。
  • scripts/model_diagnostics.py：自动诊断检查与报告生成。
  • scripts/model_comparison.py：基于 LOO/WAIC 的模型比较工具。
  • assets/linear_regression_template.py：贝叶斯线性回归完整模板（数据准备、先验检查、拟合、诊断、预测）。
  • assets/hierarchical_model_template.py：层次/多级模型完整模板（非中心化参数化与组级分析）。

其他说明

  • PyMC 与 ArviZ 集成，用于可视化和诊断。
  • 可用 pm.model_to_graphviz(model) 可视化模型结构。
  • 结果可通过 idata.to_netcdf('results.nc') 保存，并用 az.from_netcdf('results.nc') 加载。
  • 超大模型可考虑小批量 ADVI 或数据子采样。