基于ArviZ绘制贝叶斯先验/后验预测分布的技术求助

我正在复现《贝叶斯计算手册》中图1.5的先验预测分布与后验预测分布示例，当前设置为4条链、1000次抽样、20个观测值，目标是绘制x轴为成功次数、y轴为概率的直方图，但在维度处理上出现问题——尝试沿chain和draw维度求和/均值后，结果与参考示例不符。以下是我使用PyMC和ArviZ编写的代码：

import pymc as pm
import matplotlib.pyplot as plt
import scipy.stats as stats
import arviz as az
import numpy as np

az.style.use("arviz-grayscale")
plt.rcParams['figure.dpi'] = 300
np.random.seed(521)
viridish = [(0.2823529411764706, 0.11372549019607843, 0.43529411764705883, 1.0),
            (0.1450980392156863, 0.6705882352941176, 0.5098039215686274, 1.0),
            (0.6901960784313725, 0.8666666666666667, 0.1843137254901961, 1.0)]

Y = stats.bernoulli(0.7).rvs(20)
# Declare a model in PyMC3
with pm.Model() as model:
    # Specify the prior distribution of the unknown parameter
    θ = pm.Beta("θ", alpha=1, beta=1)

    # Specify the likelihood distribution and condition on the observed data
    y_obs = pm.Binomial("y_obs", n=1, p=θ, observed=Y)

    # Sample from the posterior distribution
    idata = pm.sample(1000, return_inferencedata=True)

pred_dists = (pm.sample_prior_predictive(1000, model=model),
              pm.sample_posterior_predictive(idata,model=model))

# Prior
prior_samples = pred_dists[0]['prior']['θ'].values # Prior observed
prior_pred_samples = pred_dists[0]['prior_predictive'] # Prior predictions

# Posterior
post_distribution = idata.posterior["θ"] # Posterior observed
posterior_pred_samples = pred_dists[1]['posterior_predictive'] # Posterior Predictions

posterior_pred_graph = posterior_pred_samples['y_obs'].sum(dim=['draw','chain'])
prior_pred_samples_graph = prior_pred_samples['y_obs'].sum(dim=['draw','chain'])

fig,axes =plt.subplots(4,1,gridspec_kw={'hspace': 0.1})

az.plot_dist(prior_samples, plot_kwargs={"color":"0.5"},
             fill_kwargs={'alpha':1}, ax=axes[0])
axes[0].set_title("Prior distribution", fontweight='bold',fontsize=10)
axes[0].set_xlim(0, 1)
axes[0].set_ylim(0, 4)
axes[0].tick_params(axis='both', pad=7)
axes[0].set_xlabel("θ")


az.plot_dist(prior_pred_samples_graph, plot_kwargs={"color":"0.5"},
             fill_kwargs={'alpha':1}, ax=axes[1])
axes[1].set_title("Prior predictive distribution", fontweight='bold',fontsize=10)
# axes[1].set_xlim(-1, 21)
# axes[1].set_ylim(0, 0.15)
axes[1].tick_params(axis='both', pad=7)
axes[1].set_xlabel("number of success")

az.plot_dist(post_distribution, plot_kwargs={"color":"0.5"},
             fill_kwargs={'alpha':1},ax=axes[2])
axes[2].set_title("Posterior distribution", fontweight='bold',fontsize=10)
axes[2].set_xlim(0, 1)
axes[2].set_ylim(0, 5)
axes[2].tick_params(axis='both', pad=7)
axes[2].set_xlabel("θ")

az.plot_dist(posterior_pred_graph, plot_kwargs={"color":"0.5"},
             fill_kwargs={'alpha':1}, ax=axes[3])
axes[3].set_title("Posterior predictive distribution", fontweight='bold',fontsize=10)
# axes[3].set_xlim(-1, 21)
# axes[3].set_ylim(0, 0.15)
axes[3].tick_params(axis='both', pad=7)
axes[3].set_xlabel("number of success")

问题根源 #

你错误地对chain和draw维度求和，而正确的操作应该是对每个抽样（每个chain+draw组合）对应的20个观测值求和，得到该抽样下的成功次数。每个抽样对应一组20次伯努利试验，我们需要收集所有抽样的成功次数，再基于这些样本绘制分布。

修改步骤 #

1. 调整预测样本的求和维度：对预测样本的观测值维度求和，得到每个抽样的成功次数，再展平为一维数组用于绘图。
2. 使用概率归一化的直方图：设置直方图的density=True，让y轴表示概率密度（总和为1），与手册示例一致。

修改后的关键代码段 #

替换原代码中处理预测样本的部分：

# 查看观测值维度名称（如果不确定）
# print(prior_pred_samples['y_obs'].dims)

# 处理先验预测样本：每个抽样对应20次试验，求和得到成功次数
prior_pred_success = prior_pred_samples['y_obs'].sum(dim="y_obs_dim_0").values.flatten()
# 处理后验预测样本同理
posterior_pred_success = posterior_pred_samples['y_obs'].sum(dim="y_obs_dim_0").values.flatten()

完整修改后的绘图代码 #

fig,axes =plt.subplots(4,1,gridspec_kw={'hspace': 0.1})

# 先验分布
az.plot_dist(prior_samples, plot_kwargs={"color":"0.5"},
             fill_kwargs={'alpha':1}, ax=axes[0])
axes[0].set_title("Prior distribution", fontweight='bold',fontsize=10)
axes[0].set_xlim(0, 1)
axes[0].set_ylim(0, 4)
axes[0].tick_params(axis='both', pad=7)
axes[0].set_xlabel("θ")

# 先验预测分布（概率直方图）
az.plot_dist(prior_pred_success, kind="hist", hist_kwargs={"density": True, "color":"0.5", "alpha":1},
             ax=axes[1])
axes[1].set_title("Prior predictive distribution", fontweight='bold',fontsize=10)
axes[1].set_xlim(-1, 21)
axes[1].set_ylim(0, 0.15)
axes[1].tick_params(axis='both', pad=7)
axes[1].set_xlabel("number of success")

# 后验分布
az.plot_dist(post_distribution, plot_kwargs={"color":"0.5"},
             fill_kwargs={'alpha':1},ax=axes[2])
axes[2].set_title("Posterior distribution", fontweight='bold',fontsize=10)
axes[2].set_xlim(0, 1)
axes[2].set_ylim(0, 5)
axes[2].tick_params(axis='both', pad=7)
axes[2].set_xlabel("θ")

# 后验预测分布（概率直方图）
az.plot_dist(posterior_pred_success, kind="hist", hist_kwargs={"density": True, "color":"0.5", "alpha":1},
             ax=axes[3])
axes[3].set_title("Posterior predictive distribution", fontweight='bold',fontsize=10)
axes[3].set_xlim(-1, 21)
axes[3].set_ylim(0, 0.15)
axes[3].tick_params(axis='both', pad=7)
axes[3].set_xlabel("number of success")

plt.tight_layout()
plt.show()

补充说明 #

• 最终得到的prior_pred_success和posterior_pred_success是包含4000个样本的一维数组（4链×1000抽样），每个值对应一组20次试验的成功次数。
• 若偏好Matplotlib原生绘图，可将az.plot_dist替换为：

  axes[1].hist(prior_pred_success, bins=np.arange(0,22), density=True, color="0.5", alpha=1)
  

内容的提问来源于stack exchange，提问作者JCV