PricingLoginStart free
R Academy · 课时

后验预测检验

通过比较模拟数据与观测数据的分布来验证模型拟合效果

第 4 / 4 课13 个步骤

后验预测检验 是 CoddyKit 上的免费 R Academy 课时。 这是第 4 节课，共 4 节。 你可以在下方免费阅读本课时的完整内容 — 然后在浏览器中使用内置代码编辑器和全天候 AI 导师进行实践。 这是 R Academy 学习路径的一部分，你的进度在网页和 CoddyKit 应用中同步。 R Academy 课程共包含 4 节课。

什么是后验预测检查

拟合贝叶斯模型后，您需要询问：这个模型生成的数据是否看起来像您观察到的数据？后验预测检查（PPC）通过从后验分布中模拟重复数据集 yrep，并将其与观测值 y进行图形比较来回答这个问题。

提取后验样本

extract(fit, pars='mu') 会返回一个带名称的列表；$mu 是该参数所有预热后样本组成的数值向量。使用 4 条链和每条链 1000 次预热后迭代时，您会得到 4000 个样本。

# library(rstan)
# mu_samples    <- extract(fit, pars = 'mu')$mu
# sigma_samples <- extract(fit, pars = 'sigma')$sigma
#
# cat('Samples drawn:', length(mu_samples), '
')
# cat('Posterior mean of mu:', mean(mu_samples), '
')
# cat('90% CI:', quantile(mu_samples, c(0.05, 0.95)), '
')

根据后验样本生成 yrep

对于每个后验抽样值 (mu_s, sigma_s)，模拟一个与原始数据大小相同的重复数据集。将这些数据存储在矩阵 yrep 中，其中每一行代表一个模拟数据集。

# y <- c(2.1, 1.8, 2.4, 1.9, 2.3, 2.0, 1.7, 2.2)
# n_obs <- length(y)
# S     <- length(mu_samples)  # 4000 posterior draws
#
# yrep <- matrix(NA, nrow = S, ncol = n_obs)
# for (s in seq_len(S)) {
#   yrep[s, ] <- rnorm(n_obs, mean = mu_samples[s], sd = sigma_samples[s])
# }
# dim(yrep)  # [4000, 8]

ppc_dens_overlay() — 密度比较

bayesplot::ppc_dens_overlay(y, yrep[1:50,]) 会将观测数据的核密度（深色线）与从 50 个随机选择的模拟数据集中得到的密度（浅色线）叠加显示。模型拟合良好时，深色线应位于浅色线形成的云状区域内。

# library(bayesplot)
#
# ppc_dens_overlay(y, yrep[1:50, ])
#
# Interpretation:
# - Dark line (y_obs) surrounded by light lines (yrep): good fit
# - Dark line systematically outside the cloud: model misfit
# - Light lines much wider than dark: overdispersed model
# - Light lines much narrower than dark: underdispersed model

ppc_stat() — 检查检验统计量

ppc_stat(y, yrep, stat = 'mean') 会显示在每个模拟数据集上计算得到的检验统计量（例如均值）的直方图，并在观测统计量处绘制一条竖线。如果观测值落在直方图的主体区域内，说明模型捕捉到了数据的这一特征。

# library(bayesplot)
#
# ppc_stat(y, yrep, stat = 'mean')   # does model capture the mean?
# ppc_stat(y, yrep, stat = 'sd')     # does model capture spread?
# ppc_stat(y, yrep, stat = 'max')    # does model capture extremes?
#
# If observed stat is in the tail of the histogram,
# the model fails to reproduce that statistic.

贝叶斯 p 值

贝叶斯 p 值（后验预测 p 值）是这样一部分模拟数据集所占的比例：其检验统计量比观测值更加极端。接近 0.5 的值表示校准良好；接近 0 或 1 的值表示模型在该统计量上拟合不佳。

# Bayesian p-value for the mean:
# obs_mean <- mean(y)
# rep_means <- apply(yrep, 1, mean)
# pval <- mean(rep_means >= obs_mean)
# cat('Bayesian p-value (mean):', round(pval, 3), '
')
# # 0.5 is perfect; < 0.05 or > 0.95 suggests misfit

更多 bayesplot PPC 函数

bayesplot 除了密度叠加图之外，还提供许多 PPC 可视化方法：

  • ppc_hist(y, yrep[1:8,]) — 直方图网格
  • ppc_scatter_avg(y, yrep) — 观测值与 yrep 均值的散点图
  • ppc_intervals(y, yrep) — 每个观测值周围的不确定性区间
  • ppc_rootogram(y, yrep) — 用于计数数据
# library(bayesplot)
#
# # Grid of 8 simulated histograms vs the observed
# ppc_hist(y, yrep[1:8, ])
#
# # Scatter: y_obs (x) vs mean of yrep (y) — should hug diagonal
# ppc_scatter_avg(y, yrep)
#
# # 50% and 90% posterior predictive intervals around each y_i
# ppc_intervals(y, yrep)

解读 PPC 图 — 模型拟合不佳

常见的拟合不佳模式及其原因：

  • yrep 过宽 — 先验过于分散或模型存在过度离散
  • yrep 发生偏移 — 似然族选择错误（例如对偏态数据使用正态分布）
  • yrep 未能呈现多峰性 — 需要使用混合模型
  • yrep 无法覆盖极端值 — 需要使用重尾分布
# Example: if data has a long right tail but yrep does not,
# consider switching:
# y ~ normal(mu, sigma)  =>  y ~ student_t(nu, mu, sigma)
#
# Or for count data:
# y ~ poisson(lambda)  =>  y ~ neg_binomial_2(mu, phi)  (overdispersion)
cat('PPCs guide model improvement by revealing specific failure modes
')

Stan 模型块中的 PPC

您可以直接在 Stan 中使用 generated quantities 块生成 yrep。这样可以避免在 R 中重新提取参数，计算结果也等价。

# Stan model with generated quantities:
# '
# generated quantities {
#   array[N] real y_rep;
#   for (n in 1:N) {
#     y_rep[n] = normal_rng(mu, sigma);
#   }
# }
# '
# Then extract in R:
# yrep <- extract(fit, pars = 'y_rep')$y_rep  # [S, N] matrix

留一法交叉验证

除了 PPC 之外，loo::loo(fit) 还会计算留一法交叉验证，用于比较相互竞争的模型。ELPD（预期对数预测密度）较高的模型更受推荐。使用 loo::loo_compare(loo1, loo2) 可以对模型进行排序。

# library(loo)
# loo1 <- loo(fit1)  # normal model
# loo2 <- loo(fit2)  # student-t model
#
# comparison <- loo_compare(loo1, loo2)
# print(comparison)
#
# Model with elpd_diff > 0 is preferred
# se_diff > |elpd_diff| means difference is not reliable

PPC 最佳实践

请遵循以下实践，以严谨地进行后验预测检查：

  • 始终先使用 ppc_dens_overlay() 进行全局合理性检查
  • 随后使用与分析目标相关的领域特定统计量（ppc_stat()）进行检查
  • 进行可视化检查时至少使用 50 个 yrep 抽样值；计算贝叶斯 p 值时使用全部 4000 个样本
  • PPC 未通过可以指导模型改进 — 它们是诊断信息，而不是失败
# Workflow:
# 1. Fit model -> extract() -> generate yrep matrix
# 2. ppc_dens_overlay(y, yrep[1:50,])  -- visual global check
# 3. ppc_stat(y, yrep, stat='mean')    -- check mean
# 4. ppc_stat(y, yrep, stat='sd')      -- check spread
# 5. ppc_stat(y, yrep, stat='max')     -- check tails
# 6. If misfit found -> revise model -> refit -> re-check

快速检查：贝叶斯 p 值解读

最大值统计量的贝叶斯 p 值为 0.03，这对模型意味着什么？

在交互式课程中回答

后验预测检查回顾

RStan 和 bayesplot 中的 PPC 工作流程：

  • 提取样本：extract(fit, pars='mu')$mu
  • 生成 yrep：遍历后验抽样值，调用 rnorm(n, mu_s, sigma_s)
  • 全局检查：ppc_dens_overlay(y, yrep[1:50,])
  • 统计量检查：ppc_stat(y, yrep, stat='mean')
  • 贝叶斯 p 值：mean(apply(yrep,1,stat) >= stat(y)) — 接近 0.5 表示良好
  • 比较模型时使用 loo::loo_compare()
免费开始

用 AI 导师学习 R — 免费

在浏览器中编写并运行真实代码，获得全天候 AI 导师的即时帮助，并在网页或应用中继续学习。

课程
43
课程
159
免费开始学习查看完整课程

常见问题解答

「后验预测检验」课时是免费的吗？

是的 — 「后验预测检验」的完整文本可在网页上免费阅读。要进行交互式练习（内置代码编辑器和全天候 AI 导师）并解锁 R Academy 课程的其余内容，请升级到 CoddyKit PRO。 R Academy 课程共包含 4 节课。

「后验预测检验」这节课中我会学到什么？

通过比较模拟数据与观测数据的分布来验证模型拟合效果 你通过在浏览器中直接运行的动手代码来练习 R Academy，全天候 AI 导师会在你学习这节课的过程中回答你的问题。

学习 R Academy 需要有经验吗？

无需任何先前经验。CoddyKit 上的 R Academy 课程适合初学者到高级学习者，你可以从这里开始或从头开始，按照自己的节奏学习。 这是第 4 节课，共 4 节。

「后验预测检验」课时需要多长时间？

大多数 CoddyKit 课程大约需要 5–10 分钟。每节课都很精短且互动，所以你能稳步进步，并在网页和应用中从离开的地方继续。

我能在这节 R Academy 课中编写并运行代码吗？

能。每节 R Academy 课都包含内置代码编辑器，你可以在浏览器中直接编写并运行真实代码，并获得即时 AI 反馈 — 无需本地设置。

此课程中的所有课时

  1. 贝叶斯思维入门
  2. 在 R 中编写 Stan 模型
  3. MCMC 抽样与诊断
  4. 后验预测检验
← 返回 R Academy