pymc-bayesian-modeling

PyMC 贝叶斯建模

概述

PyMC 是一个用于贝叶斯统计建模和概率编程的 Python 库。它提供了简洁灵活的语法来定义概率模型，并通过 MCMC（NUTS）和变分方法（ADVI）高效地进行推断。本技能涵盖完整的贝叶斯建模周期，从模型设定到诊断、比较和预测。

适用场景

• 在估计参数时全面量化不确定性（可信区间，而不只是点估计）
• 对分组或嵌套数据拟合层级／多层模型
• 执行先验和后验预测检验，以验证模型假设
• 使用信息准则（LOO-CV、WAIC）比较候选模型
• 在贝叶斯框架下构建回归模型（线性、逻辑、Poisson）
• 将缺失数据或测量误差作为潜在参数处理
• 使用自回归或随机游走先验对时间序列建模
• 为新观测生成带有不确定性范围的后验预测
• 对于大型模型的编译型、可扩展性更强的贝叶斯推断，请改用 Stan/PyStan；对于频率学派统计检验，请使用 statsmodels

前置条件

• Python 包：pymc >= 5.0、arviz、numpy、matplotlib
• 数据：包含数值列的 NumPy 数组或 pandas DataFrame
• 环境：大多数模型使用 CPU 即可；大型模型可通过 JAX 后端使用 GPU

pip install pymc arviz numpy matplotlib
# Optional: JAX backend for GPU acceleration
pip install pymc[jax]

快速入门

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

# Simulate data
np.random.seed(42)
X = np.random.randn(100)
y = 2.5 + 1.3 * X + np.random.randn(100) * 0.5

# Build and fit model
with pm.Model() as model:
    alpha = pm.Normal("alpha", mu=0, sigma=5)
    beta = pm.Normal("beta", mu=0, sigma=5)
    sigma = pm.HalfNormal("sigma", sigma=1)
    mu = alpha + beta * X
    y_obs = pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y)
    idata = pm.sample(1000, tune=1000, chains=4, random_seed=42)

print(az.summary(idata, var_names=["alpha", "beta", "sigma"]))
# Expected: alpha ~ 2.5, beta ~ 1.3, sigma ~ 0.5

工作流

第 1 步：准备数据

对连续预测变量进行标准化，以提高采样效率。使用命名坐标，让模型更易读，并便于与 ArviZ 集成。

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

# Load data
X = np.random.randn(200, 3)  # 200 obs, 3 predictors
y = X @ np.array([1.0, -0.5, 0.3]) + np.random.randn(200) * 0.8

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

# Define coordinates for named dimensions
coords = {
    "predictors": ["var1", "var2", "var3"],
    "obs_id": np.arange(len(y)),
}
print(f"Data shape: X={X_scaled.shape}, y={y.shape}")

第 2 步：定义模型并设置先验

在 pm.Model() 上下文中指定模型结构。使用弱信息先验、用 dims 指定命名维度，并使用 HalfNormal 或 Exponential 作为尺度参数。

with pm.Model(coords=coords) as model:
    # Priors — weakly informative, not flat
    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")

# Inspect model variables
print(model.basic_RVs)  # Lists: [alpha, beta, sigma, y_obs]

第 3 步：先验预测检验

在拟合前验证先验是否会生成合理范围的数据。如果模拟数据不合理，请调整先验。

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

# Check prior-implied data range
prior_y = prior_pred.prior_predictive["y_obs"].values.flatten()
print(f"Prior predictive range: [{prior_y.min():.1f}, {prior_y.max():.1f}]")
print(f"Observed data range:    [{y.min():.1f}, {y.max():.1f}]")

az.plot_ppc(prior_pred, group="prior", num_pp_samples=100)

第 4 步：采样后验（MCMC）

使用多个链运行 NUTS 采样。如果之后要比较模型，请包含 log_likelihood=True。

with model:
    idata = pm.sample(
        draws=2000,
        tune=1000,
        chains=4,
        target_accept=0.9,
        random_seed=42,
        idata_kwargs={"log_likelihood": True},
    )

print(f"Posterior shape: {idata.posterior['beta'].shape}")
# Expected: (4 chains, 2000 draws, 3 predictors)

第 5 步：诊断采样

在解读结果前检查收敛情况。三项诊断（R-hat、ESS、发散）都必须通过。

# Summary with convergence diagnostics
summary = az.summary(idata, var_names=["alpha", "beta", "sigma"])
print(summary[["mean", "sd", "hdi_3%", "hdi_97%", "r_hat", "ess_bulk"]])

# R-hat convergence check
bad_rhat = summary[summary["r_hat"] > 1.01]
if len(bad_rhat) > 0:
    print(f"WARNING: {len(bad_rhat)} parameters with R-hat > 1.01")
    print(bad_rhat[["r_hat"]])

# Effective sample size check
low_ess = summary[summary["ess_bulk"] < 400]
if len(low_ess) > 0:
    print(f"WARNING: {len(low_ess)} parameters with ESS < 400")

# Divergence check
n_div = idata.sample_stats.diverging.sum().item()
total = len(idata.posterior.draw) * len(idata.posterior.chain)
print(f"Divergences: {n_div}/{total} ({n_div / total * 100:.2f}%)")

# Visual diagnostics — trace plots and rank plots
az.plot_trace(idata, var_names=["alpha", "beta", "sigma"])
az.plot_rank(idata, var_names=["alpha", "beta", "sigma"])

第 6 步：后验预测检验

通过比较后验模拟数据与观测数据来验证模型拟合情况。

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

az.plot_ppc(idata, num_pp_samples=100)
# Blue = observed data, grey = posterior simulations
# Systematic deviations indicate model misspecification

第 7 步：比较模型

使用 LOO-CV 或 WAIC 比较候选模型。信息准则越低越好。

# Fit multiple models with log_likelihood=True, then compare
# Example: compare linear vs a second model
idatas = {"linear": idata}  # add more fitted models here

comparison = az.compare(idatas, ic="loo")
print(comparison[["rank", "elpd_loo", "p_loo", "d_loo", "weight"]])

# Check LOO reliability via Pareto-k diagnostics
loo_result = az.loo(idata, pointwise=True)
high_k = (loo_result.pareto_k > 0.7).sum().item()
print(f"Observations with Pareto-k > 0.7: {high_k}")
# Interpretation: Dloo < 2 = similar models; Dloo > 10 = strong evidence

az.plot_compare(comparison)

第 8 步：生成预测

为新数据生成后验预测，并全面传播不确定性。

X_new = np.array([[0.5, -1.0, 0.2]])
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
    )

y_pred = post_pred.posterior_predictive["y_obs"]
print(f"Predicted mean: {y_pred.mean().item():.3f}")
print(f"94% HDI: {az.hdi(y_pred, hdi_prob=0.94).values}")

关键参数

参数	默认值	范围／选项	作用
draws	1000	500-10000	每条链的后验样本数
tune	1000	500-5000	预热迭代次数（会被丢弃）；复杂后验需要增加此值
chains	4	2-8	独立链的数量；要可靠地计算 R-hat，至少需要 4 条链
cores	all CPUs	1-N	并行运行的链数；要实现完全并行，请设为与 chains 相同
target_accept	0.8	0.8-0.99	NUTS 接受率；提高此值可减少发散
init	"auto"	"adapt_diag"、"jitter+adapt_diag"、"advi"	采样器的初始化策略
random_seed	None	任意整数	用于确保结果可复现的种子
idata_kwargs	{}	{"log_likelihood": True}	存储对数似然，用于 LOO/WAIC 模型比较
method (pm.fit)	"advi"	"advi"、"fullrank_advi"、"svgd"	变分推断算法
n (pm.fit)	10000	5000-100000	VI 优化迭代次数
samples (先验预测)	500	100-5000	用于验证的先验预测样本数

核心概念

先验／分布选择指南

分布	适用场景	关键参数
Normal(mu, sigma)	无界实值参数（标准化数据）	mu：中心，sigma：离散程度
HalfNormal(sigma)	尺度／标准差参数（正值）	sigma：正半轴部分的离散程度
Exponential(lam)	尺度参数，是 HalfNormal 的替代方案	lam：速率（1／均值）
StudentT(nu, mu, sigma)	Normal 的稳健替代方案（抗离群值）	nu：自由度（<10 = 厚尾程度更高）
Beta(alpha, beta)	[0,1] 范围内的概率或比例	alpha=beta=2：弱信息先验
Gamma(alpha, beta)	正值参数（速率、浓度）	alpha：形状，beta：速率
LogNormal(mu, sigma)	存在乘法效应的正值参数	mu、sigma：底层 Normal 分布的参数
LKJCorr(n, eta)	相关矩阵先验	eta=1：均匀；eta>1：倾向于单位矩阵
Dirichlet(a)	概率向量（总和为 1）	a：浓度；所有值相等时为均匀分布
Bernoulli(p / logit_p)	二元结果似然	使用 logit_p 以提高数值稳定性
Poisson(mu)	计数数据（等离散）	mu：速率；过度离散时使用 NegBinomial
NegativeBinomial(mu, alpha)	过度离散的计数数据	alpha：离散度（越小表示过度离散越明显）

诊断阈值

指标	阈值	解释	未达标时的操作
R-hat	< 1.01	链已收敛	延长链运行时间；检查多峰性
ESS bulk	> 400	独立样本数量充足	增加 draws；重新参数化
ESS tail	> 400	尾部估计可靠	增加 draws
Divergences	0	NUTS 已成功探索	增加 target_accept；使用非中心参数化
Pareto-k (LOO)	< 0.7	LOO 估计可靠	使用 WAIC 或 k 折交叉验证
Max tree depth	< 10	轨迹未被截断	重新参数化或增加 max_treedepth

模型变体概览

问题类型	配方	似然	关键特征
分组／嵌套数据	层级模型	Normal（各异）	非中心参数化、部分池化
二元结果	逻辑回归	Bernoulli	logit_p 链接函数
非线性／空间数据	高斯过程	Normal	基于核的协方差、灵活的形状
计数数据	（在工作流中使用 Poisson）	Poisson / NegBinomial	对数链接；过度离散时使用 NegBinomial
时间序列	（参见参考资料）	AR / GaussianRandomWalk	自回归系数
混合／聚类	（参见参考资料）	Mixture / NormalMixture	通过 Dirichlet 设置成分权重

常见配方

配方：层级模型

适用场景：数据具有自然分组结构（医院内的患者、学校内的学生）。非中心参数化可避免漏斗几何导致的发散。

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

n_groups, n_per_group = 5, 30
group_idx = np.repeat(np.arange(n_groups), n_per_group)
group_names = [f"group_{i}" for i in range(n_groups)]

# Simulated grouped data
true_alphas = np.random.normal(3.0, 1.5, n_groups)
y_obs = np.random.normal(true_alphas[group_idx], 0.5)

with pm.Model(coords={"groups": group_names}) as hierarchical_model:
    # Hyperpriors (population level)
    mu_alpha = pm.Normal("mu_alpha", mu=0, sigma=5)
    sigma_alpha = pm.HalfNormal("sigma_alpha", sigma=2)

    # Non-centered parameterization
    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
    sigma = pm.HalfNormal("sigma", sigma=1)
    y = pm.Normal("y", mu=alpha[group_idx], sigma=sigma, observed=y_obs)

    idata_hier = pm.sample(2000, tune=1000, target_accept=0.95, random_seed=42)

print(az.summary(idata_hier, var_names=["mu_alpha", "sigma_alpha", "alpha", "sigma"]))

配方：逻辑回归

适用场景：二元结果变量（成功／失败、患病／健康）。使用 logit_p 以提高数值稳定性，而不是直接计算概率。

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

# Simulated binary outcome
np.random.seed(42)
n, p = 200, 3
X_lr = np.random.randn(n, p)
true_beta = np.array([1.0, -0.5, 0.3])
prob = 1 / (1 + np.exp(-(0.2 + X_lr @ true_beta)))
y_binary = np.random.binomial(1, prob)

with pm.Model() as logistic_model:
    alpha = pm.Normal("alpha", mu=0, sigma=2)
    beta = pm.Normal("beta", mu=0, sigma=2, shape=p)
    logit_p = alpha + pm.math.dot(X_lr, beta)
    y = pm.Bernoulli("y", logit_p=logit_p, observed=y_binary)

    idata_logit = pm.sample(2000, tune=1000, target_accept=0.9, random_seed=42)

summary = az.summary(idata_logit, var_names=["alpha", "beta"])
print(summary[["mean", "sd", "hdi_3%", "hdi_97%"]])
# Odds ratios
beta_samples = idata_logit.posterior["beta"].values.reshape(-1, p)
print(f"Odds ratios (median): {np.median(np.exp(beta_samples), axis=0)}")

配方：高斯过程

适用场景：对未知非线性函数、空间数据或平滑潜在过程建模。n > 1000 时计算开销较大；对于更大的数据集，可考虑稀疏近似。

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

# Simulated nonlinear data
np.random.seed(42)
X_gp = np.sort(np.random.uniform(0, 10, 60))[:, None]
y_gp = np.sin(X_gp[:, 0]) + np.random.randn(60) * 0.3

with pm.Model() as gp_model:
    # GP hyperparameters
    ls = pm.Gamma("lengthscale", alpha=2, beta=1)
    eta = pm.HalfNormal("amplitude", sigma=2)
    sigma_noise = pm.HalfNormal("noise", sigma=0.5)

    # Covariance function
    cov = eta**2 * pm.gp.cov.Matern52(1, ls=ls)
    gp = pm.gp.Marginal(cov_func=cov)

    # Marginal likelihood
    y_ = gp.marginal_likelihood("y", X=X_gp, y=y_gp, sigma=sigma_noise)
    idata_gp = pm.sample(1000, tune=1000, target_accept=0.9, random_seed=42)

# Predict on new points
X_new_gp = np.linspace(0, 10, 100)[:, None]
with gp_model:
    f_pred = gp.conditional("f_pred", X_new_gp)
    pred_samples = pm.sample_posterior_predictive(idata_gp, var_names=["f_pred"])

print(f"GP predictions shape: {pred_samples.posterior_predictive['f_pred'].shape}")

预期输出

• 轨迹图（az.plot_trace）——链混合可视化；健康的链看起来像“毛茸茸的毛毛虫”
• 后验汇总（az.summary）——每个参数的均值、标准差、HDI、R-hat、ESS
• 先验／后验预测图（az.plot_ppc）——模拟数据与观测数据的分布
• 森林图（az.plot_forest）——系数估计值及其可信区间
• 模型比较表（az.compare）——按 LOO／WAIC 排序的模型、权重和警告
• 秩图（az.plot_rank）——均匀直方图表示混合良好
• 能量图（az.plot_energy）——HMC 能量转换诊断
• InferenceData 文件（idata.to_netcdf("results.nc")）——序列化结果，供后续分析

故障排除

问题	原因	解决方案
转移发散	NUTS 难以处理后验几何结构（漏斗、脊线）	将 target_accept=0.95 调高；使用非中心化参数化；添加更强的先验
ESS 较低（< 400）	样本之间存在高度自相关	将 draws=5000 调高；重新参数化以降低相关性；对于相关预测变量，尝试 QR 分解
R-hat > 1.01	链尚未收敛	将 tune=2000, draws=5000 调高；检查是否存在多峰；使用 ADVI 初始化
采样缓慢	后验复杂或数据集较大	使用 ADVI 初始化；降低模型复杂度；增加 cores；尝试 JAX 后端
LOO 中 Pareto-k > 0.7	有影响力的观测值影响 LOO 估计	改用 WAIC；调查有影响力的点；考虑 k 折交叉验证
达到最大树深度	模型几何结构需要较长的轨迹	重新参数化模型；调高 max_treedepth 参数
启动时出现 SamplingError	初始化不佳	尝试 init="adapt_diag" 或 init="jitter+adapt_diag"；通过 MAP 提供 initvals
后验有偏（先验占主导）	相对于数据而言，先验过强	减弱先验（增大 sigma）；检查先验预测；确认数据传递正确
ValueError: logp = -inf	参数进入不可能区域	检查数据是否含 NaN/Inf；确保正值参数使用 HalfNormal/Gamma；确认似然与数据类型匹配

最佳实践

1. 使用弱信息先验，而非平坦先验。平坦先验（大范围内的 Uniform）会导致采样困难，且很少适用。
2. 拟合前务必进行先验预测检查。如果先验隐含的数据不合理，请在采样前调整先验。
3. 将连续预测变量标准化为均值=0、标准差=1，以改善 NUTS 几何结构并加快采样。
4. 层级模型使用非中心化参数化：使用 offset ~ Normal(0,1); param = mu + sigma * offset，而非 param ~ Normal(mu, sigma)。
5. 至少运行 4 条链，以获得可靠的 R-hat 和 ESS 诊断。切勿相信单条链的结果。
6. 检查全部三项诊断指标（R-hat < 1.01、ESS > 400、发散次数为 0），再解读后验。
7. 如果计划通过 LOO 或 WAIC 比较模型，拟合时请包含 log_likelihood=True。
8. 对于命名维度，使用 dims 而非 shape——这样可以与 ArviZ 集成，支持带标签的汇总和子集选取。
9. 报告可信区间（HDI），而不只是点估计。结果是后验分布，而不是单个数字。
10. 使用 idata.to_netcdf("results.nc") 保存结果，以便复现和后续重新分析。

捆绑资源

两个参考文件提供更深入的详情，供按需查阅：

• references/distributions_inference.md —— 按类别整理的分布目录（连续、离散、多元、混合、时间序列、特殊/修饰符），包含完整的参数签名、支持情况和常见用途。采样与推断方法（NUTS、Metropolis、Slice、CompoundStep、SMC、ADVI、fullrank ADVI、SVGD、MAP），附有代码示例和方法选择指南表。重参数化技巧（非中心化、QR 分解）附有代码。 内容涵盖：由原始的 references/distributions.md（320 行）和 references/sampling_inference.md（424 行）合并整理而成。 内联迁移：核心分布选择指南表（关键概念）、诊断阈值表（关键概念）、ADVI 基础用法（原文件中的配方：变分推断，现已在此处的参考资料中介绍）。 未纳入：形状广播示例（简单的 NumPy 风格 shape=5）、dims 基础用法（已在工作流步骤 2 中介绍）、先验预测采样（已在工作流步骤 3 中介绍）。

• references/advanced_workflows.md —— 模型比较工作流（LOO、WAIC、Pareto-k 可靠性检查、加权预测的模型平均）、全面的诊断报告生成（轨迹图、秩图、自相关、能量、ESS 演变）、先验与后验比较模式、数据准备最佳实践（标准化、中心化、缺失数据插补）、先验选择指南（弱信息先验与结合领域知识的有信息先验）、使用 xarray 子集选择的命名维度（dims）、保存/加载模式（NetCDF、pickle）。混合模型模式。 内容涵盖：由原始的 references/workflows.md（526 行）合并整理而成。 内联迁移：核心诊断检查逻辑（工作流步骤 5）、模型比较基础知识（工作流步骤 7）、先验预测检查（工作流步骤 3）、线性/逻辑/层级模型代码（配方）。 未纳入：完整的单体工作流模板（与 8 步工作流部分重复）。

各参考文件的处置情况（所有原文件）

原始 references/（3 个文件）：

1. distributions.md（320 行）-> 与 sampling_inference.md 合并到新文件 references/distributions_inference.md
2. sampling_inference.md（424 行）-> 与 distributions.md 合并到新文件 references/distributions_inference.md
3. workflows.md（526 行）-> 迁移为新文件 references/advanced_workflows.md

原始 scripts/（2 个文件）：

1. model_diagnostics.py（350 行）-> check_diagnostics() 逻辑（R-hat、ESS、发散、树深度检查）内联至工作流步骤 5；create_diagnostic_report() 的绘图生成模式在“预期输出”中描述，并在 references/advanced_workflows.md 中引用；compare_prior_posterior() 工具函数在 references/advanced_workflows.md 中介绍
2. model_comparison.py（387 行）-> compare_models() 和 check_loo_reliability() 模式内联至工作流步骤 7（直接使用 az.compare 和 az.loo）；model_averaging() 模式在 references/advanced_workflows.md 中介绍；cross_validation_comparison() 指南在 references/advanced_workflows.md 中介绍；plot_model_comparison() 是 az.plot_compare 的轻量封装（已内联至步骤 7）

相关技能

• arviz（计划中） —— 后验可视化和诊断；PyMC 返回 ArviZ InferenceData 对象
• bambi（计划中） —— PyMC 的公式接口（R 风格的 y ~ x1 + x2 语法）；为标准 GLM 提供更简单的 API
• numpyro（计划中） —— 基于 JAX 的概率编程；通过硬件加速实现更快的 NUTS
• statsmodels-statistical-modeling —— 频率学派回归和 GLM；在不需要贝叶斯框架时使用
• scikit-learn-machine-learning —— ML 分类/回归；在预测准确度比不确定性量化更重要时使用

参考资料

• PyMC 文档 —— 官方 API 参考和教程
• PyMC GitHub —— 源代码和问题跟踪器
• ArviZ 文档 —— 后验分析和可视化
• 《Python 中的贝叶斯建模与计算》 —— 使用 PyMC 的教科书
• PyMC 示例集 —— 精选笔记本合集