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

前几章介绍了网格近似、Laplace 近似、重要性抽样和手写 MCMC。这些内容有助于理解贝叶斯计算的原理。但在真实项目中，复杂模型通常依赖成熟软件完成，例如 Stan、JAGS、brms 和 rstanarm。本章的目标不是把这些工具的语法全部讲完，而是帮助读者理解它们各自适合什么场景、输出如何解释、诊断为什么重要，以及本科课程项目中应该如何组织一个可靠的贝叶斯分析。

本章代码块主要作为工作流示意，默认不在编译教材时运行。若要实际运行，需要先在本机安装相应软件和 R 包。

20.1 从算法到建模工具

手写 MCMC 的好处是透明：我们清楚知道每一步如何提议、接受或拒绝。但手写算法也有明显限制：

1. 高维参数的提议分布很难调；
2. 层次模型和广义线性模型的后验几何可能很复杂；
3. 诊断、可视化和模型比较需要大量配套代码；
4. 稍有编程错误，就可能得到看似合理但实际错误的后验样本。

现代贝叶斯工具把许多通用计算细节封装起来。用户更关注模型本身：响应变量如何生成、参数代表什么、先验如何设定、预测如何定义。软件负责把模型翻译成可计算的形式，运行抽样算法，并输出诊断指标。

不过，软件不会自动保证模型正确。贝叶斯建模仍然需要回答以下问题：

1. 似然是否对应数据类型？
2. 参数尺度是否合理？
3. 先验是否过强、过弱或与领域知识冲突？
4. 抽样诊断是否通过？
5. 后验预测数据是否像真实数据？
6. 结论是否服务于原来的经济统计问题？

20.2 工具选择

不同工具的抽象层次不同。越接近底层，模型表达越灵活，但需要写更多代码；越接近高层，使用越方便，但模型细节更容易被公式接口隐藏。

工具	典型接口	主要优势	需要注意
Stan / CmdStanR	建模语言	灵活、HMC/NUTS 效率高、诊断丰富	需要学习 Stan 语法，离散参数通常要边际化
JAGS / rjags	BUGS 风格图模型	语法接近概率图，适合教学和某些共轭层次模型	高相关连续参数模型中可能混合较慢
brms	R 公式接口	可拟合复杂回归、层次、非线性模型，底层调用 Stan	模型会编译，首次运行较慢，仍需检查先验和诊断
rstanarm	类似 lm() / glm()	常见回归模型上手快，Stan 程序预编译	灵活性低于 brms 和手写 Stan
posterior / bayesplot / loo	后处理工具	总结抽样、画诊断图、后验预测检查、LOO	不能替代模型理解

一个实用建议是：如果目标是教学和理解算法，可以从 JAGS 或手写 MCMC 入手；如果目标是应用研究中的连续参数模型，优先考虑 Stan、brms 或 rstanarm；如果模型已经超出公式接口，就写 Stan 模型。

20.3 Stan 与 Hamiltonian Monte Carlo

Stan 是现代贝叶斯建模中非常常用的工具。用户写出数据、参数和模型，Stan 使用 Hamiltonian MonteCarlo（HMC）及其自适应版本 NUTS（No-U-Turn sampler）进行后验抽样。随机游走 Metropolis 像是在参数空间中试探性移动，维度高时效率可能很低。HMC 利用后验密度的梯度信息，沿着更有方向感的路径移动，因此在连续参数模型中通常更高效。

对本科阶段而言，只需把 HMC 理解为“利用梯度信息改进 MCMC 移动方式”的算法即可，不需要掌握物理Hamiltonian 系统的细节。更重要的是知道：Stan 的诊断指标是模型分析的一部分，不是附属输出。

一个 Stan 程序通常由若干代码块组成：

代码块	作用
data	声明外部传入的数据
parameters	声明未知参数及其约束
transformed parameters	定义由参数变换得到的中间量
model	写先验和似然，定义未归一化后验
generated quantities	生成后验预测、对数似然和派生量

一个简单线性回归 Stan 模型大致如下：

data {
  int<lower=1> N;
  vector[N] y;
  vector[N] income;
}
parameters {
  real alpha;
  real beta;
  real<lower=0> sigma;
}
model {
  alpha ~ normal(0, 10);
  beta ~ normal(0, 10);
  sigma ~ exponential(1);
  y ~ normal(alpha + beta * income, sigma);
}
generated quantities {
  vector[N] y_rep;
  vector[N] log_lik;

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

其中 generated quantities 常用于两个目的：一是生成后验预测样本 y_rep，用于后验预测检查；二是保存逐点对数似然 log_lik，用于 LOO 或 WAIC 等预测比较。

使用 CmdStanR 时，典型工作流如下：

library(cmdstanr)
library(posterior)

mod <- cmdstan_model("linear_regression.stan")

stan_data <- list(
  N = nrow(survey),
  y = survey$consumption,
  income = survey$income
)

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

fit$summary(c("alpha", "beta", "sigma"))
fit$diagnostic_summary()
draws <- fit$draws()

Stan 的优势是灵活、诊断丰富、适合连续参数模型。限制是模型语法需要学习；离散未知参数不能直接作为 parameters 中的参数抽样，通常需要通过边际化或其他方式改写模型。

20.4 JAGS 与 BUGS 风格模型

JAGS（Just Another Gibbs Sampler）使用 BUGS 风格的图模型语言描述贝叶斯模型，常与 Gibbs 抽样和Metropolis 步结合。它的模型写法接近“数据如何由参数生成”的概率图，因此适合教学，也适合许多共轭结构清楚的层次模型。

一个 JAGS 线性回归模型示意如下：

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

  alpha ~ dnorm(0, 0.01)
  beta ~ dnorm(0, 0.01)
  tau ~ dgamma(1, 1)
  sigma <- 1 / sqrt(tau)
}

JAGS 中正态分布使用精度参数 tau，即方差的倒数。这一点和 R 中常用标准差参数的写法不同。上面dnorm(0, 0.01) 表示均值为 0、精度为 0.01，也就是方差为 100。

通过 rjags 调用 JAGS 的工作流大致如下：

library(rjags)
library(coda)

jags_data <- list(
  N = nrow(survey),
  y = survey$consumption,
  income = survey$income
)

model <- jags.model(
  file = "linear_regression.bug",
  data = jags_data,
  n.chains = 4
)

update(model, n.iter = 1000)

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

summary(draws)
plot(draws)

JAGS 的优点是模型表达接近概率图，很多教材和旧代码使用 BUGS/JAGS 风格；缺点是在某些强相关连续参数模型中效率可能不如 Stan。若使用 JAGS，应特别关注 trace plot、自相关和有效样本量。

20.5 brms 与 rstanarm

brms 和 rstanarm 是 R 中更高层的贝叶斯建模接口。它们允许用户使用类似 lm() 或 glm() 的公式语法拟合贝叶斯模型，底层通常调用 Stan。

brms 的特点是公式语法非常灵活，可以表达多层模型、非线性模型、分布回归、零膨胀模型等。它会根据公式自动生成 Stan 代码。使用 brms 拟合消费支出回归可以写成：

library(brms)

fit_brms <- brm(
  consumption ~ income + employed,
  data = survey,
  family = gaussian(),
  prior = c(
    prior(normal(0, 10), class = "Intercept"),
    prior(normal(0, 2), class = "b"),
    prior(exponential(1), class = "sigma")
  ),
  chains = 4,
  iter = 2000,
  seed = 2026
)

summary(fit_brms)
pp_check(fit_brms)
posterior_predict(fit_brms)

rstanarm 的目标是让常见回归模型更容易迁移到贝叶斯框架。很多函数是在熟悉的 R 函数名前加上stan_，例如 stan_glm()、stan_lmer() 等：

library(rstanarm)

fit_arm <- stan_glm(
  consumption ~ income + employed,
  data = survey,
  family = gaussian(),
  prior = normal(0, 2),
  prior_intercept = normal(0, 10),
  prior_aux = exponential(1),
  chains = 4,
  iter = 2000,
  seed = 2026
)

summary(fit_arm)
pp_check(fit_arm)
posterior_interval(fit_arm)

这类公式接口非常适合应用研究，但也容易让使用者忽略模型细节。因此，即使使用高级接口，也要清楚：响应变量分布是什么、连接函数是什么、先验是什么、抽样诊断是否通过。公式接口降低的是编程门槛，不是统计判断门槛。

20.6 抽样结果的后处理

现代贝叶斯工作流通常把“拟合模型”和“分析后验样本”分开。拟合之后，应把后验样本转成统一格式，再做摘要、诊断、可视化、后验预测和模型比较。

posterior 包用于处理后验样本，bayesplot 包用于可视化 MCMC 和后验预测检查，loo 包用于PSIS-LOO 等预测比较。

library(posterior)
library(bayesplot)
library(loo)

draws <- fit$draws()
summarize_draws(draws, "mean", "median", "sd", "q5", "q95", "rhat", "ess_bulk")

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

y_rep <- fit$draws("y_rep", format = "matrix")
ppc_dens_overlay(y = survey$consumption, yrep = y_rep[1:100, ])

log_lik <- fit$draws("log_lik", format = "matrix")
loo_fit <- loo(log_lik)
print(loo_fit)

这里有三个层次的检查：

1. 后验样本本身是否可靠，例如 R-hat、ESS、trace plot；
2. 模型是否能生成类似观测数据的复制数据，例如后验预测检查；
3. 不同模型的预测表现如何，例如 LOO、WAIC 或交叉验证。

第 19 章已经介绍后验预测和 WAIC 的基本思想。本章强调的是：现代工具通常已经把这些步骤封装成函数，但用户仍需要理解每一步在检查什么。

20.7 常见诊断指标

现代贝叶斯软件通常会输出以下诊断：

诊断	看什么	常见问题信号	可能处理方式
trace plot	多条链是否混合	链停在不同区域、趋势未稳定	增加迭代、改初始值、重新参数化
R-hat	链间与链内变异是否一致	明显大于 1	延长抽样、检查模型几何、多链初始值
bulk ESS	后验主体区域的信息量	ESS 很低	改进参数化、增加迭代
tail ESS	后验尾部的信息量	尾部 ESS 很低	谨慎解释极端分位数
divergence	HMC 是否错过困难几何	warmup 后仍有发散	提高 adapt_delta、重新参数化、检查先验
treedepth	NUTS 是否达到最大树深	频繁达到上限	增加 max_treedepth 或改进模型效率
E-BFMI	HMC 能量探索是否充分	数值过低	重新参数化、调整尺度
Pareto k	PSIS-LOO 是否可靠	个别观测 k 过大	使用 K-fold CV、检查影响点、改模型

这些指标不是软件附带的装饰，而是贝叶斯结果可信度的重要组成部分。报告后验均值之前，应先确认抽样过程本身没有明显问题。

R-hat 接近 1 只说明多条链在当前摘要上看起来一致，并不说明模型正确；ESS 足够大也不说明先验合理；后验预测检查通过也不说明模型是真实数据生成机制。诊断只能降低某些风险，不能替代建模判断。

20.8 诊断失败时怎么办

诊断失败并不罕见。比“把警告消掉”更重要的是理解警告的含义。下面给出一些常见处理思路。

如果 R-hat 偏高或 trace plot 混合差，先检查模型是否有明显错误：数据尺度是否极端、参数是否不可识别、先验是否过宽。然后可以增加 warmup 和 sampling 迭代，或者重新参数化模型。

如果 Stan 出现 divergence，不能简单忽略。少量 divergence 也可能意味着后验某些区域没有被正确探索。可以先提高 adapt_delta，例如从 0.8 调到 0.9 或 0.95；若仍有问题，应考虑重新参数化、缩放变量、使用更合理的先验，尤其是层次模型中的中心化与非中心化参数化选择。

如果 ESS 低，需要区分是整体混合慢，还是尾部估计不稳定。若目标只是后验均值，bulk ESS 更关键；若要报告 2.5% 和 97.5% 分位数，tail ESS 更重要。

如果后验预测检查显示模型生成的数据与观测数据系统性不同，说明问题可能不在抽样，而在模型本身。例如线性正态模型无法捕捉偏态、厚尾、零膨胀或异方差，此时应修改似然或加入关键预测变量。

20.9 一个推荐工作流

贝叶斯建模可以按以下流程组织：

1. 明确问题：估计、预测、比较模型，还是量化不确定性？
2. 整理数据与尺度：检查缺失、异常值、变量单位，必要时中心化或标准化。
3. 写出生成模型：包括似然、参数、层次结构、先验和预测量。
4. 做先验预测检查：不用真实响应变量拟合，先看先验能生成什么样的数据。
5. 先做小模型：用简单模型检查代码、变量尺度和先验是否合理。
6. 运行多条链：使用不同初始值，设置足够 warmup 和 sampling 迭代。
7. 检查抽样诊断：R-hat、ESS、trace plot、divergence、treedepth、E-BFMI。
8. 做后验预测检查：确认模型能再现关键数据特征。
9. 比较候选模型：使用理论判断、后验预测、LOO 或交叉验证，而不是只看拟合优度。
10. 报告结论：说明数据、模型、先验、计算设置、诊断、预测表现和经济含义。

这个流程看起来比“调用一个函数”麻烦，但它能避免许多常见错误。贝叶斯分析的可信度来自完整流程，不是来自某个软件名称。

20.10 进一步阅读

建议优先查阅官方文档，因为这些工具更新较快：

1. Stan User’s Guide: https://mc-stan.org/docs/stan-users-guide/
2. CmdStanR reference: https://mc-stan.org/cmdstanr/
3. brms reference: https://paulbuerkner.com/brms/
4. rstanarm reference: https://mc-stan.org/rstanarm/
5. JAGS manuals: https://sourceforge.net/projects/mcmc-jags/files/Manuals/4.x/
6. bayesplot reference: https://mc-stan.org/bayesplot/
7. loo reference: https://mc-stan.org/loo/

20.11 本章小结

Stan、JAGS、brms 和 rstanarm 让复杂贝叶斯模型变得可操作，但软件不能替代统计思考。Stan 适合灵活连续参数模型，brms 和 rstanarm 适合用公式接口快速建立回归模型，JAGS 适合理解 BUGS 风格图模型和某些共轭层次模型。无论使用哪种工具，都应形成完整工作流：写清模型，检查先验，运行可靠计算，查看诊断，做后验预测检查，再解释结果。

20.12 练习

1. 查阅 Stan 或 brms 的线性回归示例，指出其中的似然、先验和后验预测量分别是什么。
2. 解释 R-hat 接近 1 的含义，以及它为什么不能单独保证模型正确。
3. 思考为什么使用 brms 公式接口时仍然需要报告先验和连接函数。
4. 比较 brms 和 rstanarm：如果你要拟合一个带学校随机截距的成绩模型，会如何选择工具？
5. 在一个模拟回归问题中，写出你会放入 generated quantities 的两个量，并说明用途。
6. 设计一个贝叶斯分析报告目录，主题可以是消费支出、违约风险、就业状态或政策评估。