flowchart TB boost(boost) --> math(math) eigen(eigen) --> math(math) opencl(opencl) --> math(math) sundials(sundials) --> math(math) onetbb(onetbb) --> math(math) math(math) --> Stan(Stan) Stan(Stan) --> CmdStan(CmdStan) Stan(Stan) --> RStan(RStan) RStan --> rstanarm(rstanarm) RStan --> brms(brms) RStan --> prophet(prophet) CmdStan --> CmdStanR(CmdStanR) CmdStan --> CmdStanPy(CmdStanPy) CmdStan --> MathematicaStan(MathematicaStan)
31 广义线性模型
Stan 是一个贝叶斯统计建模和计算的概率推理框架，也是一门用于贝叶斯推断和优化的概率编程语言 (Gelman, Lee, 和 Guo 2015; Carpenter 等 2017)。它使用汉密尔顿蒙特卡罗算法（Hamiltonian Monte Carlo algorithm ，简称 HMC 算法）抽样，内置一种可以自适应调整采样步长的 No-U-Turn sampler （简称 NUTS 采样器） 。Stan 还提供自动微分变分推断（Automatic Differentiation Variational Inference algorithm 简称 ADVI 算法）算法做近似贝叶斯推断获取参数的后验分布，以及拟牛顿法（the limited memory Broyden-Fletcher-Goldfarb-Shanno algorithm 简称 L-BFGS 算法）优化算法获取参数的惩罚极大似然估计。
经过 10 多年的发展，Stan 已经形成一个相对成熟的生态，它提供统计建模、数据分析和预测能力，广泛应用于社会、生物、物理、工程、商业等领域，在学术界和工业界的影响力也不小。下 图 31.1 是 Stan 生态中各组件依赖架构图，math 库是 Stan 框架最核心的组件，它基于 boost 、eigen 、opencl 、sundials 和 oneTBB 等诸多 C++ 库，提供概率推理、自动微分、矩阵计算、并行计算、GPU 计算和求解代数微分方程等功能。
CmdStan 是 Stan 的命令行接口，可在 MacOS / Linux 的终端软件，Windows 的命令行窗口或 PowerShell 软件中使用。CmdStanR、CmdStanPy 和 MathematicaStan 分别是 CmdStan 的 R 语言、Python 语言和 Mathematica 语言接口。每次当 Stan 发布新版本时，CmdStan 也会随之发布新版，只需指定新的 CmdStan 安装路径，CmdStanR 就可以使用上，CmdStanR 包与 Stan 是相互独立的更新机制。 CmdStanR 负责处理 CmdStan 运行的结果，而编译代码，生成模型和模拟采样等都是由 CmdStan 完成。入门 CmdStanR 后，可以快速转入对 Stan 底层原理的学习，有利于编码符合实际需要的复杂模型，有利于掌握常用的炼丹技巧，提高科研和工作的效率。
rstan 是 Stan 的 R 语言接口，该接口依赖 Rcpp、RcppEigen、BH、RcppParallel 和 StanHeaders 等 R 包，由于存在众多上游 R 包依赖，RStan 的更新通常滞后于 Stan 的更新，而且滞后很多，不利于及时地使用最新的学术研究成果。 因此，相比于 rstan 包，CmdStanR 更加轻量，可以更快地将 CmdStan 的新功能融入进来，cmdstanr 和 CmdStan 是分离的，方便用户升级和维护。
rstanarm 和 brms 是 RStan 的扩展包，各自提供了一套用于表示统计模型的公式语法。它们都支持丰富的统计模型，比如线性模型、广义线性模型、线性混合效应模型、广义线性混合效应模型等。相比于 rstan， 它们使用起来更加方便，因为它内置了大量统计模型的 Stan 实现，即将公式语法翻译成 Stan 编码的模型，然后调用 rstan 或 cmdstanr 翻译成 C++，最后编译成动态链接库。除了依赖 rstan 包，rstanarm 和 brms 还依赖大量其它 R 包，因此，安装、更新都比较麻烦。
顺便一提，类似的用于概率推理和统计分析的框架，还有 Python 社区的 PyMC3 和 TensorFlow Probability。
31.1 概率编程语言
Stan 是一门用于贝叶斯推断和优化的概率编程语言。下面以一个简单示例介绍 Stan 的用法，包括 Stan 的基本用法、变量类型、代码结构等，
31.1.1 Stan 的基本用法
考虑一个已知方差的正态分布，设 \(-3, -2, -1, 0, 1, 2, 3\) 是取自正态分布 \(\mathcal{N}(\mu,1)\) 的一个样本，也是取自该正态分布的一组随机数。现在的问题是估计该正态分布的均值参数 \(\mu\) 。Stan 编码的正态分布模型如下：
transformed data代码块是一组已知的数据，这部分数据是不需要从外部传递进来的。这个样本是以向量存储的，需要声明向量的长度和类型（默认类型是实数），每一行以分号结尾，这与 C++ 的语法一样。parameters代码块是未知的参数，需要声明各个参数的类型。这里只有一个参数，且只是一个未知的实数，声明类型即可。model代码块是抽样语句表示的模型结构，符号~表示服从的意思，函数y ~ normal(mu, 1)是正态分布的抽样语句。
接下来，编译 Stan 代码，准备参数初值，配置采样的参数。首先加载 cmdstanr 包，设置 2 条迭代链，给每条链设置相同的参数初始值。代码编译后，生成一个模型对象 mod_gaussian，接着，调用方法 sample() ，传递迭代初值 init，初始化阶段的迭代次数 iter_warmup ，采样阶段的迭代次数 iter_sampling，采样的链条数 chains 及并行时 分配的 CPU 核心数 parallel_chains ，随机数种子 seed 。
library(cmdstanr)
nchains <- 2 # 2 条迭代链
# 给每条链设置相同的参数初始值
inits_data_gaussian <- lapply(1:nchains, function(i) {
list(
mu = 1
)
})
fit_gaussian <- mod_gaussian$sample(
init = inits_data_gaussian, # 迭代初值
iter_warmup = 200, # 每条链初始化迭代次数
iter_sampling = 200, # 每条链采样迭代次数
chains = nchains, # 马尔科夫链的数目
parallel_chains = nchains,# 指定 CPU 核心数，可以给每条链分配一个
seed = 20232023 # 设置随机数种子，不要使用 set.seed() 函数
)Running MCMC with 2 parallel chains...
Chain 1 Iteration: 1 / 400 [ 0%] (Warmup)
Chain 1 Iteration: 100 / 400 [ 25%] (Warmup)
Chain 1 Iteration: 200 / 400 [ 50%] (Warmup)
Chain 1 Iteration: 201 / 400 [ 50%] (Sampling)
Chain 1 Iteration: 300 / 400 [ 75%] (Sampling)
Chain 1 Iteration: 400 / 400 [100%] (Sampling)
Chain 2 Iteration: 1 / 400 [ 0%] (Warmup)
Chain 2 Iteration: 100 / 400 [ 25%] (Warmup)
Chain 2 Iteration: 200 / 400 [ 50%] (Warmup)
Chain 2 Iteration: 201 / 400 [ 50%] (Sampling)
Chain 2 Iteration: 300 / 400 [ 75%] (Sampling)
Chain 2 Iteration: 400 / 400 [100%] (Sampling)
Chain 1 finished in 0.0 seconds.
Chain 2 finished in 0.0 seconds.
Both chains finished successfully.
Mean chain execution time: 0.0 seconds.
Total execution time: 0.3 seconds.默认情况下，采样过程中会输出一些信息，以上是 2 条链并行采样的过程，给出百分比进度及时间消耗。采样完成后，调用方法 summary() 汇总和展示采样结果。
# A tibble: 2 × 10
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <num> <num> <num> <num> <num> <num> <num> <num> <num>
1 lp__ -14.4 -14.2 0.708 0.212 -15.6 -14.0 1.01 177. 220.
2 mu 0.0365 0.0307 0.348 0.298 -0.511 0.572 1.01 229. 173.输出模型中各个参数的后验分布的一些统计量，如均值（mean）、中位数（median）、标准差（sd），0.05 分位点（q5），0.95 分位点（q95）等。此外，还有 lp__ 后验对数概率密度值，每个模型都会有该值。summary() 方法有一些参数可以控制数字的显示方式和精度。下面展示的是保留 4 位有效数字的结果。
# A tibble: 2 × 10
variable mean median sd mad q5 q95 rhat ess_bulk
<chr> <dec:4> <dec:4> <dec:4> <dec:4> <dec:4> <dec:4> <dec:> <dec:4>
1 lp__ -14.43 -14.15 0.7083 0.2118 -15.60 -14.00 1.007 177.2
2 mu 0.03647 0.03073 0.3476 0.2978 -0.5115 0.5722 1.007 229.0
# ℹ 1 more variable: ess_tail <dec:4>接下来，要介绍 Stan 代码中的保留字 target 的含义，因为它在 Stan 代码中很常见，与输出结果中的 lp__ 一行紧密相关。
-
lp__表示后验概率密度函数的对数。 - target 累加一些和参数无关的数不影响参数的估计，但影响
lp__的值。 - 抽样语句表示模型会扔掉后验概率密度函数的对数的常数项。
为此，不妨在之前的 Stan 代码的基础上添加两行，新的 Stan 代码如下：
接着，再次编译代码、采样，为了节约篇幅，设置两个参数 show_messages 和 refresh ，不显示中间过程和采样进度。其它参数设置不变，代码如下：
fit_gaussian <- mod_gaussian_target$sample(
init = inits_data_gaussian,
iter_warmup = 200,
iter_sampling = 200,
chains = nchains,
parallel_chains = nchains,
show_messages = FALSE, # 不显示中间过程
refresh = 0, # 不显示采样进度
seed = 20232023
)
fit_gaussian$summary(.num_args = list(sigfig = 4, notation = "dec"))# A tibble: 2 × 10
variable mean median sd mad q5 q95 rhat
<chr> <dec:4> <dec:4> <dec:4> <dec:4> <dec:4> <dec:4> <dec:4>
1 lp__ 12335. 12335. 0.7074 0.1483 12334. 12336. 1.008
2 mu 0.03647 0.03073 0.3476 0.2978 -0.5115 0.5722 1.007
# ℹ 2 more variables: ess_bulk <dec:4>, ess_tail <dec:4>可以清楚地看到 lp__ 的值发生了变化，而参数 mu 的值没有变化。这是因为抽样语句 y ~ normal(mu, 1); 隐含一个 lp__ ，target 指代 lp__ 的值，符号 += 表示累加。两次累加后得到 12335.09。
下面从概率密度函数出发，用 R 语言来计算逐点对数似然函数值。一般地，不妨设 \(x_1,x_2,\cdots,x_n\) 是来自正态总体 \(\mathcal{N}(\mu,1)\) 的一个样本。则正态分布的概率密度函数 \(f(x)\) 的对数如下：
\[ \log f(x) = \log \frac{1}{\sqrt{2\pi}} - \frac{(x - \mu)^2}{2} \]
已知参数 \(\mu\) 是一个非常接近 0 的数，不妨将 \(\mu = 0\) 代入计算。
去掉常数项后，计算概率密度函数值的对数和。
这就比较接近原 lp__ 的值了，所以，lp__ 表示后验概率密度函数的对数，扔掉了与参数无关的常数项。若以概率密度函数的对数 normal_lpdf 替代抽样语句，则常数项是保留的。normal_lpdf 是 Stan 内置的函数，输入值为随机变量的取值 y 、位置参数 mu 和尺度参数 sigma，返回值为 real 实数。
real normal_lpdf(reals y | reals mu, reals sigma)
接着，编译上述代码以及重复采样的步骤，参数设置也一样。
fit_gaussian <- mod_gaussian_lpdf$sample(
init = inits_data_gaussian,
iter_warmup = 200,
iter_sampling = 200,
chains = nchains,
parallel_chains = nchains,
show_messages = FALSE,
refresh = 0,
seed = 20232023
)
fit_gaussian$summary(.num_args = list(sigfig = 4, notation = "dec"))# A tibble: 2 × 10
variable mean median sd mad q5 q95 rhat ess_bulk
<chr> <dec:4> <dec:4> <dec:4> <dec:4> <dec:4> <dec:4> <dec:> <dec:4>
1 lp__ -20.86 -20.58 0.7083 0.2119 -22.03 -20.43 1.007 176.7
2 mu 0.03647 0.03073 0.3476 0.2978 -0.5115 0.5722 1.007 229.0
# ℹ 1 more variable: ess_tail <dec:4>可以看到，此时 lp__ 的值包含常数项，两种表示方式对参数的计算结果没有影响。
31.1.2 Stan 的变量类型
变量的声明没有太多的内涵，就是 C++ 和 Stan 定义的语法，比如整型用 int 声明。建模过程中，时常需要将 R 语言环境中的数据传递给 Stan 代码编译出来的模型，而 Stan 是基于 C++ 语言，在变量类型方面有继承有发展。下表给出 Stan 与 R 语言中的变量类型对应关系。值得注意， R 语言的类型检查是不严格的，使用变量也不需要提前声明和初始化。Stan 语言中向量、矩阵的类型都是实数，下标也从 1 开始，元组类型和 R 语言中的列表类似，所有向量默认都是列向量。
| 变量类型 | Stan 语言 | R 语言 |
|---|---|---|
| 整型 | int x = 1; |
x = 1L |
| 实数 | real x = 3.14; |
x = 3.14 |
| 向量 | vector[3] x = [1, 2, 3]'; |
x = c(1L, 2L, 3L) |
| 矩阵 | matrix[3,1] x; |
matrix(data = c(1L, 2L, 3L)) |
| 数组 | array[3] int x; |
array(data = c(1L, 2L, 3L), dim = c(3, 1, 1)) |
| 元组 | tuple(vector[3]) x; |
list(x = c(1L, 2L, 3L)) |
31.1.3 Stan 的代码结构
functions {
// ... function declarations and definitions ...
}
data {
// ... declarations ...
}
transformed data {
// ... declarations ... statements ...
}
parameters {
// ... declarations ...
}
transformed parameters {
// ... declarations ... statements ...
}
model {
// ... declarations ... statements ...
}
generated quantities {
// ... declarations ... statements ...
}31.2 广义线性模型
31.2.1 生成模拟数据
先介绍泊松广义线性模型，包括模拟和计算，并和 Stan 实现的结果比较。
泊松广义线性模型如下：
\[ \begin{aligned} \log(\lambda) &= \beta_0 + \beta_1 x_1 + \beta_2 x_2 \\ Y &\sim \mathrm{Poisson}(u\lambda) \end{aligned} \]
设定参数向量 \(\beta = (\beta_0, \beta_1, \beta_2) = (0.5, 0.3, 0.2)\)，观测变量 \(X_1\) 和 \(X_2\) 的均值都为 0，协方差矩阵 \(\Sigma\) 为
\[ \left[ \begin{matrix} 1.0 & 0.8 \\ 0.8 & 1.0 \end{matrix} \right] \]
模拟观测到的响应变量值和协变量值，添加漂移项
31.2.2 拟合泊松模型
拟合泊松回归模型
fit_poisson_glm <- glm(y ~ X, family = poisson(link = "log"), offset = log(u))
summary(fit_poisson_glm)
Call:
glm(formula = y ~ X, family = poisson(link = "log"), offset = log(u))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.488932 0.009427 51.86 <2e-16 ***
X1 0.289984 0.014298 20.28 <2e-16 ***
X2 0.214846 0.014420 14.90 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 6052.9 on 2499 degrees of freedom
Residual deviance: 2675.5 on 2497 degrees of freedom
AIC: 10773
Number of Fisher Scoring iterations: 4# 对数似然函数值
log_poisson_lik <- logLik(fit_poisson_glm)
# 计算 AIC AIC(fit_poisson_glm)
-2 * c(log_poisson_lik) + 2 * attr(log_poisson_lik, "df")[1] 10772.79下面用 Stan 编码泊松回归模型，模型代码如下：
data {
int<lower=1> k;
int<lower=0> n;
matrix[n, k] X;
array[n] int<lower=0> y;
vector[n] log_offset;
}
parameters {
vector[k] beta;
real alpha;
}
model {
target += std_normal_lpdf(beta);
target += std_normal_lpdf(alpha);
target += poisson_log_glm_lpmf(y | X, alpha + log_offset, beta);
}
generated quantities {
vector[n] log_lik; // pointwise log-likelihood for LOO
vector[n] y_rep; // replications from posterior predictive dist
for (i in 1 : n) {
real y_hat_i = alpha + X[i] * beta + log_offset[i];
log_lik[i] = poisson_log_lpmf(y[i] | y_hat_i);
y_rep[i] = poisson_log_rng(y_hat_i);
}
}Stan 代码主要分三部分：
数据部分
data：声明模型的输入数据，数据类型、大小、约束。参数部分
parameters：类似数据部分，声明模型的参数，参数类型、大小。模型部分
model：指定模型参数的先验分布。生成量
generated quantities：拟合模型获得参数估计值后，计算一些统计量。
下面准备数据
编译模型，抽样获取参数的后验分布
# 加载 cmdstanr 包
library(cmdstanr)
# 编译模型
mod_poisson <- cmdstan_model(
stan_file = "code/poisson_log_glm.stan",
compile = TRUE,
cpp_options = list(stan_threads = TRUE)
)
# 采样拟合模型
fit_poisson_stan <- mod_poisson$sample(
data = poisson_d, # 观测数据
init = inits_data, # 迭代初值
iter_warmup = 1000, # 每条链预处理迭代次数
iter_sampling = 2000, # 每条链总迭代次数
chains = nchains, # 马尔科夫链的数目
parallel_chains = 1, # 指定 CPU 核心数，可以给每条链分配一个
threads_per_chain = 1, # 每条链设置一个线程
show_messages = FALSE, # 不显示迭代的中间过程
refresh = 0, # 不显示采样的进度
seed = 20222022 # 设置随机数种子，不要使用 set.seed() 函数
)
# 迭代诊断
fit_poisson_stan$diagnostic_summary()$num_divergent
[1] 0 0 0 0
$num_max_treedepth
[1] 0 0 0 0
$ebfmi
[1] 1.172526 1.114983 1.010774 1.125754# A tibble: 4 × 10
variable mean median sd mad q5 q95 rhat ess_bulk
<chr> <num> <num> <num> <num> <num> <num> <num> <num>
1 alpha 0.489 0.489 0.00948 0.00953 0.473 5.04e-1 1.00 4448.
2 beta[1] 0.290 0.290 0.0143 0.0144 0.267 3.14e-1 1.00 3144.
3 beta[2] 0.214 0.214 0.0145 0.0147 0.191 2.38e-1 1.00 3133.
4 lp__ -5388. -5388. 1.20 1.02 -5390. -5.39e+3 1.00 3232.
# ℹ 1 more variable: ess_tail <num>31.2.3 参数后验分布
加载 bayesplot 包，bayesplot 包提供一系列描述数据分布的绘图函数，比如绘制散点图 mcmc_scatter() 。\(\beta_1\) 和 \(\beta_2\) 的联合分布
library(ggplot2)
library(bayesplot)
mcmc_scatter(fit_poisson_stan$draws(c("beta[1]", "beta[2]"))) +
theme_classic() +
labs(x = expression(beta[1]), y = expression(beta[2]))
如果提取采样的数据，也可使用 ggplot2 包绘图，不局限于 bayesplot 设定的风格。
beta_df <- fit_poisson_stan$draws(c("beta[1]", "beta[2]"), format = "draws_df")
ggplot(data = beta_df, aes(x = `beta[1]`, y = `beta[2]`)) +
geom_density_2d_filled() +
facet_wrap(~.chain, ncol = 2) +
theme_classic() +
labs(x = expression(beta[1]), y = expression(beta[2]))
\(\beta_1\) 和 \(\beta_2\) 的热力图
mcmc_hex(fit_poisson_stan$draws(c("beta[1]", "beta[2]"))) +
theme_classic() +
labs(x = expression(beta[1]), y = expression(beta[2]))
各个参数的轨迹图
mcmc_trace(fit_poisson_stan$draws(c("beta[1]", "beta[2]")),
facet_args = list(
labeller = ggplot2::label_parsed, strip.position = "top", ncol = 1
)
) +
theme_classic()
可以将模型参数的后验分布图展示出来
mcmc_dens(fit_poisson_stan$draws(c("beta[1]", "beta[2]")),
facet_args = list(
labeller = ggplot2::label_parsed, strip.position = "top", ncol = 1
)
) +
theme_classic()
岭线图就是将各个参数的后验分布图放在一起。
mcmc_areas_ridges(x = fit_poisson_stan$draws(), pars = c("beta[1]", "beta[2]")) +
scale_y_discrete(labels = scales::parse_format()) +
theme_classic()
参数的 \(\hat{R}\) 潜在尺度收缩因子
后验预测诊断的想法是检查根据拟合模型生成的随机数 \(y^{rep}\) 与真实观测数据 \(y\) 的接近程度。为直观起见，可以用一系列描述数据分布的图来可视化检验。
y 是真实数据，yrep 是根据贝叶斯拟合模型生成的数据。下图是真实数据的密度图和50组生成数据的密度图。
31.2.4 模型评估指标
loo 包可以计算 WAIC
fit_poisson_waic <- loo::waic(fit_poisson_stan$draws(variables = "log_lik"))
print(fit_poisson_waic)
Computed from 8000 by 2500 log-likelihood matrix
Estimate SE
elpd_waic -5386.4 37.7
p_waic 2.9 0.1
waic 10772.8 75.5loo 包推荐使用 LOO-CV ，它还提供诊断信息、有效样本量和蒙特卡罗估计。
Computed from 8000 by 2500 log-likelihood matrix
Estimate SE
elpd_loo -5386.4 37.7
p_loo 2.9 0.1
looic 10772.8 75.5
------
Monte Carlo SE of elpd_loo is 0.0.
All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.31.2.5 可选替代实现
对于常见的统计模型，rstanarm 和 brms 包都内置了预编译的 Stan 程序，下面用 brms 包的函数 brm() 拟合带上述漂移项的泊松广义线性模型，参数估计结果和 Base R 函数 glm() 的几乎一致，因编译和抽样的过程比较花费时间，速度不及 Base R。
# brms
dat <- data.frame(y = y, X = X, u = u)
colnames(dat) <- c("y", "x1", "x2", "u")
fit_poisson_brm <- brms::brm(y ~ x1 + x2 + offset(log(u)),
data = dat, family = poisson(link = "log")
)
fit_poisson_brm Family: poisson
Links: mu = log
Formula: y ~ x1 + x2 + offset(log(u))
Data: dat (Number of observations: 2500)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.49 0.01 0.47 0.51 1.00 2509 2171
x1 0.29 0.01 0.26 0.32 1.00 1771 1645
x2 0.21 0.01 0.19 0.24 1.00 1727 1847
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).调用函数 brm() 拟合模型后返回一个 brmsfit 对象 fit_poisson_brm，brms 包提供很多函数处理该数据对象，比如 brms::loo() 计算 LOO-CV
Computed from 4000 by 2500 log-likelihood matrix
Estimate SE
elpd_loo -5386.3 37.8
p_loo 2.9 0.1
looic 10772.6 75.5
------
Monte Carlo SE of elpd_loo is 0.0.
All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.输出结果中， LOO IC 信息准则 Loo information criterion，looic 指标的作用类似频率派模型中的 AIC 指标，所以也几乎相同的。
31.3 选择先验分布
响应变量服从伯努利分布的广义线性模型有 11 个协变量（含截距），实际有用的是 3 个变量，其它变量应该被收缩掉。贝叶斯收缩 （bayesian shrinkage）与变量选择是有关系的，先验分布影响收缩的力度。下面模拟生成 2500 个样本，非 0 的回归系数分别是 \(\alpha = 1,\beta_1 = 3,\beta_2 = -2\) 。
31.3.1 正态先验
data {
int<lower=1> k;
int<lower=0> n;
matrix[n, k] X;
array[n] int<lower=0, upper=1> y;
}
parameters {
vector[k] beta;
real alpha;
}
model {
target += normal_lpdf(beta | 0, 1000);
target += normal_lpdf(alpha | 0, 1000);
target += bernoulli_logit_glm_lpmf(y | X, alpha, beta);
}
generated quantities {
vector[n] log_lik;
for (i in 1 : n) {
log_lik[i] = bernoulli_logit_lpmf(y[i] | alpha + X[i] * beta);
}
}mod_logit_normal <- cmdstan_model(
stan_file = "code/bernoulli_logit_glm_normal.stan",
compile = TRUE, cpp_options = list(stan_threads = TRUE)
)
fit_logit_normal <- mod_logit_normal$sample(
data = mdata,
chains = 2,
parallel_chains = 2,
iter_warmup = 1000,
iter_sampling = 1000,
threads_per_chain = 2,
seed = 20232023,
show_messages = FALSE,
refresh = 0
)
# 输出结果
fit_logit_normal$summary(c("alpha", "beta", "lp__"))# A tibble: 12 × 10
variable mean median sd mad q5 q95 rhat ess_bulk
<chr> <num> <num> <num> <num> <num> <num> <num> <num>
1 alpha 1.02 1.02 0.0710 0.0699 0.902 1.13e+0 1.00 2648.
2 beta[1] 3.14 3.13 0.136 0.139 2.92 3.37e+0 1.00 2054.
3 beta[2] -2.03 -2.03 0.0994 0.0992 -2.19 -1.87e+0 1.00 2041.
4 beta[3] 0.0604 0.0609 0.0663 0.0667 -0.0455 1.70e-1 1.00 3339.
5 beta[4] -0.0266 -0.0244 0.0665 0.0655 -0.140 8.30e-2 0.999 3548.
6 beta[5] 0.0146 0.0174 0.0674 0.0689 -0.0969 1.24e-1 1.00 4352.
7 beta[6] 0.0215 0.0231 0.0655 0.0652 -0.0862 1.26e-1 1.00 3111.
8 beta[7] 0.103 0.105 0.0611 0.0600 0.00131 2.02e-1 1.00 3377.
9 beta[8] -0.0303 -0.0312 0.0649 0.0600 -0.138 7.59e-2 1.00 3577.
10 beta[9] -0.0882 -0.0873 0.0660 0.0666 -0.194 2.09e-2 1.00 3487.
11 beta[10] 0.0814 0.0803 0.0647 0.0653 -0.0237 1.86e-1 1.00 4377.
12 lp__ -843. -842. 2.38 2.25 -847. -8.39e+2 1.00 800.
# ℹ 1 more variable: ess_tail <num>31.3.2 Lasso 先验
data {
int<lower=1> k;
int<lower=0> n;
matrix[n, k] X;
array[n] int<lower=0, upper=1> y;
}
parameters {
vector[k] beta;
real alpha;
real<lower=0> lambda;
}
model {
target += double_exponential_lpdf(beta | 0, lambda);
target += double_exponential_lpdf(alpha | 0, lambda);
target += cauchy_lpdf(lambda | 0, 0.01);
target += bernoulli_logit_glm_lpmf(y | X, alpha, beta);
}
generated quantities {
vector[n] log_lik;
for (i in 1 : n) {
log_lik[i] = bernoulli_logit_lpmf(y[i] | alpha + X[i] * beta);
}
}mod_logit_lasso <- cmdstan_model(
stan_file = "code/bernoulli_logit_glm_lasso.stan",
compile = TRUE, cpp_options = list(stan_threads = TRUE)
)
fit_logit_lasso <- mod_logit_lasso$sample(
data = mdata,
chains = 2,
parallel_chains = 2,
iter_warmup = 1000,
iter_sampling = 1000,
threads_per_chain = 2,
seed = 20232023,
show_messages = FALSE,
refresh = 0
)
# 输出结果
fit_logit_lasso$summary(c("alpha", "beta", "lambda", "lp__"))# A tibble: 13 × 10
variable mean median sd mad q5 q95 rhat ess_bulk
<chr> <num> <num> <num> <num> <num> <num> <num> <num>
1 alpha 0.993 0.994 0.0756 0.0782 0.871 1.12e+0 1.00 2117.
2 beta[1] 3.08 3.08 0.135 0.137 2.87 3.31e+0 1.00 1910.
3 beta[2] -1.99 -1.99 0.0975 0.0949 -2.15 -1.82e+0 1.00 1926.
4 beta[3] 0.0539 0.0508 0.0623 0.0645 -0.0430 1.59e-1 1.00 3336.
5 beta[4] -0.0223 -0.0212 0.0615 0.0638 -0.122 7.95e-2 1.00 3525.
6 beta[5] 0.0108 0.0105 0.0651 0.0647 -0.0977 1.20e-1 1.00 4587.
7 beta[6] 0.0191 0.0183 0.0592 0.0623 -0.0767 1.20e-1 1.00 3167.
8 beta[7] 0.0940 0.0937 0.0629 0.0621 -0.00795 1.95e-1 1.00 2785.
9 beta[8] -0.0256 -0.0251 0.0612 0.0629 -0.129 7.35e-2 1.00 3390.
10 beta[9] -0.0801 -0.0793 0.0621 0.0635 -0.182 2.06e-2 1.00 2988.
11 beta[10] 0.0743 0.0734 0.0622 0.0632 -0.0241 1.78e-1 1.00 3208.
12 lambda 0.599 0.564 0.201 0.162 0.355 9.73e-1 1.00 2921.
13 lp__ -775. -775. 2.44 2.42 -779. -7.72e+2 1.00 795.
# ℹ 1 more variable: ess_tail <num>计算 LOO-CV 比较正态先验和 Lasso 先验
fit_logit_normal_loo <- fit_logit_normal$loo(variables = "log_lik", cores = 2)
print(fit_logit_normal_loo)
Computed from 2000 by 2500 log-likelihood matrix
Estimate SE
elpd_loo -762.2 27.3
p_loo 11.2 0.5
looic 1524.4 54.5
------
Monte Carlo SE of elpd_loo is 0.1.
All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.fit_logit_lasso_loo <- fit_logit_lasso$loo(variables = "log_lik", cores = 2)
print(fit_logit_lasso_loo)
Computed from 2000 by 2500 log-likelihood matrix
Estimate SE
elpd_loo -761.6 26.8
p_loo 10.5 0.5
looic 1523.3 53.7
------
Monte Carlo SE of elpd_loo is 0.1.
All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.loo 包的函数 loo_compare() 比较两个模型
elpd_diff se_diff
model1 0.0 0.0
model0 -0.6 0.5 输出结果中最好的模型放在第一行。looic 越小越好，所以，Lasso 先验更好。
31.3.3 Horseshoe 先验
# horseshoe 先验
mod_logit_horseshoe <- cmdstan_model(
stan_file = "code/bernoulli_logit_glm_horseshoe.stan",
compile = TRUE, cpp_options = list(stan_threads = TRUE)
)
fit_logit_horseshoe <- mod_logit_horseshoe$sample(
data = mdata,
chains = 2,
parallel_chains = 2,
iter_warmup = 1000,
iter_sampling = 1000,
threads_per_chain = 2,
seed = 20232023,
show_messages = FALSE,
refresh = 0
)
fit_logit_horseshoe$summary(c("alpha", "beta", "lambda", "lp__"))
注释
brms 包在模型中全局性地使用 horseshoe 先验信息。
31.3.4 SpikeSlab 先验
31.4 选择推理算法
开篇提及 Stan 内置了多种推理算法，不同的算法获得的结果是存在差异的。
- full Bayesian statistical inference with MCMC sampling (NUTS, HMC)
- approximate Bayesian inference with variational inference (ADVI)
- penalized maximum likelihood estimation with optimization (L-BFGS)
31.4.1 惩罚极大似然算法
L-BFGS 算法拟合模型，速度非常快。
# L-BFGS 算法拟合模型
fit_optim_logit <- mod_logit_lasso$optimize(
data = mdata, # 观测数据
init = 0, # 所有参数初值设为 0
refresh = 0, # 不显示迭代进程
algorithm = "lbfgs", # 优化器
threads = 1, # 单线程
seed = 20232023 # 随机数种子
)Finished in 0.2 seconds.# A tibble: 13 × 2
variable estimate
<chr> <num>
1 alpha 0.981
2 beta[1] 3.05
3 beta[2] -1.96
4 beta[3] 0.0488
5 beta[4] -0.0166
6 beta[5] 0.00528
7 beta[6] 0.0126
8 beta[7] 0.0923
9 beta[8] -0.0204
10 beta[9] -0.0777
11 beta[10] 0.0721
12 lambda 0.488
13 lp__ -768. 31.4.2 变分近似推断算法
ADVI 算法拟合模型，可选的优化器有 meanfield 和 fullrank ，相比于 L-BFGS 稍慢
# ADVI 算法拟合模型
fit_advi_logit <- mod_logit_lasso$variational(
data = mdata, # 观测数据
init = 0, # 所有参数初值设为 0
refresh = 0, # 不显示迭代进程
algorithm = "meanfield", # 优化器
threads = 1, # 单线程
seed = 20232023 # 随机数种子
)Finished in 1.3 seconds.# A tibble: 13 × 7
variable mean median sd mad q5 q95
<chr> <num> <num> <num> <num> <num> <num>
1 alpha 1.02 1.02 0.0615 0.0630 0.914 1.11
2 beta[1] 3.07 3.07 0.0899 0.0870 2.93 3.22
3 beta[2] -1.98 -1.98 0.0675 0.0666 -2.08 -1.86
4 beta[3] 0.0161 0.0159 0.0678 0.0670 -0.0945 0.129
5 beta[4] -0.0199 -0.0221 0.0639 0.0621 -0.121 0.0857
6 beta[5] -0.00128 -0.00202 0.0722 0.0713 -0.116 0.121
7 beta[6] -0.0423 -0.0446 0.0705 0.0689 -0.156 0.0754
8 beta[7] 0.0760 0.0750 0.0490 0.0489 -0.00517 0.152
9 beta[8] -0.0742 -0.0752 0.0659 0.0637 -0.181 0.0347
10 beta[9] -0.0495 -0.0501 0.0802 0.0818 -0.185 0.0805
11 beta[10] 0.0444 0.0440 0.0520 0.0522 -0.0444 0.128
12 lambda 0.698 0.662 0.243 0.228 0.379 1.13
13 lp__ -777. -777. 3.39 3.22 -783. -773. 31.4.3 拉普拉斯近似算法
CmdStan 已经集成，但 cmdstanr 包对该算法的支持还在开发中。
31.5 习题
-
分析挑战者号航天飞机 O 型环数据。DAAG 包的 orings 数据集记录美国挑战者号航天飞机 O 型环在不同温度下发生 Erosion 腐蚀和 Blowby 串气的失效数量。 图 31.10 展示航天飞机 O 型环在不同温度下失效的分布图（条件密度图）：随着温度升高，O 型环越来越不容易失效。请分别用 Base R 函数
glm()和 cmdstanr 包建模分析 O 型环数据。代码
# data(orings, package = "DAAG") orings <- readRDS(file = "data/orings.rds") ggplot(orings, aes(x = Temperature, y = after_stat(count))) + geom_density(aes(fill = Total > 0), position = "fill", bw = 2) + scale_y_continuous(labels = scales::label_percent()) + scale_fill_grey(labels = c("TRUE" = "是", "FALSE" = "否")) + theme_classic() + labs(x = "温度", y = "比例", fill = "失效")
图 31.10: 航天飞机 O 型环在不同温度下失效的条件密度图 根据 章节 27 的数据，建立贝叶斯空间广义线性混合模型，用 Stan 预测核辐射强度的分布。