Stan 概率编程语言

【摘要】Stan 是一种灵活的概率建模语言，可以使用贝叶斯技术直接估计多种类型的概率模型。
【原文】https://nbviewer.org/github/QuantEcon/QuantEcon.notebooks/blob/master/IntroToStan_basics_workflow.ipynb
【时间】 2016
【作者】Jim Savage, Lendable Inc.

1 概述

1.1 简介

Stan 是一种灵活的建模语言，可以直接使用贝叶斯技术来估计非常广泛的概率模型。有几个原因，使人们可能想花些时间来学习 Stan ：

Stan 使用 Hamilton Monte Carlo、Variational Inference 和 Penalised Maximum Likelihood 实现了对大数据集概率模型的有效估计。
Stan 可以从用户可能用于数据准备的许多环境中调用，包括 R、Python、Stata、Julia 和 Matlab。
Stan 允许用户通过将 Stan 函数（已编译的 C++函数）输出到 R 中来加速他们的 R 代码。这对于希望获得与 Julia 类似的性能，但又被 R 的大型库生态系统所束缚的用户特别有用。
Stan 也许是最容易学习和使用的贝叶斯库，它有直接的语法和配套的库，便于模型检查和比较。
Stan 社区所鼓励的工作流程，迫使用户必须比其他方式更仔细地思考模型。学习 Stan 会发挥很多统计学知识的优势，而其中很多知识可能是你不曾重视的；即使你不使用 Stan，它也应当能够帮助你更仔细地思考建模问题。
在 Stan 中，使用非共轭先验分布没有任何难度，支持与 MCMC 方法相当的丰富模型。

在本教程结束时，您应该对以下内容感到满意：

（1）了解贝叶斯工作流程；

（2）知道如何写出 Stan 模型；

（3）知道如何在 R 中估计和检查模型；

（4）知道从哪里获得帮助。

下面给出的示例将使用 Stan 的 R 接口（ rstan ）。但 Stan 也可以从 Python（使用 PyStan）、Julia（使用 Stan.jl）和其他环境中调用，但有一些非常有用的配套库可能在 R 中开发得更好。

1.2 关于符号的简单说明

下面的教程使用如下概率符号：

在概率建模中，每个随机变量都与概率分布函数（对于连续随机变量）或概率质量函数（对于离散随机变量）相关联，一般符号记为 $p()$。
给定点 $y_{i}$ 的概率分布为 $p(y_{i})$ 。
随机变量 $y$ 的某个特定值 $y_{i}$ 落于 $lower$ 和 $upper$ 两点之间的概率为 $\int_{lower}^{upper}p(y)dy$ ，即这两点之间概率分布函数曲线下的面积。因此，概率分布函数下的总面积应当等于 1 。
不依赖于任何其他随机变量的数据 $y$ 的分布，可以记为 $p(y)$ 。由于其分布不以任何其他随机变量的值为条件，我们习惯称之为无条件分布。无条件分布可以有参数，但其参数应被视为固定值。
在贝叶斯分析中，一个模型的参数通常被汇聚在一个向量 $\theta$ 中，参数的无条件分布则被称为先验分布 ，记作 $p(\theta)$。
如果一个随机变量的分布取决于另一个随机变量的值，则我们将概率分布写为条件概率。例如： $p(y \mid \theta)$ 表示给定参数向量 $\theta$ 的值时 $y$ 的概率分布。
我们可以从概率分布函数中随机抽取样本，而抽中某个给定值的概率，应当与该值的分布成正比。
在数据 $y$ 和参数 $\theta$ 固定的情况下，符号 $L(\theta) = p(y \mid \theta)$ 也用于表示似然函数的值，它是所有数据点 $y_{i} \in y$ 分布的乘积。
令 $\theta = (\mu, \sigma)$ ，其中 $\mu$ 为期望值（平均值），$\sigma$ 为尺度参数。如果我们说 $y \mid \theta \sim \text{Normal}(\mu, \sigma)$ ，则会被读作：变量 $y$ 服从期望值为 $\mu$ 、标准差为 $\sigma$ 的正态分布。
在下面的教程中，我们将去掉左边的条件，将表达式缩写为 $y \sim \text{Normal}(\mu, \sigma)$ ，因为很明显 $y$ 的分布取决于 $ \mu$ 和 $\sigma$。
通常将数据元素 $y_{i} \in y$ 建模为服从其自身特定参数定义的分布。例如，在正常情况下， $y_{i} \sim \text{Normal}(\mu_{i}, \sigma_{i})$ 。许多统计模型只是通过简单地为 $\mu_{i}$ 和 $\sigma_{i}$ 提供函数形式来建模；例如，在高斯线性模型中， $\mu_{i} = X_{i}\beta$ ，并且 $\sigma_{i} = \sigma$ ，其中参数 $\theta = (\beta, \sigma)$ 被建模为随机变量。
高斯线性模型的似然通常写为 $p(y \mid \beta, \sigma)$ 。但 $y$ 的分布显然取决于条件期望值 $X \beta$ ，因此有时你会看到其完整形式： $p(y \mid \beta, \sigma, X)$ 。应该提醒的是，条件信息 $X$ 并不是随机变量，因此它通常不会出现在管道符号的右侧，不过为了清楚起见，下文中我们将包含它。

2. 贝叶斯工作流

围绕着促进高质量建模的工作流问题，贝叶斯社区内存在着非常强大的文化。经济学家应该熟悉其中许多步骤，但有一些步骤与贝叶斯建模不同。

（1）写出完整的概率模型，列出所有观测变量和参数（隐变量）的联合概率模型。
（2）用已知的参数值来模拟该模型，生成模拟数据。
（3）利用模拟数据来估计模型，试着恢复参数值。
（4）根据真实的观测数据估计模型。
（5）检查估计是否正常运行。
（6）运行后验预测检查或时间序列交叉验证，评估模型的拟合效果。
（7）利用拟合结果执行推断和预测任务。
（8）模型比较。
（9）迭代！

通常，我们从最简单的模型开始上述工作流程。只有当其能够在整个工作流程中没问题地运行时，才会逐步增加模型的复杂性并进入下一迭代周期。从简单开始逐层增加复杂度可以有效减少每次迭代的错误范围；如果模型出现了问题，它能够提供一个可用的工作模型，让你快速确定可能发生问题的地方。

在实际工作中，为每一个复杂度创建不同的 git 分支是有价值的。这可以最大限度地减少模型之间污染的可能性。一旦更复杂的模型工作正常，就可以将其合并回主分支中。

最初严格遵守此工作流程可能会让你感觉有些麻烦，但一段时间后你会感觉到轻车熟路。它会减少您工作中犯错的数量，并帮助你聚焦思考建模任务的细节。

下面让我们基于 Stan 概率编程系统来逐个介绍每一步。

2.1 以概率形式写出模型

2.1.1 写出模型

Stan 是以概率模型为核心执行估计的。概率模型假设所有未知参数 $\theta$ 和结果变量 $y$ 均来自某个（条件）概率分布。

我们选取线性回归模型作为第一个示例。该模型有 $N$ 个样本点，因此，有一个长度为 $N$ 的结果向量 $\mathbf{y}$、一个 $N \times P$ 的协变量矩阵 $\mathbf{X}$ （ $P$ 为协变量的数量）、一个长度为 $P$ 的未知系数向量 $\beta$ ，以及一个长度为 $N$ 的残差向量 $\epsilon$ 。该模型通常写成矩阵形式为：

$$
y = X \beta + \epsilon
$$

这目前还不是一个概率模型。真正的概率模型还需要对不确定性做假设，此处需要对残差 $\epsilon$ 和参数 $\beta$ 的概率分布做出假设。我们首先假设残差服从正态分布分布 $\epsilon \sim \text{Normal}(0, \sigma)$ ，其中 $\sigma$ 是另一个模型参数 —— 残差 $\epsilon$ 的标准差。

一旦做出了该分布假设，我们就可以用概率符号来表达上面的线性模型：

$$
y = X \beta + \epsilon \text{ and } \epsilon \sim \text{Normal}(0, \sigma) \implies y \sim \text{Normal}(X\beta, \sigma)
$$

（这是因为如果将 $X \beta$ 和 $\epsilon$ 相加，就会得到一个期望值为 $X \beta$ 的新分布。）

注意正态分布 $\text{Normal}(\mu, \sigma)$ 是一个双参数分布，它有期望值 $\mu$ 和标准差 $\sigma$。但这并不妨碍我们估计一个包含各参数函数的模型，而且参数可能还有自身的参数（就像上面的 $\epsilon$ 一样），此处我们使用了一个线性函数 $\mu = X\beta$ 作为均值的模型，其中 $\beta$ 向量中含有 $P$ 个回归系数。

模型 $y \sim \text{Normal}(X\beta, \sigma)$ 是数据模型。如果我们确实地知道参数 $\beta$ 和 $\sigma$ 的值，则在给定协变量 $X_{i}$ 时，可以通过简单地从 $\text{Normal}(X_{i}\beta,\sigma)$ 中抽取的随机样本，模拟输出 $y_{i}^{\mathrm{sim}}$ 的可信值。这涉及到了贝叶斯模型的一个重要方面，即贝叶斯模型没有预测值，只有预测分布 。但是，实际情况是我们不确定 $\beta$ 和 $\sigma$，因此将为它们生成概率估计。

在贝叶斯统计中，我们希望能够掌握后验分布，在上例中，参数向量 $\theta = (\beta,\sigma)$，其后验分布为 $p(\theta | y, X)$ ，是在给定数据后模型参数的联合分布。根据贝叶斯定理，该分布由下式给出：

$$
p(\theta | y, X) = \frac{p(y | \theta, X)\times p(\theta)}{p(y)}
$$

因为分布 $p(y)$ 不依赖于 $\theta$ 的值，所以通常忽略该表达式的分母项，并将其写成一个比例常数。也就是：

$$
p(\theta | y, X) \propto p(y | \theta, X)\times p(\theta)
$$

如果假设 $\beta$ 和 $\sigma$ 的先验分布相互独立，那么 $p(\theta) = p(\beta) \times p(\sigma)$ ，上式就变成：

$$
p(\beta, \sigma | y, X) \propto p(y | \beta, \sigma, X)\times p(\beta)\times p(\sigma)
$$

现在观测 $p(y|\beta, \sigma, X)$ 项，即给定参数和协变量时 $y$ 的概率分布。我们已经对该分布做了建模假设，它就是我们的数据模型 $\text{Normal}(X_{i}\beta,\sigma)$ 。

2.1.2 选择先验

为了完成（和估计）我们的概率模型，需要为参数 $\beta$ 和 $\sigma$ 指定先验。

首先，先验是用来做什么的？非正式地：

先验将后验分布估计从似然方向拉向先验方向，而且数据中包含的信息越少，这种影响就越明显。
先验有助于在数学上（非因果上）识别模型。例如，假设回归 $y=α+βx+ϵ$，其中对于特定样本 $x=0$ 时的每个观测值， $β$ 无法被识别 — 它可以取任何值而不影响模型的合理性。添加有关 $β$ 的先验信息有助于识别该模型。在此情况下，我们对 $β$ 的后验估计只是先验。
先验可以帮助我们约束对后验分布的估计。

对于组件不变化的回归系数，我们经常使用单变量的正态先验；如果系数因组而异，我们通常指定多元正态先验。先验分布的期望值和尺度可用于包括在观测数据之前已知的信息，例如来自元研究的参数估计，或理论上强加的值。使用所谓的收缩先验或正则化先验也很常见。这些是将估计值缩小到 0（或组均值）的先验。在许多预测任务中，收缩先验有助于防止过度拟合。

关于先验分布的一个重要说明：
如果先验分布不对参数的特定值施加权重，则后验分布也不能对该值施加权重。我们可以通过选择将 $β$ 的估计值限制为具有经济意义的值的先验分布来使用此属性。例如，如果我们估计一个非常简单的成本函数 $\text{costs} = α + β \text{quantity_sold} + ϵ$ ，我们可能想要排除零或负固定成本的可能性（即，我们认为 $α> 0$ ）。在这种情况下，我们可以给 $α$ 一个仅针对正值定义的先验分布，例如截断正态分布、对数正态分布或伽马分布。
类似地，$σ$ 的先验应该被限制为正数。标准偏差的一个方便的先验是半柯西分布（仅限于正数），它提供了一些先验信息，但允许潜在的大标准偏差。

2.1.3 将它们放在一起

我们现在有了概率模型。它是由两个参数的先验分布（一个是 $β$，一个是 $σ$）和一个给定参数时的数据模型组成的。其中 $μ_p$、$σ_p$、$x_0$ 和 $γ$ 为固定值而非待估计的参数，它们汇总了一些先验知识并由研究人员指定。我们有先验：

$$
\text{for }p\in [1 \dots P] \text{ }\beta_{p} \sim \text{Normal}(\mu_{p}, \sigma_{p})
$$

$$
\sigma \sim \text{Cauchy}{+}(x{0}, \gamma)
$$

而数据模型为：

$$
y \sim \text{Normal}(X\beta, \sigma)
$$

2.1.4 可以让模型更丰富一些

上图显示了如何组合一个简单的概率模型，但可能并非您想要的目标。概率建模的真正力量在于，定义更丰富的模型非常简单。

例如，假设我们要对上面的模型进行两个更改。（1）我们愿意接受 $y$ 可能来自所谓的 “重尾” 分布，例如学生的 t 分布。（2）我们要将 $β$ 的第一个元素限制为非负。则我们可以很容易地定义一个稍微丰富的模型，其中：

先验为：

$$
\beta_{1} \sim\text{Normal}{+}(\mu{1}, \sigma_{1})
$$

$$
\text{for }p\in [2 \dots P] \text{ }\beta_{p} \sim\text{Normal}(\mu_{p}, \sigma_{p})
$$

$$
\sigma \sim\text{Cauchy}{+}(x{0}, \gamma)
$$

数据模型（似然）为：

$$
y \sim\text{Student’s t}(\nu, X\beta, \sigma)
$$

其中 $\nu$ 是学生 $t$ 分布的自由度。我们还需要给这个参数设定一个先验，限制在 $0$ 以下。鉴于学生 $t$ 分布在 $\nu$ 大于 $20$ 时接近正态分布，我们可能希望该分布以 $7$ 为中心，并具有相当宽的延展。

$$
\nu \sim \text{Cauchy}_{+}(7, 5)
$$

2.2 用已知参数模拟模型（Stan 简介）

我们现在已经指定了两个概率模型，接下来要做的是模拟来自第二个模型的一些数据，然后通过估计上面正确和错误指定的模型，来检查是否可以恢复已知的模型参数。模拟和恢复已知参数是模型构建中的一个重要检查程序；它通常有助于捕获模型中的错误，并在建模者脑海中形成清晰的模型。本课程以 Stan 作为主要统计软件工具箱。

Stan 模型由代码块组成，每个代码块都执行特定的任务。下文中粗体部分是必须出现在所有 Stan 程序中的内容（即使它们不包含任何参数）：

（1）函数（functions），这个代码块定义了在下面的块中使用的函数。我们将在这儿写一个从假设模型中抽取样本的随机数生成器。

（2）数据，声明要用于模型的数据。

（3）转换数据，对上面传入的数据进行转换。

（4）参数，定义待估计的未知数，包括对其值的任何。

（5）转换后的参数，通常最好使用上面声明的参数和数据转换。

（6）模型，其中定义了全概率模型。

（7）生成量，从模型中生成一系列输出（后验预测、预报、损失函数值等）

下面是执行任务的 R 和 Stan 脚本。

首先我们需要加载一些库。

我通常使用 dplyr 进行数据处理，使用 ggplot2 绘制参数，使用 rstan 建模，使用 reshape2来操控参数的绘制。

# In R:
# Load necessary libraries and set up multi-core processing for Stan
options(warn=-1, message =-1)
library(dplyr); library(ggplot2); library(rstan); library(reshape2)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

接下来我们定义数据生成过程。

请注意，您可以将其定义为 R 中定义的字符串，就像我在下面所做的那样，或者定义为使用 .stan 后缀保存的文本文件。下面的脚本有注释。需要注意的几点：

Stan 编写 C++ 程序，因此需要静态类型（即，按类型定义数据和变量）。在我们下面声明的函数中，参数包括使用向量的向量、使用整数的整数、一些实数和使用矩阵的矩阵。
Stan 中的所有随机数生成器都是分布名称后跟 _rng。这包括您编写的任何函数。
Stan 建模语言用户指南和参考手册中提供了对 Stan 语言的非常完整的说明。

# In R, or you could save the contents of the string in a file with .stan file type

dgp_string <- "

functions {
    /**
   * Return draws from a linear regression with data matrix X,
   * coefficients beta, and student-t noise with degrees of freedom nu
   * and scale sigma.
   *
   * @param X Data matrix (N x P)
   * @param beta Coefficient vector (P x 1)
   * @param nu Residual distribution degrees of freedom.
   * @param sigma Residual distribution scale.
   * @return Return an N-vector of draws from the model.
   */

    vector dgp_rng(matrix X, vector beta, real nu, real sigma) {
      vector[rows(X)] y; // define the output vector to be as long as the number of rows in X
      
      // Now fill it in
      for (n in 1:rows(X))
        y[n] <- student_t_rng(nu, X[n] * beta, sigma);
      return y;
   }
}
data {
 // If we were estimating a model, we'd define the data inputs here
}
parameters {
  // ... and the parameters we want to estimate would go in here
}
model {
  // This is where the probability model we want to estimate would go
}
"

现在我们已经写出了数据生成模型，让我们生成一些已知的参数和协变量，并模拟模型。第一：为数据和参数生成一些值。

# Generate a matrix of random numbers, and values for beta, nu and sigma

set.seed(42) # Set the random number generator seed so that we get the same parameters
N <- 1000 # Number of observations
P <- 10 # Number of covariates
X <- matrix(rnorm(N*P), N, P) # generate an N*P covariate matrix of random data
nu <- 5 # Set degrees of freedom
sigma <- 5 # And scale parameter
beta <- rnorm(10) # Generate some random coefficients that we'll try to recover
# Make sure the first element of beta is positive as in our chosen DGP
beta[1] <- abs(beta[1])

现在我们有了模型的输入，我们应该编译上面的脚本。 Stan 使用模板将您的脚本转换为 C++ 脚本，然后对其进行编译（这将需要一些时间）。我们将使用stan_model() 编译脚本，然后使用expose_stan_functions() 使我们声明的函数对R 可用。

这是一个有用的功能：特别是当您使用嵌套循环时，Stan 函数可以比 R 函数快几个数量级。

# Compile the script
compiled_function <- stan_model(model_code = dgp_string) # you could use file = "path/to/yourfile.stan" if you have saved it as so
# And make the function available to the user in R
expose_stan_functions(compiled_function)

现在我们的 dgp_rng() 在 R 中可用。让我们将它与我们上面声明的数据一起使用来模拟一些假数据。同样——我们想要这样做的原因是为了确保我们可以从我们的数据中恢复已知参数。

# Draw a vector of random numbers for known Xs and parameters
y_sim <- dgp_rng(nu = nu, X = X, sigma = sigma, beta = beta)

# Plot the data
data_frame(y_sim = y_sim) %>% # Declare a data frame and pipe it into a ggplot
    ggplot(aes( x = y_sim)) + # Where we state the x-axis aesthetic (our simulated values)
        geom_histogram(binwidth = 3) # And tell ggplot what sort of chart to build

2.3 估计已知参数的模型

现在我们有 yy 和 XX，我们想要估计 ββ、σσ 和取决于模型的 νν。我们有两个候选概率模型要估计并检查哪个是更合理的数据模型。为此，我们需要在 Stan 中指定两个模型，然后对其进行估计。

让我们直接进入并定义错误指定的模型——下面有大量评论。这是不正确的，因为它假设平局来自正态分布。它也会稍微低效，因为它认为 β1β1 可以是任何值（即使我们可能有一些理由相信 β1β1 是正的）。

# In R, or in your .stan file (contents from within the quotes only)

incorrect_model <- "
data {
  // In this section, we define the data that must be passed to Stan (from whichever environment you are using)

  int N; // number of observations
  int P; // number of covariates
  matrix[N, P] X; //covariate matrix
  vector[N] y; //outcome vector
}
parameters {
  // Define the parameters that we will estimate, as well as any restrictions on the parameter values (standard deviations can't be negative...)

  vector[P] beta; // the regression coefficients
  real<lower = 0> sigma; // the residual standard deviation (note that it's restricted to be non-negative)
}
model {
  // This is where we write out the probability model, in very similar form to how we would using paper and pen

  // Define the priors
  beta ~ normal(0, 5); // same prior for all betas; we could define a different one for each, or use a multivariate prior
  sigma ~ cauchy(0, 2.5);
  
  // The likelihood
  y ~ normal(X*beta, sigma);
}
generated quantities {
  // For model comparison, we'll want to keep the likelihood contribution of each point
  // We will also generate posterior predictive draws (yhat) for each data point. These will be elaborated on below. 
  
  vector[N] log_lik;
  vector[N] y_sim;
  for(i in 1:N){
    log_lik[i] <- normal_log(y[i], X[i,]*beta, sigma);
    y_sim[i] <- normal_rng(X[i,]*beta, sigma);
  }
}
"

现在我们定义正确指定的模型。它与上面相同，但有一些变化：

我们单独定义受约束的参数，然后将其加入到不受约束的参数中
我们定义一个自由度参数 nu，并给它一个先验。
似然模型假设数据来自学生 t 分布。

# In R, or in your .stan file (contents from within the quotes only)

correct_model <- "
data {
  int N; // number of observations
  int P; // number of covariates
  matrix[N, P] X; //covariate matrix
  vector[N] y; //outcome vector
}
parameters {
  // We need to define two betas--the first is the restricted value, the next are the others. We'll join these in the next block
  real<lower = 0> beta_1;
  vector[P-1] beta_2; // the regression coefficients
  real<lower = 0> sigma; // the residual scale (note that it's restricted to be non-negative)
  real<lower = 0> nu; 
}
transformed parameters {
  vector[P] beta;
  beta <- append_row(rep_vector(beta_1, 1), beta_2);
}
model {
  // Define the priors
  beta ~ normal(0, 5); // same prior for all betas; we could define a different one for each, or use a multivariate prior. The first beta will have a prior of the N+(0, 5)
  sigma ~ cauchy(0, 2.5);
  nu ~ cauchy(7, 5);

  // The likelihood
  y ~ student_t(nu, X*beta, sigma);
}
generated quantities {
  // For model comparison, we'll want to keep the likelihood contribution of each point
  vector[N] log_lik;
  vector[N] y_sim;
  for(i in 1:N){
    log_lik[i] <- student_t_log(y[i], nu, X[i,]*beta, sigma);
    y_sim[i] <- student_t_rng(nu, X[i,]*beta, sigma);
  }
}
"

现在我们已经指定了两个模型，让我们用上面生成的 yy 和 XX 来估计它们。

# In R

# Specify the data list that we will pass to Stan. This gives Stan everything declared in the data{} block. 
data_list_2 <- list(X = X, N = N, y = y_sim, P = P)

# Call Stan. You'll need to give it either model_code (like the ones we defined above), a file (.stan file), 
# or a fitted Stan object (fit)
# You should also pass Stan a data list, number of cores to estimate on (jupyter only has access to one), 
# the number of Markov chains to run (4 by default)
# and number of iterations (2000 by default). 
# We use multiple chains to make sure that the posterior distribution that we converge on 
# is stable, and not affected by starting values. 

# The first time you run the models, they will take some time to compile before sampling. 
# On subsequent runs, it will only re-compile if you change the model code. 

incorrect_fit <- stan(model_code = incorrect_model, data = data_list_2, cores = 1, chains = 2, iter = 2000)
correct_fit <- stan(model_code = correct_model, data = data_list_2, cores = 1, chains = 2, iter = 2000)

我们现在已经将我们的两个竞争模型拟合到数据中。估计了什么？

2.3.1 这些拟合的对象包含什么？

如果您习惯于使用最小二乘法（ OLS ）、最大似然法（ MLE ）或高斯混合模型（ GMM ）来估计模型，那么您可能会期望参数的点估计：回归结果中包含参数的点估计以及一些标准误差。贝叶斯倾向于不使用点估计，而是估计参数的分布。除少数模型外，大部分模型都无法获得后验分布的封闭形式解，因此我们使用蒙特卡罗估计。

蒙特卡罗近似非常简单。假设参数 $θ$ 服从 SomeDistribution() 分布，我们无法得到该分布的解析表达形式，但可以从中生成样本。我们的目的是对该分布进行推断；获得该分布的期望值、中位数、标准差、分位数等特性。蒙特卡罗方法只需简单地从分布中抽取足够多的独立样本，然后从这些样本中计算出感兴趣的统计数据，进而做出推断。

重要提示：这些样本来自我们感兴趣的分布，并且趋向于从分布的高概率区域中抽取。

例如，我们要估计 SomeDistribution() 的平均值，则只需要从 $θ_k \sim SomeDistribution()$ 中独立地抽取 $N$ 个样本，然后估计 $E[\theta]=\frac{1}{N}\sum_{n=1}^{N} \theta_n$ 。该估计的标准差为 $O(1/\sqrt{N})$ ，随着样本数量的增加会逐步减小。

一个拟合后的 Stan 对象包含了每个参数的样本。如果模型拟合正确的话，这些样本应当来自模型的后验分布。这些样本抽取自所有参数的联合后验分布，该分布中蕴含参数之间的相关性，即使在先验分布中可能没有。

在上述两个模型的 生成量 代码块中，我们还创建了另外两个变量类型。

第一个 log_lik 是对数似然，主要用于模型比较。针对每个 $y_i$ 值，我们为参数样本计算该值。因此，如果您有 $N$ 个观测值和 $iter$ 次参数抽样，这将包含 $N \times \text{iter}$ 个对数似然值（如果在许多数据点上做模型估计，需要注意一点）。
第二个 yhat 是后验预测样本。对于每次参数的抽样，我们都从数据模型中为 $y$ 观测抽取一个对应的可信值。因此，贝叶斯框架并不为每个观测给出一个确切的“预测值”，而是给出一个预测分布，该分布考虑了参数回归估计中的残差和不确定性。

2.3.2 一些遗留的问题

现在利用我们拟合后的模型生成推断或预测还为时过早。还存在几个问题：

我们对后验分布的估计收敛了吗？
模型拟合中是否还有其他问题？
哪个模型更好？

2.4 使用 `Shinystan` 对模型拟合做检查

为了解决上面的问题 1 和 2，我们需要分析从模型中抽取的参数样本，以检查一些常见问题：

缺乏混合。 我们可以使用不同的初始值运行多个马尔可夫链。这是确定我们是否收敛于同一个后验分布的好方法。如果几个链不“混合”，那么我们不太可能从一个明确指定的后验中采样。此错误最常见的原因是模型指定不当。
平稳性。 马尔可夫链应该是协方差平稳的（即链的均值和方差与抽取观测样本的时间无关）。非平稳性通常是由于模型指定不当或迭代次数不足造成的。
自相关。 特别是在指定不明确或识别不充分的模型中，给定的参数样本可能高度依赖于先前的参数样本。这导致蒙特卡罗估计不可靠。如果您已经检查了模型，那么瘦化（每 n 次抽样保留一个样本）是解决问题的一种方法 —— 尽管并不一定理想。重参数化模型通常是解决此问题的最佳方法（参见手册的第 21 节，关于优化 Stan 代码）。
Divergent Transitions 在后验分布非常扭曲或不规则的模型中，我们经常会遇到 “Divergent Transitions” 的问题。这意味着模型可能没有很好地采样，并且可能需要重新指定或更改采样程序。解决此问题的最简单方法是使用 $\text{control} = \text{list} (\text{adapt_delta} = 0.99)$ 或其他一些接近 1 的数字。这会减慢模型拟合的速度，但有助于 Divergent Transitions 。

为了检查所有这些潜在问题，我们在 R 中使用了一个名为 Shinystan 的工具。你可以用 install.packages("shinystan", dependencies = T) 安装，然后用 Shinystan::launch_shinystan(correct_fit) 运行它。它将在您的网络浏览器中打开一个交互式会话，您可以在其中探索估计的参数。可在此处获取有关 Shinystan 的更多信息。

2.5 比较模型

让我们从查看模型的输出开始。从每个参数中得出的摘要信息可以用 print 简洁地打印出来：

# In R:

print(incorrect_fit, pars = c("beta", "sigma"))
# Notice that we specify which parameters we want; else we'd get values for `log_lik` and `yhat` also

# Some things to note: 

# - mean is the mean of the draws for each observation
# - se_mean is the Monte Carlo error (standard error of the Monte Carlo estimate from the true mean)
# - sd is the standard deviation of the parameter's draws
# - the quantiles are self-explanatory
# - n_eff is the effective number of independent draws. If there is serial correlation between sequential draws, 
#   the draws cannot be considered independent. In Stan, high serial correlation is typically a problem in 
#   poorly specified models
# - Rhat: this is the Gelman Rubin convergence diagnostic. Values close to 1 indicate that the multiple chains
#   that you estimated have converged to the same distribution and are "mixing" well.

也可以指定更多参数：

1
2
3

# In R

print(correct_fit, pars = c("beta", "sigma", "nu"))

绘制已知值和参数估计值通常很有用。如果一切顺利，您估计的参数分布应该对真实值具有合理的权重：

# In R

# Declare a data frame that contains the known parameter names in one column `variable` and their known values
known_parameters <- data_frame(variable = c(paste0("beta[",1:P,"]"),"sigma", "nu"), real_value = c(beta, sigma, nu))


# Extract params as a (draws * number of chains * number of params) array
extract(correct_fit, permuted = F, pars = c("beta", "sigma", "nu")) %>% 
# Stack the chains on top of one another and drop the chains label
    plyr::adply(2) %>% 
    dplyr::select(-chains) %>% 
# Convert from wide form to long form (stack the columns on one another)
    melt() %>% 
# Perform a left join with the known parameters
    left_join(known_parameters, by = "variable") %>%
# Generate the plot
    ggplot(aes(x = value)) + 
    geom_density(fill = "orange", alpha = 0.5) + # Make it pretty
    facet_wrap(~ variable, scales = "free") +
    geom_vline(aes(xintercept = real_value), colour = "red") +
    ggtitle("Actual parameters and estimates\ncorrectly specified model\n")

extract(incorrect_fit, permuted = F, pars = c("beta", "sigma")) %>% 
# Extract params as a (draws * number of chains * number of params) array
    plyr::adply(2) %>% 
    dplyr::select(-chains) %>% 
# Stack the chains on top of one another and drop the chains label
    melt() %>% 
    left_join(known_parameters, by = "variable") %>% # Join the known parameter table
# Convert from wide form to long form (stack the columns on one another)
# Write out the plot
    ggplot(aes(x = value)) + 
    geom_density(fill = "orange", alpha = 0.5) + # Make it pretty
    facet_wrap(~ variable, scales = "free") + # small sub-plots of each variable
    geom_vline(aes(xintercept = real_value), colour = "red") + # red vertical lines for the known parameters
    ggtitle("Actual parameters and estimates\nincorrectly specified model\n") # A title

目前，似乎两个模型在估计回归系数 $β$ 方面做得一样好，但是错误指定的模型严重高估了 $σ$ 。其原因是：采用 $\nu=5$ 的 Student-t 分布会有重尾，正态分布会尝试通过扩大方差来复制极值。

我们还能如何比较这两个模型？

一种方法是使用留一法 loo 软件包。该包的想法是估计每个模型的留一法 (LOO) 交叉验证误差，允许通过 “LOO 信息标准” 进行模型比较。该包估计 $\sum_{n=1}^N \log p(y_n|y_1,…,y_{n−1},y_{n+1} ,…,y_N)$ 而无需重新估计模型 N 次。这种方法的一大优点是，它使我们能够产生有关“模型倾向于生成最佳预测的程度”的概率性估计。我们这样使用 loo ：

# in R

library(loo) # Load the library

# Extract the log likelihoods of both models. Note that we need to declare log_lik in the generated quantities {} block
llik_incorrect <- extract_log_lik(incorrect_fit, parameter_name = "log_lik")
llik_correct <- extract_log_lik(correct_fit, parameter_name = "log_lik")

# estimate the leave-one-out cross validation error 
loo_incorrect <- loo(llik_incorrect)
loo_correct <- loo(llik_correct)

# Print the LOO statistics
print("Incorrect model")
print(loo_incorrect)
sprintf("\n\nCorrect model")
print(loo_correct)

# Print the comparison between the two models
print(compare(loo_incorrect, loo_correct), digits = 2)

统计量 elpd_loo 是预期的对数逐点预测分布。这是对模型对数似然的粗略类比。 p_loo 为我们提供了有效的参数数量。我们可以将 elpd_loo 乘以 $−2$ 来计算 looic，您可以将其想象为来自贝叶斯框架的 AIC 或 BIC（在偏差尺度上）。有关这些统计数据的更多详细信息，请参阅本文。

elpd_diff 为我们提供了两个（或更多）模型之间对数后验分布的预期差异。 elpd_diff 的值大于 $0$ 表示第二个模型比第一个模型生成了更合理的预测 —— 这正是我们所期望的。

2.5 为模型生成后验预测

现在我们已经确定了一个用于推断和/或预测的模型，我们可以开始进行一些预测。

如上所述，在贝叶斯分析中没有预测值，只有预测分布。

回想一下我们的数据模型：$y_{i} \sim \text{Student-t}(\nu, X_{i}\beta, \sigma)$ 。在这个模型下，假设参数的值是固定的。但在我们刚估计的模型中，存在许多合理的 $ν$、$β$ 和 $σ$ 值。这些可信值均来自后验分布，进而将我们带到了后验预测分布。

后验预测由以下过程构建：

（1）从后验分布中抽取参数集 $θ^{\text{draw}}$

（2）对于每个观测 $i$，从 $p(y|θ^{\text{draw}},X_i)$ 中抽取一个 $y_i^{sim}$ 值

（3）重复为所有参数抽取样本

对于每个数据点，最终得到的可信结果与从后验中抽取的一样多。这些样本既考虑了数据模型中的预期变化，也考虑了后验的不确定性。 这允许我们使用蒙特卡罗方法来计算预测的统计量：

# In R

known_y <- data_frame(variable = paste0("y_sim[",1:N,"]"), real_y = y_sim)


# Extract params as a (draws * number of chains * number of params) array
plot_data <- extract(correct_fit, permuted = F, pars = c("y_sim")) %>%
    plyr::adply(2) %>% 
    dplyr::select(-chains) %>% 
# Stack the chains on top of one another and drop the chains label
    melt() %>% 
    left_join(known_y, by = "variable") %>% # Join the known parameter table
# Convert from wide form to long form (stack the columns on one another)
# Write out the plot
    group_by(variable) %>%
    summarise(median = median(value),
              lower = quantile(value, 0.025),
              upper = quantile(value, 0.975),
              actual = first(real_y)) 

plot_data %>%
    ggplot(aes(x = median)) + 
    geom_ribbon(aes(ymin = lower, ymax = upper), fill = "orange", alpha = 0.5) + 
    geom_line(aes(y = median)) +
    geom_point(aes(y = actual)) +
    ggtitle("Actual outcomes and 95% posterior predictive interval\n") # A title