共计 5138 个字符，预计需要花费 13 分钟才能阅读完成。
原文链接：http://tecdat.cn/?p=25453
最近咱们被客户要求撰写对于贝叶斯回归的钻研报告，包含一些图形和统计输入。
当初有了对贝叶斯办法的概念了解，咱们将理论钻研应用它的回归模型
为了简略起见，咱们从回归的规范线性模型开始。而后增加对采样散布或先验的更改。咱们将通过 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 ~ normal(0, 1) 示意 y 正态分布，均值为 0，标准差为 1。
此外，要装置 rstan，须要通过 CRAN 或 GitHub。它不须要独自装置 Stan 自身，但它的确须要几个步骤并且须要 C++ 编译器。一旦你装置了 rstan，它就会像任何其余 R 包一样被调用。
# 应用语法创立模型对象
stanmodelcode = "
data { // 数据块
int<lower=1> N; // 样本大小
int<lower=1> K; // 模型矩阵的尺寸
/*
转换后的数据 {// 转换后的数据块。目前不应用。}
*/
参数 // 参数块
real<lower=0> sigma; // 误差比例
}
模型 // 模型块
mu = X * beta; // 创立线性预测器
// 先验指标
sigma ~ cauchy(0, 5); // 因为 sigma 被限定在 0
// 似然
y ~ normal(mu, sigma);
/*
生成量 {// 生成量块Stan 代码
第一局部是数据块，咱们通知 Stan 它应该从数据列表中取得的数据。设置边界作为对数据输出的查看，这就是 < >。申明的前两个变量是 N 和 K，都是整数。接下来代码别离申明模型矩阵和指标向量。咱们申明变量的类型和维度，而后申明其名称。在 Stan 中，在一个块中申明的所有内容都可用于后续块，但在一个块中申明的内容不会在更早的块中应用，例如申明 N 和 K, 而后能够随后应用，就像咱们指定模型矩阵的维度一样 X。
作为参考，以下内容来注明了感兴趣的变量以及将在其中申明它们的相干块。
| 品种 | 申明块 |
|---|---|
| 建模的，未建模的数据 | 数据，转换后的数据 |
| 建模参数，缺失数据 | 参数，转换参数 |
| 未建模参数 | 数据，转换后的数据 |
| 产生量 | 转换数据、转换参数、生成量 |
| 循环索引 | 循环语句 |
转换后的数据块是您能够执行对数或中心化等相似操作的中央，即您能够依据输出数据或个别状况下创立新数据。然而，如果您应用的是 R，那么首先在 R 中执行这些操作并将它们蕴含在数据列表中。您还能够在此处申明任何未建模的参数，例如您心愿固定为某个值的参数。
点击题目查阅往期内容
R 语言贝叶斯 MCMC：用 rstan 建设线性回归模型剖析汽车数据和可视化诊断
左右滑动查看更多
01
02
03
04
要预计的次要感兴趣的参数位于参数块中。与数据块一样，您只能申明这些变量，不能进行任何赋值。在这里，咱们留神到要预计的 β 和 σ，后者的上限为零。在实践中，如果截距或其余系数在显着不同的尺度上，您可能更违心将它们离开建模。
转换后的参数块是可能蕴含可选参数的中央。为了提高效率，您通常只想搁置依赖于参数块的特定趣味的货色。
模型块是指定您的先验和可能性以及任何必要变量的申明的中央。例如，此处蕴含线性预测器，因为它将趋向于似然. 请留神，咱们能够将线性预测器放在转换后的参数局部，但这会减慢过程，而且咱们对这些特定值不太感兴趣。
我对系数应用的是正态先验，平均值为零，标准差很大。对于 σ 的预计，我应用的是 Cauchy 散布。许多应用 BUGS 的回归例子都会应用反伽马先验，这对这个模型来说是齐全能够的，只管它对其余方差参数的成果并不现实。如果咱们没有为参数的先验散布指定任何货色，均匀分布是默认的。
最初，你想计算的任何货色都能够放在这里 – 对新数据的预测、参数的比率、参数大于 x 的次数、为报告目标的参数转换等等。咱们稍后将对此进行演示。
运行模型
当初咱们对代码的作用有了一个概念。贝叶斯预计，像最大似然法一样，以初始猜想为终点，而后以迭代的形式运行，每一步都从后验散布中产生模仿抽样，而后纠正这些抽样，直到最初达到某个指标，或安稳散布。这一部分是要害，与经典的统计学不同。咱们的指标是一个散布，而不是一个点估计。
这个模仿过程被称为马尔科夫链蒙特卡洛，简称 MCMC。这个过程的具体细节使许多贝叶斯编程语言 / 办法不同凡响。在 MCMC 中，所有来自后验的模仿抽样都是基于以前的抽样并与之相干的，因为这个过程是沿着走向安稳散布的路线后退的。咱们通常会让这个过程预烧，或者说从最后的终点开始有点稳定下来，这可能会有很大的偏差，因而在最后的几次迭代中，后续的估计值也会有很大的偏差。然而，作为进一步的查看，咱们将屡次运行整个过程，也就是说，有一个以上的链。因为这些链将从不同的中央开始，如果多个链最初都达到了同一个后果，咱们就能够对后果更有信念。
你会留神到 Stan 将其代码编译为 C ++ 的工夫可能比运行模型的工夫要长，而在我的电脑上，每条链只须要一秒钟多一点的工夫。然而，贝叶斯办法已经须要很长的工夫，即便是像这样的规范回归，这兴许是贝叶斯剖析在过来几十年里才流行起来的次要起因；咱们基本没有机器来无效地做这件事。即便是当初，对于高度简单的模型和大数据集来说，它依然须要很长的工夫来运行，只管通常不会太长。
在上面的代码中，咱们留神到蕴含 stan 模型代码的对象，数据列表，咱们想要多少次迭代（5000），咱们想要这个过程在开始保留任何估计值之前运行多长时间（warmup=2500），咱们想要保留多少次后验的抽取（thin=10 意味着每十次抽取），以及链的数量（chains=4）。最初，咱们将有四条链，从参数的后验散布中抽取 1000 次。即便在 verbose = FALSE 的状况下，Stan 也会向 R 控制台输入大量的输入，我在这里省略，然而你会看到一些对于编译过程的初始信息，当每条链通过 iter 参数中指定的 10% 的迭代时的更新，以及最初对耗时的预计。
library(rstan)
### 运行模型并查看后果
fit (
warmup = 2500,
chains = 4)随着模型的运行，咱们当初能够查看后果。在下文中，咱们指定要显示的数字精度，咱们想要哪些参数，以及咱们想要哪些后验抽样的量级，在本例中是中位数和那些会产生 95% 区间预计的参数。
# 摘要
print(fit到目前为止还不错。均匀估计值反映了感兴趣的参数的后验后果的平均值，是规范回归剖析中报告的典型系数。值得注意的是 95% 的概率或置信区间，因为它们不是你所晓得的置信区间。这里没有反复抽样的解释。概率区间是更直观的。它的意思很简略，依据这个模型的后果，实在值有 95% 的可能性会落在这两点之间。
将这些后果与 R 的 lm 函数的后果相比拟，咱们能够看到咱们失去了相似的估计值，因为它们在小数点后两位是雷同的。事实上，如果咱们应用对立先验，模型与用规范最大似然预计所做的基本相同。
summary然而咱们怎么晓得咱们的模型是否运作良好呢？有几个规范的诊断办法，但让咱们看一下目前的一些状况。在摘要中，se\_mean 是蒙特卡洛误差，是对只有无限数量的后验抽样所带来的不确定性的预计。n\_eff 是给定所有链的无效样本量，基本上占了链中的自相干，即当咱们从一次抽样到下一次抽样时预计的相关性。它实际上不须要很大，但如果它绝对于所需的总抽样数来说很小，那就可能引起关注了。Rhat 是掂量链的混合水平的指标，当链被容许运行有限次抽样时，它就会变成 1。在这种状况下，n_eff 和 Rhat 表明咱们有很好的收敛性，但咱们也能够用跟踪图来直观地查看。
# 可视化
srace(fit'))咱们能够看到跟踪图中，每条链的估计值都能很快地从终点找到一个或多或少的稳固状态（灰色的初始预烧迭代）。此外，所有三条链（每条链都用不同的色彩示意）都混合得很好，并在同一论断左近跳动。
stan 开发团队通过 shinystan 包使交互式摸索诊断变得很容易。此外，coda 包中还有其余诊断办法，Stan 模型的后果能够很容易地转换为与之配合。上面的代码演示了如何开始。
bets = extract$beta除了制作数据列表和产生特定语言的模型代码的初始设置之外，绝对于规范模型，运行贝叶斯回归模型并不一定须要太多的工夫。最次要的兴许是采纳不同的思维形式，并从这个新的角度来解释后果。对于你所相熟的规范模型，可能不须要太长的工夫就能像应用那些模型一样自若。
本文摘选 《 R 语言 MCMC 的 rstan 贝叶斯回归模型和规范线性回归模型比拟 》，点击“ 浏览原文”获取全文残缺材料。
点击题目查阅往期内容
R 语言 RSTAN MCMC：NUTS 采样算法用 LASSO 构建贝叶斯线性回归模型剖析职业声望数据
R 语言 STAN 贝叶斯线性回归模型剖析气候变化影响北半球海冰范畴和可视化查看模型收敛性
R 语言贝叶斯 MCMC：用 rstan 建设线性回归模型剖析汽车数据和可视化诊断
R 语言贝叶斯 MCMC：GLM 逻辑回归、Rstan 线性回归、Metropolis Hastings 与 Gibbs 采样算法实例
R 语言贝叶斯 Poisson 泊松 - 正态分布模型剖析职业足球比赛进球数
R 语言用 Rcpp 减速 Metropolis-Hastings 抽样预计贝叶斯逻辑回归模型的参数
R 语言逻辑回归、Naive Bayes 贝叶斯、决策树、随机森林算法预测心脏病
R 语言中贝叶斯网络（BN）、动静贝叶斯网络、线性模型剖析错颌畸形数据
R 语言中的 block Gibbs 吉布斯采样贝叶斯多元线性回归
Python 贝叶斯回归剖析住房累赘能力数据集
R 语言实现贝叶斯分位数回归、lasso 和自适应 lasso 贝叶斯分位数回归剖析
Python 用 PyMC3 实现贝叶斯线性回归模型
R 语言用 WinBUGS 软件对学术能力测验建设档次（分层）贝叶斯模型
R 语言 Gibbs 抽样的贝叶斯简略线性回归仿真剖析
R 语言和 STAN,JAGS：用 RSTAN,RJAG 建设贝叶斯多元线性回归预测选举数据
R 语言基于 copula 的贝叶斯分层混合模型的诊断准确性钻研
R 语言贝叶斯线性回归和多元线性回归构建工资预测模型
R 语言贝叶斯推断与 MCMC：实现 Metropolis-Hastings 采样算法示例
R 语言 stan 进行基于贝叶斯推断的回归模型
R 语言中 RStan 贝叶斯层次模型剖析示例
R 语言应用 Metropolis-Hastings 采样算法自适应贝叶斯预计与可视化
R 语言随机搜寻变量抉择 SSVS 预计贝叶斯向量自回归（BVAR）模型
WinBUGS 对多元随机稳定率模型：贝叶斯预计与模型比拟
R 语言实现 MCMC 中的 Metropolis–Hastings 算法与吉布斯采样
R 语言贝叶斯推断与 MCMC：实现 Metropolis-Hastings 采样算法示例
R 语言应用 Metropolis-Hastings 采样算法自适应贝叶斯预计与可视化
视频：R 语言中的 Stan 概率编程 MCMC 采样的贝叶斯模型
R 语言 MCMC:Metropolis-Hastings 采样用于回归的贝叶斯预计
