以下是使用贝叶斯框架和arviz库绘制先验预测和后验预测分布的示例代码:
import numpy as np
import pymc3 as pm
import arviz as az
import matplotlib.pyplot as plt
# 生成数据
np.random.seed(0)
n = 100
x = np.linspace(-5, 5, n)
y = 2 * x + np.random.normal(0, 1, n)
# 定义贝叶斯模型
with pm.Model() as model:
# 设置先验分布
alpha = pm.Normal('alpha', mu=0, sd=10)
beta = pm.Normal('beta', mu=0, sd=10)
sigma = pm.HalfCauchy('sigma', 5)
# 定义线性回归模型
mu = alpha + beta * x
# 定义似然函数
likelihood = pm.Normal('y', mu=mu, sd=sigma, observed=y)
# 采样参数后验分布
trace = pm.sample(1000, tune=1000)
# 绘制先验预测分布
with model:
prior = pm.sample_prior_predictive(1000)
plt.figure(figsize=(10, 6))
for i in range(100):
plt.plot(x, prior['alpha'][i] + prior['beta'][i] * x, color='lightgray', alpha=0.1)
plt.plot(x, y, 'ko', label='Observations')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Prior Predictive Distribution')
plt.legend()
plt.show()
# 绘制后验预测分布
with model:
posterior = pm.sample_posterior_predictive(trace)
plt.figure(figsize=(10, 6))
for i in range(100):
plt.plot(x, posterior['y'][i], color='lightgray', alpha=0.1)
plt.plot(x, y, 'ko', label='Observations')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Posterior Predictive Distribution')
plt.legend()
plt.show()
# 使用arviz绘制先验预测和后验预测分布
az.plot_ppc(az.from_pymc3(posterior), mean=False, figsize=(10, 6))
plt.plot(x, y, 'ko', label='Observations')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Posterior Predictive Distribution (arviz)')
plt.legend()
plt.show()
这段代码首先使用贝叶斯框架定义了一个线性回归模型,然后使用arviz库绘制了先验预测分布和后验预测分布。先验预测分布表示在观测数据之前的模型预测,后验预测分布表示在观测数据之后的模型预测。通过比较观测数据和预测数据,可以对模型的准确性进行评估。
在代码的最后,使用arviz的plot_ppc
函数绘制了后验预测分布的直方图,并将观测数据的散点图叠加在上面,以便直观地比较模型预测和观测数据。
要运行这段代码,需要安装pymc3、arviz和matplotlib库。