数据集来自 Gelfand et al (1990) [^Alan-Gelfand]，研究 30 只小鼠的体重变化，每周测量其体重，连续测量了5周。用$Y_{ij}$表示第$i$只鼠在时间$x_j$测量的体重。部分数据集见下表，其中 Week 1 表示小鼠一周大 $x_1 = 8$，Week 2 表示小鼠两周大 $x_2 = 15$ ，依次类推。
N <- 30
T <- 5
y <- structure(c(
151, 145, 147, 155, 135, 159, 141, 159, 177, 134,
160, 143, 154, 171, 163, 160, 142, 156, 157, 152, 154, 139, 146,
157, 132, 160, 169, 157, 137, 153, 199, 199, 214, 200, 188, 210,
189, 201, 236, 182, 208, 188, 200, 221, 216, 207, 187, 203, 212,
203, 205, 190, 191, 211, 185, 207, 216, 205, 180, 200, 246, 249,
263, 237, 230, 252, 231, 248, 285, 220, 261, 220, 244, 270, 242,
248, 234, 243, 259, 246, 253, 225, 229, 250, 237, 257, 261, 248,
219, 244, 283, 293, 312, 272, 280, 298, 275, 297, 350, 260, 313,
273, 289, 326, 281, 288, 280, 283, 307, 286, 298, 267, 272, 285,
286, 303, 295, 289, 258, 286, 320, 354, 328, 297, 323, 331, 305,
338, 376, 296, 352, 314, 325, 358, 312, 324, 316, 317, 336, 321,
334, 302, 302, 323, 331, 345, 333, 316, 291, 324
), .Dim = c(30, 5))
x <- c(8.0, 15.0, 22.0, 29.0, 36.0)
xbar <- 22.0 # x 的平均值
# 第一张表格
rownames(y) <- paste("Rat", seq(30))
knitr::kable(y,
format = "markdown", caption = "Rats 数据集",
row.names = TRUE,
col.names = paste("Week", seq(5))
)
| Week 1 | Week 2 | Week 3 | Week 4 | Week 5 | |
|---|---|---|---|---|---|
| Rat 1 | 151 | 199 | 246 | 283 | 320 |
| Rat 2 | 145 | 199 | 249 | 293 | 354 |
| Rat 3 | 147 | 214 | 263 | 312 | 328 |
| … | … | … | … | … | … |
| Rat 28 | 157 | 205 | 248 | 289 | 316 |
| Rat 29 | 137 | 180 | 219 | 258 | 291 |
| Rat 30 | 153 | 200 | 244 | 286 | 324 |
为了比较直观地观察数据集，我们绘制如下图形
# 第一张图
png(filename = "rats.png", res = 200, width = 480*200/72, height = 480*200/72,
type = "cairo", bg = "transparent")
par(mar = c(4, 4, 0.5, 0.5))
matplot(t(y),
col = hcl.colors(30), type = "b",
xlab = "Weeks", ylab = "Weight", lty = 2, pch = 1
)
dev.off()
根据图中的显示，我们建立带随机效应的线性模型
$$ Y_{ij} \sim \mathsf{Normal}(\alpha_i + \beta_i(X_{j} - X_{bar}), \tau_{c})\\ \alpha_i \sim \mathsf{Normal}(\alpha_{c}, \tau_{\alpha}) \\ \beta_j \sim \mathsf{Normal}(\beta_{c},\tau_{\beta}) \\ i = 1, 2, \ldots, 30;\quad j = 1, 2, \ldots, 5. $$
根据的先验选择，我们设置的先验如下
$$ \alpha_{c} \sim \mathsf{Normal}(0, 100) \\ \beta_{c} \sim \mathsf{Normal}(0, 100) \\ \tau^2_{c} \sim \mathsf{inv\_{gamma}}(0.001, 0.001) \\ \tau^2_{\alpha} \sim \mathsf{inv\_{gamma}}(0.001, 0.001) \\ \tau^2_{\beta} \sim \mathsf{inv\_{gamma}}(0.001, 0.001). $$
- 为什么选择这些先验？后验分布的具体形式能推导吗？参数的似然函数能推导吗？[^hnm-with-mcmc] 见安德鲁的书籍 BDA 第三版
- 它们是如何表示成 Stan 代码中的模型块？见 Stan 文档
- Stan 是如何编码计算这些统计模型的，自动微分，变分推断？见 Stan 文档
- 分层正态模型的 JAGS nimble Stan 实现，参考 bridgesampling 包的文档，该 R 包基于 Xiao-Li Meng and Wing Hung Wong (1996) Simulating ratios of normallizing constants via A simple identity: a theoretical exploration http://www3.stat.sinica.edu.tw/statistica/j6n4/j6n43/j6n43.htm 提出的 Bridge Sampling 比较贝叶斯模型，还可以计算边际似然、贝叶斯因子、后验模型概率，bfrms 和 BayesFactor 可用于计算贝叶斯因子，前者适用范围更广
加载 rstan 包后，我们需要设置一些模型运行时和模型结果展示阶段的参数
# BUGS rats example (Vol 1, Example 1)
# http://www.openbugs.net/Examples/Rats.html
library(rstan, quietly = TRUE)
# 将编译后的模型写入磁盘，避免重新编译
rstan_options(auto_write = TRUE)
# 如果CPU和内存足够，设置成与马尔科夫链一样多
options(mc.cores = 2)
# Windows 上可能需要指定
# Sys.setenv(USE_Cxx14 = 1)
custom_colors <- c(
"#4285f4", # GoogleBlue
"#34A853", # GoogleGreen
"#FBBC05", # GoogleYellow
"#EA4335" # GoogleRed
)
rstan_ggtheme_options(
panel.background = element_rect(fill = "white"),
legend.position = "top"
)
rstan_gg_options(
fill = "#4285f4", color = "white",
pt_color = "#EA4335", chain_colors = custom_colors
)
后续以Stan编写的模型中，每个参数迭代的链条数为4，即从后验分布中抽样的次数为4，每次抽样的数量为 1000
chains = 4 # 链条数
iter = 1000 # 样本量
接着设置模型参数初始值，每条链开始迭代的初值值都一样，这里选用 WinBUGS 给的初始值，后续我们还考虑随机选择的初始值与之比较
# 复制 list 的姿势 https://d.cosx.org/d/420668
init <- rep(list(list(
alpha = rep(250, 30), beta = rep(6, 30),
alpha_c = 150, beta_c = 10,
tausq_c = 1, tausq_alpha = 1,
tausq_beta = 1
)), chains)
调用 rstan 包中的关键函数 stan 编译运行 Stan 编写的分层正态模型
file- 指定模型文件的路径
model_name- Stan 代码编译成功后，指定的模型名称
model_code- 传递一段字符串，字符串的内容即是 Stan 代码编写的模型
data- 参数传递模型的观测数据集，注意数据结构是 R 语言的列表 list，其每个元素的属性需和 Stan 代码编写的模型中设定的参数属性保持一致
algorithm- Stan 可选的抽样算法有
c("NUTS", "HMC", "Fixed_param")默认和推荐的参数值是"NUTS"，这三者的算法细节及使用情况以后介绍。 chains- 参数迭代时设置的马氏链条数目
init- 模型中各个参数的初始值，也是用列表来存储，列表中的元素属性和 Stan 代码中的参数也要保持一致，参数名称也要一致，默认值
random表示随机选择 iter- 模型每条链迭代的次数，在这里设置了四条链，迭代总计4000次，抽样间隔默认是1，预处理阶段
thin- 抽样间隔，即每迭代多少次抽一组迭代值，默认是1，就是迭代多少次就抽取多少次
warmup- 预处理阶段迭代的次数，默认设置
floor(iter/2)就是每条链迭代总次数的一半，向下取整 seed- 模型运行时设置的随机数种子，Stan 内置的算法都是基于 MCMC 的，为了复现结果，设置相同的随机数种子有时是必要的
cores- 设置模型运行时可使用的系统 CPU 的核心数，在内存和 CPU 资源充足的情况下，越多越好。注意，可用的 CPU 核心，可用
parallel::detectCores()检测出来，不要超过这个，此外，数目也不要超过参数chains设置的数值，超过这个意义也不大。其默认值是getOption("mc.cores", 1L)即获取R语言环境中全局参数列表中的mc.cores，若无设置就取1 verbose- 是否显示模型编译和运行过程中的状态，默认是 FALSE 即不显示中间的迭代过程
control- 传递一个列表，包含算法的细节调控，调整采样的策略，在采样效率和准确率中寻找一个平衡，此处不详细介绍，留待详细介绍的算法细节的时候介绍
fit- 参数传递一个类型为 S4 的
stanfit对象实例，默认是NA，如果参数值不是NA而是前一次相关模型拟合的结果，那么模型拟合就会重用之前的结果，可以节省下重编译 Stan 代码的时间 boost_lib = NULL, eigen_lib = NULL- Stan 语言大量使用了 boost 和 eigen 两个 C++ 库，模型编译过程允许使用用户指定的版本
save_dso- 参数默认值是
TRUE，该参数只有在参数fit = NA时才有作用。意义是从 C++ 编译的动态分享对象 DSO (dynamic shared object) 是否需要保留。如果保留，那么从另一个 R 控制台运行模型时，不再重复编译 C++ 代码，而是直接使用保留的 DSO，这在并行的时候是可以节省编译时间的。 ......这种参数在 R 语言当中是比较高级的传参方式了，其实际是传递可变长度的参数向量。其中，参数refresh控制采样进度的报告频率，默认值是max(iter/10, 1)，取值 0 表示不显示采样的进度；save_warmup是否保留预处理阶段的采样过程，默认参数值TRUE表示保留。还有参数chain_id和init_r、test_grad等，这里不一一介绍了，其它参数可见帮助文档?rstan::stan
rats_fit <- stan(
# file = "code/bugs/vol1/rats/rats.stan",
model_name = "rats",
model_code = "
data {
int<lower=0> N;
int<lower=0> T;
real x[T];
real y[N,T];
real xbar;
}
parameters {
real alpha[N];
real beta[N];
real alpha_c;
real beta_c; // beta.c in original bugs model
real<lower=0> tausq_c;
real<lower=0> tausq_alpha;
real<lower=0> tausq_beta;
}
transformed parameters {
real<lower=0> tau_c; // sigma in original bugs model
real<lower=0> tau_alpha;
real<lower=0> tau_beta;
tau_c = sqrt(tausq_c);
tau_alpha = sqrt(tausq_alpha);
tau_beta = sqrt(tausq_beta);
}
model {
alpha_c ~ normal(0, 100);
beta_c ~ normal(0, 100);
tausq_c ~ inv_gamma(0.001, 0.001);
tausq_alpha ~ inv_gamma(0.001, 0.001);
tausq_beta ~ inv_gamma(0.001, 0.001);
alpha ~ normal(alpha_c, tau_alpha); // vectorized
beta ~ normal(beta_c, tau_beta); // vectorized
for (n in 1:N)
for (t in 1:T)
y[n,t] ~ normal(alpha[n] + beta[n] * (x[t] - xbar), tau_c);
}
generated quantities {
real alpha0;
alpha0 = alpha_c - xbar * beta_c;
}
",
data = list(N = N, T = T, y = y, x = x, xbar = xbar), # 观测数据集
chains = chains, init = init, iter = iter)
分层正态模型 Stan 拟合结果
print(rats_fit, digits = 1)
Inference for Stan model: rats.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
alpha[1] 240.0 0.1 2.7 234.7 238.2 240.0 241.8 245.1 2379 1
alpha[2] 247.7 0.1 2.7 242.3 245.9 247.7 249.4 253.2 2293 1
alpha[3] 252.4 0.1 2.8 246.7 250.6 252.4 254.3 257.8 2296 1
alpha[4] 232.6 0.1 2.7 227.4 230.9 232.7 234.5 237.7 2553 1
alpha[5] 231.6 0.1 2.7 226.2 229.8 231.6 233.4 237.0 2443 1
alpha[6] 249.8 0.1 2.6 244.7 248.0 249.7 251.6 254.7 2171 1
alpha[7] 228.7 0.1 2.6 223.5 226.9 228.6 230.4 234.0 1980 1
alpha[8] 248.3 0.1 2.7 242.7 246.6 248.2 250.0 253.7 2162 1
alpha[9] 283.3 0.1 2.7 277.9 281.5 283.4 285.2 288.5 2060 1
alpha[10] 219.2 0.0 2.5 214.3 217.5 219.2 220.9 224.2 2708 1
alpha[11] 258.3 0.0 2.7 253.1 256.5 258.2 260.1 263.7 2995 1
alpha[12] 228.1 0.1 2.6 223.0 226.3 228.1 229.8 233.3 1985 1
alpha[13] 242.5 0.1 2.7 237.4 240.7 242.4 244.2 247.9 1926 1
alpha[14] 268.3 0.1 2.7 262.9 266.5 268.3 270.1 273.4 2332 1
alpha[15] 242.7 0.1 2.8 237.4 240.8 242.7 244.6 248.3 2225 1
alpha[16] 245.3 0.1 2.6 240.3 243.6 245.3 247.0 250.4 2365 1
alpha[17] 232.2 0.1 2.6 226.8 230.5 232.2 233.9 237.2 2353 1
alpha[18] 240.5 0.1 2.7 235.2 238.7 240.5 242.3 245.7 2023 1
alpha[19] 253.8 0.1 2.7 248.5 252.0 253.8 255.6 258.8 1920 1
alpha[20] 241.6 0.1 2.7 236.2 239.8 241.7 243.5 246.9 2394 1
alpha[21] 248.5 0.1 2.6 243.2 246.8 248.6 250.3 253.5 2314 1
alpha[22] 225.3 0.1 2.7 220.1 223.4 225.2 227.2 230.7 1805 1
alpha[23] 228.5 0.1 2.7 223.6 226.7 228.4 230.2 233.7 2334 1
alpha[24] 245.1 0.1 2.7 239.7 243.3 245.0 246.7 250.5 2420 1
alpha[25] 234.5 0.1 2.8 229.1 232.6 234.6 236.3 239.9 1869 1
alpha[26] 253.9 0.1 2.8 248.5 252.0 254.0 255.8 259.4 1879 1
alpha[27] 254.3 0.1 2.6 249.4 252.5 254.3 256.0 259.3 2400 1
alpha[28] 242.9 0.1 2.7 237.8 241.1 242.9 244.7 248.3 2405 1
alpha[29] 217.9 0.1 2.7 212.5 216.0 217.9 219.6 223.2 2335 1
alpha[30] 241.4 0.1 2.8 235.9 239.6 241.4 243.1 246.9 2132 1
beta[1] 6.1 0.0 0.2 5.6 5.9 6.1 6.2 6.5 2285 1
beta[2] 7.1 0.0 0.2 6.6 6.9 7.0 7.2 7.5 1882 1
beta[3] 6.5 0.0 0.3 6.0 6.3 6.5 6.6 7.0 2247 1
beta[4] 5.3 0.0 0.3 4.8 5.2 5.3 5.5 5.9 1729 1
beta[5] 6.6 0.0 0.2 6.1 6.4 6.6 6.7 7.0 2280 1
beta[6] 6.2 0.0 0.2 5.7 6.0 6.2 6.3 6.7 2230 1
beta[7] 6.0 0.0 0.2 5.5 5.8 6.0 6.1 6.4 2666 1
beta[8] 6.4 0.0 0.2 5.9 6.3 6.4 6.6 6.9 1885 1
beta[9] 7.1 0.0 0.3 6.5 6.9 7.0 7.2 7.6 2122 1
beta[10] 5.8 0.0 0.2 5.4 5.7 5.8 6.0 6.3 2389 1
beta[11] 6.8 0.0 0.2 6.3 6.6 6.8 7.0 7.3 2002 1
beta[12] 6.1 0.0 0.2 5.6 6.0 6.1 6.3 6.6 2128 1
beta[13] 6.2 0.0 0.2 5.7 6.0 6.2 6.3 6.6 2698 1
beta[14] 6.7 0.0 0.2 6.2 6.5 6.7 6.9 7.2 1970 1
beta[15] 5.4 0.0 0.3 4.9 5.2 5.4 5.6 5.9 1975 1
beta[16] 5.9 0.0 0.2 5.4 5.8 5.9 6.1 6.4 2187 1
beta[17] 6.3 0.0 0.2 5.8 6.1 6.3 6.4 6.7 2442 1
beta[18] 5.8 0.0 0.2 5.4 5.7 5.8 6.0 6.3 2299 1
beta[19] 6.4 0.0 0.2 5.9 6.2 6.4 6.6 6.9 2126 1
beta[20] 6.1 0.0 0.2 5.6 5.9 6.1 6.2 6.5 2023 1
beta[21] 6.4 0.0 0.2 5.9 6.2 6.4 6.6 6.9 1801 1
beta[22] 5.9 0.0 0.2 5.4 5.7 5.9 6.0 6.3 2136 1
beta[23] 5.7 0.0 0.3 5.3 5.6 5.8 5.9 6.2 1862 1
beta[24] 5.9 0.0 0.2 5.4 5.7 5.9 6.1 6.4 2427 1
beta[25] 6.9 0.0 0.3 6.4 6.7 6.9 7.1 7.4 2053 1
beta[26] 6.5 0.0 0.2 6.1 6.4 6.5 6.7 7.0 2407 1
beta[27] 5.9 0.0 0.2 5.4 5.7 5.9 6.1 6.4 2528 1
beta[28] 5.8 0.0 0.2 5.4 5.7 5.8 6.0 6.3 2228 1
beta[29] 5.7 0.0 0.2 5.2 5.5 5.7 5.8 6.2 2038 1
beta[30] 6.1 0.0 0.2 5.7 6.0 6.1 6.3 6.6 2378 1
alpha_c 242.5 0.1 2.8 236.4 240.6 242.5 244.3 248.0 1936 1
beta_c 6.2 0.0 0.1 6.0 6.1 6.2 6.3 6.4 1565 1
tausq_c 37.3 0.2 5.7 27.8 33.2 36.8 41.0 49.4 1178 1
tausq_alpha 220.2 1.6 66.6 125.7 173.9 207.8 251.5 378.8 1812 1
tausq_beta 0.3 0.0 0.1 0.1 0.2 0.3 0.3 0.5 1207 1
tau_c 6.1 0.0 0.5 5.3 5.8 6.1 6.4 7.0 1171 1
tau_alpha 14.7 0.0 2.1 11.2 13.2 14.4 15.9 19.5 1919 1
tau_beta 0.5 0.0 0.1 0.4 0.5 0.5 0.6 0.7 1161 1
alpha0 106.5 0.1 3.7 98.9 104.1 106.6 109.0 113.6 1540 1
lp__ -438.4 0.3 6.7 -452.3 -442.8 -438.1 -433.7 -426.0 567 1
Samples were drawn using NUTS(diag_e) at Mon May 13 18:53:46 2019.
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).
用图表示参数向量 $\alpha,\beta$ 各个分量的后验分布
# 保存图形
png(filename = "rats-alpha.png", res = 200, width = 480*200/72, height = 480*200/72,
type = "cairo", bg = "transparent")
plot(rats_fit, pars="alpha")
dev.off()
png(filename = "rats-beta.png", res = 200, width = 480*200/72, height = 480*200/72,
type = "cairo", bg = "transparent")
plot(rats_fit, pars="beta")
dev.off()
各个参数的迭代轨迹图全部画出，实在太多，但是作为重要的环节，一定要逐一验证，每个参数的迭代过程都是平稳的，只有基于平稳的马氏链才可以做出合理的后验估计。下面以参数 $\alpha_1$ 为例
png(filename = "rats-traceplot-alpha1.png", res = 200, width = 480*200/72, height = 480*200/72,
type = "cairo", bg = "transparent")
traceplot(rats_fit, pars = c("alpha[1]"), inc_warmup = FALSE) +
theme_minimal() +
labs(x = "Iteration", y = expression(alpha[1]))
dev.off()
考虑使用随机生成的初始值，为了节省版面，常常将 Stan 代码保存到磁盘文件，这里就不妨将前面的 model_code 参数值保存到文件 rats.stan
# 使用随机生成的初始值
rats <- stan(
#file = "code/bugs/vol1/rats/rats.stan",
model_code = "
data {
int<lower=0> N;
int<lower=0> T;
real x[T];
real y[N,T];
real xbar;
}
parameters {
real alpha[N];
real beta[N];
real alpha_c;
real beta_c; // beta.c in original bugs model
real<lower=0> tausq_c;
real<lower=0> tausq_alpha;
real<lower=0> tausq_beta;
}
transformed parameters {
real<lower=0> tau_c; // sigma in original bugs model
real<lower=0> tau_alpha;
real<lower=0> tau_beta;
tau_c = sqrt(tausq_c);
tau_alpha = sqrt(tausq_alpha);
tau_beta = sqrt(tausq_beta);
}
model {
alpha_c ~ normal(0, 100);
beta_c ~ normal(0, 100);
tausq_c ~ inv_gamma(0.001, 0.001);
tausq_alpha ~ inv_gamma(0.001, 0.001);
tausq_beta ~ inv_gamma(0.001, 0.001);
alpha ~ normal(alpha_c, tau_alpha); // vectorized
beta ~ normal(beta_c, tau_beta); // vectorized
for (n in 1:N)
for (t in 1:T)
y[n,t] ~ normal(alpha[n] + beta[n] * (x[t] - xbar), tau_c);
}
generated quantities {
real alpha0;
alpha0 = alpha_c - xbar * beta_c;
}
",
data = list(N = N, T = T, y = y, x = x, xbar = xbar),
chains = chains, iter = iter, init = "random"
)
Stan 模型在如前所述的设置中，运行上面的代码块是非常快的，因为之前编译的模型没有变，只是改变了算法迭代的初始值，可以直接使用
print(rats, digits = 1)
Inference for Stan model: rats.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
alpha[1] 239.9 0.1 2.7 234.7 238.1 239.8 241.6 245.3 1628 1
alpha[2] 247.7 0.1 2.6 242.6 245.9 247.7 249.4 252.6 2227 1
alpha[3] 252.4 0.1 2.7 247.0 250.6 252.3 254.2 257.7 2812 1
alpha[4] 232.6 0.1 2.7 227.5 230.7 232.5 234.3 237.9 2284 1
alpha[5] 231.7 0.1 2.6 226.5 229.8 231.7 233.5 236.8 2229 1
alpha[6] 249.7 0.1 2.7 244.6 247.9 249.7 251.5 255.0 2098 1
alpha[7] 228.7 0.1 2.7 223.1 226.9 228.7 230.5 234.0 1874 1
alpha[8] 248.4 0.1 2.7 243.2 246.6 248.4 250.2 253.6 2736 1
alpha[9] 283.3 0.1 2.7 278.0 281.6 283.4 285.2 288.3 2257 1
alpha[10] 219.3 0.1 2.7 214.1 217.5 219.3 221.1 224.4 2117 1
alpha[11] 258.2 0.1 2.8 252.6 256.4 258.1 260.0 263.8 1893 1
alpha[12] 228.1 0.1 2.7 222.8 226.4 228.1 229.9 233.6 1975 1
alpha[13] 242.4 0.1 2.6 237.4 240.8 242.4 244.1 247.6 2328 1
alpha[14] 268.2 0.1 2.6 263.4 266.5 268.2 270.0 273.2 2181 1
alpha[15] 242.8 0.1 2.7 237.6 240.9 242.8 244.7 248.2 1262 1
alpha[16] 245.3 0.1 2.7 240.1 243.4 245.3 247.1 250.5 2137 1
alpha[17] 232.2 0.1 2.7 227.0 230.4 232.2 234.1 237.4 2480 1
alpha[18] 240.3 0.1 2.6 234.9 238.7 240.3 242.1 245.5 2476 1
alpha[19] 253.8 0.1 2.6 248.6 252.0 253.8 255.5 259.0 2621 1
alpha[20] 241.5 0.1 2.7 236.3 239.7 241.6 243.4 247.0 2086 1
alpha[21] 248.6 0.1 2.7 243.2 246.8 248.5 250.4 253.9 1916 1
alpha[22] 225.2 0.1 2.6 219.9 223.5 225.2 226.9 230.3 2170 1
alpha[23] 228.6 0.1 2.5 223.7 226.9 228.6 230.2 233.6 2365 1
alpha[24] 245.1 0.1 2.6 240.0 243.3 245.1 246.8 249.9 2269 1
alpha[25] 234.5 0.1 2.7 229.4 232.6 234.4 236.3 240.0 2272 1
alpha[26] 253.9 0.1 2.7 248.6 252.0 253.9 255.6 259.4 2608 1
alpha[27] 254.4 0.1 2.7 249.2 252.6 254.3 256.1 259.8 2714 1
alpha[28] 242.9 0.1 2.6 237.8 241.2 242.9 244.7 248.0 2290 1
alpha[29] 217.9 0.1 2.7 212.8 216.1 217.9 219.8 223.4 1949 1
alpha[30] 241.4 0.1 2.8 235.8 239.5 241.4 243.3 246.7 2101 1
beta[1] 6.1 0.0 0.2 5.6 5.9 6.1 6.2 6.5 2438 1
beta[2] 7.1 0.0 0.3 6.6 6.9 7.0 7.2 7.5 1905 1
beta[3] 6.5 0.0 0.2 6.0 6.3 6.5 6.6 6.9 2662 1
beta[4] 5.3 0.0 0.3 4.8 5.2 5.3 5.5 5.8 1703 1
beta[5] 6.6 0.0 0.2 6.1 6.4 6.6 6.7 7.0 2734 1
beta[6] 6.2 0.0 0.2 5.7 6.0 6.2 6.3 6.6 2530 1
beta[7] 6.0 0.0 0.2 5.5 5.8 6.0 6.1 6.5 2470 1
beta[8] 6.4 0.0 0.2 5.9 6.3 6.4 6.6 6.9 2131 1
beta[9] 7.0 0.0 0.3 6.5 6.9 7.0 7.2 7.5 2072 1
beta[10] 5.8 0.0 0.2 5.4 5.7 5.8 6.0 6.3 2436 1
beta[11] 6.8 0.0 0.2 6.3 6.6 6.8 7.0 7.3 2030 1
beta[12] 6.1 0.0 0.2 5.7 6.0 6.1 6.3 6.6 2384 1
beta[13] 6.2 0.0 0.2 5.7 6.0 6.2 6.3 6.6 2385 1
beta[14] 6.7 0.0 0.2 6.2 6.5 6.7 6.9 7.2 2611 1
beta[15] 5.4 0.0 0.3 4.9 5.3 5.4 5.6 5.9 2452 1
beta[16] 5.9 0.0 0.2 5.4 5.8 5.9 6.1 6.4 2055 1
beta[17] 6.3 0.0 0.2 5.8 6.1 6.3 6.4 6.7 2314 1
beta[18] 5.8 0.0 0.2 5.4 5.7 5.8 6.0 6.3 2368 1
beta[19] 6.4 0.0 0.3 5.9 6.2 6.4 6.6 6.9 2431 1
beta[20] 6.1 0.0 0.2 5.6 5.9 6.1 6.2 6.5 2612 1
beta[21] 6.4 0.0 0.2 5.9 6.2 6.4 6.6 6.9 2163 1
beta[22] 5.9 0.0 0.2 5.4 5.7 5.9 6.0 6.4 2000 1
beta[23] 5.7 0.0 0.2 5.3 5.6 5.7 5.9 6.2 2699 1
beta[24] 5.9 0.0 0.2 5.4 5.7 5.9 6.1 6.4 3197 1
beta[25] 6.9 0.0 0.2 6.4 6.7 6.9 7.1 7.4 2084 1
beta[26] 6.5 0.0 0.2 6.1 6.4 6.5 6.7 7.1 2199 1
beta[27] 5.9 0.0 0.2 5.4 5.7 5.9 6.1 6.4 3171 1
beta[28] 5.9 0.0 0.2 5.4 5.7 5.9 6.0 6.3 2536 1
beta[29] 5.7 0.0 0.2 5.2 5.5 5.7 5.8 6.2 2137 1
beta[30] 6.1 0.0 0.2 5.6 6.0 6.1 6.3 6.6 2700 1
alpha_c 242.5 0.1 2.7 237.3 240.8 242.4 244.3 247.7 1826 1
beta_c 6.2 0.0 0.1 6.0 6.1 6.2 6.3 6.4 1860 1
tausq_c 37.1 0.2 5.7 27.6 33.0 36.6 40.5 50.1 955 1
tausq_alpha 218.4 2.1 71.2 122.8 169.9 206.5 251.2 391.7 1166 1
tausq_beta 0.3 0.0 0.1 0.1 0.2 0.3 0.3 0.5 1531 1
tau_c 6.1 0.0 0.5 5.3 5.7 6.0 6.4 7.1 954 1
tau_alpha 14.6 0.1 2.2 11.1 13.0 14.4 15.8 19.8 1389 1
tau_beta 0.5 0.0 0.1 0.4 0.4 0.5 0.6 0.7 1449 1
alpha0 106.3 0.1 3.6 98.9 104.0 106.3 108.6 113.4 1666 1
lp__ -438.0 0.3 7.3 -453.9 -442.7 -437.5 -432.9 -425.4 532 1
Samples were drawn using NUTS(diag_e) at Mon May 13 19:03:38 2019.
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).
从上面的结果来看，随机选择的初始值对结果没有影响，基于 NUTS 采样器实现的算法对初值不敏感，Stan 也具有非常高的采样效率。
Stan 代码尽可能地使用内置的向量化函数，可以极大地加快编译速度，运行速度，达到更好的优化性能。
# 使用随机生成的初始值
rats_vec <- stan(
# file = "code/bugs/vol1/rats/rats_vec.stan",
model_code = "
data {
int<lower=0> N;
int<lower=0> T;
real x[T];
real y[N,T];
real xbar;
}
transformed data {
real x_minus_xbar[T];
real y_linear[N*T];
for (t in 1:T)
x_minus_xbar[t] = x[t] - xbar;
for (n in 1:N)
for (t in 1:T)
y_linear[(n-1)*T + t] = y[n, t];
}
parameters {
real alpha[N];
real beta[N];
real alpha_c;
real beta_c;
real<lower=0> tausq_c;
real<lower=0> tausq_alpha;
real<lower=0> tausq_beta;
}
transformed parameters {
real<lower=0> tau_c;
real<lower=0> tau_alpha;
real<lower=0> tau_beta;
tau_c = sqrt(tausq_c);
tau_alpha = sqrt(tausq_alpha);
tau_beta = sqrt(tausq_beta);
}
model {
real pred[N*T];
for (n in 1:N)
for (t in 1:T)
pred[(n-1)*T + t] = fma(beta[n], x_minus_xbar[t], alpha[n]);
alpha_c ~ normal(0, 100);
beta_c ~ normal(0, 100);
tausq_c ~ inv_gamma(0.001, 0.001);
tausq_alpha ~ inv_gamma(0.001, 0.001);
tausq_beta ~ inv_gamma(0.001, 0.001);
alpha ~ normal(alpha_c, tau_alpha); // vectorized
beta ~ normal(beta_c, tau_beta); // vectorized
y_linear ~ normal(pred, tau_c); // vectorized
}",
data = list(N = N, T = T, y = y, x = x, xbar = xbar),
chains = chains, iter = iter, init = "random"
)
print(rats_vec, digits = 1)
Inference for Stan model: 83c0b6cb8d216e4723b1edc83099db35.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
alpha[1] 239.9 0.1 2.7 234.8 238.0 239.9 241.8 245.2 1278 1.0
alpha[2] 247.6 0.1 2.7 242.1 245.8 247.6 249.5 252.8 361 1.0
alpha[3] 252.2 0.3 3.1 244.3 250.4 252.4 254.2 258.0 124 1.0
alpha[4] 232.8 0.3 3.1 227.2 230.7 232.6 234.6 240.4 116 1.0
alpha[5] 232.0 0.3 3.1 226.6 230.0 231.7 233.7 240.7 95 1.0
alpha[6] 249.5 0.2 3.0 243.2 247.6 249.6 251.5 255.3 239 1.0
alpha[7] 229.2 0.4 3.4 223.7 227.1 228.9 230.9 239.8 64 1.0
alpha[8] 248.2 0.1 2.9 242.4 246.4 248.3 250.1 254.2 385 1.0
alpha[9] 282.0 1.3 6.9 246.0 281.1 283.1 284.8 288.1 27 1.1
alpha[10] 219.9 0.7 4.6 214.0 217.5 219.3 221.4 239.2 38 1.1
alpha[11] 257.8 0.5 3.7 245.2 256.1 258.1 260.1 263.5 63 1.1
alpha[12] 228.5 0.5 3.4 223.0 226.4 228.2 230.0 239.7 58 1.1
alpha[13] 242.4 0.1 2.8 236.9 240.5 242.4 244.3 247.6 1890 1.0
alpha[14] 267.6 0.9 4.9 246.1 266.1 268.2 270.1 273.6 32 1.1
alpha[15] 242.8 0.1 2.6 237.6 241.0 242.8 244.5 247.9 2598 1.0
alpha[16] 245.2 0.1 2.6 240.1 243.4 245.2 247.0 250.5 1877 1.0
alpha[17] 232.6 0.3 3.1 227.1 230.6 232.4 234.3 240.8 91 1.0
alpha[18] 240.4 0.1 2.8 234.9 238.6 240.4 242.3 245.8 1618 1.0
alpha[19] 253.6 0.3 3.3 245.2 251.8 253.8 255.6 259.3 118 1.0
alpha[20] 241.6 0.1 2.8 236.2 239.7 241.5 243.4 247.1 2047 1.0
alpha[21] 248.4 0.2 2.9 242.7 246.6 248.5 250.4 254.0 355 1.0
alpha[22] 225.7 0.5 3.8 220.0 223.6 225.4 227.1 240.3 60 1.1
alpha[23] 228.8 0.4 3.4 223.4 226.8 228.5 230.3 240.1 71 1.0
alpha[24] 245.1 0.1 2.7 239.7 243.3 245.1 246.9 250.4 1982 1.0
alpha[25] 234.7 0.3 3.1 229.1 232.6 234.5 236.5 241.7 143 1.0
alpha[26] 253.7 0.3 3.2 244.4 252.0 253.9 255.7 259.2 98 1.0
alpha[27] 254.0 0.3 3.3 244.3 252.5 254.2 256.0 259.7 106 1.0
alpha[28] 242.9 0.1 2.7 237.7 241.3 242.9 244.8 248.1 2227 1.0
alpha[29] 218.6 0.7 4.8 212.6 216.1 218.1 220.0 240.0 41 1.1
alpha[30] 241.5 0.1 2.7 236.0 239.6 241.5 243.3 246.8 1777 1.0
beta[1] 6.1 0.0 0.2 5.6 5.9 6.1 6.2 6.5 785 1.0
beta[2] 7.0 0.0 0.3 6.4 6.8 7.0 7.2 7.6 91 1.0
beta[3] 6.5 0.0 0.2 6.0 6.3 6.5 6.6 6.9 1648 1.0
beta[4] 5.4 0.0 0.3 4.8 5.2 5.4 5.5 6.2 73 1.0
beta[5] 6.6 0.0 0.2 6.1 6.4 6.6 6.7 7.0 1072 1.0
beta[6] 6.2 0.0 0.2 5.7 6.0 6.2 6.3 6.6 2051 1.0
beta[7] 6.0 0.0 0.2 5.5 5.8 6.0 6.1 6.5 466 1.0
beta[8] 6.4 0.0 0.2 6.0 6.2 6.4 6.6 6.9 1865 1.0
beta[9] 7.0 0.0 0.3 6.3 6.8 7.0 7.2 7.6 108 1.0
beta[10] 5.9 0.0 0.3 5.4 5.7 5.9 6.0 6.4 237 1.0
beta[11] 6.8 0.0 0.3 6.2 6.6 6.8 7.0 7.3 185 1.0
beta[12] 6.1 0.0 0.2 5.7 6.0 6.1 6.3 6.6 1571 1.0
beta[13] 6.2 0.0 0.2 5.7 6.0 6.2 6.3 6.7 2226 1.0
beta[14] 6.7 0.0 0.2 6.2 6.5 6.7 6.8 7.2 337 1.0
beta[15] 5.4 0.0 0.3 4.9 5.2 5.4 5.6 6.2 56 1.1
beta[16] 5.9 0.0 0.2 5.5 5.8 5.9 6.1 6.4 314 1.0
beta[17] 6.3 0.0 0.2 5.8 6.1 6.3 6.4 6.7 2943 1.0
beta[18] 5.9 0.0 0.2 5.4 5.7 5.9 6.0 6.4 236 1.0
beta[19] 6.4 0.0 0.2 5.9 6.2 6.4 6.6 6.9 2360 1.0
beta[20] 6.1 0.0 0.2 5.6 5.9 6.1 6.2 6.6 711 1.0
beta[21] 6.4 0.0 0.2 5.9 6.2 6.4 6.6 6.9 2454 1.0
beta[22] 5.9 0.0 0.3 5.4 5.7 5.9 6.0 6.4 253 1.0
beta[23] 5.8 0.0 0.3 5.3 5.6 5.8 5.9 6.3 183 1.0
beta[24] 5.9 0.0 0.3 5.4 5.7 5.9 6.1 6.4 318 1.0
beta[25] 6.9 0.0 0.3 6.3 6.7 6.9 7.1 7.4 133 1.0
beta[26] 6.5 0.0 0.2 6.1 6.4 6.5 6.7 7.0 1081 1.0
beta[27] 5.9 0.0 0.3 5.4 5.7 5.9 6.1 6.4 288 1.0
beta[28] 5.9 0.0 0.2 5.4 5.7 5.9 6.0 6.3 254 1.0
beta[29] 5.7 0.0 0.3 5.2 5.5 5.7 5.9 6.3 116 1.0
beta[30] 6.1 0.0 0.2 5.7 6.0 6.1 6.3 6.6 1706 1.0
alpha_c 242.4 0.1 2.6 237.0 240.8 242.4 244.1 247.7 2163 1.0
beta_c 6.2 0.0 0.1 6.0 6.1 6.2 6.3 6.4 555 1.0
tausq_c 43.8 7.6 37.2 27.7 33.5 37.1 41.6 198.1 24 1.1
tausq_alpha 209.0 6.9 68.9 5.5 170.5 203.1 245.2 353.5 100 1.0
tausq_beta 0.3 0.0 0.1 0.0 0.2 0.3 0.3 0.5 111 1.0
tau_c 6.4 0.3 1.7 5.3 5.8 6.1 6.4 14.1 25 1.1
tau_alpha 14.2 0.4 2.9 2.3 13.1 14.3 15.7 18.8 43 1.1
tau_beta 0.5 0.0 0.1 0.0 0.4 0.5 0.6 0.7 47 1.1
lp__ -437.9 0.6 8.0 -453.3 -442.5 -437.8 -433.5 -424.1 178 1.0
Samples were drawn using NUTS(diag_e) at Mon May 13 19:38:13 2019.
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).
毫无疑问，向量化只是为了提高计算效率，是编程技巧，不会显著影响结果的准确度。