pymc-bayesian-modeling

PyMC 贝叶斯建模

概述

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

何时使用此技能

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

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

标准贝叶斯工作流

按照此工作流构建和验证贝叶斯模型：

1. 数据准备

import pymc as pm
import arviz as az
import numpy as np

# Load and prepare data
X = ...  # Predictors
y = ...  # Outcomes

# Standardize predictors for better sampling
X_mean = X.mean(axis=0)
X_std = X.std(axis=0)
X_scaled = (X - X_mean) / X_std

关键实践：

• 标准化连续预测变量（提高采样效率）
• 尽可能对结果变量进行中心化
• 显式处理缺失数据（将其视为参数）
• 使用带 coords 的命名维度以提高清晰度

2. 模型构建

coords = {
    'predictors': ['var1', 'var2', 'var3'],
    'obs_id': np.arange(len(y))
}

with pm.Model(coords=coords) as model:
    # Priors
    alpha = pm.Normal('alpha', mu=0, sigma=1)
    beta = pm.Normal('beta', mu=0, sigma=1, dims='predictors')
    sigma = pm.HalfNormal('sigma', sigma=1)

    # Linear predictor
    mu = alpha + pm.math.dot(X_scaled, beta)

    # Likelihood
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y, dims='obs_id')

关键实践：

• 使用弱信息先验（而非平坦先验）
• 对尺度参数使用 HalfNormal 或 Exponential
• 尽可能使用命名维度（dims）而非 shape
• 对将用于预测更新的值使用 pm.Data()

3. 先验预测检查

在拟合之前始终验证先验：

with model:
    prior_pred = pm.sample_prior_predictive(samples=1000, random_seed=42)

# Visualize
az.plot_ppc(prior_pred, group='prior')

检查：

• 先验预测是否覆盖合理的值范围？
• 根据领域知识，极端值是否合理？
• 如果先验生成不合理的数据，请调整并重新检查

4. 拟合模型

with model:
    # Optional: Quick exploration with ADVI
    # approx = pm.fit(n=20000)

    # Full MCMC inference
    idata = pm.sample(
        draws=2000,
        tune=1000,
        chains=4,
        target_accept=0.9,
        random_seed=42,
        idata_kwargs={'log_likelihood': True}  # For model comparison
    )

关键参数：

• draws=2000：每条链的样本数
• tune=1000：预热样本数（丢弃）
• chains=4：运行 4 条链以检查收敛性
• target_accept=0.9：对于困难的后验使用更高值（0.95-0.99）
• 包含 log_likelihood=True 以进行模型比较

5. 检查诊断

使用诊断脚本：

from scripts.model_diagnostics import check_diagnostics

results = check_diagnostics(idata, var_names=['alpha', 'beta', 'sigma'])

检查：

• R-hat < 1.01：链已收敛
• ESS > 400：有效样本充足
• 无散度：NUTS 采样成功
• 轨迹图：链应混合良好（毛虫状）

如果出现问题：

• 散度 → 增加 target_accept=0.95，使用非中心化参数化
• ESS 低 → 采样更多 draws，重新参数化以减少相关性
• R-hat 高 → 运行更长时间，检查是否存在多峰性

6. 后验预测检查

验证模型拟合：

with model:
    pm.sample_posterior_predictive(idata, extend_inferencedata=True, random_seed=42)

# Visualize
az.plot_ppc(idata)

检查：

• 后验预测是否捕捉到观测数据的模式？
• 是否存在系统性偏差（模型设定错误）？
• 如果拟合不佳，考虑替代模型

7. 分析结果

# Summary statistics
print(az.summary(idata, var_names=['alpha', 'beta', 'sigma']))

# Posterior distributions
az.plot_posterior(idata, var_names=['alpha', 'beta', 'sigma'])

# Coefficient estimates
az.plot_forest(idata, var_names=['beta'], combined=True)

8. 进行预测

X_new = ...  # New predictor values
X_new_scaled = (X_new - X_mean) / X_std

with model:
    pm.set_data({'X_scaled': X_new_scaled})
    post_pred = pm.sample_posterior_predictive(
        idata.posterior,
        var_names=['y_obs'],
        random_seed=42
    )

# Extract prediction intervals
y_pred_mean = post_pred.posterior_predictive['y_obs'].mean(dim=['chain', 'draw'])
y_pred_hdi = az.hdi(post_pred.posterior_predictive, var_names=['y_obs'])

常见模型模式

线性回归

适用于具有线性关系的连续结果变量：

with pm.Model() as linear_model:
    alpha = pm.Normal('alpha', mu=0, sigma=10)
    beta = pm.Normal('beta', mu=0, sigma=10, shape=n_predictors)
    sigma = pm.HalfNormal('sigma', sigma=1)

    mu = alpha + pm.math.dot(X, beta)
    y = pm.Normal('y', mu=mu, sigma=sigma, observed=y_obs)

使用模板： assets/linear_regression_template.py

逻辑回归

适用于二值结果变量：

with pm.Model() as logistic_model:
    alpha = pm.Normal('alpha', mu=0, sigma=10)
    beta = pm.Normal('beta', mu=0, sigma=10, shape=n_predictors)

    logit_p = alpha + pm.math.dot(X, beta)
    y = pm.Bernoulli('y', logit_p=logit_p, observed=y_obs)

层次模型

适用于分组数据（使用非中心化参数化）：

with pm.Model(coords={'groups': group_names}) as hierarchical_model:
    # Hyperpriors
    mu_alpha = pm.Normal('mu_alpha', mu=0, sigma=10)
    sigma_alpha = pm.HalfNormal('sigma_alpha', sigma=1)

    # Group-level (non-centered)
    alpha_offset = pm.Normal('alpha_offset', mu=0, sigma=1, dims='groups')
    alpha = pm.Deterministic('alpha', mu_alpha + sigma_alpha * alpha_offset, dims='groups')

    # Observation-level
    mu = alpha[group_idx]
    sigma = pm.HalfNormal('sigma', sigma=1)
    y = pm.Normal('y', mu=mu, sigma=sigma, observed=y_obs)

使用模板： assets/hierarchical_model_template.py

关键： 对层次模型始终使用非中心化参数化以避免散度。

泊松回归

适用于计数数据：

with pm.Model() as poisson_model:
    alpha = pm.Normal('alpha', mu=0, sigma=10)
    beta = pm.Normal('beta', mu=0, sigma=10, shape=n_predictors)

    log_lambda = alpha + pm.math.dot(X, beta)
    y = pm.Poisson('y', mu=pm.math.exp(log_lambda), observed=y_obs)

对于过度离散的计数数据，改用 NegativeBinomial。

时间序列

适用于自回归过程：

with pm.Model() as ar_model:
    sigma = pm.HalfNormal('sigma', sigma=1)
    rho = pm.Normal('rho', mu=0, sigma=0.5, shape=ar_order)
    init_dist = pm.Normal.dist(mu=0, sigma=sigma)

    y = pm.AR('y', rho=rho, sigma=sigma, init_dist=init_dist, observed=y_obs)

模型比较

比较模型

使用 LOO 或 WAIC 进行模型比较：

from scripts.model_comparison import compare_models, check_loo_reliability

# Fit models with log_likelihood
models = {
    'Model1': idata1,
    'Model2': idata2,
    'Model3': idata3
}

# Compare using LOO
comparison = compare_models(models, ic='loo')

# Check reliability
check_loo_reliability(models)

解释：

• Δloo < 2：模型相似，选择更简单的模型
• 2 < Δloo < 4：对更好模型的弱证据
• 4 < Δloo < 10：中等证据
• Δloo > 10：对更好模型的强证据

检查 Pareto-k 值：

• k < 0.7：LOO 可靠
• k > 0.7：考虑 WAIC 或 k 折交叉验证

模型平均

当模型相似时，对预测结果取平均：

from scripts.model_comparison import model_averaging

averaged_pred, weights = model_averaging(models, var_name='y_obs')

分布选择指南

先验

尺度参数（σ、τ）：

• 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.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.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
    start = approx.sample(return_inferencedata=False)[0]
    idata = pm.sample(start=start)

权衡：

• 比 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：使用 LOO/WAIC 的模型比较工具。函数：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': X_new})
    pred = pm.sample_posterior_predictive(idata.posterior)

附加说明

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