5

从stan 网页启动标准示例时，如下所示：

schools_code <- '
  data {
   int<lower=0> J; // number of schools 
   real y[J]; // estimated treatment effects
   real<lower=0> sigma[J]; // s.e. of effect estimates 
 } 
 parameters {
   real theta[J]; 
   real mu; 
   real<lower=0> tau; 
 } 
 model {
   theta ~ normal(mu, tau); 
   y ~ normal(theta, sigma);
 } 
 '

  schools_dat <- list(J = 8, 
                 y = c(28,  8, -3,  7, -1,  1, 18, 12),
                 sigma = c(15, 10, 16, 11,  9, 11, 10, 18))

 fit <- stan(model_code = schools_code, data = schools_dat, 
           iter = 1000, n_chains = 4)

（这是从这里获得的）

然而，这只为我提供了参数后验的分位数。所以我的问题是：如何获得其他百分位数？我想它应该类似于错误（？）

备注：我试图介绍标签stan，但是我的名声太小了；）抱歉

4

2 回答 2

9

从rstan v1.0.3（尚未发布）开始，您将能够apply()直接在stanfit class由stan() function. 如果 fit 是从 获得的对象stan()，那么例如，

apply(fit, MARGIN = "parameters", FUN = quantile, probs = (1:100) / 100)

或者

apply(as.matrix(fit), MARGIN = 2, FUN = quantile, probs = (1:100) / 100)

前者将 FUN 应用于每个链中的每个参数，而后者在将 FUN 应用于每个参数之前将链组合起来。如果您只对一个参数感兴趣，那么类似

beta <- extract(fit, pars = "beta", inc_warmup = FALSE, permuted = TRUE)[[1]]
quantile(beta, probs = (1:100) / 100)

是一种选择。

于 2012-10-21T10:55:21.060 回答
0

这是我的尝试，希望这是正确的：

假设fit是从 获得的对象stan(...)。那么任何百分位数的后验来自：

quantile(fit@sim$sample[[1]]$beta, probs=c((1:100)/100))

方括号中的数字是我猜的链。万一这还不清楚：我用rstan

于 2012-09-04T19:41:07.423 回答