MCMC的rstan贝叶斯回归模型和标准线性回归模型比较

为了简单起见，我们从回归的标准线性模型开始。然后添加对采样分布或先验的更改。我们将通过 R 和相关的 R 包 rstan 使用编程语言 Stan。

相关视频

示例：线性回归模型

在下文中，我们将设置一些初始数据，并使用标准 lm 函数运行模型比较。

设置

首先，我们需要创建在此处使用的数据。

# 设置可复制种子

set.seed(8675309)

# 运行 lm 以供稍后比较； 但如果需要，请立即检查

modlm = lm(y~., data=data.frame)

此时我们有三个协变量和一个 y，它是正态分布线性函数，标准差等于 2。系数的总体值包括截距分别为 5、0.2、-1.5 和 0.9，尽管添加了噪声，但样本的实际估计值略有不同。现在我们准备好为输入到 Stan 的数据设置一个 R 列表对象，以及对这些数据进行建模的相应 Stan 代码。

我将展示在 R 中通过单个字符串实现的所有 Stan 代码，然后提供每个相应模型块的一些细节。但是，这里的目标不是专注于工具，而是专注于概念。

Stan 的数据列表应包括 Stan 代码中可能使用的任何矩阵、向量或值。例如，与数据一起，可以包括样本大小、组指标（例如混合模型）等。在这里，我们可以只使用样本大小 (N)、模型矩阵中的列数 (K)、目标变量 (y) 和模型矩阵 (X)。

# 为stan输入创建数据列表对象

dat = list

接下来是 Stan 代码。在 R2OpenBugs 或 rjags 中，可以使用代码调用单独的文本文件，并且可以对 rstan 执行相同操作，但出于我们的目的，我们在 R 代码中显示它。首先要注意的是模型代码。接下来，Stan 有必须按顺序调用的编程块。我将在代码中列出所有块来记录它们的顺序并依次讨论每个块。// 或 # 之后或 // 之间的任何内容都是与代码相关的注释。而分布用 ∼∼ 指定，例如，  表示 y 正态分布，均值为 0，标准差为 1。

此外，要安装 rstan，需要通过 CRAN 或 GitHub。它不需要单独安装 Stan 本身，但它确实需要几个步骤并且需要 C++ 编译器。一旦你安装了 rstan，它就会像任何其他 R 包一样被调用。

# 使用语法创建模型对象

stanmodelcode = "

data { // 数据块

int N; // 样本大小

int K; // 模型矩阵的尺寸

/*

转换后的数据 { // 转换后的数据块。目前不使用。

}

*/

参数// 参数块

real sigma; // 误差比例

}

模型 // 模型块

mu = X * beta; // 创建线性预测器

// 先验指标

sigma ~ cauchy(0, 5); // 由于sigma被限定在0

// 似然

y ~ normal(mu, sigma);

/*

生成量 { // 生成量块Stan代码

第一部分是数据块，我们告诉 Stan 它应该从数据列表中获得的数据。设置边界作为对数据输入的检查，这就是  。声明的前两个变量是  和 ，都是整数。接下来代码分别声明模型矩阵和目标向量。我们声明变量的类型和维度，然后声明其名称。在 Stan 中，在一个块中声明的所有内容都可用于后续块，但在一个块中声明的内容不会在更早的块中使用，例如声明  和 , 然后可以随后使用，就像我们指定模型矩阵的维度一样 。

作为参考，以下内容来注明了感兴趣的变量以及将在其中声明它们的相关块。

转换后的数据块是您可以执行对数或中心化等类似操作的地方，即您可以根据输入数据或一般情况下创建新数据。但是，如果您使用的是 R，那么首先在 R 中执行这些操作并将它们包含在数据列表中。您还可以在此处声明任何未建模的参数，例如您希望固定为某个值的参数。