probabilistic-reasoning-framework.qmd

# 概率推理框架 {#sec-probabilistic-reasoning-framework}

```{r}
#| echo: false

Sys.setenv(CMDSTANR_NO_VER_CHECK = TRUE)
source("_common.R")
```

::: hidden
$$
 \def\bm#1{{\boldsymbol #1}}
$$
:::

本章的目的是让读者快速熟悉和上手，主要分为以下几个部分。

1.  Stan 的概览，介绍 Stan 是什么，怎么样。
2.  Stan 的入门，以推理一个正态分布均值参数为例，从基础语法、类型声明和代码结构三个方面介绍 Stan 的使用。
3.  选择先验分布，先验分布在贝叶斯推理中的重要性不言而喻，本节以一个简单的广义线性模型为例，介绍常见的几个先验分布对模型的影响。
4.  选择推理算法，接上一节的例子，围绕怎么用、效果如何介绍 Stan 内置的几个推理算法。

## Stan 概览 {#sec-stan-overview}

[Stan](https://github.com/stan-dev) 是一个贝叶斯统计建模和计算的概率推理框架，也是一门用于贝叶斯推断和优化的概率编程语言 [@Gelman2015; @Carpenter2017]。它使用汉密尔顿蒙特卡罗算法（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 已经形成一个相对成熟的生态，它提供统计建模、数据分析和预测能力，广泛应用于社会、生物、物理、工程、商业等领域，在学术界和工业界的影响力也不小。下 @fig-stan-api 是 Stan 生态中各组件依赖架构图，[math](https://github.com/stan-dev/math) 库[@Carpenter2015]是 Stan 框架最核心的组件，它基于 [Boost](https://github.com/boostorg/boost) 、[Eigen](https://gitlab.com/libeigen/eigen) 、[OpenCL](https://www.khronos.org/opencl/) 、[SUNDIALS](https://github.com/LLNL/sundials) 和 [oneTBB](https://github.com/oneapi-src/oneTBB) 等诸多 C++ 库，提供概率推理、自动微分、矩阵计算、并行计算、GPU 计算和求解代数微分方程等功能。

```{mermaid}
%%| label: fig-stan-api
%%| fig-width: 6.5
%%| fig-cap: Stan、CmdStan 和 CmdStanR 等的依赖关系图

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)
  CmdStanR --> bayesplot(bayesplot)
  CmdStanR --> loo(loo)
  CmdStanR --> posterior(posterior)
  CmdStanR --> projpred(projpred)
```

CmdStan 是 Stan 的命令行接口，可在 MacOS / Linux 的终端软件，Windows 的命令行窗口或 PowerShell 软件中使用。**CmdStanR** [@Gabry2023]**、**CmdStanPy 和 MathematicaStan 分别是 CmdStan 的 R 语言、Python 语言和 Mathematica 语言接口。每次当 Stan 发布新版本时，CmdStan 也会随之发布新版，只需指定新的 CmdStan 安装路径，**CmdStanR** 就可以使用上，**CmdStanR** 包与 Stan 是相互独立的更新机制。 **CmdStanR** 负责处理 CmdStan 运行的结果，而编译代码，生成模型和模拟采样等都是由 **CmdStan 完成**。入门 **CmdStanR** 后，可以快速转入对 Stan 底层原理的学习，有利于编码符合实际需要的复杂模型，有利于掌握常用的炼丹技巧，提高科研和工作的效率。

此外，**bayesplot** 包 [@Gabry2019] 针对 cmdstanr 包生成的拟合模型对象提供一系列可视化图形，用于诊断采样过程、展示后验分布等。**loo** 包[@Vehtari2017]计算 LOOIC （留一交叉验证信息准则）和 WAIC （通用信息准则）等指标，用于模型评估与比较。**posterior** 包 [@Vehtari2021] 对采样数据提供统一的操作方法和类型转化，计算常用的后验分布的统计量等。**projpred** 包 [@Piironen2017b; @Piironen2020] 实现投影预测推断用于模型预测和特征选择。

[**rstan**](https://github.com/stan-dev/rstan) 包[@RStan2023]是 Stan 的 R 语言接口，该接口依赖 **Rcpp** [@Rcpp2011; @Rcpp2018]、**RcppEigen** [@Bates2013]、**BH** [@BH2023]、**RcppParallel** [@RcppParallel2023]和 **StanHeaders** [@StanHeaders2023]等 R 包，由于存在众多上游 R 包依赖和兼容性问题，尤其在 Windows 系统环境中，因此，**RStan 的**安装、更新都比较麻烦。**RStan** 的更新通常严重滞后于 Stan 的更新，不利于及时地使用最新的学术研究成果。 而相比于 **rstan** 包，**CmdStanR** 更加轻量，可以更快地将 CmdStan 的新功能融入进来，而且 **cmdstanr** 和 CmdStan 是分离的，方便用户升级和维护。

[**rstanarm**](https://github.com/stan-dev/rstanarm) [@Goodrich2023] 和 [**brms**](https://github.com/paul-buerkner/brms) [@brms2017] 是 **RStan** 的扩展包，各自提供了一套用于表示统计模型的公式语法。它们都支持丰富的统计模型，比如线性模型、广义线性模型、线性混合效应模型、广义线性混合效应模型等。相比于 **rstan**， 它们使用起来更加方便，因为它内置了大量统计模型的 Stan 实现，即将公式语法翻译成 Stan 编码的模型，然后调用 **rstan** 或 **cmdstanr** 翻译成 C++，最后编译成动态链接库。除了依赖 **rstan** 包，**rstanarm** 和 **brms** 还依赖大量其它 R 包。

顺便一提，类似的用于概率推理和统计分析的框架，还有 Python 社区的 [PyMC](https://github.com/pymc-devs/pymc) [@pymc2023]和 [TensorFlow Probability](https://github.com/tensorflow/probability) [@Dillon2017]，它们采用的 MCMC 采样算法也是基于 NUTS 的 HMC 算法。

## Stan 入门 {#sec-getting-started}

### Stan 的基础语法 {#sec-stan-syntax}

下面以一个简单示例介绍 Stan 的用法，包括 Stan 的基本用法、变量类型、代码结构等，

```{r}
#| echo: false

# 注册 Stan 引擎替换 Quarto 文档中默认的 Stan 块
# 原 Stan 块的编译采用 rstan 包
# eng_cmdstan 不支持传递函数 cmdstan_model 的其他参数选项
knitr::knit_engines$set(stan = cmdstanr::eng_cmdstan)
```

考虑一个已知方差的正态分布，设 $-3, -2, -1, 0, 1, 2, 3$ 是取自正态分布 $\mathcal{N}(\mu,1)$ 的一个样本，也是取自该正态分布的一组随机数。现在的问题是估计该正态分布的均值参数 $\mu$ 。Stan 编码的正态分布模型如下：

```{stan output.var="mod_gaussian"}
transformed data {
  vector[7] y = [-3, -2, -1, 0, 1, 2, 3]';
}
parameters {
  real mu;
}
model {
  y ~ normal(mu, 1);
}
```

-   `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` 。

```{r}
#| message: false

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() 函数
)
```

默认情况下，采样过程中会输出一些信息，以上是 2 条链并行采样的过程，给出百分比进度及时间消耗。采样完成后，调用方法 `summary()` 汇总和展示采样结果。

```{r}
fit_gaussian$summary()
```

输出模型中各个参数的后验分布的一些统计量，如均值（mean）、中位数（median）、标准差（sd），0.05 分位点（q5），0.95 分位点（q95）等。此外，还有 `lp__` 后验对数概率密度值，每个模型都会有该值。`summary()` 方法有一些参数可以控制数字的显示方式和精度。下面展示的是保留 4 位有效数字的结果。

```{r}
fit_gaussian$summary(.num_args = list(sigfig = 4, notation = "dec"))
```

接下来，要介绍 Stan 代码中的保留字 target 的含义，因为它在 Stan 代码中很常见，与输出结果中的 `lp__` 一行紧密相关。

-   `lp__` 表示后验概率密度函数的对数。
-   target 累加一些和参数无关的数不影响参数的估计，但影响 `lp__` 的值。
-   抽样语句表示模型会扔掉后验概率密度函数的对数的常数项。

```{r}
#| label: fig-stan-lp
#| fig-cap: lp__ 的后验分布
#| fig-showtext: true
#| fig-width: 5
#| fig-height: 4
#| message: false

library(ggplot2)
library(bayesplot)
mcmc_hist(fit_gaussian$draws("lp__")) +
  theme_classic()
```

为此，不妨在之前的 Stan 代码的基础上添加两行，新的 Stan 代码如下：

```{stan output.var="mod_gaussian_target"}
transformed data {
  vector[7] y = [-3, -2, -1, 0, 1, 2, 3]';
}
parameters {
  real mu;
}
model {
  y ~ normal(mu, 1);
  target += 12345;
  target += mean(exp(y));
}
```

接着，再次编译代码、采样，为了节约篇幅，设置两个参数 `show_messages` 和 `refresh` ，不显示中间过程和采样进度。其它参数设置不变，代码如下：

```{r}
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"))
```

可以清楚地看到 `lp__` 的值发生了变化，而参数 `mu` 的值没有变化。这是因为抽样语句 `y ~ normal(mu, 1);` 隐含一个 `lp__` ，target 指代 `lp__` 的值，符号 `+=` 表示累加。两次累加后得到 12335.09。

``` stan
model {
  y ~ normal(mu, 1);
  target += 12345;
  target += mean(exp(y));
}
```

```{r}
y <- c(-3, -2, -1, 0, 1, 2, 3)
12345 + mean(exp(y)) - 14.45 
```

下面从概率密度函数出发，用 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$ 代入计算。

```{r}
sum(dnorm(x = y, mean = 0, sd = 1, log = TRUE))
```

去掉常数项后，计算概率密度函数值的对数和。

```{r}
# 扔掉常数
f <- function(y, mu) {
  return(-0.5 * (y - mu)^2)
}
sum(f(-3:3, 0))
```

这就比较接近原 `lp__` 的值了，所以，`lp__` 表示后验概率密度函数的对数，扔掉了与参数无关的常数项。若以概率密度函数的对数 `normal_lpdf` 替代抽样语句，则常数项是保留的。`normal_lpdf` 是 Stan 内置的函数，输入值为随机变量的取值 `y` 、位置参数 `mu` 和尺度参数 `sigma`，返回值为 `real` 实数。

`real` **`normal_lpdf`**`(reals y | reals mu, reals sigma)`

```{stan output.var="mod_gaussian_lpdf"}
transformed data {
  vector[7] y = [-3, -2, -1, 0, 1, 2, 3]';
}
parameters {
  real mu;
}
model {
  target += normal_lpdf(y | mu, 1);
}
```

接着，编译上述代码以及重复采样的步骤，参数设置也一样。

```{r}
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"))
```

可以看到，此时 `lp__` 的值包含常数项，两种表示方式对参数的计算结果没有影响。

### Stan 的变量类型 {#sec-stan-variables}

Stan 语言和 C/C++ 语言比较类似，变量需要先声明再使用，函数需要用 `return` 返回值，总而言之，类型声明比较严格。变量的声明没有太多的内涵，就是 C++ 和 Stan 定义的语法，比如整型用 `int` 声明。建模过程中，时常需要将 R 语言环境中的数据传递给 Stan 代码编译出来的模型，而 Stan 是基于 C++ 语言，在变量类型方面有继承有发展。下表给出 Stan 与 R 语言中的变量类型对应关系。值得注意， R 语言的类型检查是不严格的，使用变量也不需要提前声明和初始化。Stan 语言中向量、矩阵的类型都是实数，下标也从 1 开始，元组类型和 R 语言中的列表类似，所有向量默认都是列向量。

下表第一列表示 Stan 语言的变量类型，第二列给出使用该变量的声明示例，第三列给出 R 语言中构造该类型变量的示例。

| 类型 | Stan 语言                       | R 语言                                          |
|--------|-------------------------------|-------------------------------------|
| 整型 | `int x = 1;`                    | `x = 1L`                                        |
| 实数 | `real x = 3.14;`                | `x = 3.14`                                      |
| 向量 | `vector[3] x = [1, 2, 3]';`     | `x = c(1, 2, 3)`                                |
| 矩阵 | `matrix[3,1] x;`                | `matrix(data = c(1, 2, 3), nrow = 3)`           |
| 数组 | `array[3] int x;`               | `array(data = c(1L, 2L, 3L), dim = c(3, 1, 1))` |
| 元组 | `tuple(vector[3],vector[3]) x;` | `list(x = c(1, 2, 3), y = c(4, 5, 6))`          |

: Stan 变量类型和 R 语言中的对应 {#tbl-stan-var-dec}

### Stan 的代码结构 {#sec-stan-code}

一般地，Stan 代码文件包含数据、参数和模型三块内容，一个简单的示例如 @sec-stan-syntax 所示。Stan 代码文件最多有如下 7 块内容，函数块 `functions` 放一些自定义的函数，数据变换块 `transformed data` 对输入数据做一些变换，预计算，以便放入后续模型，参数变换块 `transformed parameters` 作用类似数据变换块，方便在模型中使用，它们也会作为参数在输出结果中显示。生成量块 `generated quantities` 计算一些统计量，概率分布的随机数、分位数等。模拟、拟合和预测模型会用到其中的一部分或全部。

``` 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 ...
}
```

### Stan 的函数使用 {#sec-stan-function}

Stan 有大量的内建函数，然而，有时候，Stan 内建的函数不足以满足需求，需要自己创建函数。下面以函数 `cholesky_decompose` 为例介绍 Stan 内置/一般函数的调用，在该函数的基础上自定义函数 `cholesky_decompose2` ，这不过是对它改个名字，其它内容只要符合 Stan 语言即可，不甚重要。

根据 Stan 官网函数 `cholesky_decompose` 帮助文档，Cholesky 分解的形式（Cholesky 分解有多种形式）如下：

$$
M = LL^{\top}
$$

$M$ 是一个对称正定的矩阵，而 $L$ 是一个下三角矩阵。函数 `cholesky_decompose` 有一个参数 A， A 需要传递一个对称正定的矩阵。不妨设这个对称正定的矩阵为

$$
M = \begin{bmatrix}
4 & 1 \\
1 & 1
\end{bmatrix}
$$

```{r}
# 准备函数
stan_file <- write_stan_file("
functions {
 matrix cholesky_decompose2(matrix A) {
   return cholesky_decompose(A);
 }
}
parameters {
 real x;
}
model {
 x ~ std_normal();
}
")
```

接着，将以上 Stan 代码编译

```{r}
#| message: false

mod_cholesky_decompose <- cmdstan_model(stan_file = stan_file, compile = TRUE)
```

准备测试数据，只要是一个对称正定的矩阵都可以做 cholesky 分解。

```{r}
# 测试矩阵
M <- rbind(c(4, 1), c(1, 1))
```

**cmdstanr** 包导出函数的方法将以上 Stan 代码中的函数部分独立导出。

```{r}
#| message: false

# 编译独立的函数
mod_cholesky_decompose$expose_functions()
```

现在，可以直接调用导出的函数 `cholesky_decompose2` 。

```{r}
# cholesky 分解
mod_cholesky_decompose$functions$cholesky_decompose2(A = M)
```

最后，将 Stan 函数计算的结果与 R 语言内置的 cholesky 分解函数的结果比较。发现，函数 `chol()` 的结果正好是 `cholesky_decompose2` 的转置。

```{r}
chol(M)
```

查看帮助文档，可知 R 软件对 Cholesky 分解的定义如下：

$$
M = L^{\top}L
$$

根据数学表达式，感觉上都是差不多的，但还是有差异。R 与 Stan 混合编程就需要注意这些表达上不同的，不然，排错会很麻烦。

::: callout-tip
StanHeaders 可以编译和调用 Stan 的内置的数学函数，比如 Cholesky 分解函数 `cholesky_decompose` 。

```{r}
library(StanHeaders)
stanFunction("cholesky_decompose", A = M)
```

可以看到，结果和前面一样。
:::

## 先验分布 {#sec-choose-prior}

考虑一个响应变量服从伯努利分布的广义线性模型。

$$
\begin{aligned}
&\boldsymbol{y} \sim \mathrm{Bernoulli}(\boldsymbol{p}) \\
&\mathrm{logit}(\boldsymbol{p}) = \log (\frac{\boldsymbol{p}}{1-\boldsymbol{p}})=  \alpha + X \boldsymbol{\beta}
\end{aligned}
$$

下面模拟生成 2500 个样本，其中 10 个正态协变量，非 0 的回归系数是截距 $\alpha = 1$ 和向量 $\boldsymbol{\beta}$ 中的 $\beta_1 = 3,\beta_2 = -2$ 。对模型实际有用的是 3 个变量，采用贝叶斯建模，其它变量应该被收缩掉。贝叶斯收缩 （Bayesian shrinkage）与变量选择 （Variable selection） 是有关系的，先验分布影响收缩的力度。

```{r}
set.seed(2023)
n <- 2500
k <- 10
X <- matrix(rnorm(n * k), ncol = k)
y <- rbinom(n, size = 1, prob = plogis(1 + 3 * X[, 1] - 2 * X[, 2]))
# 准备数据
mdata <- list(k = k, n = n, y = y, X = X)
```

在贝叶斯先验分布中，有几个常用的概率分布，分别是正态分布、拉普拉斯分布（双指数分布）、柯西分布，下图集中展示了这几个的标准分布。

```{r}
#| label: fig-prior
#| fig-cap: 几个常用的概率分布
#| fig-width: 4.5
#| fig-height: 3.5
#| dev: 'tikz'
#| fig-process: !expr to_png
#| code-fold: true
#| echo: !expr knitr::is_html_output()

dlaplace <- function(x, mu = 0, sigma = 1) {
  1 / (2*sigma) * exp(- abs(x - mu) / sigma)
}

ggplot() +
  geom_function(
    fun = dnorm, args = list(mean = 0, sd = 1),
    aes(colour = "正态分布"), linewidth = 1.2, xlim = c(-6, 6)
  ) +
  geom_function(
    fun = dlaplace, args = list(mu = 0, sigma = 1),
    aes(colour = "双指数分布"), linewidth = 1.2, xlim = c(-6, 6)
  ) +
  geom_function(
    fun = dcauchy, args = list(location = 0, scale = 0.5),
    aes(colour = "柯西分布"), linewidth = 1.2, xlim = c(-6, 6)
  ) +
  theme_classic() +
  theme(legend.position = "inside", legend.position.inside = c(0.8, 0.8)) +
  labs(x = "$x$", y = "$f(x)$", colour = "先验分布")
```

接下来，考虑几种常见的先验设置。

### 正态先验 {#sec-prior-normal}

指定回归系数 $\alpha,\beta$ 的先验分布如下

$$
\begin{aligned}
\alpha &\sim \mathcal{N}(0, 1000) \\
\beta &\sim \mathcal{N}(0, 1000)
\end{aligned}
$$

正态分布中设置相当大的方差意味着分布相当扁平， $\alpha,\beta$ 的取值在区间 $(-\infty,+\infty)$ 上比较均匀。

```{verbatim, file="code/bernoulli_logit_glm_normal.stan", lang="stan"}
```

```{r}
#| message: false

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__"))
```

### Lasso 先验 {#sec-prior-lasso}

指定回归系数 $\alpha,\beta$ 的先验分布如下

$$
\begin{aligned}
\lambda &\sim \mathrm{Half\_Cauchy}(0,0.01) \\
\alpha  &\sim \mathrm{Double\_exponential}(0, \lambda) \\
\beta   &\sim \mathrm{Double\_exponential}(0, \lambda)
\end{aligned}
$$

其中， $\alpha,\beta$ 服从双指数分布，惩罚因子 $\lambda$ 服从柯西分布。顺便一提，若把双指数分布改为正态分布，则 Lasso 先验变为岭先验。相比于岭先验，Lasso 先验有意将回归系数往 0 上收缩，这非常类似于频率派中的岭回归与 Lasso 回归的关系 [@Bhadra2019]。

```{verbatim, file="code/bernoulli_logit_glm_lasso.stan", lang="stan"}
```

```{r}
#| message: false
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__"))
```

计算 LOO-CV 比较正态先验和 Lasso 先验

```{r}
fit_logit_normal_loo <- fit_logit_normal$loo(variables = "log_lik", cores = 1)
print(fit_logit_normal_loo)

fit_logit_lasso_loo <- fit_logit_lasso$loo(variables = "log_lik", cores = 1)
print(fit_logit_lasso_loo)
```

loo 包的函数 `loo_compare()` 比较两个模型

```{r}
loo::loo_compare(list(model0 = fit_logit_normal_loo, 
                      model1 = fit_logit_lasso_loo))
```

输出结果中最好的模型放在第一行。LOOIC 越小越好，所以，Lasso 先验更好。

### Horseshoe 先验 {#sec-prior-horseshoe}

Horseshoe 先验（Horse shoe）[@Piironen2017a] 指定回归系数 $\alpha,\bm{\beta}$ 的先验分布如下

$$
\begin{aligned}
\lambda_i &\sim \mathrm{Half\_Cauchy}(0,1) \\
\alpha | \lambda_0,\tau  &\sim \mathcal{N}(0, \tau^2\lambda_0^2) \\
\beta_i | \lambda_i,\tau  &\sim \mathcal{N}(0, \tau^2\lambda_i^2),\quad i = 1,2,\cdots,10
\end{aligned}
$$

其中，$\tau$ 称之为全局超参数，它将所有的回归系数朝着 0 收缩。而作用在局部超参数 $\lambda_i$ 上的重尾柯西先验允许某些回归系数逃脱收缩。

```{verbatim, file="code/bernoulli_logit_glm_horseshoe.stan", lang="stan"}
```

```{r}
#| message: false

# 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", "tau", "lambda", "lp__")) 
```

可以看到回归系数小的压缩效果很明显，而回归系数大的几乎没有压缩。

```{r}
fit_logit_horseshoe_loo <- fit_logit_horseshoe$loo(variables = "log_lik", cores = 1)
print(fit_logit_horseshoe_loo)
```

LOOIC 比之 Lasso 先验的情况更小了。

::: callout-note
```{r}
#| eval: false

library(rstanarm)
# set up the prior, use hyperprior tau ∼ half-Cauchy(0,tau0^2) 
D <- ncol(X) # 10 变量
n <- nrow(X) # 2500 样本量
p0 <- 5 # prior guess for the number of relevant variables 
sigma <- 1 / sqrt(mean(y)*(1-mean(y))) # pseudo sigma
tau0 <- p0 / (D - p0) * sigma / sqrt(n)
# hs() 函数指定层次收缩先验 Hierarchical shrinkage
# 拟合模型
fit <- stan_glm(
  y ~ X, family = binomial(), data = data.frame(I(X), y), 
  # horseshoe 先验
  prior = hs(df = 1, global_df = 1, global_scale = tau0)
)
# 输出结果
summary(fit, digits = 4)
```

模型输出如下：

``` markdown
Model Info:
 function:     stan_glm
 family:       binomial [logit]
 formula:      y ~ X
 algorithm:    sampling
 sample:       4000 (posterior sample size)
 priors:       see help('prior_summary')
 observations: 2500
 predictors:   11

Estimates:
              mean    sd      10%     50%     90%  
(Intercept)  1.0016  0.0732  0.9080  1.0007  1.0947
X1           3.0921  0.1343  2.9219  3.0888  3.2660
X2          -1.9907  0.1002 -2.1200 -1.9903 -1.8631
X3           0.0205  0.0429 -0.0180  0.0084  0.0804
X4          -0.0069  0.0364 -0.0534 -0.0018  0.0297
X5           0.0045  0.0367 -0.0336  0.0008  0.0474
X6           0.0062  0.0364 -0.0323  0.0014  0.0519
X7           0.0469  0.0559 -0.0068  0.0327  0.1258
X8          -0.0082  0.0376 -0.0545 -0.0021  0.0296
X9          -0.0342  0.0492 -0.1042 -0.0196  0.0109
X10          0.0310  0.0472 -0.0125  0.0180  0.0971

Fit Diagnostics:
           mean   sd     10%    50%    90% 
mean_PPD 0.5915 0.0087 0.5804 0.5916 0.6028

The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).

MCMC diagnostics
              mcse   Rhat   n_eff
(Intercept)   0.0010 0.9994 5422 
X1            0.0021 0.9996 3994 
X2            0.0016 1.0000 3817 
X3            0.0006 0.9997 4531 
X4            0.0005 0.9998 4652 
X5            0.0005 0.9993 5052 
X6            0.0005 0.9994 4795 
X7            0.0010 1.0002 3045 
X8            0.0006 1.0000 4397 
X9            0.0009 1.0002 3034 
X10           0.0008 1.0003 3292 
mean_PPD      0.0001 0.9994 4742 
log-posterior 0.1367 1.0012 1206 

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
```

**rstanarm** 包可以获得与前面一致的结果，甚至收缩效果比手写的 Stan 代码好一点点。
:::

### SpikeSlab 先验 {#sec-prior-spikeslab}

SpikeSlab 先验（Spike Slab）放在非 0 协变量的个数上，是离散的先验。回归系数的先验分布的有限混合，常用有限混合多元正态分布。参考文章 [Discrete Mixture Models](https://betanalpha.github.io/assets/case_studies/modeling_sparsity.html#221_Discrete_Mixture_Models)

::: callout-note
**BoomSpikeSlab** 包是 **Boom** 包的扩展，提供基于 SpikeSlab 先验的贝叶斯变量选择功能。

```{r}
#| eval: false

set.seed(2023)
n <- 2500
k <- 10
X <- matrix(rnorm(n * k), ncol = k)
y <- rbinom(n, size = 1, prob = plogis(1 + 3 * X[, 1] - 2 * X[, 2]))
# 加载 BoomSpikeSlab
library(BoomSpikeSlab)
fit_logit_spike <- logit.spike(y ~ X, niter = 500)
# 模型输出
summary(fit_logit_spike)
```

``` markdown

null log likelihood:            -1690.677 
posterior mean log likelihood:  -766.5283 
posterior max log likelihood:   -754.8686 
mean deviance R-sq:             0.5466147 

predicted vs observed success rates, by decile:
                  predicted    observed
(0.00596,0.0279] 0.01388670 0.008032129
(0.0279,0.108]   0.06371528 0.060000000
(0.108,0.273]    0.17839881 0.176000000
(0.273,0.496]    0.39146661 0.404000000
(0.496,0.734]    0.61807048 0.608000000
(0.734,0.865]    0.80694458 0.764000000
(0.865,0.942]    0.90690322 0.928000000
(0.942,0.979]    0.96436544 0.976000000
(0.979,0.992]    0.98636201 0.992000000
(0.992,0.996]    0.99441504 1.000000000

summary of coefficients:
             mean    sd mean.inc sd.inc inc.prob
(Intercept)  1.02 0.105     1.02  0.105     1.00
X2          -2.00 0.232    -2.02  0.118     0.99
X1           3.10 0.354     3.13  0.170     0.99
X10          0.00 0.000     0.00  0.000     0.00
X9           0.00 0.000     0.00  0.000     0.00
X8           0.00 0.000     0.00  0.000     0.00
X7           0.00 0.000     0.00  0.000     0.00
X6           0.00 0.000     0.00  0.000     0.00
X5           0.00 0.000     0.00  0.000     0.00
X4           0.00 0.000     0.00  0.000     0.00
X3           0.00 0.000     0.00  0.000     0.00
```
:::

## 推理算法 {#sec-choose-inference}

开篇提及 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)

### 惩罚极大似然算法 {#sec-optimization-algorithms}

L-BFGS 算法拟合模型，速度非常快。

```{r}
# L-BFGS 算法拟合模型
fit_optim_logit <- mod_logit_lasso$optimize(
  data = mdata, # 观测数据
  init = 0,     # 所有参数初值设为 0
  refresh = 0,  # 不显示迭代进程
  algorithm = "lbfgs", # 优化器
  threads = 1,    # 单线程
  seed = 20232023 # 随机数种子
)

fit_optim_logit$summary(c("alpha", "beta", "lambda", "lp__"))
```

### 变分近似推断算法 {#sec-variational-approximation-algorithms}

ADVI 算法拟合模型，可选的优化器有 `meanfield` 和 `fullrank` ，相比于 L-BFGS 稍慢

```{r}
# ADVI 算法拟合模型
fit_advi_logit <- mod_logit_lasso$variational(
  data = mdata, # 观测数据
  init = 0,     # 所有参数初值设为 0
  refresh = 0,  # 不显示迭代进程
  algorithm = "meanfield", # 优化器
  threads = 1,    # 单线程
  seed = 20232023 # 随机数种子
)

fit_advi_logit$summary(c("alpha", "beta", "lambda", "lp__"))
```

### 拉普拉斯近似算法 {#sec-laplace-approximation-algorithms}

Stan 内置的 Laplace 近似算法是对后验分布的 Laplace 正态近似，再从近似的后验分布中采样获得样本，最后，对样本进行统计分析获得参数的后验估计。详见 Stan 语言参考手册的[Laplace Approximation 一章](https://mc-stan.org/docs/reference-manual/laplace-approximation.html)。

```{r}
# Laplace 算法
fit_laplace_logit <- mod_logit_lasso$laplace(
  data = mdata, # 观测数据
  init = 0,     # 所有参数初值设为 0
  refresh = 0,  # 不显示迭代进程
  threads = 1,    # 单线程
  seed = 20232023 # 随机数种子
)

fit_laplace_logit$summary(c("alpha", "beta", "lambda", "lp__"))
```

### 探路者变分算法 {#sec-pathfinder-algorithms}

探路者算法 Pathfinder 属于变分法，针对可微的对数目标密度函数，沿着逆牛顿优化算法的迭代路径，获得目标密度函数的正态近似。正态近似中的局部协方差的估计采用 LBFGS 计算的负逆 Hessian 矩阵。探路者算法的优势是可以极大地减少对数密度函数和梯度的计算次数，缓解迭代陷入局部最优点和鞍点（何为鞍点，一个可视化示例详见 @sec-bayesian-gaussian-processes ）。

```{r}
# Pathfinder 算法
fit_pathfinder_logit <- mod_logit_lasso$pathfinder(
  data = mdata, # 观测数据
  init = 0,     # 所有参数初值设为 0
  refresh = 0,  # 不显示迭代进程
  num_threads = 1,    # 单线程
  seed = 20232023 # 随机数种子
)

fit_pathfinder_logit$summary(c("alpha", "beta", "lambda", "lp__"))
```

## 习题

1.  在 @sec-choose-prior 的基础上，比较 Stan 实现的贝叶斯 Lasso 和 R 包 **glmnet** 的结果，发现 **glmnet** 包是很有竞争力的。在选择 Lasso 先验的情况下，收缩效果比 Stan 还好，运行速度也很快。Stan 的优势在于不限于先验分布的选取，当选择 Horseshoe 先验时，Stan 的收缩又比 **glmnet** 包更好。Stan 的优势还在于不限于 **glmnet** 包支持的常见分布族，如高斯、二项、泊松、多项、Cox 等。本质上，这两点都是 Stan 作为一门概率编程语言的优势，只要知道概率分布的数学表达式，总是可以用 Stan 编码出来的。

    ```{r}
    #| message: false

    library(glmnet)
    # 10 折交叉验证 Lasso 回归
    fit_lasso <- cv.glmnet(x = X, y = y, family = "binomial", alpha = 1, nfolds = 10)
    # 回归系数
    coef(fit_lasso, s = fit_lasso$lambda.min)
    ```

2.  基于德国信用卡评分数据，建立逻辑回归模型，分析 20 个协变量对响应变量的贡献，采用合适的先验分布选择适当的变量数。

    ```{r}
    german_credit_data <- readRDS(file = "data/german_credit_data.rds")
    str(german_credit_data)
    ```

3.  下图是美国黄石公园老忠实间歇泉喷发时间和等待时间的分布规律，请建立合适的正态混合模型，用 Stan 拟合模型，并对结果做出说明。（提示：参考 Stan 用户手册的[有限混合章节](https://mc-stan.org/docs/stan-users-guide/mixture-modeling.html)）

    ```{r}
    #| label: fig-faithful-mixture
    #| fig-cap: 黄石公园老忠实间歇泉
    #| fig-width: 7
    #| fig-height: 3
    #| fig-showtext: true
    #| code-fold: true
    #| echo: !expr knitr::is_html_output()

    p1 <- ggplot(data = faithful, aes(x = eruptions)) +
      geom_histogram(aes(y = after_stat(density)), 
                     bins = 30, fill = "white", color = "gray") +
      geom_density() +
      theme_classic() +
      labs(x = "喷发时间", y = "概率密度值")

    p2 <- ggplot(data = faithful, aes(x = waiting)) +
      geom_histogram(aes(y = after_stat(density)), 
                     bins = 30, fill = "white", color = "gray") +
      geom_density() +
      theme_classic() +
      labs(x = "等待时间", y = "概率密度值")

    library(patchwork)
    p1 | p2
    ```

    ```{verbatim, file="code/faithful_finite_mixtures.stan", lang="stan"}
    ```

    ```{r}
    #| code-fold: true
    #| echo: !expr knitr::is_html_output()
    #| eval: false
    #| message: false

    library(cmdstanr)

    faithful_d <- list(
      K = 2, # 几个正态分布混合
      N = 272, # 样本量
      # y = faithful$waiting,
      y = faithful$eruptions
    )

    mod_faithful_normal <- cmdstan_model(
      stan_file = "code/faithful_finite_mixtures.stan",
      compile = TRUE, cpp_options = list(stan_threads = TRUE)
    )

    fit_faithful_normal <- mod_faithful_normal$sample(
      data = faithful_d,
      chains = 2,
      parallel_chains = 2,
      iter_warmup = 1000, 
      iter_sampling = 1000, 
      threads_per_chain = 2, 
      seed = 20232023,
      show_messages = FALSE,
      refresh = 0
    )

    # 输出结果
    fit_faithful_normal$summary(c("theta", "mu", "sigma", "lp__"))
    # theta[1] = 0.350 混合比例
    # theta[2] = 0.650
    # mu[1] = 2.02 均值
    # mu[2] = 4.27
    # sigma[1] = 0.243 标准差
    # sigma[2] = 0.437
    ```

4.  对美国黄石公园的老忠实泉喷发规律，建立二项分布和二维正态分布的混合模型，请用 Stan 编码估计模型中的参数。

    $$
    f(\bm{x};p,\bm{\mu_1},\Sigma_1,\bm{\mu_2},\Sigma_2) = p\mathcal{N}(\bm{x};\bm{\mu_1},\Sigma_1) + (1-p) \mathcal{N}(\bm{x};\bm{\mu_2},\Sigma_2) 
    $$

    其中，参数 $p$ 是一个介于 0 到 1 之间的常数，参数 $\bm{\mu_1} = (\mu_{11},\mu_{12})^\top,\bm{\mu_2}=(\mu_{21},\mu_{22})^\top$ 是二维的列向量，参数 $\Sigma_1 = (\sigma_{ij}),\Sigma_2 = (\delta_{ij}),i=1,2,j=1,2$ 是二阶的协方差矩阵。（提示：因有限混合模型存在可识别性问题，简单起见，考虑各个多元正态分布的协方差矩阵相同的情况。）

    ```{verbatim, file="code/faithful_2d_finite_mixtures.stan", lang="stan"}
    ```

    ```{r}
    #| code-fold: true
    #| echo: !expr knitr::is_html_output()
    #| eval: false
    #| message: false

    data("faithful")
    library(cmdstanr)
    # 准备数据
    faithful_2d <- list(
      K = 2,    # 2 个分布混合
      N = 272,  # 样本量 nrow(faithful)
      D = 2,    # 二维正态分布
      y = faithful # 数据集
    )
    # 编译模型
    mod_faithful_normal_2d <- cmdstan_model(
      stan_file = "code/faithful_2d_finite_mixtures.stan",
      compile = TRUE, cpp_options = list(stan_threads = TRUE)
    )
    # 采样
    fit_faithful_normal_2d <- mod_faithful_normal_2d$sample(
      data = faithful_2d,
      chains = 2,
      parallel_chains = 2,
      iter_warmup = 1000, 
      iter_sampling = 1000, 
      threads_per_chain = 2, 
      seed = 20232023,
      show_messages = FALSE,
      refresh = 0
    )
    # 输出结果
    fit_faithful_normal_2d$summary(c("theta", "mu", "lp__"))
    ```