技能简介
贝叶斯统计提供了一种直观的参数不确定性量化方法，相比传统频率学派方法， 它能够直接给出参数的概率分布解释，并允许将先验知识纳入模型。
本教程将介绍如何使用 PyMC 进行概率编程， 以及 statsmodels 进行传统统计建模， 并对比两种方法的结果差异。
📊 贝叶斯 vs 频率学派对比
| 方面 | 贝叶斯 | 频率学派 |
|---|---|---|
| 参数解释 | 概率分布 | 点估计 |
| 不确定性 | 可信区间 | 置信区间 |
| 先验信息 | ✓ 可以纳入 | ✗ 不考虑 |
| 小样本 | 表现更好 | 可能不稳定 |
| 计算方式 | MCMC采样 | 解析解/优化 |
💡 使用场景
小样本临床研究
当样本量较小时，贝叶斯方法可以通过先验信息提高估计稳定性
不确定性量化
需要给出预测的完整概率分布，而非单一置信区间
临床试验设计
适应性试验设计需要实时更新后验概率
复杂层次模型
多中心研究、纵向数据等需要层次结构的数据
🛠️ Step 1: 模型定义
首先使用PyMC定义一个贝叶斯逻辑回归模型，用于预测心衰患者的30天死亡率。 我们使用弱信息先验（weakly informative priors）避免对结果产生过大影响。
import pymc as pm import numpy as np import pandas as pd # 准备数据 # X: 特征矩阵 (n_samples, n_features) # y: 二分类结局 (0/1) with pm.Model() as bayesian_model: # 1. 定义先验分布 # 使用正态分布作为系数的弱信息先验 intercept = pm.Normal('intercept', mu=0, sigma=10) beta = pm.Normal('beta', mu=0, sigma=10, shape=X.shape[1]) # 2. 定义线性预测器 logit_p = intercept + pm.math.dot(X, beta) # 3. 定义似然函数 # 使用Bernoulli分布（二分类） y_obs = pm.Bernoulli('y_obs', logit_p=logit_p, observed=y) # 4. 查看模型结构 pm.model_to_graphviz(bayesian_model)
💡 先验选择建议：
- • 弱信息先验：N(0, 10) 或 N(0, 5)，影响小但正则化
- • 无信息先验：N(0, 100) 或更大，让数据说话
- • 正则化先验：如Laplace(0, 1)，类似LASSO
- • 避免使用过于强的先验，除非有充分的先验知识
🛠️ Step 2: MCMC采样
使用马尔可夫链蒙特卡洛（MCMC）方法从后验分布中采样。 PyMC默认使用NUTS（No-U-Turn Sampler）算法，这是HMC的高效变体。
with bayesian_model: # 1. MCMC采样 trace = pm.sample( draws=2000, # 每条链的采样次数 tune=1000, # 调谐步数（burn-in） chains=4, # 链的数量 cores=4, # 并行计算 return_inferencedata=True ) # 2. 检查收敛性 # R-hat < 1.05 表示收敛良好 print(pm.summary(trace, round_to=2)) # 3. 可视化采样轨迹 pm.plot_trace(trace) pm.plot_posterior(trace)
📌 收敛性诊断指标：
- • R-hat：< 1.05 表示收敛，< 1.01 更佳
- • ESS：有效样本量，越大越好（>400通常可接受）
- • 轨迹图：应像"毛毛虫"，无规律波动
- • 自相关图：快速下降到0表示采样效率高
🛠️ Step 3: 后验分析
后验分布包含了参数的所有信息。我们可以计算各种统计量，并进行后验预测检验。
import arviz as az # 1. 提取后验摘要 summary = pm.summary(trace) print(summary) # 输出示例: # mean sd hdi_3% hdi_97% ... # intercept -2.1 0.5 -3.0 -1.2 # beta[0] 0.8 0.3 0.3 1.3 # 2. 计算95%可信区间 (HDI) # HDI: Highest Density Interval hdi = pm.hdi(trace, hdi_prob=0.95) # 3. 后验预测检验 with bayesian_model: # 生成后验预测样本 ppc = pm.sample_posterior_predictive(trace) # 比较观测数据和预测数据 az.plot_ppc(az.from_pymc3(posterior_predictive=ppc, model=bayesian_model))
⚠️ 可信区间 vs 置信区间：
- • 可信区间 (Credible Interval)：参数有95%概率落在该区间内（贝叶斯）
- • 置信区间 (Confidence Interval)：重复抽样95%的区间会包含真值（频率学派）
- • 贝叶斯的解释更直观，符合大多数研究者的直觉
🛠️ Step 4: 与频率学派模型比较
使用statsmodels拟合传统的逻辑回归模型，并与贝叶斯结果进行比较。
import statsmodels.api as sm import statsmodels.formula.api as smf # 1. 使用statsmodels拟合逻辑回归 X_sm = sm.add_constant(X) # 添加截距项 model_freq = sm.Logit(y, X_sm) result_freq = model_freq.fit() # 2. 查看结果 print(result_freq.summary()) # 3. 比较系数估计 comparison = pd.DataFrame({ '贝叶斯后验均值': summary['mean'], '频率学派估计': result_freq.params, '差异': summary['mean'] - result_freq.params }) print(comparison)
💡 结果解读：
- • 当样本量较大时，两种方法的点估计通常很接近
- • 贝叶斯方法提供完整的后验分布，而非单一估计值
- • 小样本时，先验信息会使贝叶斯估计更稳定
- • 贝叶斯方法可以自然处理复杂模型结构
🛠️ Step 5: 预测与决策支持
使用拟合好的贝叶斯模型进行新样本预测，并量化预测的不确定性。
with bayesian_model: # 1. 对新数据进行预测 X_new = ... # 新样本特征 # 生成后验预测分布 ppc_new = pm.sample_posterior_predictive( trace, predictions=True ) # 2. 计算预测概率 pred_prob = ppc_new['y_obs'].mean(axis=0) # 3. 预测不确定性 pred_std = ppc_new['y_obs'].std(axis=0) # 4. 风险分层示例 risk_categories = [] for prob, uncertainty in zip(pred_prob, pred_std): if prob > 0.7: risk_categories.append('高风险') elif prob > 0.3: risk_categories.append('中风险') else: risk_categories.append('低风险')
⚠️ 注意事项
先验选择很重要
过强的先验可能主导结果，过弱的先验可能导致收敛问题。建议使用弱信息先验作为起点。
检查模型收敛
不要忽略R-hat和ESS等诊断指标。未收敛的模型结果不可信。
计算成本较高
贝叶斯分析通常比频率学派方法需要更多计算资源，特别是复杂模型。
模型验证必不可少
使用后验预测检验、交叉验证等方法验证模型预测能力。