[翻译] RStan: the R interface to Stan

作者：Stan开发团队 翻译：Rhenium Yuan

2023-11-07

原文链接https://cran.r-project.org/web/packages/rstan/vignettes/rstan.html

这篇文章简单介绍了RStan，Stan的R语言接口。Stan是一个用于贝叶斯推断的C++库，使用No-U-Turn采样（一种Hamiltonian Monte Carlo的变式）或（频率学派的）最优化推断。本文使用Gelman等（2003）的例子演示RStan。

简介

Stan是一个用于贝叶斯建模和推断的C++库，使用No-U-Turn采样（NUTS）（Hoffman & Gelman, 2012）获得后验分布模拟。另外，Stan还可以利用LBFGS优化算法最大化目标函数，如对数似然（log-likelihood）。

R语言程序包rstan提供了Stan的R语言接口。rstan包可以让我们在R中便捷地拟合Stan模型并读取结果，包括后验推断和一些中间量，比如对数后验密度及其梯度。

在本文中我们将对rstan包的功能进行简要介绍。Stan的官网mc-stan.org介绍了更多关于如何使用Stan及其各种接口的详细内容，比如RStan Getting Started（Stan开发团队，2014）。

前置条件

Stan是一种建模语言，该语言和贝叶斯图形建模程序BUGS（Lunn et al., 2000）类似，但不完全相同。解析器将Stan语言表达的模型翻译成C++代码，然后编译成可执行程序，并在R中加载为动态共享对象（Dynamic Shared Object，DSO），用户可以调用该动态共享对象。

这个过程中需要用到C++编译器，比如g++或clang++。关于如何安装RStan可以使用的C++编译器可以看RStan-Getting-Started。

rstan包依赖下列R包：

在安装rstan包时这些包会自动下载。

典型的工作流

以下是通过RStan接口使用Stan进行贝叶斯推断的工作流。

  1. 使用Stan语言写出统计模型的对数后验密度（到一个不依赖于模型中未知参数的归一化常数），以此来表示模型。我们建议使用扩展名为 .stan的单独文件，当然也可以使用R中的字符串。

  2. 使用stanc函数将Stan程序翻译为C++代码。

  3. 将C++代码编译并生成可以被R加载的DSO（也叫dynamic link library, DLL）。

  4. 运行DSO并从后验分布中抽样。

  5. 诊断MCMC链（马尔科夫链蒙特卡洛）的收敛性。

  6. 根据后验样本（后验分布的MCMC抽样得到）进行推断。

一般来说第2,3,4步都在stan函数中进行了。

示例

在接下来的部分中，我们将使用Gelman等（2003）的研究中的一个分层元分析模型（hierarchical meta-analysis model）作为示例进行演示。案例中，分层模型被用于研究辅导计划对大学入学考试的影响。数据如下表所示，总结了在8所高中进行的实验结果和估计的标准误。这些数据和模型作为完全贝叶斯推理的示例具有历史意义（Rubin, 1981）。为了方便期间，我们简称“八校”示例（Eight Schools examples）。

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

我们在这里使用“八校”示例是因为它很简单同时代表了一种不无意义的马尔科夫链模拟问题，因为研究问题的原始参数——八所学校的效应——以及代表这些效应在模型总体中变化的超参数之间存在某种互相依赖的关系。在这个例子中，某些Gibbs采样或Hamiltonian蒙特卡洛采样可能收敛较慢。

我们的统计模型可以指定为

\[\begin{aligned} y_j &\sim \mathsf{Normal}(\theta_j, \sigma_j), \quad j=1,\ldots,8 \\ \theta_j &\sim \mathsf{Normal}(\mu, \tau), \quad j=1,\ldots,8 \\ p(\mu, \tau) &\propto 1, \end{aligned}\]

其中\(\sigma_j\)未知。

编写Stan程序

RStan允许将Stan程序编码为文本文件（通常后缀为.stan）或R字符串向量（长度为1）。我们将以下八校模型的代码写在schools.stan文件中：

data {
  int<lower=0> J;          // number of schools 
  real y[J];               // estimated treatment effects
  real<lower=0> sigma[J];  // s.e. of effect estimates 
}
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\)，估计值向量\((y_1,\dots,y_J)\)，以及估计标准误向量\((\sigma_1,\dots,\sigma_J)\)。数据声明为整型或实型，如果指定了维数，则可以是向量（或更一般的数组）。数据也可以被约束，比如在上面的模型里\(J\)被限制为不小于1而\(\sigma_y\)必须为正数。

parameter模块声明了后验分布的参数。包括效应的均值\(\mu\)，标准差\(\tau\)，以及标准化的学校水平的效应\(\eta\)。在模型中，我们让非标准化的学校水平的效应\(\theta\)作为一个转换后的参数而非直接作为模型参数，其构建方法为：用\(\tau\)对效应进行标准化，并减去均值\(\mu\)。我们这样对模型设置参数是为了让采样的运行效率更高，因为这样产生的多个变量的分布更适合Hamiltonian蒙特卡洛采样（Neal, 2011）。

最后，model模块看起来和一般的统计记号很像。（注意Stan里的正态分布的第二个参数是标准差，而非一般统计学记号里的方差）。我们用向量来编写模型，这样Stan就能利用更高效的算法微分（algorithmic differentiation , AD）。我们也可以用一个循环来代替normal_lpdf(y|θ,sigma)来编写模型，但效率较低。

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

Stan中有很多R中非常好用的统计建模的函数，包括概率分布、矩阵运算还有各种特殊函数。但是这些Stan函数的名字可能和R里不一样，而且需要注意的是，Stan中概率分布函数的参数可能与R中不同。为了解决这个问题，我们可以使用lookup函数查找一个R函数在Stan中对应的函数及其参数，以及在Stan开发团队的文档（2016）中该函数的页码。

lookup('dnorm')
##           StanFunction
## 374 normal_id_glm_lpmf
## 375 normal_id_glm_lpmf
## 376      normal_id_glm
## 379        normal_lpdf
## 380             normal
##                                                        Arguments ReturnType
## 374   (vector y , matrix x, real alpha, vector beta, real sigma)       real
## 375 (vector y , matrix x, vector alpha, vector beta, real sigma)       real
## 376                                                            ~       real
## 379                            (reals y , reals mu, reals sigma)       real
## 380                                                            ~       real
lookup('dwilcox') #Stan中没有Wilcoxon秩和分布的对应函数
## [1] "no matching Stan functions"

如果lookup没能查找到对应的Stan函数，则会将其参数视为正则表达式，并尝试查找与Stan函数名称匹配的参数。

准备数据

stan函数接受的数据包括命名列表、对象名称字符向量或环境（environment）。也可以省略data参数，此时R将搜索与Stan程序data模块中声明的名称相同的对象。以下是八校示例的数据：

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             # 不显示运行过程
  )
## Warning in readLines(file, warn = TRUE):
## 读'E:\贝叶斯统计2023\schools.stan'时最后一行未遂
## Warning: There were 3 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems

stan函数包括下面三步：

只需调用一次stan就可以完成所有三个步骤，但也可以逐个执行（参见stanc、stan_model和sampling的帮助页面），这对调试非常有用。此外，Stan 还会保存DSO，以便再次拟合同一模型（可能使用新数据和设置）时避免重新编译。如果在模型编译后、采样前发生错误（如数据和初始值等输入出现问题），我们仍可重复使用编译后的模型。

stan函数返回一个stanfit对象，即类型为”stanfit“的S4对象。如果您不熟悉R中类和S4类的概念，请参阅Chambers(2008)。S4类由一些属性（数据）和一些方法组成，前者用于对对象进行建模，后者用于对对象的行为进行建模。从用户的角度来看，一旦创建了stanfit对象，我们主要关心的是定义了哪些方法。

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

“stanfit”类定义了许多方法，如print和plot，用于处理MCMC样本。例如，下面使用print显示了八校模型的参数摘要：

print(fit1, pars=c("theta", "mu", "tau", "lp__"), probs=c(.1,.5,.9))
## Inference for Stan model: anon_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##            mean se_mean   sd    10%    50%    90% n_eff Rhat
## theta[1]  11.50    0.16 8.37   2.53  10.40  21.94  2752    1
## theta[2]   7.69    0.10 6.25  -0.08   7.68  15.51  4060    1
## theta[3]   6.29    0.13 7.83  -3.34   6.70  15.14  3696    1
## theta[4]   7.71    0.10 6.61  -0.50   7.70  15.88  4021    1
## theta[5]   5.23    0.11 6.32  -2.92   5.62  12.69  3559    1
## theta[6]   6.34    0.11 6.69  -1.89   6.72  14.18  3715    1
## theta[7]  10.79    0.11 6.59   3.09  10.14  19.52  3619    1
## theta[8]   8.40    0.14 8.02  -0.52   8.18  17.72  3283    1
## mu         8.02    0.12 5.05   1.69   8.02  14.25  1862    1
## tau        6.64    0.15 5.60   0.96   5.24  14.01  1432    1
## lp__     -39.50    0.07 2.65 -43.00 -39.24 -36.37  1385    1
## 
## Samples were drawn using NUTS(diag_e) at Tue Nov  7 23:08:04 2023.
## For each parameter, 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).

输出的最后一行lp__是Stan计算的（未归一化的）后验密度的对数。这个后验密度对数可用于模型评估和比较（参见Vehtari & Ojanen (2012)）。

stan函数的参数

采样的主要参数（在函数stan和sampling中）包括数据、初始值和采样方法的选项，如chains、iter和warmup。其中，warmup指定了NUTS采样在采样开始前的预热阶段所使用的迭代次数。预热后，采样会关闭预热，并继续迭代，直到完成总迭代次数（包括预热）。理论上无法保证在预热过程中获得的抽样来自后验分布，因此预热抽样只能用于诊断，而不能用于推断。print方法显示的参数摘要仅使用预热后的抽样进行计算。

可选的init参数可用于指定马尔科夫链的初始值。有几种方法可以指定初始值，详情请参见stan函数的文档。绝大多数情况下，让Stan随机生成自己的初始值就足够了。不过有时最好至少为Stan程序parameter模块中声明的对象的子集指定初始值。

Stan使用的随机数生成器（random number generator, RNG）支持并行性。RNG的初始化由参数seed和chain_id决定。即使一次调用 stan函数运行多条链，我们也只需指定一个种子，如果没有指定，种子将由R随机生成。

数据预处理和传递

传递给stan的数据将经过预处理。有关预处理的详细信息，请参阅stan函数的文档。这里我们强调几个重要步骤。首先，RStan允许用户传递比data模块中声明的更多对象作为数据（静默省略任何不必要的对象）。一般来说，从R传递给Stan的数据列表中的元素应为数值型，其维度应与模型data模块中的声明相匹配。例如，R中的factor类型不支持作为RStan的数据元素，必须通过as.integer转换为整型。Stan语言区分整型和双精度类型（在Stan语言中分别为 int和real类型）。如果可能，stan函数会将一些R数据（通常是双精度数据）转换为整型。

Stan语言有标量和其他标量集合类型，如向量、矩阵和数组。由于R没有真正的标量，RStan将长度为一的向量视为标量。然而，考虑一个data模块定义为

data {                
  int<lower=1> N;      
  real y[N];
}

其中N可以为1，即个案。所以如果我们指导N始终大于1，我们可以使用R中长度为N的向量作为y的数据输入（比如，y <- rnorm(N)）。如果我们想要防止RStan在N为1时将输入y的数据作为标量，我们需要将y转换为数组，R代码如下：

y <- as.array(y)

Stan不能自动处理数据中的缺失值，所以数据中的任何元素不能含有NA值。RStan数据预处理的一个重要步骤是检查缺失值，如果发现任何缺失值，就会出错。不过，有多种方法可以编写可考虑缺失数据的Stan程序（参见Stan开发团队(2016)）。

用于stanfit类型的方法

rstan包中的另一个简介更详细地讨论了stanfit对象，并举例说明了如何访问对象中包含的最重要内容（如后验绘图、诊断摘要）。此外，“stanfit”类的全部可用方法列表可在help("stanfit", "rstan")的帮助文档中找到。这里我们仅举几个例子。

plot(fit1)
## 'pars' not specified. Showing first 10 parameters by default.
## ci_level: 0.8 (80% intervals)
## outer_level: 0.95 (95% intervals)

可选的plotfun参数可用于从各种可用绘图中进行选择。请参阅help("plot,stanfit-method").

使用traceplot方法绘制后验抽样的轨迹图。如果我们通过设置inc_warmup=TRUE来包含预热抽样，预热区域的背景颜色就会与预热后阶段不同：

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

为了评估马尔科夫链的收敛性，除了检查轨迹图我们还可以计算半分\(\hat{R}\)统计量（split R statistic）。半分\(\hat{R}\)是Gelman和Rubin（1992）提出的\(\hat{R}\)的改进版，它基于将马尔科夫链分成两半。更多细节详见Stan手册。每个参数估计的\(\hat{R}\)在summary和print输出中有一列显示。

print(fit1, pars = c("mu", "tau"))
## Inference for Stan model: anon_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##     mean se_mean   sd  2.5%  25%  50%   75% 97.5% n_eff Rhat
## mu  8.02    0.12 5.05 -2.02 4.76 8.02 11.20 18.08  1862    1
## tau 6.64    0.15 5.60  0.24 2.50 5.24  9.27 20.51  1432    1
## 
## Samples were drawn using NUTS(diag_e) at Tue Nov  7 23:08:04 2023.
## For each parameter, 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).

更多详情，请再次参阅有关stanfi对象的附加简介。

采样中的困难

可视化模型输出的最佳方式是通过ShinyStan界面，该界面可通过shinystan R包访问。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)
##  accept_stat__    stepsize__      treedepth__   n_leapfrog__  divergent__   
##  Min.   :0.00   Min.   : 0.027   Min.   :0.0   Min.   :  1   Min.   :0.000  
##  1st Qu.:0.76   1st Qu.: 0.279   1st Qu.:3.0   1st Qu.:  7   1st Qu.:0.000  
##  Median :0.95   Median : 0.349   Median :3.0   Median : 15   Median :0.000  
##  Mean   :0.82   Mean   : 0.395   Mean   :3.4   Mean   : 12   Mean   :0.011  
##  3rd Qu.:0.99   3rd Qu.: 0.422   3rd Qu.:4.0   3rd Qu.: 15   3rd Qu.:0.000  
##  Max.   :1.00   Max.   :11.097   Max.   :7.0   Max.   :223   Max.   :1.000  
##     energy__ 
##  Min.   :35  
##  1st Qu.:42  
##  Median :44  
##  Mean   :45  
##  3rd Qu.:47  
##  Max.   :63
# 每条链分开
lapply(sampler_params, summary, digits = 2)
## [[1]]
##  accept_stat__    stepsize__     treedepth__   n_leapfrog__  divergent__   
##  Min.   :0.00   Min.   :0.027   Min.   :0.0   Min.   :  1   Min.   :0.000  
##  1st Qu.:0.73   1st Qu.:0.349   1st Qu.:3.0   1st Qu.:  7   1st Qu.:0.000  
##  Median :0.94   Median :0.349   Median :3.0   Median : 15   Median :0.000  
##  Mean   :0.81   Mean   :0.400   Mean   :3.3   Mean   : 12   Mean   :0.012  
##  3rd Qu.:0.99   3rd Qu.:0.398   3rd Qu.:4.0   3rd Qu.: 15   3rd Qu.:0.000  
##  Max.   :1.00   Max.   :6.975   Max.   :7.0   Max.   :223   Max.   :1.000  
##     energy__ 
##  Min.   :35  
##  1st Qu.:42  
##  Median :44  
##  Mean   :45  
##  3rd Qu.:47  
##  Max.   :59  
## 
## [[2]]
##  accept_stat__    stepsize__     treedepth__   n_leapfrog__  divergent__   
##  Min.   :0.00   Min.   :0.033   Min.   :0.0   Min.   :  1   Min.   :0.000  
##  1st Qu.:0.76   1st Qu.:0.324   1st Qu.:3.0   1st Qu.:  7   1st Qu.:0.000  
##  Median :0.95   Median :0.324   Median :3.0   Median : 15   Median :0.000  
##  Mean   :0.82   Mean   :0.378   Mean   :3.4   Mean   : 13   Mean   :0.011  
##  3rd Qu.:0.99   3rd Qu.:0.379   3rd Qu.:4.0   3rd Qu.: 15   3rd Qu.:0.000  
##  Max.   :1.00   Max.   :7.442   Max.   :7.0   Max.   :127   Max.   :1.000  
##     energy__ 
##  Min.   :35  
##  1st Qu.:42  
##  Median :44  
##  Mean   :45  
##  3rd Qu.:47  
##  Max.   :63  
## 
## [[3]]
##  accept_stat__    stepsize__     treedepth__   n_leapfrog__  divergent__   
##  Min.   :0.00   Min.   :0.032   Min.   :0.0   Min.   :  1   Min.   :0.000  
##  1st Qu.:0.66   1st Qu.:0.395   1st Qu.:3.0   1st Qu.:  7   1st Qu.:0.000  
##  Median :0.91   Median :0.422   Median :3.0   Median :  7   Median :0.000  
##  Mean   :0.78   Mean   :0.435   Mean   :3.1   Mean   : 11   Mean   :0.013  
##  3rd Qu.:0.98   3rd Qu.:0.422   3rd Qu.:3.0   3rd Qu.: 15   3rd Qu.:0.000  
##  Max.   :1.00   Max.   :6.859   Max.   :7.0   Max.   :191   Max.   :1.000  
##     energy__ 
##  Min.   :35  
##  1st Qu.:42  
##  Median :44  
##  Mean   :44  
##  3rd Qu.:47  
##  Max.   :61  
## 
## [[4]]
##  accept_stat__    stepsize__     treedepth__   n_leapfrog__  divergent__   
##  Min.   :0.00   Min.   : 0.04   Min.   :0.0   Min.   : 1    Min.   :0.000  
##  1st Qu.:0.89   1st Qu.: 0.25   1st Qu.:3.0   1st Qu.: 7    1st Qu.:0.000  
##  Median :0.98   Median : 0.25   Median :4.0   Median :15    Median :0.000  
##  Mean   :0.88   Mean   : 0.37   Mean   :3.6   Mean   :13    Mean   :0.007  
##  3rd Qu.:0.99   3rd Qu.: 0.44   3rd Qu.:4.0   3rd Qu.:15    3rd Qu.:0.000  
##  Max.   :1.00   Max.   :11.10   Max.   :6.0   Max.   :95    Max.   :1.000  
##     energy__ 
##  Min.   :35  
##  1st Qu.:42  
##  Median :44  
##  Mean   :44  
##  3rd Qu.:47  
##  Max.   :59

在这里我们可以看到有少量的发散转换（divergent transitions），这些分歧转换的标识是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可以了解尾部或图案附近是否存在采样困难的问题：

pairs(fit1, pars = c("mu", "tau", "lp__"), las = 1)
## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

在上图中，每个选定参数的边缘分布都以直方图的形式沿对角线标出。默认情况下，低于中位数accept_stat__（MCMC建议接受率）的抽样绘制在对角线下方，高于中位数accept_stat__的抽样绘制在对角线上方（可以使用条件参数进行更改）。每个对角线外的正方形代表行变量和列变量交叉点抽样的二元分布。理想情况下，同两个变量的对角线以下交集和对角线以上交集的分布应该互为镜像。任何黄点都表示达到treedepth__的过渡，红点表示发散转换。

其他问题

用户定义的Stan函数

Stan还允许用户定义自己的函数，这些函数可在整个Stan程序中使用。这些函数在function模板中定义。function模板是可选的，但如果存在，则必须放在任何其他模块之前。通过这种机制，用户可以实现统计分布或Stan目前不具备的其他功能。不过，即使用户的函数只是封装了对现有Stan函数的调用，如果将完成一个（或两个）任务的几行Stan代码替换为对用户定义函数的调用，model模块中的代码也会变得更加易读。

使用用户自定义函数的另一个原因是，RStan提供了expose_stan_functions函数，用于将此类函数导出到R全局环境中，以便在R中对其进行测试，确保其工作正常。例如

model_code <-
'
functions {
  real standard_normal_rng() {
    return normal_rng(0,1);
  }
}
model {}
'
expose_stan_functions(stanc(model_code = model_code))
standard_normal_rng()
## [1] -0.9529876

对数后验（函数和梯度）

Stan定义了后验分布的概率密度函数的对数，但不包括一个加的未知常数。我们用lp__来表示每次迭代时该对数核的实现情况（lp__在摘要和半分\(\hat{R}\)和有效样本量的计算中被视为未知量。

rstan包的一个特点是，它提供了用于计算给定stanfit对象的lp__及其梯度的函数。这两个函数分别是log_prob和grad_log_prob。这两个函数都使用无约束空间上的参数，即使参数不在整个实数轴上定义。Stan手册（Stan开发团队，2016）详细介绍了Stan用于从整个实轴映射到其中某个子空间（反之亦然）的特定变换。

无约束参数的数量可能少于参数总数。例如，对于K个单纯形参数，实际上只有K-1个无约束参数，因为简单方程的所有元素都必须是非负数且总和为1。get_num_upars方法用于获取非约束参数的数量，而unconstrain_pars和constrain_pars方法则分别用于计算参数的非约束值和约束值。前者将参数列表作为输入，并将其转换为无约束向量，后者则相反。利用这些函数，我们可以实现贝叶斯模型的最大后验估计等其他算法。

Stan中的优化

RStan还提供了Stan优化器的接口，可以通过最大化Stan程序定义的似然函数（可能是受惩罚的）来获得点估计值。我们将用一个非常简单的例子来说明这一功能：从假定取自标准差已知的正态分布的样本中估计平均值。也就是说，我们假设

\[y_n \sim \mathsf{Normal}(\mu,1), \quad n = 1, \ldots, N.\]

通过确定先验\(p(\mu)\propto 1\)，\(\mu\)的最大后验估计就是样本均值。由于\(p(\mu)\propto 1\)是默认先验，我们不需要额外的代码设置\(\mu\)的先验。

我们首先创建一个"stanmodel"类对象，然后使用优化方法，向该方法输入数据和其他参数。

ocode <- "
  data {
    int<lower=1> N;
    real y[N];
  } 
  parameters {
    real mu;
  } 
  model {
    target += normal_lpdf(y | mu, 1);
  } 
"

sm <- stan_model(model_code = ocode)
y2 <- rnorm(20)
mean(y2)
## [1] 0.1563328
optimizing(sm, data = list(y = y2, N = length(y2)), hessian = TRUE)
## $par
##        mu 
## 0.1563328 
## 
## $value
## [1] -28.55551
## 
## $return_code
## [1] 0
## 
## $hessian
##     mu
## mu -20
## 
## $theta_tilde
##             mu
## [1,] 0.1563328

模型编译

正如前面提到的，Stan程序是用Stan语言编写的，翻译成C++代码后编译成动态共享对象（DSO）。然后由R加载DSO并执行以获得后验样本。将C++代码编译成DSO的过程有时需要一段时间。当模型相同时，我们可以重复使用之前运行的DSO。stan函数接受一个可选参数fit，它可以用来传递一个现有的拟合模型对象，从而重复使用编译后的模型。在重复使用以前的拟合模型时，我们仍然可以为stan的其他参数指定不同的值，包括为data参数传递不同的数据。

此外，如果使用save和save.image等函数保存拟合模型，RStan还可以保存 DSO，以便在不同的R会话中使用。要避免保存DSO，请在调用stan函数时指定save_dso=FALSE。

如果用户执行rstan_options(auto_write=TRUE)，那么编译后模型的序列化版本将自动保存到硬盘中与.stan文件位于同一目录下，如果Stan程序以字符串形式表示，则保存到R的临时目录下。虽然根据CRAN的规定，默认情况下不启用该选项，但用户通常应该指定该选项，以消除冗余编译。

当代码以最高优化级别编译时，Stan运行速度会快很多，大多数C++编译器的优化级别为-O3。然而，R的默认值是-O2，这对大多数R包来说都是合适的，但对Stan来说却会略微减慢运行速度。你可以按照CRAN - Customizing-package-compilation中的说明在本地更改默认值。不过，需要注意的是，设置CXXFLAGS = -O3可能会对其他R包产生不利的副作用。

有关Stan程序解析和编译的更多详情，请参阅stanc和stan_model函数的文档。

并行运算多条链

运行马尔可夫链的数量可以通过stan或采样函数的chains参数来指定。默认情况下，链会使用父R进程串行执行（即一次执行一个）。还有一个可选的核心参数，可以设置为链的数量（如果硬件有足够的处理器和内存），这在大多数笔记本电脑上都是合适的。我们通常建议每个R会话先调用一次options(mc.cores=parallel::detectCores())，这样就可以使用所有可用的内核，而无需手动指定内核参数。

对于使用不同并行化方案（可能是远程集群）的用户，rstan软件包提供了一个名为sflist2stanfit的函数，用于将多个stanfit对象（由相同的Stan程序创建，并使用相同的预热和采样迭代次数）合并为一个stanfit对象。为所有链指定相同的种子非常重要，使用不同的链ID（参数chain_id）也同样重要，两者结合可确保Stan中为所有链生成的随机数基本上是独立的。当cores大于1时这些会自动进行。

使用CmdStan

rstan软件包提供了一些函数，用于为CmdStan（Stan的命令行接口）创建数据和读取输出。

首先，当Stan读取数据或初始值时，它只支持R dump数据格式的语法子集。因此，如果我们使用基础R中的dump函数来准备数据，Stan可能无法读取其中的内容。rstan中的stan_rdump函数旨在将R中的数据转储为Stan支持的格式，其语义与dump函数非常相似。

其次，read_stan_csv函数通过读取CmdStan生成的CSV文件创建stanfit对象。生成的stanfit对象与各种诊断和后验分析方法兼容。

参见

Chambers, John M. 2008. Software for Data Analysis : Programming with R. New York: Springer.

Gelman, Andrew, J. B. Carlin, Hal S. Stern, and Donald B. Rubin. 2003. Bayesian Data Analysis. 2nd ed. London: CRC Press.

Gelman, Andrew, and Donald B. Rubin. 1992. “Inference from Iterative Simulation Using Multiple Sequences.” Statistical Science 7 (4): 457–72.

Hoffman, Matthew D., and Andrew Gelman. 2012. “The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo.” Journal of Machine Learning Research.

Lunn, D.J., A. Thomas, N. Best, and D. Spiegelhalter. 2000. “WinBUGS — a Bayesian Modelling Framework: Concepts, Structure, and Extensibility.” Statistics and Computing, 325–37.

Neal, Radford. 2011. “MCMC Using Hamiltonian Dynamics.” In Handbook of Markov Chain Monte Carlo, edited by Steve Brooks, Andrew Gelman, Galin L. Jones, and Xiao-Li Meng, 116–62. Chapman; Hall/CRC.

R Core Team. 2014. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Rubin, Donald B. 1981. “Estimation in Parallel Randomized Experiments.” Journal of Educational and Behavioral Statistics 6 (4): 377–401.

The Stan Development Team. 2014. “RStan Getting Started.” https://mc-stan.org/.

———. 2016. Stan Modeling Language: User’s Guide and Reference Manual. https://mc-stan.org.

Vehtari, A., and J. Ojanen. 2012. “A Survey of Bayesian Predictive Methods for Model Assessment, Selection and Comparison.” Statistics Surveys 6: 142–228.