贝叶斯 RStan 包入门教程-阿里云开发者社区

小编有话说

这是小编在学习 stan 时发现的一篇相关文章，很适合初学者。小编花了写时间对该文章进行了翻译，分享给大家。

本文将介绍 R 与 Stan 的接口（RStan）并通过 Gelman 等人著作^[1] 中的例子来说明 RStan 的特性。

介绍

Stan 是一个用于贝叶斯建模和推断的 C++ 库，主要使用 No-U-Turn 采样器^[2] (NUTS) 在给定指定模型和数据的情况下获得后验样本。或者，Stan 可以利用 LBFGS 优化算法来最大化目标函数（例如，对数似然函数）。rstan 是一个提供 R 接口到 Stan的 R 包。它允许人们方便地拟合 R 中的 Stan 模型并访问输出，包括后验推断和其他中间量，例如对数后验密度及其梯度的评估。

本文将简要介绍 rstan 包中的主要功能。如果读者想进一步了解该包功能，可见 Stan 的官方网站^[3]。该网站提供有关如何操作 Stan 及与其他接口相关的最新信息。关于 RStan 的教程也可见 RStan 入门^[4]。

预备知识

Stan 是一种建模语言，它与贝叶斯图形建模包 BUGS^[5]的语言相似但不完全相同。解析器将模型（ Stan 语言表示）翻译成 C++ 代码，然后将其编译为可执行的程序并作为 R 中的动态共享对象 (Dynamic Shared Object，DSO) 进行加载，然后可以用于用户调用。

注意：此过程需要 C++ 编译器，例如 g++ 或 clang++。安装用于 RStan 的 C++ 编译器的使用教程，请参阅RStan-Getting-Started^[6]。

rstan 包还依赖于其他几个 R 包：

StanHeaders（Stan C++ headers）
BH（Boost C++ headers）
RcppEigen（Eigen C++ headers）
Rcpp（在 R 中使用 C++）
inline（编译 C++ 以与 R 一起使用）

读者在安装 rstan 包时，默认情况下这些依赖项会自动安装。

工作流程

以下是通过 RStan 使用 Stan 进行贝叶斯分析的常用工作流程。

使用 Stan 语言来编写对数后验密度来表示统计模型。建议使用带有 .stan 扩展名的单独文件。
使用 stanc（） 函数将 Stan 程序转换为 C++ 代码。
编译 C++ 代码来创建可由 R 加载的 DSO（也称为：动态链接库 (dynamic link library，DLL)）。
运行 DSO 并从后验分布中采样。
诊断 MCMC 链的收敛性。
基于后验样本进行统计推断。

注意：通过调用 stan() 函数即可隐式执行步骤 2、3 和 4 。

例子

这里以《贝叶斯数据分析》第 5.5 节中的分层模型作为例子。利用分层模型来刻画辅导项目对大学入学考试的影响。下表显示了八所高中得到的实验结果，并为每所高中估计了标准误差。方便起见，该数据简称为“八所学校”数据。

学校名称	估计值()	标准差 ()
A	28	15
B	8	10
C	-3	16
D	7	11
E	-1	9
F	1	11
G	18	10
H	12	18

假设统计模型指定为

其中，假定每个是已知的。的表示不同高中的不同影响。

Stan 程序

RStan 允许将 Stan 程序编码为文本文件（通常带有后缀.stan）或 R 字符向量（长度为 1）。这里，我们推荐使用第一种方式，将所提模型代码放入到文件名为 schools.stan 的文件中：

data {
  int<lower=0> J;          // 学校数量
  real y[J];               // 估计实验效果
  real<lower=0> sigma[J];  // 效应估计的标准差 
}
parameters {
  real mu; 
  real<lower=0> tau;
  vector[J] eta;
}
transformed parameters {
  vector[J] theta;
  theta = mu + tau * eta;
}
model {
  target += normal_lpdf(eta | 0, 1);
  target += normal_lpdf(y | theta, sigma);
}

现在对上面的 Stan 程序进行说明。

第一部分（data），指定了所需数据：学校数量（J），估计向量（）和估计值的标准误差向量（）。数据被声明为整数或实数，并且如果指定了维度，则可以是向量（或更一般地，数组）。
第二部分（parameters），声明了寻求后验分布的参数。这些是学校效应的平均值和标准差，加上标准化的学校水平效应。

注意：在这个模型中，我们让非标准化的学校级效应成为一个转换参数，通过将标准化效应缩放并将它们移动而不是直接声明作为参数。通过这种方式对模型进行参数化，采样器的运行效率更高，因为所产生的多变量几何学更适合于哈密顿蒙特卡罗^[7]。

第三部分（model）看起来类似于标准统计符号。（注意：Stan中的分布的第二个参数是标准差，而不是统计符号中的方差）。我们用矢量符号来写模型，这使得 Stan 可以利用更有效的算法微分。- 第三部分（model）看起来类似于标准统计符号。（注意：Stan中的分布的第二个参数是标准差，而不是统计符号中的方差）。我们用矢量符号来写模型，这使得 Stan 可以利用更有效的算法微分。我们也可以用一个循环来代替normal_lpdf(y | theta,sigma)来写这个模型，但效率较低。

for (j in 1:J) 
  target += normal_lpdf(y[j] | theta[j],sigma[j]);

Stan 有许多对统计建模最有用的 R 函数的版本，包括概率分布、矩阵运算和各种特殊函数。然而，Stan 函数的名称可能与 R 函数的名称不同，更微妙的是，Stan 中的概率分布的参数化可能与 R 中相同分布的参数化不同。为了缓解这个问题，可以向 lookup 函数传递一个 R 函数或命名R函数的字符串，RStan 将尝试查找相应的 Stan 函数，显示其参数，并给出The Stan Development Team（2016）中讨论该函数的页码。

lookup("dnorm")

lookup(dwilcox)   # no corresponding Stan function
# [1] "no matching Stan functions"

如果该 lookup 函数找不到与 Stan 函数对应的 R 函数，它将把它的参数视为正则表达式并尝试查找与 Stan 函数名称的匹配项。

准备数据

stan() 函数接受的数据是一个命名的列表、一个对象名称的字符向量或一个environment。另外，读者也可以省略数据参数，R 将搜索与 Stan 程序的数据块中声明的对象具有相同的名称。下面是八校例子的数据。

schools_data <- list(
  J = 8,
  y = c(28,  8, -3,  7, -1,  1, 18, 12),
  sigma = c(15, 10, 16, 11,  9, 11, 10, 18)
)

实际上，也可以从文件中读取数据，而不是直接在 R 脚本中输入数字。

来自后验分布的样本

接下来，调用 stan() 函数来抽取后验样本：

library(rstan)
fit1 <- stan(
  file = "schools.stan",  # Stan 程序
  data = schools_data,    # 数据
  chains = 4,             # 马尔可夫链数量
  warmup = 1000,          # 每条链的预热迭代次数
  iter = 2000,            # 每条链的总迭代次数
  cores = 1,              # 电脑核数(每条链可以使用一个)
  refresh = 0             # 不显示进展
  )

执行 stan() 函数即可完成以下三个步骤：

将 Stan 代码中的模型转换为 C++ 代码
将 C++ 代码编译为动态共享对象 (DSO) 并加载 DSO
给定一些用户指定的数据和其他设置的样本

对stan的一次调用就可以完成这三个步骤，但也可以逐一执行（见 stanc、stan_model 和 sampling 的帮助页面）。这种方式对调试很有用此外。Stan 保存了 DSO，以便再次拟合相同模型时（可能使用新数据和设置）可以避免重新编译。如果在模型编译之后但在采样之前发生错误（例如，数据和初始值等输入的问题），读者仍然可以重用编译后的模型。

stan() 函数返回一个 stanfit 对象，它是一个 "stanfit" 类的 S4 对象。对于那些不熟悉 R 语言中类和 S4 类概念的人，可以参考Chambers (2008)^[8]。一个 S4 类由一些属性（数据）组成，用来对一个对象进行建模，还有一些方法用来对该对象的行为进行建模。从用户的角度来看，一旦一个 stanfit 对象被创建，我们主要关注的是定义了哪些方法。

如果没有发生错误，返回的 stanfit 对象包括从模型参数的后验分布中抽取的样本和模型中定义的其他数量。如果出现错误（例如， Stan 程序中的语法错误），stan 将退出或返回一个不包含后验抽样的 stanfit 对象。

"stanfit" 类定义了许多方法，如 print 和 plot，用于处理 MCMC 的样本。例如，下面显示了使用 print 方法对八校模型的参数进行的总结。

print(fit1, pars=c("theta", "mu", "tau", "lp__"), probs=c(.1,.5,.9))

此输出的最后一行 lp__ 是 Stan 计算的（未归一化的）后验密度的对数。该对数密度可以多种方式用于模型评估和比较（参见Vehtari 和 Ojanen (2012)^[9]）。

stan函数的参数

采样的主要参数包括：数据、初始值和采样器的选项，例如 chains、iter 和 warmup。特别是，warmup 指定 NUTS 采样器在采样开始之前用于适应阶段的迭代次数。预热之后，采样器关闭适应，继续进行，直到完成总共的迭代次数（包括预热）。理论上不能保证预热期间得到的抽样是来自后验分布，所以热身抽样只能用于诊断而不能用于推断。print方法所显示的参数摘要只使用热身后的抽样计算。

可选 init 参数可用于指定马尔可夫链的初始值。有几种指定初始值的方法，详细信息可以在 stan 函数的文档中找到。绝大多数情况下，让 Stan 随机生成自己的初始值就足够了。然而，有时最好至少为 Stan 程序的参数块中声明的对象的一个子集指定初始值。

Stan 使用支持并行性的随机数生成器 (random number generator，RNG)。RNG 的初始化是由参数 seed 和 chain_id 决定的。即使我们从对 stan() 函数的一次调用中运行多个链，我们也只需要指定一个种子，如果不指定的话，它是由 R 随机生成的。

`"stanfit"`类的方法

rstan包中包含的另一个小册子更详细地讨论了 stanfit 对象，并给出了访问对象中包含的最重要内容的例子（例如，后验图、诊断摘要）。另外，可用方法的完整列表可以在 "stanfit" 类的文档中找到：help("stanfit", "rstan")。这里我们只举几个例子。

stanfit 对象的绘图方法提供了各种输出的图形概览。默认的绘图显示了所有参数的后验不确定性区间（默认为 80%（内部）和 95%（外部））和后验中值，以及 lp__（后验密度函数的对数加上一个常数）。

stanfit 对象的 plot 方法提供了输出的各种图形概览。默认展示后验不确定性区间（默认为 80%（内部）和 95%（外部））和所有参数的后验中值以及 lp__ 的结果。

plot(fit1)

可选 plotfun 参数可用于在各种可用图中进行选择。详细介绍见看 help("plot,stanfit-method")。

traceplot 方法用于绘制后验图的时间序列。如果我们通过设置包含热身绘制，inc_warmup=TRUE 热身区域的背景颜色与热身后阶段不同：

traceplot(fit1, pars = c("mu", "tau"), inc_warmup = TRUE, nrow = 2)

为了评估马尔可夫链的收敛性，除了目测轨迹图之外，我们还可以计算 Split 统计量。Split 是统计的更新版本。由 Gelman 和 Rubin (1992) 提出。有关详细信息，请参阅 Stan 手册。每个参数的包含在 summary 和 print 中。

print(fit1, pars = c("mu", "tau"))

同样，更多的细节请参见关于 stanfit 对象的附加小册子。

抽样困难问题

模型输出的可视化最好方法是通过 ShinyStan 接口，可以通过 shinystan^[10] 包访问。ShinyStan 既有利于参数分布的可视化，也有利于诊断采样器的问题。shinystan 包的文档提供了关于使用 stanfit 对象接口的说明。

除了使用 ShinyStan 之外，还可以使用 rstan 包中的函数来诊断一些采样问题。get_sampler_params()函数返回与采样器性能相关的参数信息：

# 所有链一起运行
sampler_params <- get_sampler_params(fit1, inc_warmup = TRUE)
summary(do.call(rbind, sampler_params), digits = 2)

# 每个链分别运行
lapply(sampler_params, summary, digits = 2)

在这里，我们看到有少量的发散转移，标识为 divergent__ = 1。理想情况下，预热阶段后不应有链发散的情况。解决该问题的最佳方法是增加目标接受概率，默认情况下为 0.8。在这种情况下，所有链的均值 accept_stat__ 接近 0.8，但由于中位数接近0.95，所以分布非常偏斜。我们可以回去再次调用 stan()，并指定可选的参数 control=list(adapt_delta=0.9) 来尝试消除发散的过渡。然而，有时当目标接受率很高时，步长非常小，采样器在每次迭代中都会遇到它的跃迁步数的限制。在这种情况下，这不是问题，因为每个链 treedepth__ 最多为 7，而而默认的是 10。但是，如果任何 treedepth__ 是 11，那么通过向 stan 传递 control=list(max_treedepth=12) 来增加限制是明智的选择。关于 get_sampler_params 返回的对象的结构，请参见 stanfit 对象的小节。

我们还可以使用对（大部分）相同信息进行图形表示pairs。“成对”图可用于了解尾部或模式附近是否出现任何采样困难：

pairs(fit1, pars = c("mu", "tau", "lp__"), las = 1)

在上图中，每个选定参数的边际分布都作为沿对角线的直方图包含在内。默认情况下，低于中位数 accept_stat__（MCMC 建议接受率）的抽样被绘制在对角线下方，而高于中位数 accept_stat__ 的抽样被绘制在对角线上方（可以使用 condition 参数更改）。每个非对角线的正方形代表了行变量和列变量交汇处的抽样的双变量分布。理想情况下，同一两个变量的下方对角线交点和上方对角线交点的分布应该是彼此的镜像。任何黄色的点都表示达到最大 treedepth__ 的过渡，红色的点表示有分歧的过渡。

贝叶斯 RStan 包入门教程

小编有话说

目录

介绍

预备知识

工作流程

例子

Stan 程序

准备数据

来自后验分布的样本

`"stanfit"`类的方法

抽样困难问题

相关资料

参考资料

热门文章

最新文章

相关课程

相关电子书

热门

活动广场

任务中心

开发者评测

高校计划

乘风者计划

训练营

阿里云MVP

话题

直播

下载

镜像站

技术资料

插件

贝叶斯 RStan 包入门教程

小编有话说

目录

介绍

预备知识

工作流程

例子

Stan 程序

准备数据

来自后验分布的样本

"stanfit"类的方法

抽样困难问题

相关资料

参考资料

热门文章

最新文章

相关课程

相关电子书

`"stanfit"`类的方法