第 20 章 现代贝叶斯计算工具

第 17 章通过手写算法解释 MCMC，第 18 章把这些算法用于回归，第 19 章又说明后验样本怎样进入预测和模型评价。到这里，读者还需要把这些环节连接到真实建模工具：模型怎样翻译成程序，抽样器怎样运行，输出又怎样接受诊断？

现代贝叶斯软件使复杂模型的计算变得可行，但也容易制造一种错觉：只要模型成功运行并打印出一张系数表，分析就已经结束。事实上，软件只能对“写进程序的模型”做计算。似然选错、先验尺度失当、参数不可识别或预测目标定义含混，都可能得到数值整齐但没有意义的结果。本章因此不把 Stan、JAGS、brms 和 rstanarm 当作若干互不相干的软件来介绍，而是围绕一条完整主线展开：把统计模型翻译成程序，获得后验样本，检查计算质量，再用这些样本回答预测和决策问题。

本章代码用于说明工作流，均设置为不执行。读者实际运行时需要另行安装 R 包、CmdStan 或 JAGS，并根据本机环境配置编译工具链。

20.1 学习目标

学完本章后，读者应当能够：

  1. 区分统计模型、抽样算法、后验样本和模型评价四个层次；
  2. 根据模型结构在 Stan、JAGS、brms 和 rstanarm 之间作出有理由的选择；
  3. 读懂一个基本 Stan 程序，并说明预热（warmup）、NUTS 和重参数化各自解决什么问题；
  4. 从多条链的输出中检查秩标准化 \(\widehat R\)、有效样本量、Monte Carlo 标准误和 HMC 警告；
  5. 保存后验预测量与逐点对数似然，并正确交给 posterior、bayesplot 和 loo；
  6. 把计算诊断、模型检查和样本外评价组织成可复现的分析流程。

20.2 从统计模型到计算工具

20.2.1 统计分析的四个层次

贝叶斯分析至少包含四个层次。混淆这些层次，是实际工作中许多错误的来源。

层次 核心问题 典型对象
统计模型 数据被假定如何生成？参数和先验表示什么？ 似然、先验、层次结构
计算方法 如何近似后验分布？ Gibbs、Metropolis、HMC、NUTS
计算结果 抽样是否可靠，误差是否足够小？ 多条链、\(\widehat R\)、ESS、MCSE
模型评价 模型能否再现和预测关心的现象？ 复制数据、LOO、留出预测

例如，Stan 使用 NUTS 得到很大的有效样本量，只能说明某些后验期望被较精确地计算出来；它不能证明正态似然适合数据。反过来，一张后验预测图看起来不错，也不能抵消大量发散转移带来的计算疑问。稳妥的顺序是先确认程序表达了预期模型，再检查抽样，最后讨论模型表现。

现代工具替我们完成自动微分、抽样器调节、样本存储和许多诊断计算。它们没有替我们回答以下问题：

  1. 观测单位是什么，似然的条件独立假设是否合理？
  2. 参数采用什么量纲，先验在数据尺度上意味着什么？
  3. 模型是否可识别，层次结构是否符合数据收集过程？
  4. 预测的是原样本复制、新个体，还是未来时期？
  5. 结果怎样进入研究问题或实际决策？

所以，现代贝叶斯计算不是“选一个包并调用函数”，而是统计建模与数值计算共同组成的工作流。

20.2.2 工具的抽象层次与选择

不同工具提供不同程度的抽象。选择时首先看模型结构和分析目的，而不是只比较代码长短。

工具 适合的任务 主要优势 主要限制
Stan / CmdStanR 连续参数的回归、层次和自定义模型 模型表达灵活，NUTS 效率高，诊断完整 需要学习 Stan；离散未知参数须边际化
JAGS / rjags BUGS 风格模型、共轭结构、教学演示 生成式写法直观，可包含离散随机节点 强相关或复杂连续后验可能混合较慢
brms 多层、非线性、分布回归等公式模型 公式丰富，可查看自动生成的 Stan 代码 抽象层较高，仍须理解先验和参数化
rstanarm 常见回归和多层模型 接口接近 lm()、glm()，常用模型易上手 自定义结构的自由度较低
posterior、bayesplot、loo 统一后处理、绘图和预测评价 可跨接口处理抽样结果 不负责修正模型或抽样问题

如果课程目标是理解 Gibbs 抽样和概率图，JAGS 很合适；如果要写一个含特殊似然或参数约束的连续模型，Stan 通常更直接；如果模型属于回归族且公式能够清楚表达，先用 brms 或 rstanarm 往往更省力。高级接口并不比手写 Stan “不严谨”，严谨性取决于模型、检查和报告，而不是代码行数。

20.3 贯穿本章的回归例子

20.3.1 数据、尺度与先验

设某调查记录家庭消费支出 consumption 与收入 income。为使本章代码能够独立复现，下面先生成一份具有缺失记录的模拟调查数据；真实分析应把这一步替换为有来源说明的数据导入。为突出计算过程，暂时采用正态线性模型：

\[ Y_i\mid x_i,\alpha,\beta,\sigma \overset{\mathrm{ind}}{\sim}N(\alpha+\beta x_i,\sigma^2), \qquad i=1,\ldots,n. \]

其中 \(Y_i\) 是随机消费响应，其观测实现记为 \(y_i\)；\(x_i\) 是作为给定条件进入模型的收入协变量。

直接在人民币原值上给 \(\alpha\)、\(\beta\) 和 \(\sigma\) 指定诸如 \(N(0,10^2)\)、\(N(0,2^2)\) 和\(\operatorname{Exp}(1)\) 的先验，通常没有清楚含义。变量的数值单位一变，这些先验的实际约束也会完全改变。这里把变换中心和尺度作为分析方案的一部分预先固定，而不是从即将评价的响应数据反推：

set.seed(2001)
n_household <- 320
survey <- data.frame(
  income = exp(rnorm(n_household, log(9000), 0.45))
)
survey$consumption <- 1500 + 0.70 * survey$income +
  rnorm(n_household, sd = 1800)
survey$income[c(37, 144)] <- NA_real_
survey$consumption[c(82, 219)] <- NA_real_

vars <- c("consumption", "income")
survey_model <- survey[complete.cases(survey[vars]), vars]

stopifnot(
  nrow(survey_model) > 1,
  all(is.finite(survey_model$consumption)),
  all(is.finite(survey_model$income)),
  sd(survey_model$consumption) > 0,
  sd(survey_model$income) > 0
)

y_center <- 8000
y_scale  <- 5000
x_center <- 9000
x_scale  <- 5000

survey_model$y_std <-
  (survey_model$consumption - y_center) / y_scale
survey_model$x_std <-
  (survey_model$income - x_center) / x_scale

这段代码不只是数据清洗。它明确了进入模型的观测、缺失值处理和参数尺度。这里删除缺失记录只是为了简化例子；正式分析必须说明完整案例为何成立，必要时应建立缺失机制或插补模型。固定的四个换算常数来自事先约定的业务量纲，因此不会让留出观测进入预处理。如果改用样本均值和标准差，训练—测试划分与 K 折验证都必须在各自训练部分重新估计，再原样应用到对应留出部分。

在标准化尺度上，可先考虑

\[ \alpha\sim N(0,1),\qquad \beta\sim N(0,1),\qquad \sigma\sim N^+(0,1), \]

其中 \(N^+\) 表示限制在正半轴的正态分布。这些先验不是自动正确，但含义容易检查：在先验下，响应均值大多位于若干个响应标准差以内，收入每增加一个标准差，条件均值通常不会移动许多个响应标准差。是否合理还要通过先验预测检查判断。

标准化不会使解释消失。原单位下的斜率和截距分别为

\[ \beta_{\mathrm{raw}}=\beta\frac{s_y}{s_x},\qquad \alpha_{\mathrm{raw}} =c_y+s_y\alpha-\beta\frac{s_yc_x}{s_x}. \]

其中 \(c_y,s_y,c_x,s_x\) 分别是预先固定的响应中心、响应尺度、自变量中心和自变量尺度。对每个后验抽样值作同样变换，就得到原单位下的后验分布，而不是只变换一个后验均值。

20.3.2 Stan 模型文件

Stan 程序由若干功能明确的代码块组成。并非每个模型都需要全部代码块。

代码块 用途
data 声明从外部传入且在抽样中固定的数据
parameters 声明未知的连续参数及其约束
transformed parameters 定义抽样过程中反复使用的参数变换
model 累加先验和似然的对数密度
generated quantities 在每次抽样后生成预测值、逐点对数似然和派生量

标准化回归可以写成下面的 Stan 文件。这里使用 stan 代码块而不是 R 代码块，因为这段程序应保存为例如随书提供的 rcode/consumption_regression.stan 文件。

data {
  int<lower=1> N;
  vector[N] y;
  vector[N] x;
}
parameters {
  real alpha;
  real beta;
  real<lower=0> sigma;
}
model {
  alpha ~ normal(0, 1);
  beta  ~ normal(0, 1);
  sigma ~ normal(0, 1);

  y ~ normal(alpha + beta * x, sigma);
}
generated quantities {
  vector[N] y_rep;
  vector[N] log_lik;

  for (i in 1:N) {
    real mu_i = alpha + beta * x[i];
    y_rep[i] = normal_rng(mu_i, sigma);
    log_lik[i] = normal_lpdf(y[i] | mu_i, sigma);
  }
}

parameters 中的下界保证 sigma 为正；在这个约束下，sigma ~ normal(0, 1) 对应半正态形状。model 中的最后一行采用向量化写法，与逐个观测累加同一个似然等价。

generated quantities 不改变后验分布。y_rep 表示在原协变量条件下生成的复制数据\(\mathbf y^{\mathrm{rep}}\)，用于后验预测检查；log_lik[i] 是第 \(i\) 个观测的对数似然贡献，用于 PSIS–LOO 等逐点预测评价。两者都应在建模时预先设计，而不是得到结果后才临时拼凑。

20.4 NUTS、后验几何与抽样

20.4.1 HMC、NUTS 与预热

随机游走 Metropolis 每次在当前位置附近试探，参数维数增加或后验高度相关时，常出现移动缓慢的问题。HMC 利用对数后验的梯度构造一条数值轨迹，可以在保持较高接受概率的同时移动较远。数值积分并不精确，因此算法仍需接受—拒绝校正。

NUTS 在 HMC 基础上自动决定每次轨迹的大致长度，避免使用者手工指定固定步数。预热期间还会调节两个重要部分：

  1. 步长决定每次数值积分前进多远。步长过大容易产生较大积分误差，过小则计算昂贵；
  2. 质量矩阵用来适应不同参数的尺度和相关性，使轨迹更适合后验几何。

自动调节降低了操作门槛，却没有消除困难后验。若后验包含狭窄弯曲的区域、漏斗结构或几乎不可识别的方向，有限步长的数值轨迹仍可能失败。Stan 报告的 divergence 正是在提醒我们：某些轨迹出现了异常大的数值误差，后验的一部分可能没有被可靠探索。它不是普通的“提议被拒绝”，也不应仅按占比小就忽略。

20.4.2 参数化与模型限制

层次模型最能说明“统计上等价的写法，计算上未必等价”。设组参数满足

\[ \theta_j\sim N(\mu,\tau^2),\qquad j=1,\ldots,J. \]

直接把 \(\theta_j\) 作为参数的中心化写法为：

下面两段只展示参数化中不同的部分，并不是可独立编译的完整模型；共同的 \(\mu\)、\(\tau\) 先验以及数据似然被省略了。

parameters {
  real mu;
  real<lower=0> tau;
  vector[J] theta;
}
model {
  theta ~ normal(mu, tau);
}

非中心化写法先引入标准正态变量，再作确定性变换：

parameters {
  real mu;
  real<lower=0> tau;
  vector[J] z;
}
transformed parameters {
  vector[J] theta = mu + tau * z;
}
model {
  z ~ std_normal();
}

两种写法定义同一个条件分布，但参数空间的几何不同。当每组信息较少、\(\tau\) 可能接近 0 时，中心化参数常形成漏斗，非中心化写法往往更易抽样；每组数据非常充足时，中心化写法有时反而更有效率。不存在对所有数据都占优的参数化，应结合数据强度、发散位置、链图和有效样本量判断。

变量缩放、更有信息的合理先验，以及把强相关参数改写成较弱相关参数，也都在改变计算所见的几何。这说明重参数化不是为了讨好软件，而是在不改变研究问题的前提下，让同一个概率模型更容易计算。

参数化之外，还要注意 HMC 对参数类型的限制。Stan 的 parameters 块不能包含未知离散参数，因为HMC 需要梯度。有限个离散状态常可通过求和从后验中边际化，再在 generated quantities 中恢复状态概率。JAGS 可以直接表示许多离散随机节点，但这不意味着混合、多峰或标签交换问题会自动消失。

20.4.3 从模型文件到后验样本

理解抽样器的任务以后，再把前面的回归程序交给 CmdStanR。随书文件rcode/consumption_regression.stan 与正文中的 Stan 程序一致。一个最小工作流如下：

library(cmdstanr)
library(posterior)

stan_data <- list(
  N = nrow(survey_model),
  y = survey_model$y_std,
  x = survey_model$x_std
)

mod <- cmdstan_model("rcode/consumption_regression.stan")

fit <- mod$sample(
  data = stan_data,
  chains = 4,
  parallel_chains = 4,
  iter_warmup = 1000,
  iter_sampling = 1000,
  seed = 2026
)

fit$summary(
  variables = c("alpha", "beta", "sigma"),
  posterior::default_summary_measures(),
  posterior::default_convergence_measures(),
  posterior::default_mcse_measures()
)
fit$diagnostic_summary()

draws <- fit$draws(
  variables = c("alpha", "beta", "sigma"),
  format = "draws_array"
)

四条链是四次独立运行，而不是同一条链的复制。CmdStanR 会在给定总种子后管理各链的随机数流。iter_warmup 对应前面所述的步长和质量矩阵调节，预热样本不进入正式后验汇总。抽样结束后，fit$summary() 汇总参数后验以及秩标准化 \(\widehat R\)、有效样本量和 MCSE；fit$diagnostic_summary() 汇总发散、最大树深和能量等 HMC 诊断。两者检查的对象不同。

draws_array 保留“迭代 × 链 × 变量”的结构，适合链间诊断；矩阵格式通常合并各链，适合某些后续矩阵运算。转换格式前应先判断后续任务是否仍需链信息。

20.5 其他建模接口

20.5.1 JAGS：基于条件分布的节点更新

JAGS 使用 BUGS 风格的生成式语言描述模型，并自动为不同节点安排 Gibbs 或 Metropolis 类更新。同一个标准化回归可写成：

model {
  for (i in 1:N) {
    y[i] ~ dnorm(mu[i], tau_y)
    mu[i] <- alpha + beta * x[i]
    y_rep[i] ~ dnorm(mu[i], tau_y)
  }

  alpha ~ dnorm(0, 1)
  beta  ~ dnorm(0, 1)
  sigma ~ dnorm(0, 1) T(0, )
  tau_y <- pow(sigma, -2)
}

JAGS 的 dnorm(mean, precision) 使用精度而不是标准差，\(\text{precision}=1/\sigma^2\)。因此dnorm(0, 1) 在这里表示方差为 1。把 R 的 rnorm(mean, sd) 习惯直接照搬到 JAGS，是非常常见的参数化错误。

通过 rjags 运行时，应把适应、烧入和正式留样明确分开：

library(rjags)
library(coda)
library(posterior)

jags_data <- list(
  N = nrow(survey_model),
  y = survey_model$y_std,
  x = survey_model$x_std
)

inits <- list(
  list(alpha = -0.5, beta = -0.5, sigma = 0.8,
       .RNG.name = "base::Mersenne-Twister", .RNG.seed = 20261),
  list(alpha =  0.0, beta =  0.5, sigma = 1.0,
       .RNG.name = "base::Mersenne-Twister", .RNG.seed = 20262),
  list(alpha =  0.5, beta = -0.2, sigma = 1.2,
       .RNG.name = "base::Mersenne-Twister", .RNG.seed = 20263),
  list(alpha =  0.2, beta =  0.2, sigma = 1.5,
       .RNG.name = "base::Mersenne-Twister", .RNG.seed = 20264)
)

jags_fit <- jags.model(
  file = "rcode/consumption_regression.bug",
  data = jags_data,
  inits = inits,
  n.chains = 4,
  n.adapt = 1000
)

# 转移核已经适应完成；这一步作为固定核的烧入期
update(jags_fit, n.iter = 1000)

jags_draws <- coda.samples(
  model = jags_fit,
  variable.names = c("alpha", "beta", "sigma"),
  n.iter = 2000
)

summary(jags_draws)
plot(jags_draws)

jags_draws_array <- posterior::as_draws_array(jags_draws)
posterior::summarize_draws(
  jags_draws_array,
  "rhat", "ess_bulk", "ess_tail", "mcse_mean"
)

不同初始值的目的不是人为制造差异，而是观察链能否从不同位置到达同一后验区域。正式分析还应记录每条链的随机数生成器与种子。coda 的经典诊断和 posterior 的现代诊断定义并不完全相同；这里把mcmc.list 转成保留链结构的 draws 数组，再计算秩标准化 \(\widehat R\)、bulk ESS、tail ESS 和 MCSE。若自相关高、链移动慢，只增加留样数虽然可能提高有效样本量，却不会解决不可识别、多峰或参数化不良；先找出混合慢的原因更重要。

Stan 与 JAGS 的结果也不能只凭“谁每秒生成的样本更多”比较。应比较目标估计量的有效样本量、MCSE、诊断警告和总计算成本，同时确认两个程序表达的是同一个模型，尤其要核对标准差与精度、链接函数、截断分布和先验常数。

20.5.2 brms 与 rstanarm：从公式生成模型

许多应用问题并不需要手写建模语言。brms 和 rstanarm 使用熟悉的 R 公式描述回归模型，底层仍可利用 Stan 的计算。它们降低的是程序表达成本，不是模型检查要求。

以 brms 为例，先查看参数类别，再明确指定先验：

library(brms)

form <- bf(y_std ~ 1 + x_std)

get_prior(
  formula = form,
  data = survey_model,
  family = gaussian()
)

priors <- c(
  prior(normal(0, 1), class = "Intercept"),
  prior(normal(0, 1), class = "b"),
  prior(normal(0, 1), class = "sigma")
)

fit_brms_prior <- brm(
  formula = form,
  data = survey_model,
  family = gaussian(),
  prior = priors,
  backend = "cmdstanr",
  chains = 4,
  iter = 1000,
  seed = 2026,
  sample_prior = "only"
)
pp_check(fit_brms_prior)

fit_brms <- brm(
  formula = form,
  data = survey_model,
  family = gaussian(),
  prior = priors,
  backend = "cmdstanr",
  chains = 4,
  iter = 2000,
  seed = 2026,
  sample_prior = "yes"
)

summary(fit_brms)
pp_check(fit_brms)

这里把 sigma 也设为正半轴上的标准正态形状，以便与前面的 Stan 和 JAGS 模型保持一致。get_prior() 很重要，因为公式中的截距、一般系数、组间标准差、残差标准差和相关参数属于不同类别，同一个 prior() 不会自动作用于所有参数。sample_prior = "yes" 可额外保存多数参数的先验抽样值；总体截距的 prior draws 有专门的参数化要求，应按 brms 文档核对。真正的先验预测检查应在使用似然更新之前进行；上面的 fit_brms_prior 正是用 sample_prior = "only" 生成先验预测。此时所有参数都必须有适当先验，再检查模型生成的响应范围和关键统计量，确认无误后才拟合 fit_brms。

brms 可以用 stancode() 查看自动生成的 Stan 程序，用 make_standata() 查看传入 Stan 的数据。当公式变复杂时，这两个函数有助于确认软件实际拟合了什么，而不是只根据公式名称猜测。

rstanarm 更接近常用的 lm()、glm() 和多层模型接口，适合从经典回归平稳过渡到贝叶斯分析。使用时仍要检查 prior_summary()、链接函数、辅助参数先验、标准化约定以及 pp_check()。如果模型逐渐需要自定义似然、特殊缺失机制或复杂参数变换，应转向 brms 的扩展功能或手写 Stan，而不是在不合适的公式接口中勉强拼接。

20.6 从后验样本到模型评价

20.6.1 参数摘要与链的可视化

模型拟合以后，首先保留链结构做计算诊断，再根据任务转换为矩阵或数据框。posterior 包提供统一的draws 格式，使不同接口的后验样本可以使用相近的处理方式。

library(posterior)
library(bayesplot)

draws <- fit$draws(
  variables = c("alpha", "beta", "sigma"),
  format = "draws_array"
)

draw_summary <- summarize_draws(
  draws,
  "mean", "median", "sd", "q5", "q95",
  "rhat", "ess_bulk", "ess_tail", "mcse_mean"
)

draw_summary

mcmc_trace(draws, pars = c("alpha", "beta", "sigma"))
mcmc_rank_overlay(draws, pars = c("alpha", "beta", "sigma"))

均值、标准差和分位数描述后验；\(\widehat R\)、ESS 和 MCSE 描述这些数值由有限相关样本估计得多稳定。二者不能混写成同一类“不确定性”。例如，参数后验标准差很大可能是真实的信息不足，MCSE 很大则表示目前这次计算还不足以精确汇总该后验。

20.6.2 把结果变回研究量纲

如果最终要报告原单位下的收入效应，应逐抽样变换：

draws_df <- as_draws_df(draws)

draws_df$beta_raw <- draws_df$beta * y_scale / x_scale
draws_df$alpha_raw <- y_center + y_scale * draws_df$alpha -
  draws_df$beta * y_scale * x_center / x_scale
draws_df$sigma_raw <- y_scale * draws_df$sigma

raw_draws <- subset_draws(
  draws_df,
  variable = c("alpha_raw", "beta_raw", "sigma_raw")
)
summarize_draws(raw_draws)

斜率、截距和误差标准差都要回到原量纲。先把每个抽样值变换后再取分位数，可以保留联合后验的不对称性和相关性。把一张系数表中的均值、上下限分别代入非线性公式，通常不是同一件事。

20.6.3 后验预测与逐点对数似然

计算诊断通过以后，才进入模型检查。Stan 文件中已经保存了 y_rep，可以把它变回消费支出的原单位：

y_rep_std <- fit$draws("y_rep", format = "matrix")
y_rep_raw <- y_center + y_scale * y_rep_std

set.seed(2026)
keep <- sample(seq_len(nrow(y_rep_raw)), size = 100)

ppc_dens_overlay(
  y = survey_model$consumption,
  yrep = y_rep_raw[keep, ]
)

ppc_stat(
  y = survey_model$consumption,
  yrep = y_rep_raw,
  stat = "sd"
)

密度图检查整体形状，标准差检查离散程度。实际分析应根据研究问题继续检查偏度、极端值、零值比例、组间差异或时间连续性。只选择模型最容易通过的统计量，会使后验预测检查失去作用。

逐点对数似然可交给 loo：

library(loo)

log_lik <- fit$draws("log_lik", format = "draws_array")
r_eff <- loo::relative_eff(exp(log_lik))
waic_fit <- loo::waic(log_lik)
loo_fit <- loo::loo(log_lik, r_eff = r_eff)

print(waic_fit)
print(loo_fit)
pareto_k_table(loo_fit)

三维 draws_array 保留了各条链，relative_eff() 因而能把 MCMC 自相关纳入 PSIS 的有效样本比例。waic() 与 loo() 使用同一逐点对数似然，前者计算 WAIC，后者提供带 Pareto \(k\) 诊断的PSIS–LOO；两者的评价目标和差异解释见第 19 章。log_lik 的每一列必须对应预先定义的预测单位。独立横截面数据中通常是一条观测；同一人的重复测量、时间序列或空间数据中，按单行留出可能破坏依赖结构，此时应按个体、时间块或空间块重定义留出单位。预测一个新群组时，通常还要对尚未观察的群组效应积分；预测未来时期时，应使用留未来法或滚动预测起点。矩阵形状正确并不保证预测问题定义正确。

本例使用事先固定的中心和尺度，所以 LOO 条件于同一变换。若缩放、缺失插补或变量筛选由样本估计，逐点 LOO 只重算似然并不能评价完整建模流程；这些步骤也应在每个训练折内重做，或改用能包住全部预处理的 K 折/滚动验证。

PSIS–LOO 的 Pareto \(k\) 诊断反映重要性权重是否稳定。出现较大的 \(k\) 时，应先检查相应观测为何对后验有强影响，以及模型是否遗漏厚尾、非线性或分组结构。根据模型接口和任务，可使用矩匹配、精确逐点重拟合或 K 折交叉验证。不能仅为消除警告而删除观测。有关 ELPD、差异标准误和模型权重的解释，见第 19 章。

20.7 从诊断到修正

20.7.1 计算、模型与预测诊断

第 17 章已经解释轨迹、秩标准化 \(\widehat R\)、ESS 和 MCSE 的含义。本节不再重复定义，而是说明现代工具给出的警告分别属于哪一层。三类检查的对象不同，不能由其中一个替代其余指标。

诊断 它检查什么 需要警惕的现象
\(\widehat R\)、bulk/tail ESS、MCSE 多条链是否一致，目标量的数值误差是否足够小 链间不一致、尾部 ESS 过低或 MCSE 影响报告位数
发散转移（divergence） HMC 是否可靠穿过困难几何 预热后仍有发散，且集中在漏斗或边界附近
最大树深（maximum treedepth） NUTS 轨迹是否频繁触及计算上限 大量迭代达到上限、单位时间 ESS 很低
能量诊断（E-BFMI） 动量更新能否有效探索能量分布 软件报告能量探索警告
后验预测检查 模型能否再现任务相关的数据特征 系统遗漏偏态、异方差、极端值或分组结构
Pareto \(k\) PSIS–LOO 的重要性近似是否稳定 个别观测的 \(k\) 较大，近似可能失效

\(\widehat R\) 接近 1 只是必要检查之一。若所有链都困在同一个局部区域，它也可能看起来良好。ESS 也没有统一的“够用”数值：估计均值与估计 0.995 分位数需要的有效信息不同。更直接的原则是让目标估计量的MCSE 明显小于报告精度。例如一个概率只报告到小数点后两位，MCSE 却为 0.02，继续解释第二位小数就没有计算依据。

表中后两项不是 MCMC 收敛诊断：Pareto \(k\) 检查 PSIS–LOO 近似，后验预测检查模型与数据的相容性。把计算、模型与预测诊断分开，处理警告时才不会找错方向。

诊断失败后，警告不应被当成需要逐项“调参消除”的软件提示。较稳妥的处理顺序如下。

20.7.2 诊断失败后的处理顺序

第一步是排除程序和数据错误。核对样本数、缺失值、变量单位、因子编码和参数约束；用很小的模拟数据检查参数是否能够恢复；比较不同接口时，逐项核对先验参数化。若模型程序写错，增加一百万次迭代只会更精确地计算错误模型。

第二步是检查后验几何。若 \(\widehat R\) 偏高、链混合慢或出现发散，先看变量尺度、参数可识别性、先验是否允许极端区域，以及层次模型能否重参数化。成对图中发散点集中在漏斗或边界附近，通常比单一警告计数更有信息。

适当提高 adapt_delta 会使 NUTS 使用更小步长，有时能减少发散；但它会增加计算量，也不能修复错误模型或根本性的几何问题。提高 max_treedepth 只允许轨迹走得更久，适用于确认模型几何合理后仍偶尔触及上限的情形。二者都不是看到警告后的第一反应。

第三步才判断是否需要更多样本。如果链已经稳定、没有结构性警告，只是目标量的 MCSE 偏大，那么增加iter_sampling 才是对症处理。延长抽样会增加 ESS，却不会自动降低 \(\widehat R\)、消除多峰或修复发散。是否增加预热长度，要看适应是否充分，而不是把它和正式留样数机械地一起翻倍。

计算诊断通过以后，才根据模型检查修改统计结构。若复制数据仍不能再现偏态、异方差、极端值或组间结构，问题在统计模型。此时应考虑改变似然、均值结构、方差结构或层次结构，而不是继续调 NUTS。若 PSIS–LOO 的高 \(k\) 来自某些强影响观测，也要判断它们是录入错误、真实但稀有的现象，还是暴露了模型缺陷。

下面只保留现代工具中特有的处置映射；\(\widehat R\)、ESS、MCSE 与稀疏化（thinning）的一般误区见第17 章。

现象 不充分的做法 更有信息的处理
有发散 只把 adapt_delta 调到很高 定位发散区域，检查尺度、先验和漏斗结构
达到最大树深 只提高上限 同时看发散、参数相关和每秒有效样本量
后验预测失败 调整抽样器 修改似然或结构，再重新检查
Pareto \(k\) 高 删除对应观测 调查影响来源，采用稳健模型或更可靠的交叉验证

20.8 组织一次完整分析

20.8.1 从问题定义到结果解释

本章的工具可以组织成五个有先后关系的关口，而不是一张需要机械打勾的软件菜单。

阶段 核心工作 进入下一阶段前应确认
研究任务 界定观测单位、目标总体、预测对象和使用场景 参数、预测或决策目标已经写清
模型与数据 写出生成模型，固定纳入、缺失、编码、尺度与验证边界 概率程序与数据契约一致
先验与程序 做先验预测；用小型模拟数据检查输出维度和参数恢复 模型能生成合理数据，程序通过基本核验
后验计算 从分散初值运行多条链，检查 \(\widehat R\)、ESS、MCSE 与 HMC 警告 关键参数和派生量均有可信数值精度
模型使用 做后验预测检查与符合依赖结构的样本外评价，再解释和决策 结论附带预测误差、模型限制和决策条件

前一关未通过，不能靠后一关的漂亮图表掩盖问题。程序尚未通过模拟核验时，增加迭代没有意义；抽样存在发散时，后验预测图也不可靠；候选模型共同遗漏关键结构时，最高的 ELPD 仍不能把模型变成可信的决策依据。

20.8.2 可复现性与结果报告

随机种子只是可复现计算的一部分。软件、编译器和硬件环境改变后，浮点轨迹可能不同；真正需要保证的是分析过程可以被审查和重新运行，而不是要求所有平台逐位生成相同样本。报告内容和应保存的计算材料可以放在同一张表中核对。

环节 报告中说明 应保存的材料
数据 来源、观测单位、纳入规则、缺失处理和尺度变换 原始数据或获取说明、完整预处理代码
模型与先验 似然、链接函数、层次结构、先验参数化及先验预测检查 Stan/JAGS 文件或公式、先验配置
计算 软件版本、链数、预热、正式留样数、种子和非默认设置 环境信息、运行配置、未经合并的后验样本
诊断 关键量的 \(\widehat R\)、bulk/tail ESS、MCSE 及 HMC 警告 诊断输出、链图和警告处理记录
模型评价 任务相关的复制统计量、留出单位、ELPD 差和 Pareto \(k\) y_rep、逐点 log_lik、数据分折
解释 原量纲下的效应或预测、敏感性、限制和决策条件 生成表图与派生量的代码

只保存后验均值表会丢失最重要的核查依据：没有链信息就无法复查计算，没有 y_rep 就难以重做模型检查，没有逐点 log_lik 就无法定位不稳定的留一近似。“所有 \(\widehat R\) 小于 1.01”也不是完整报告；还要说明诊断覆盖了哪些参数和派生量、是否出现 HMC 警告，以及结论是否依赖先验和模型结构。

20.9 本章小结

现代贝叶斯计算工具连接了概率模型与可用的后验样本。Stan 通过自动微分和 NUTS 高效处理许多连续参数模型；JAGS 采用 BUGS 风格的节点更新，便于表达概率图并理解 Gibbs 类计算；brms 和rstanarm 用公式接口降低常见回归模型的程序成本。选择工具时应根据模型结构和分析目标，而不是把任何软件当作默认答案。

可信结果需要三类检查。\(\widehat R\)、ESS、MCSE、发散和能量诊断回答“后验是否被可靠计算”；后验预测回答“模型能否生成关心的数据特征”；PSIS–LOO 和留出实验回答“模型对目标预测单位表现如何”。三者不能互相替代。调高 adapt_delta、增加迭代或提高树深也不是通用修复；先识别问题属于程序、后验几何、Monte Carlo 精度还是统计模型，才能采取合适措施。

最终，贝叶斯计算的产物不应只是一张系数表，而应包括模型程序、可复查的后验样本、计算诊断、复制数据、逐点预测信息和清楚的分析记录。这些材料共同说明结果是怎样计算出来的，以及为什么值得相信。

20.10 思考题

  1. 为什么“Stan 没有报错”不能推出模型正确？请分别从统计模型和数值计算两方面回答。
  2. 标准化变量以后，参数更容易抽样是否意味着结果更难解释？写出恢复原单位斜率的方法。
  3. 预热与正式抽样有什么区别？为什么不宜把预热只理解成需要丢弃的前若干次迭代？
  4. \(\widehat R\) 接近 1、ESS 很大而后验预测检查失败时，应修改抽样器还是统计模型？为什么？
  5. 发散转移与普通的 Metropolis 拒绝有什么不同？为什么不能只根据发散比例很小就忽略它？
  6. 在什么数据情形下，层次模型的非中心化参数化可能比中心化参数化更合适？它们是否定义不同模型？
  7. JAGS 的 dnorm(0, 0.01) 与 R 的 rnorm(n, 0, 0.01) 有何区别？这种混淆会怎样改变先验？
  8. 为什么 Pareto \(k\) 不属于 MCMC 收敛诊断？高 \(k\) 与高 \(\widehat R\) 分别指向什么问题？
  9. 对同一人的多期观测逐行计算 LOO，可能违背什么预测目标？请提出更合适的留出单位。
  10. 为什么在没有存储限制时，不应把稀疏化当作处理低 ESS 的首选方法？

20.11 上机实验

  1. 先验预测。 为标准化线性回归分别设置 \(\beta\sim N(0,1)\) 和 \(\beta\sim N(0,10^2)\)，生成先验预测数据，比较响应范围和回归线。说明哪一个先验更符合你的案例背景。
  2. 参数恢复。 自定 \(\alpha\)、\(\beta\) 和 \(\sigma\) 生成数据，分别用 Stan 或 brms 拟合。重复改变样本量，比较后验区间覆盖、ESS 和 MCSE。
  3. 尺度实验。 分别用收入原值和标准化收入拟合同一模型，比较后验相关、有效样本量、树深和运行时间。确保两种写法经过变换后表达相同的先验，再讨论计算差异。
  4. 层次参数化。 模拟每组样本很少的随机截距数据，分别写中心化和非中心化 Stan 模型，比较发散、\(\widehat R\)、ESS 和每秒有效样本量；再增大每组样本量重复实验。
  5. 跨软件核对。 在 Stan 与 JAGS 中拟合同一个标准化正态回归。逐项核对先验参数化，再比较后验摘要和计算效率。不得只比较两张默认 summary() 表。
  6. 目标量与 MCSE。 同时估计后验均值、0.95 分位数和 \(P(\beta>0)\)，考察增加留样数时三者 MCSE如何变化，并据此决定报告位数。
  7. 后验预测。 在数据中制造异方差或一个厚尾观测，用同方差正态模型拟合。设计至少三个复制统计量，说明每个统计量发现了模型的哪项不足。
  8. LOO 诊断。 保存逐点 log_lik 并运行 PSIS–LOO，定位最大的 Pareto \(k\)。比较删除观测、稳健似然和 K 折交叉验证三种处理思路，说明为什么“删除警告点”不是默认方案。
  9. 计算报告。 按本章模板整理一次完整分析，附模型文件、数据处理、版本信息、诊断图、后验预测图和预测评价。请另一位同学只根据这些材料复现主要结论。

20.12 延伸阅读

Stan 的模型语言、重参数化与离散参数边际化可查阅Stan User’s Guide；CmdStanR 的抽样、draws 格式和诊断接口可查阅 CmdStanR Reference。公式建模可参考brms 文档和rstanarm 文档。JAGS 的分布参数化与采样设置应以JAGS manuals为准。后验样本、诊断图与预测评价分别可参考 posterior、bayesplot和 loo 的官方文档。