共计 9246 个字符，预计需要花费 24 分钟才能阅读完成。
原文链接：http://tecdat.cn/?p=11161
最近咱们被客户要求撰写对于贝叶斯模型的钻研报告，包含一些图形和统计输入。
概率编程使咱们可能实现统计模型，而不用放心技术细节。这对于基于 MCMC 采样的贝叶斯模型特地有用
R 语言中 RStan 贝叶斯层次模型剖析示例
stan 简介
Stan 是用于贝叶斯推理的 C ++ 库。它基于 No-U-Turn 采样器（NUTS），该采样器用于依据用户指定的模型和数据预计后验散布。应用 Stan 执行剖析波及以下步骤：
- 应用 Stan 建模语言指定统计模型。通过专用的 \_.stan\_ 文件实现此操作。
- 筹备要提供给模型的数据。
- 应用该
stan函数从后验散布中采样。 - 剖析后果。
在本文中，我将通过两个层次模型展现 Stan 的用法。我将应用第一个模型探讨 Stan 的基本功能，并应用第二个示例演示更高级的利用。
学校数据集
咱们要应用的第一个数据集是 学校的数据集。该数据集掂量了教练打算对大学入学考试（在美国应用的学业能力测验（SAT））的影响。数据集如下所示：
正如咱们所看到的：对于八所学校中的大多数，短期教练打算确实进步了 SAT 分数。对于此数据集，咱们有趣味估算与每所学校相干的实在教练打算成果大小。咱们思考两种代替办法。首先，咱们能够假如所有学校彼此独立。然而，这将难以解释，因为学校的后验区间因为高标准差而在很大水平上重叠。第二，假如所有学校的实在成果都雷同，则能够汇总所有学校的数据。然而，这也是不合理的，因为该打算有针对学校的不同成果（例如，不同的老师和学生应该有不同的打算）。
因而，须要另一个模型。分层模型的长处是能够合并来自所有八所学校的信息，而无需假设它们具备独特的实在成果。咱们能够通过以下形式指定档次贝叶斯模型：
依据该模型，教练的成果遵循正态分布，其均值是实在成果 θj，其标准偏差为 σj（从数据中得悉）。真正的影响 θj 遵循参数 μ 和 τ 的正态分布。
定义 Stan 模型文件
在指定了要应用的模型之后，咱们当初能够探讨如何在 Stan 中指定此模型。在为上述模型定义 Stan 程序之前，让咱们看一下 Stan 建模语言的构造。
变量
在 Stan 中，能够通过以下形式定义变量：
int<lower=0> n; # 下界是 0int<upper=5> n; # 下限是 5int<lower=0,upper=5> n; # n 的范畴是 [0,5]
留神，如果先验已知变量，则应指定变量的高低边界。
多维数据能够通过方括号指定：
vector[n] numbers; // 长度为 n 的向量
real[n] numbers; // 长度为 n 的浮点数组
matrix[n,n] matrix; // n 乘 n 矩阵
程序
Stan 中应用以下程序：
- data：用于指定以贝叶斯规定为条件的数据
- 转换后的数据 ：用于预处理数据
- 参数 （必填）：用于指定模型的参数
- 转换后的参数 ：用于计算后验之前的参数解决
- 模型 （必填）：用于指定模型
- 生成数量 ：用于对后果进行后处理
点击题目查阅往期内容
MCMC 的 rstan 贝叶斯回归模型和规范线性回归模型比拟
左右滑动查看更多
01
02
03
04
对于 模型 程序块，能够两种等效形式指定散布。第一个，应用以下统计符号：
y ~ normal(mu, sigma); # y 遵从正态分布
第二种办法应用基于对数概率密度函数（lpdf）的程序化表示法：
target += normal_lpdf(y | mu, sigma); # 减少正态对数密度
Stan 反对大量的概率分布。通过 Stan 指定模型时，该 lookup 函数会派上用场：它提供从 R 函数到 Stan 函数的映射。思考以下示例：
library(rstan) # 加载 stan 包 lookup(rnorm)
<!—->
## StanFunction Arguments ReturnType Page
## 355 normal_rng (real mu, real sigma) real 494
在这里，咱们看到 R 中的 rnorm 等价于 Stan 的 normal_rng。
模型
当初，咱们理解了 Stan 建模语言的基础知识，咱们能够定义模型，并将其存储在一个名为的文件中 schools.stan：
留神，θ 永远不会呈现在参数中。这是因为咱们没有显式地对 θ 进行建模，而是对 η（各个学校的标准化成果）进行了建模。而后，依据 μ，τ 和 η 在 \_变换后的参数 \_局部结构 θ。此参数化使采样器更高效。
筹备数据进行建模
在拟合模型之前，咱们须要将输出数据编码为一个列表，其参数应与 Stan 模型的数据局部绝对应。对于学校数据，数据如下：
schools.data <- list(
n = 8,
y = c(28, 8, -3, 7, -1, 1, 18, 12),
sigma = c(15, 10, 16, 11, 9, 11, 10, 18)
)
从后验散布抽样
咱们能够应用 stan 函数从后验散布中采样，函数执行以下三个步骤：
- 它将模型标准转换为 C ++ 代码。
- 它将 C ++ 代码编译为共享对象。
- 它依据指定的模型，数据和设置从后验散布中采样。
如果 rstan_options(auto_write = TRUE)，则雷同模型的后续调用将比第一次调用快得多，因为该 stan 函数随后跳过了前两个步骤（转换和编译模型）。此外，咱们将设置要应用的内核数：
options(mc.cores = parallel::detectCores()) # 并行化 rstan_options(auto_write = TRUE) # 存储编译的 stan 模型
当初，咱们能够从后验中编译模型和样本。
模型解释
咱们将首先对模型进行根本解释，而后钻研 MCMC 程序。
根本模型解释
要应用拟合模型执行推断，咱们能够应用 print 函数。
print(fit1) # 可选参数：pars，probs
<!—->
## Inference for Stan model: schools.
## 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 7.67 0.15 5.14 -2.69 4.42 7.83 10.93 17.87 1185 1## tau 6.54 0.16 5.40 0.31 2.52 5.28 9.05 20.30 1157 1## eta[1] 0.42 0.01 0.92 -1.47 -0.18 0.44 1.03 2.18 4000 1## eta[2] 0.03 0.01 0.87 -1.74 -0.54 0.03 0.58 1.72 4000 1## eta[3] -0.18 0.02 0.92 -1.95 -0.81 -0.20 0.45 1.65 3690 1## eta[4] -0.03 0.01 0.92 -1.85 -0.64 -0.02 0.57 1.81 4000 1## eta[5] -0.33 0.01 0.86 -2.05 -0.89 -0.34 0.22 1.43 3318 1## eta[6] -0.20 0.01 0.87 -1.91 -0.80 -0.21 0.36 1.51 4000 1## eta[7] 0.37 0.02 0.87 -1.37 -0.23 0.37 0.96 2.02 3017 1## eta[8] 0.05 0.01 0.92 -1.77 -0.55 0.05 0.69 1.88 4000 1## theta[1] 11.39 0.15 8.09 -2.21 6.14 10.30 15.56 30.22 2759 1## theta[2] 7.92 0.10 6.25 -4.75 4.04 8.03 11.83 20.05 4000 1## theta[3] 6.22 0.14 7.83 -11.41 2.03 6.64 10.80 20.97 3043 1## theta[4] 7.58 0.10 6.54 -5.93 3.54 7.60 11.66 20.90 4000 1## theta[5] 5.14 0.10 6.30 -8.68 1.40 5.63 9.50 16.12 4000 1## theta[6] 6.08 0.10 6.62 -8.06 2.21 6.45 10.35 18.53 4000 1## theta[7] 10.60 0.11 6.70 -0.94 6.15 10.01 14.48 25.75 4000 1## theta[8] 8.19 0.14 8.18 -8.13 3.59 8.01 12.48 25.84 3361 1## lp__ -39.47 0.07 2.58 -45.21 -41.01 -39.28 -37.70 -34.99 1251 1##
## Samples were drawn using NUTS(diag_e) at Thu Nov 29 11:17:50 2018.## 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).
在此，行名称示意预计的参数：mu 是后验散布的平均值，而 tau 是其标准偏差。eta 和 theta 的条目别离示意矢量 η 和 θ 的估计值。这些列示意计算值。百分比示意置信区间。例如，教练打算的总体成果的 95%可信区间 μ 为 [-1.27,18.26]。因为咱们不确定平均值，因而 θj 的 95%置信区间也很宽。例如，对于第一所学校，95%置信区间为 [−2.19,32.33]。
咱们能够应用以下 plot 函数来可视化预计中的不确定性：
黑线示意 95%的距离，而红线示意 80%的距离。圆圈示意平均值的预计。
咱们能够应用以下 extract 函数获取生成的样本：
# 获取样本 samples <- extract(fit1, permuted = TRUE) # 每个参数 1000 个样本
MCMC 诊断
通过绘制采样过程的轨迹图，咱们能够确定采样期间是否出了问题。例如，链条在一个地位停留的工夫过长或在一个方向上走了太多步，就会有问题。咱们能够应用 traceplot 函数绘制模型中应用的四个链的轨迹：
# 诊断:
要从各个马尔可夫链中获取样本，咱们能够 extract 再次应用函数：
## parameters
## chains mu tau eta[1] eta[2] eta[3] eta[4]
## chain:1 1.111120 2.729124 -0.1581242 -0.8498898 0.5025965 -1.9874554## chain:2 3.633421 2.588945 1.2058772 -1.1173221 1.4830778 0.4838649## chain:3 13.793056 3.144159 0.6023924 -1.1188243 -1.2393491 -0.6118482## chain:4 3.673380 13.889267 -0.0869434 1.1900236 -0.0378830 -0.2687284## parameters
## chains eta[5] eta[6] eta[7] eta[8] theta[1]
## chain:1 0.3367602 -1.1940843 0.5834020 -0.08371249 0.6795797## chain:2 -1.8057252 0.7429594 0.9517675 0.55907356 6.7553706## chain:3 -1.5867789 0.6334288 -0.4613463 -1.44533007 15.6870727## chain:4 0.1028605 0.3481214 0.9264762 0.45331024 2.4657999## parameters
## chains theta[2] theta[3] theta[4] theta[5] theta[6] theta[7]
## chain:1 -1.208335 2.482769 -4.31289292 2.030181 -2.147684 2.703297## chain:2 0.740736 7.473028 4.88612054 -1.041502 5.556902 6.097494## chain:3 10.275294 9.896345 11.86930758 8.803971 15.784656 12.342510## chain:4 20.201935 3.147213 -0.05906019 5.102037 8.508530 16.541455## parameters
## chains theta[8] lp__
## chain:1 0.8826584 -41.21499## chain:2 5.0808317 -41.17178## chain:3 9.2487083 -40.35351## chain:4 9.9695268 -36.34043
为了对采样过程进行更高级的剖析，咱们能够应用该 shinystan 软件包。应用该软件包，能够通过以下形式启动 Shiny 应用程序来剖析拟合模型：
library(shinystan)launch_shinystan(fit1)
档次回归
当初，咱们对 Stan 有了根本的理解，咱们能够深入研究更高级的应用程序：让咱们尝试一下档次回归。在惯例回归中，咱们对以下模式的关系进行建模
此示意假如所有样本都具备雷同的散布。如果只存在一组样本，那么咱们就会遇到问题，因为将疏忽组内和组之间的潜在差别。
另一种抉择是为每个组建设一个回归模型。然而，在这种状况下，预计单个模型时，小样本量会带来问题。
档次回归是两个极其之间的折衷。该模型假如组是类似的，但存在差别。
假如每个样本都属于 K 组之一。而后，档次回归指定如下：
其中 Yk 是第 k 组的后果，αk 是截距，Xk 是特色，β（k）示意权重。层次模型不同于其中 Yk 别离拟合每个组的模型，因为假设参数 αk 和 β（k）源自独特的散布。
数据集
分层回归的经典示例是 老鼠数据集。该数据集蕴含 5 周内测得的 鼠体重。让咱们加载数据：
## day8 day15 day22 day29 day36
## 1 151 199 246 283 320## 2 145 199 249 293 354## 3 147 214 263 312 328## 4 155 200 237 272 297## 5 135 188 230 280 323## 6 159 210 252 298 331
让咱们考察数据：
library(ggplot2)ggplot(ddf, aes(x = variable, y = value, group = Group)) + geom_line() + geom_point()
数据显示线性增长趋势对于不同的大鼠十分类似。然而，咱们还看到，大鼠的初始体重不同，须要不同的截距，并且成长速度也须要不同的斜率。因而，分层模型仿佛是适当的。
档次回归模型的标准
该模型能够指定如下：
第 i 个大鼠的截距由 αi 示意，斜率由 βi 示意。留神，测量工夫的核心是 x = 22，它是工夫序列数据的中值测量值（第 22 天）。
当初，咱们能够指定模型并将其存储在名为 rats.stan 的文件中：
请留神，模型代码估算的是方差（sigmasq 变量）而不是标准差。
材料筹备
为了筹备模型数据，咱们首先将测量点提取为数值，而后将所有内容编码为列表构造：
data <- list(N = nrow(df), T = ncol(df), x = days,
y = df, xbar = median(days))
拟合回归模型
当初，咱们能够为老鼠体重数据集拟合贝叶斯档次回归模型：
# 模型蕴含截距（alpha）和斜率（beta）的预计
档次回归模型的预测
在确定了每只大鼠的 α 和 β 之后，咱们当初能够预计任意工夫点单个大鼠的体重。在这里，咱们寻找从第 0 天到第 100 天的大鼠体重。
ggplot(pred.df[pred.df$Rat %in% sel.rats,],
aes(x = Day, y = Weight, group = Rat,
geom_line() +
与原始数据相比，该模型的预计是平滑的，因为每条曲线都遵循线性模型。钻研最初一个图中所示的置信区间，咱们能够看到方差预计是正当的。咱们对采样时（第 8 至 36 天）的老鼠体重充满信心，然而随着来到采样区域，不确定性会减少。
点击文末 “浏览原文”
获取全文残缺材料。
本文选自《R 语言中的 Stan 概率编程 MCMC 采样的贝叶斯模型》。
点击题目查阅往期内容
R 语言贝叶斯 MCMC：用 rstan 建设线性回归模型剖析汽车数据和可视化诊断 \
【视频】马尔可夫链蒙特卡罗办法 MCMC 原理与 R 语言实现 | 数据分享 \
R 语言实现 MCMC 中的 Metropolis–Hastings 算法与吉布斯采样 \
R 语言贝叶斯 METROPOLIS-HASTINGS GIBBS 吉布斯采样器预计变点指数分布剖析泊松过程车站等待时间 \
R 语言马尔可夫 MCMC 中的 METROPOLIS HASTINGS，MH 算法抽样（采样）法可视化实例 \
python 贝叶斯随机过程：马尔可夫链 Markov-Chain，MC 和 Metropolis-Hastings，MH 采样算法可视化 \
Python 贝叶斯推断 Metropolis-Hastings（M-H）MCMC 采样算法的实现 \
Metropolis Hastings 采样和贝叶斯泊松回归 Poisson 模型 \
Matlab 用 BUGS 马尔可夫区制转换 Markov switching 随机稳定率模型、序列蒙特卡罗 SMC、M H 采样剖析工夫序列 R 语言 RSTAN MCMC：NUTS 采样算法用 LASSO 构建贝叶斯线性回归模型剖析职业声望数据 \
R 语言 BUGS 序列蒙特卡罗 SMC、马尔可夫转换随机稳定率 SV 模型、粒子滤波、Metropolis Hasting 采样工夫序列剖析 \
R 语言 Metropolis Hastings 采样和贝叶斯泊松回归 Poisson 模型 \
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 采样用于回归的贝叶斯预计
