如何用PyMC的后验预测分布实现模型检验并生成最大值直方图？

1. 提取后验预测样本 #

假设你已经通过pm.sample_posterior_predictive()得到后验预测结果ppc，先从中提取所有预测样本。假设你的观测变量名为obs：

# 提取所有后验预测样本，形状为(链数, 抽样数, 样本量)
predicted_samples = ppc.posterior_predictive["obs"].values
# 重塑为(总抽样组数, 单组样本量)，方便后续计算
predicted_samples_flat = predicted_samples.reshape(-1, predicted_samples.shape[-1])

2. 计算最大值分布 #

分别计算每组后验预测样本的最大值，以及原始样本的最大值：

import numpy as np

# 每组后验预测样本的最大值
predicted_maxes = np.max(predicted_samples_flat, axis=1)
# 原始样本的最大值
original_max = np.max(your_original_data)

替换your_original_data为你手动生成的20个样本数组。

3. 绘制对比直方图 #

用Matplotlib绘制直方图，标记原始样本最大值的位置：

import matplotlib.pyplot as plt

plt.hist(predicted_maxes, bins=30, alpha=0.7, label='后验预测样本最大值')
plt.axvline(original_max, color='red', linestyle='--', linewidth=2, label='原始样本最大值')
plt.xlabel('样本最大值')
plt.ylabel('频数')
plt.legend()
plt.title('后验预测最大值分布 vs 原始样本最大值')
plt.show()

额外提示 #

• 若要获取更多后验预测样本，在调用sample_posterior_predictive时增加draws参数，比如draws=5000，这样直方图会更平滑。
• 观察原始最大值在预测分布中的位置：如果它落在分布的极端尾部甚至外侧，说明模型对极端值的拟合可能存在问题。

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