Pyro:由Uber开发的Python库,用于深度概率建模。Pyro提供了灵活的API,并且与PyTorch紧密集成,后者是另一个流行的深度学习库。Pyro支持多种推理算法,包括变分推理和马尔可夫链蒙特卡洛(MCMC)。
import pyro
import pyro.distributions as dist
import torch
import numpy as np
import matplotlib.pyplot as plt
from pyro.infer import MCMC, NUTS
plt.rcParams['font.sans-serif'] = ['SimHei'] # 设置中文字体为黑体
plt.rcParams['axes.unicode_minus'] = False # 解决负号显示为方块的
# 定义模型
def model(x, y):
# 先验:权重和噪声的先验分布
w_prior = dist.Normal(torch.tensor(0.), torch.tensor(10.))
b_prior = dist.Normal(torch.tensor(0.), torch.tensor(10.))
sigma_prior = dist.Uniform(torch.tensor(0.), torch.tensor(10.))
# 随机样本
w = pyro.sample("weights", w_prior)
b = pyro.sample("bias", b_prior)
sigma = pyro.sample("sigma", sigma_prior)
# 似然函数
with pyro.plate("data", x.shape[0]):
pyro.sample("obs", dist.Normal(w * x + b, sigma), obs=y)
# 示例数据
N = 100 # 数据点数量
x_data = torch.linspace(-8, 8, N)
y_data = torch.tensor(1.5 + x_data + torch.randn(N) * 2., dtype=torch.float)
# 设置MCMC
nuts_kernel = NUTS(model)
mcmc = MCMC(nuts_kernel, num_samples=10000, warmup_steps=200)
mcmc.run(x_data, y_data)
# 获取MCMC样本
mcmc_samples = mcmc.get_samples()
# 打印MCMC样本的统计信息
print(mcmc_samples["weights"].mean())
print(mcmc_samples["bias"].mean())
print(mcmc_samples["sigma"].mean())
# 可视化MCMC样本
w_samples = mcmc_samples["weights"].detach().numpy()
b_samples = mcmc_samples["bias"].detach().numpy()
sigma_samples = mcmc_samples["sigma"].detach().numpy()
plt.subplot(3, 1, 1)
plt.hist(w_samples, bins=50, density=True)
plt.title("权重 w")
plt.subplot(3, 1, 2)
plt.hist(b_samples, bins=50, density=True)
plt.title("偏差 b")
plt.subplot(3, 1, 3)
plt.hist(sigma_samples, bins=50, density=True)
plt.title("噪声尺度 sigma")
plt.tight_layout()
plt.show()
# 绘制线性回归线和0.95 HDI
x_line = torch.linspace(x_data.min(), x_data.max(), 100)
w_mean = mcmc_samples["weights"].mean()
b_mean = mcmc_samples["bias"].mean()
y_line = w_mean * x_line + b_mean
# 生成预测值的样本
pred_samples = []
for i in range(10000):
w = mcmc_samples["weights"][i]
b = mcmc_samples["bias"][i]
sigma = mcmc_samples["sigma"][i]
y_pred = dist.Normal(w * x_line + b, sigma).sample()
pred_samples.append(y_pred)
pred_samples = torch.stack(pred_samples)
# 计算0.95 HDI
pred_lower = torch.quantile(pred_samples, 0.025, dim=0)
pred_upper = torch.quantile(pred_samples, 0.975, dim=0)
# 绘制数据点和线性回归线
plt.scatter(x_data, y_data, label='数据点')
plt.plot(x_line, y_line, label='回归线')
# 绘制0.95 HDI
plt.fill_between(x_line, pred_lower, pred_upper, color='b', alpha=0.2, label='0.95 HDI')
plt.legend()
plt.show()