当前位置：首页 > news >正文

概率编程实战：使用Pyro/PyMC3构建贝叶斯模型

news 2025/9/24 10:33:38

点击 “AladdinEdu，同学们用得起的【H卡】算力平台”，注册即送-H卡级别算力，沉浸式云原生的集成开发环境，80G大显存多卡并行，按量弹性计费，教育用户更享超低价。

引言：从数学推导到声明式建模的革命

传统贝叶斯建模面临着一个核心矛盾：贝叶斯统计提供了强大的理论框架来处理不确定性，但实际实现却异常复杂。研究人员需要掌握复杂的采样算法、耐心调试MCMC链的收敛性、手动计算边缘概率——这些技术门槛让许多实践者对贝叶斯方法望而却步。

概率编程语言（PPL）的出现彻底改变了这一局面。正如TensorFlow和PyTorch降低了深度学习的门槛一样，Pyro和PyMC3等概率编程库让贝叶斯建模变得直观、可扩展、工业化。通过声明式的建模方式，研究者可以专注于模型本身而非计算细节，将复杂的推断过程交给库自动处理。

本文将带你全面掌握使用Pyro和PyMC3进行概率编程的实战技能，从基础模型到高级应用，展示如何快速构建复杂的贝叶斯模型并进行高效的后验推断。

第一章：概率编程基础与库的选择

1.1 什么是概率编程？

概率编程是一种编程范式，它允许用户用程序代码的形式表示概率模型，并自动进行贝叶斯推断。核心思想是：

模型定义：用编程语言描述生成过程（先验和似然）
自动推断：系统自动计算后验分布
灵活扩展：轻松构建复杂的层次模型

1.2 Pyro vs PyMC3：如何选择？

两个主流库各有特色：

特性	Pyro	PyMC3
后端框架	PyTorch	Theano/Aesara
编程风格	命令式、面向对象	声明式、上下文管理器
推断方法	变分推断、HMC、NUTS	NUTS、Metropolis、变分推断
深度学习集成	与PyTorch无缝集成	通过TensorFlow Probability间接集成
学习曲线	较陡峭，需要PyTorch基础	相对平缓，API直观
适用场景	复杂模型、贝叶斯深度学习	传统统计模型、层次模型

选择建议：

如果你是PyTorch用户或需要贝叶斯深度学习，选择Pyro
如果你是统计学背景或需要传统贝叶斯模型，选择PyMC3
初学者建议从PyMC3开始，API更加友好

1.3 环境安装与配置

# 安装PyMC3
pip install pymc3 arviz pandas matplotlib seaborn# 安装Pyro
pip install pyro-ppl torch# 安装可选依赖（用于更快的计算）
pip install pymc3[advanced]  # 包含Aesara等优化

第二章：PyMC3入门实战

2.1 第一个贝叶斯模型：简单线性回归

让我们从一个经典的线性回归开始，了解PyMC3的基本工作流程：

import pymc3 as pm
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import arviz as az# 生成模拟数据
np.random.seed(42)
n_samples = 100
true_alpha = 2.0  # 截距
true_beta = 3.5   # 斜率
true_sigma = 1.0  # 噪声标准差x = np.linspace(0, 1, n_samples)
y_true = true_alpha + true_beta * x
y_obs = y_true + np.random.normal(0, true_sigma, n_samples)# 使用PyMC3构建模型
with pm.Model() as linear_model:# 先验分布alpha = pm.Normal('alpha', mu=0, sigma=10)  # 截距先验beta = pm.Normal('beta', mu=0, sigma=10)    # 斜率先验sigma = pm.HalfNormal('sigma', sigma=1)     # 噪声标准差先验（必须为正）# 期望值（确定性变量）mu = alpha + beta * x# 似然函数likelihood = pm.Normal('y', mu=mu, sigma=sigma, observed=y_obs)# 采样推断trace = pm.sample(2000, tune=1000, cores=2, random_seed=42)print("采样完成！")
print(trace)

2.2 结果分析与可视化

# 后验分析
with linear_model:# 后验摘要summary = pm.summary(trace)print("后验统计摘要:")print(summary)# 轨迹图pm.traceplot(trace)plt.tight_layout()plt.show()# 后验分布图pm.plot_posterior(trace, var_names=['alpha', 'beta', 'sigma'])plt.tight_layout()plt.show()# 模型预测与验证
with linear_model:# 后验预测检查ppc = pm.sample_posterior_predictive(trace, samples=1000)plt.figure(figsize=(10, 6))
plt.scatter(x, y_obs, alpha=0.5, label='观测数据')
plt.plot(x, y_true, 'r-', linewidth=2, label='真实关系')# 绘制后验预测区间
y_pred_mean = ppc['y'].mean(axis=0)
y_pred_std = ppc['y'].std(axis=0)
plt.plot(x, y_pred_mean, 'g-', label='后验预测均值')
plt.fill_between(x, y_pred_mean - 2*y_pred_std, y_pred_mean + 2*y_pred_std, alpha=0.3, color='green', label='95%预测区间')plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.title('线性回归后验预测')
plt.grid(True, alpha=0.3)
plt.show()

2.3 模型比较与诊断

# 模型比较
with linear_model:# 计算WAIC和LOOcompare_dict = {linear_model: trace}comparison = pm.compare(compare_dict, ic='waic')print("模型比较（WAIC）:")print(comparison)# 收敛诊断print("\nRhat统计量（应接近1.0）:")rhat = pm.rhat(trace)print(rhat)# 有效样本数print("\n有效样本数:")ess = pm.ess(trace)print(ess)# 自相关分析
az.plot_autocorr(trace, var_names=['alpha', 'beta', 'sigma'])
plt.tight_layout()
plt.show()

第三章：高级PyMC3模型

3.1 层次贝叶斯模型：八所学校问题

八所学校问题是展示层次模型威力的经典案例：

# 八所学校数据（教育干预效果评估）
schools_data = {'J': 8,  # 学校数量'y': [28, 8, -3, 7, -1, 1, 18, 12],  # 干预效果估计'sigma': [15, 10, 16, 11, 9, 11, 10, 18]  # 估计标准误
}with pm.Model() as schools_model:# 超先验mu = pm.Normal('mu', mu=0, sigma=10)  # 总体均值tau = pm.HalfCauchy('tau', beta=5)    # 总体标准差# 学校特定效应（非中心参数化，提高采样效率）theta_tilde = pm.Normal('theta_tilde', mu=0, sigma=1, shape=schools_data['J'])theta = pm.Deterministic('theta', mu + tau * theta_tilde)# 似然函数y_obs = pm.Normal('y_obs', mu=theta, sigma=schools_data['sigma'], observed=schools_data['y'])# 采样trace_schools = pm.sample(3000, tune=1000, target_accept=0.95, random_seed=42)# 分析结果
print("八所学校模型后验摘要:")
print(pm.summary(trace_schools, var_names=['mu', 'tau', 'theta']))# 可视化学校效应收缩
theta_samples = trace_schools['theta']
raw_estimates = schools_data['y']plt.figure(figsize=(12, 6))
for j in range(schools_data['J']):# 原始估计plt.scatter([j-0.1], [raw_estimates[j]], color='red', s=80, label='原始估计' if j==0 else "")# 后验估计plt.scatter([j+0.1], [theta_samples[:, j].mean()], color='blue', s=80, label='后验估计' if j==0 else "")# 连接线plt.plot([j-0.1, j+0.1], [raw_estimates[j], theta_samples[:, j].mean()], 'k-', alpha=0.5)plt.axhline(y=trace_schools['mu'].mean(), color='green', linestyle='--', label='总体均值')
plt.xlabel('学校编号')
plt.ylabel('干预效果')
plt.title('层次模型的收缩效应')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

3.2 时间序列模型：贝叶斯结构时间序列

# 生成时间序列数据
np.random.seed(42)
n_timepoints = 100
time = np.arange(n_timepoints)# 真实成分：趋势 + 季节性 + 噪声
trend = 0.05 * time
seasonal = 2 * np.sin(2 * np.pi * time / 12)
noise = np.random.normal(0, 0.5, n_timepoints)
y_ts = trend + seasonal + noisewith pm.Model() as time_series_model:# 趋势成分（随机游走）sigma_trend = pm.HalfNormal('sigma_trend', sigma=1)trend_component = pm.GaussianRandomWalk('trend', sigma=sigma_trend, shape=n_timepoints)# 季节性成分（傅里叶级数）n_harmonics = 3for i in range(1, n_harmonics + 1):# 每个谐波的振幅beta_sin = pm.Normal(f'beta_sin_{i}', mu=0, sigma=1)beta_cos = pm.Normal(f'beta_cos_{i}', mu=0, sigma=1)harmonic = (beta_sin * np.sin(2 * np.pi * i * time / 12) + beta_cos * np.cos(2 * np.pi * i * time / 12))if i == 1:seasonal_component = harmonicelse:seasonal_component += harmonic# 观测噪声sigma_obs = pm.HalfNormal('sigma_obs', sigma=1)# 期望值mu = trend_component + seasonal_component# 似然函数likelihood = pm.Normal('y', mu=mu, sigma=sigma_obs, observed=y_ts)# 采样（使用NUTS）trace_ts = pm.sample(2000, tune=1000, target_accept=0.9, random_seed=42)# 时间序列分解可视化
trend_mean = trace_ts['trend'].mean(axis=0)
seasonal_mean = np.zeros(n_timepoints)
for i in range(1, n_harmonics + 1):beta_sin_mean = trace_ts[f'beta_sin_{i}'].mean()beta_cos_mean = trace_ts[f'beta_cos_{i}'].mean()seasonal_mean += (beta_sin_mean * np.sin(2 * np.pi * i * time / 12) + beta_cos_mean * np.cos(2 * np.pi * i * time / 12))plt.figure(figsize=(15, 10))
plt.subplot(3, 1, 1)
plt.plot(time, y_ts, 'o-', alpha=0.7, label='观测数据')
plt.plot(time, trend_mean + seasonal_mean, 'r-', label='模型拟合')
plt.legend()
plt.title('时间序列拟合')plt.subplot(3, 1, 2)
plt.plot(time, trend_mean, 'g-', label='估计趋势')
plt.plot(time, trend, 'g--', alpha=0.5, label='真实趋势')
plt.legend()
plt.title('趋势成分')plt.subplot(3, 1, 3)
plt.plot(time, seasonal_mean, 'b-', label='估计季节性')
plt.plot(time, seasonal, 'b--', alpha=0.5, label='真实季节性')
plt.legend()
plt.title('季节性成分')plt.tight_layout()
plt.show()

第四章：Pyro概率编程实战

4.1 Pyro基础：重新实现线性回归

import torch
import torch.distributions as dist
import pyro
import pyro.distributions as pyro_dist
from pyro.infer import MCMC, NUTS, Predictive
from pyro.infer import SVI, Trace_ELBO
from pyro.optim import Adam# 设置随机种子
pyro.set_rng_seed(42)
torch.manual_seed(42)# 准备数据（使用PyTorch张量）
x_tensor = torch.tensor(x, dtype=torch.float32)
y_tensor = torch.tensor(y_obs, dtype=torch.float32)def linear_model_pyro(x, y=None):"""Pyro线性回归模型"""# 先验分布alpha = pyro.sample("alpha", pyro_dist.Normal(0., 10.))beta = pyro.sample("beta", pyro_dist.Normal(0., 10.))sigma = pyro.sample("sigma", pyro_dist.HalfNormal(1.))# 线性预测mu = alpha + beta * x# 观测模型with pyro.plate("data", len(x)):pyro.sample("obs", pyro_dist.Normal(mu, sigma), obs=y)return mu# 使用NUTS进行MCMC采样
nuts_kernel = NUTS(linear_model_pyro)
mcmc = MCMC(nuts_kernel, num_samples=2000, warmup_steps=1000)
mcmc.run(x_tensor, y_tensor)# 获取后验样本
posterior_samples = mcmc.get_samples()
print("Pyro MCMC采样完成!")
print(f"后验样本形状: alpha {posterior_samples['alpha'].shape}")# 后验分析
def analyze_pyro_posterior(samples):"""分析Pyro后验样本"""for param_name, samples in samples.items():mean_val = samples.mean().item()std_val = samples.std().item()hdi_95 = torch.quantile(samples, torch.tensor([0.025, 0.975]))print(f"{param_name}: 均值={mean_val:.3f}, 标准差={std_val:.3f}, "f"95% HDI=[{hdi_95[0]:.3f}, {hdi_95[1]:.3f}]")analyze_pyro_posterior(posterior_samples)# 后验预测
predictive = Predictive(linear_model_pyro, posterior_samples, num_samples=1000)
predictions = predictive(x_tensor, None)
y_pred_pyro = predictions["obs"]# 可视化比较Pyro和PyMC3结果
plt.figure(figsize=(12, 5))plt.subplot(1, 2, 1)
plt.scatter(x, y_obs, alpha=0.5, label='观测数据')
plt.plot(x, y_true, 'r-', label='真实关系')# Pyro预测
y_pred_mean_pyro = y_pred_pyro.mean(dim=0).numpy()
y_pred_std_pyro = y_pred_pyro.std(dim=0).numpy()
plt.plot(x, y_pred_mean_pyro, 'g-', label='Pyro预测')
plt.fill_between(x, y_pred_mean_pyro - 2*y_pred_std_pyro, y_pred_mean_pyro + 2*y_pred_std_pyro, alpha=0.3, color='green')plt.title('Pyro线性回归')
plt.legend()plt.subplot(1, 2, 2)
# 比较参数后验分布
import seaborn as sns
sns.kdeplot(posterior_samples['alpha'].numpy(), label='Pyro alpha', fill=True)
sns.kdeplot(trace['alpha'], label='PyMC3 alpha', fill=True)
plt.title('参数后验分布比较')
plt.legend()plt.tight_layout()
plt.show()

4.2 变分推断：更快的近似贝叶斯方法

# 使用随机变分推断（SVI）进行快速近似
def guide_linear(x, y=None):"""变分分布（指导函数）"""# 变分参数alpha_loc = pyro.param("alpha_loc", torch.tensor(0.0))alpha_scale = pyro.param("alpha_scale", torch.tensor(1.0), constraint=pyro.distributions.constraints.positive)beta_loc = pyro.param("beta_loc", torch.tensor(0.0))beta_scale = pyro.param("beta_scale", torch.tensor(1.0),constraint=pyro.distributions.constraints.positive)sigma_loc = pyro.param("sigma_loc", torch.tensor(1.0),constraint=pyro.distributions.constraints.positive)# 变分分布pyro.sample("alpha", pyro_dist.Normal(alpha_loc, alpha_scale))pyro.sample("beta", pyro_dist.Normal(beta_loc, beta_scale))pyro.sample("sigma", pyro_dist.LogNormal(sigma_loc, 0.1))# 变分推断
pyro.clear_param_store()
svi = SVI(linear_model_pyro, guide_linear, Adam({"lr": 0.01}), Trace_ELBO())# 训练变分分布
num_iterations = 2000
for i in range(num_iterations):elbo = svi.step(x_tensor, y_tensor)if i % 500 == 0:print(f"Iteration {i}: ELBO = {elbo:.1f}")# 获取变分参数
print("\n变分参数:")
for name, value in pyro.get_param_store().items():print(f"{name}: {value.data}")# 比较MCMC和VI结果
vi_predictive = Predictive(linear_model_pyro, guide=guide_linear, num_samples=1000)
vi_predictions = vi_predictive(x_tensor, None)
y_pred_vi = vi_predictions["obs"]plt.figure(figsize=(10, 6))
plt.scatter(x, y_obs, alpha=0.5, label='观测数据')# MCMC预测
plt.plot(x, y_pred_mean_pyro, 'r-', label='MCMC预测', linewidth=2)# VI预测
y_pred_mean_vi = y_pred_vi.mean(dim=0).numpy()
y_pred_std_vi = y_pred_vi.std(dim=0).numpy()
plt.plot(x, y_pred_mean_vi, 'b--', label='VI预测', linewidth=2)
plt.fill_between(x, y_pred_mean_vi - 2*y_pred_std_vi, y_pred_mean_vi + 2*y_pred_std_vi, alpha=0.2, color='blue')plt.title('MCMC vs 变分推断预测比较')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

第五章：贝叶斯深度学习实战

5.1 贝叶斯神经网络

import torch.nn as nnclass BayesianNN(nn.Module):"""简单的贝叶斯神经网络"""def __init__(self, input_dim, hidden_dim, output_dim):super().__init__()self.input_dim = input_dimself.hidden_dim = hidden_dimself.output_dim = output_dim# 网络层self.fc1 = nn.Linear(input_dim, hidden_dim)self.fc2 = nn.Linear(hidden_dim, output_dim)self.relu = nn.ReLU()def forward(self, x):x = self.relu(self.fc1(x))x = self.fc2(x)return xdef model_bayesian_nn(x, y=None):"""贝叶斯神经网络模型"""# 创建网络实例net = BayesianNN(input_dim=1, hidden_dim=10, output_dim=1)# 为权重和偏置设置先验for name, param in net.named_parameters():if 'weight' in name:pyro.sample(f"{name}_prior", pyro_dist.Normal(0, 1).expand(param.shape).to_event(param.dim()))elif 'bias' in name:pyro.sample(f"{name}_prior", pyro_dist.Normal(0, 1).expand(param.shape).to_event(param.dim()))# 观测噪声sigma = pyro.sample("sigma", pyro_dist.HalfNormal(1.))# 前向传播mu = net(x.unsqueeze(1)).squeeze()# 观测模型with pyro.plate("data", len(x)):pyro.sample("obs", pyro_dist.Normal(mu, sigma), obs=y)return mudef guide_bayesian_nn(x, y=None):"""贝叶斯神经网络的变分指导函数"""net = BayesianNN(input_dim=1, hidden_dim=10, output_dim=1)# 为每个参数定义变分分布for name, param in net.named_parameters():param_shape = param.shapeparam_dim = param.dim()# 变分参数loc = pyro.param(f"{name}_loc", torch.zeros(param_shape))scale = pyro.param(f"{name}_scale", torch.ones(param_shape) * 0.1,constraint=pyro.distributions.constraints.positive)pyro.sample(f"{name}_prior", pyro_dist.Normal(loc, scale).to_event(param_dim))# 噪声参数的变分分布sigma_loc = pyro.param("sigma_loc", torch.tensor(0.0))sigma_scale = pyro.param("sigma_scale", torch.tensor(0.1),constraint=pyro.distributions.constraints.positive)pyro.sample("sigma", pyro_dist.LogNormal(sigma_loc, sigma_scale))# 生成非线性数据
x_nn = torch.linspace(-3, 3, 100)
y_true_nn = torch.sin(x_nn) + 0.3 * x_nn
y_obs_nn = y_true_nn + 0.2 * torch.randn_like(x_nn)# 训练贝叶斯神经网络
pyro.clear_param_store()
svi_nn = SVI(model_bayesian_nn, guide_bayesian_nn, Adam({"lr": 0.01}), Trace_ELBO())# 训练循环
print("训练贝叶斯神经网络...")
for i in range(3000):elbo = svi_nn.step(x_nn, y_obs_nn)if i % 1000 == 0:print(f"Iteration {i}: ELBO = {elbo:.1f}")# 预测
predictive_nn = Predictive(model_bayesian_nn, guide=guide_bayesian_nn, num_samples=1000)
predictions_nn = predictive_nn(x_nn, None)
y_pred_nn = predictions_nn["obs"]# 可视化贝叶斯神经网络结果
plt.figure(figsize=(12, 5))plt.subplot(1, 2, 1)
plt.scatter(x_nn, y_obs_nn, alpha=0.5, label='观测数据')
plt.plot(x_nn, y_true_nn, 'r-', label='真实函数', linewidth=2)y_pred_mean_nn = y_pred_nn.mean(dim=0).detach().numpy()
y_pred_std_nn = y_pred_nn.std(dim=0).detach().numpy()
plt.plot(x_nn, y_pred_mean_nn, 'g-', label='贝叶斯NN预测', linewidth=2)
plt.fill_between(x_nn.detach().numpy(),y_pred_mean_nn - 2*y_pred_std_nn,y_pred_mean_nn + 2*y_pred_std_nn,alpha=0.3, color='green')plt.title('贝叶斯神经网络拟合')
plt.legend()plt.subplot(1, 2, 2)
# 显示不确定性如何随输入变化
uncertainty = y_pred_std_nn
plt.plot(x_nn, uncertainty, 'b-', linewidth=2)
plt.xlabel('x')
plt.ylabel('预测标准差')
plt.title('预测不确定性')
plt.grid(True, alpha=0.3)plt.tight_layout()
plt.show()

第六章：实际应用案例

6.1 A/B测试的贝叶斯方法

# A/B测试数据：两种网页设计的转化率
ab_data = {'design_a': {'visitors': 1000, 'conversions': 120},  # 设计A：120/1000'design_b': {'visitors': 1050, 'conversions': 150}   # 设计B：150/1050
}with pm.Model() as ab_test_model:# 先验分布（Beta分布是二项分布的共轭先验）theta_a = pm.Beta('theta_a', alpha=2, beta=20)  # 预期转化率约10%theta_b = pm.Beta('theta_b', alpha=2, beta=20)# 似然函数conv_a = pm.Binomial('conv_a', n=ab_data['design_a']['visitors'], p=theta_a, observed=ab_data['design_a']['conversions'])conv_b = pm.Binomial('conv_b', n=ab_data['design_b']['visitors'], p=theta_b, observed=ab_data['design_b']['conversions'])# 计算差异diff = pm.Deterministic('diff', theta_b - theta_a)relative_improvement = pm.Deterministic('relative_improvement', (theta_b - theta_a) / theta_a)# 采样trace_ab = pm.sample(2000, tune=1000, random_seed=42)# A/B测试结果分析
print("A/B测试后验分析:")
print(pm.summary(trace_ab))# 计算设计B优于设计A的概率
prob_b_better = (trace_ab['diff'] > 0).mean()
print(f"\n设计B优于设计A的概率: {prob_b_better:.3f}")# 可视化后验分布
fig, axes = plt.subplots(2, 2, figsize=(12, 10))# 转化率后验
az.plot_posterior(trace_ab, var_names=['theta_a'], ax=axes[0,0])
axes[0,0].set_title('设计A转化率后验')az.plot_posterior(trace_ab, var_names=['theta_b'], ax=axes[0,1])
axes[0,1].set_title('设计B转化率后验')# 差异后验
az.plot_posterior(trace_ab, var_names=['diff'], ax=axes[1,0])
axes[1,0].set_title('转化率差异后验 (B - A)')# 相对改进后验
az.plot_posterior(trace_ab, var_names=['relative_improvement'], ax=axes[1,1])
axes[1,1].set_title('相对改进后验')plt.tight_layout()
plt.show()

6.2 生存分析：贝叶斯Cox比例风险模型

# 生存分析数据（模拟）
np.random.seed(42)
n_patients = 200
treatment = np.random.binomial(1, 0.5, n_patients)  # 治疗组 vs 对照组
age = np.random.normal(60, 10, n_patients)          # 年龄# 生成生存时间（Weibull分布）
true_beta_treatment = -0.5  # 治疗降低风险
true_beta_age = 0.02        # 年龄增加风险linear_predictor = true_beta_treatment * treatment + true_beta_age * (age - 60)
scale = 100  # 尺度参数
shape = 1.5  # 形状参数# 生存时间
survival_time = scale * np.random.weibull(shape, n_patients) * np.exp(-linear_predictor)# 随机删失
censoring_time = np.random.exponential(150, n_patients)
observed_time = np.minimum(survival_time, censoring_time)
event_observed = (survival_time <= censoring_time).astype(int)with pm.Model() as cox_model:# 先验分布beta_treatment = pm.Normal('beta_treatment', mu=0, sigma=1)beta_age = pm.Normal('beta_age', mu=0, sigma=0.1)# 基线风险参数alpha = pm.Normal('alpha', mu=0, sigma=1)  # 对数尺度参数rho = pm.HalfNormal('rho', sigma=1)        # 形状参数# 线性预测项lin_pred = beta_treatment * treatment + beta_age * (age - 60)# Weibull生存函数的对数似然def weibull_log_likelihood(observed_time, event, alpha, rho, lin_pred):scale = pm.math.exp(alpha - lin_pred)log_p = (event * (pm.math.log(rho) - rho * pm.math.log(scale) + (rho - 1) * pm.math.log(observed_time)) - (observed_time / scale) ** rho)return log_p# 自定义似然函数survival_likelihood = pm.Potential('survival_likelihood', weibull_log_likelihood(observed_time, event_observed, alpha, rho, lin_pred))# 采样trace_cox = pm.sample(2000, tune=1000, target_accept=0.9, random_seed=42)print("Cox模型后验分析:")
print(pm.summary(trace_cox, var_names=['beta_treatment', 'beta_age']))# 风险比分析
hr_treatment = np.exp(trace_cox['beta_treatment'])
hr_age = np.exp(trace_cox['beta_age'])plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
az.plot_posterior(hr_treatment, ref_val=1)
plt.title('治疗风险比后验 (HR < 1 表示有益)')plt.subplot(1, 2, 2)
az.plot_posterior(hr_age, ref_val=1)
plt.title('年龄风险比后验 (每增加1岁)')plt.tight_layout()
plt.show()

第七章：性能优化与最佳实践

7.1 采样效率优化

# 1. 非中心参数化（提高层次模型采样效率）
with pm.Model() as non_centered_model:mu = pm.Normal('mu', 0, 1)sigma = pm.HalfNormal('sigma', 1)# 非中心参数化offset = pm.Normal('offset', 0, 1, shape=100)theta = pm.Deterministic('theta', mu + offset * sigma)# 比较采样效率trace_nc = pm.sample(1000, tune=1000, random_seed=42)# 2. 使用ADVI快速近似
with linear_model:# 自动微分变分推断approx = pm.fit(n=50000, method='advi')trace_advi = approx.sample(1000)print("ADVI与NUTS比较:")
print("NUTS有效样本数:", pm.ess(trace).to_array().mean())
print("ADVI有效样本数:", pm.ess(trace_advi).to_array().mean())# 3. 并行采样
with linear_model:# 使用多核心并行采样trace_parallel = pm.sample(2000, tune=1000, cores=4, random_seed=42)

7.2 模型诊断与验证

def comprehensive_diagnostics(model, trace):"""综合模型诊断"""# 1. 收敛诊断max_rhat = pm.rhat(trace).max()min_ess = pm.ess(trace).min()print(f"收敛诊断: Rhat最大值 = {max_rhat:.3f} (应<1.1), ESS最小值 = {min_ess:.0f}")# 2. 后验预测检查with model:ppc = pm.sample_posterior_predictive(trace, samples=500)# 3. 留一法交叉验证（近似）with model:loo_result = pm.loo(trace, pointwise=True)print(f"LOO评估: LOO值 = {loo_result.loo:.1f}, 标准误 = {loo_result.loo_se:.1f}")# 4. 残差分析y_pred = ppc['y'].mean(axis=0)residuals = y_obs - y_predfig, axes = plt.subplots(2, 2, figsize=(12, 10))# 残差分布axes[0,0].hist(residuals, bins=30, alpha=0.7)axes[0,0].set_title('残差分布')axes[0,0].axvline(0, color='red', linestyle='--')# Q-Q图import scipy.stats as statsstats.probplot(residuals, dist="norm", plot=axes[0,1])axes[0,1].set_title('Q-Q图')# 残差vs拟合值axes[1,0].scatter(y_pred, residuals, alpha=0.5)axes[1,0].axhline(0, color='red', linestyle='--')axes[1,0].set_xlabel('预测值')axes[1,0].set_ylabel('残差')axes[1,0].set_title('残差vs拟合值')# 自相关az.plot_autocorr(trace, ax=axes[1,1])axes[1,1].set_title('参数自相关')plt.tight_layout()plt.show()return ppc, loo_result# 运行综合诊断
ppc_diag, loo_diag = comprehensive_diagnostics(linear_model, trace)

第八章：总结与展望

8.1 概率编程的核心优势

建模灵活性：可以轻松构建复杂的层次模型、时间序列模型等
不确定性量化：自动提供完整的后验分布，而不仅仅是点估计
代码简洁性：声明式语法让模型代码更接近数学公式
自动推断：无需手动实现采样算法，专注于模型设计

8.2 选择指南总结

场景	推荐工具	理由
传统统计模型	PyMC3	API直观，文档丰富，社区成熟
贝叶斯深度学习	Pyro	与PyTorch集成，支持复杂神经网络
生产环境	PyMC3	稳定性好，推断速度快
研究原型	Pyro	灵活性高，易于扩展新算法
初学者	PyMC3	学习曲线平缓，示例丰富