How to obtain the full marginal distribution of a parameter in stan

Question

when starting a standard example from the stan webpage like the following:

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)

(this has been obtained from here)

however this does only provide me with the quantiles of the posterior of the parameters. so my question is: how to obtain other percentiles? i guess it should be similar to bugs(?)

remark: i tried to introduce the tag stan however, i have too little reputation ;) sorry for that

Answer 1

here's my attempt hope this is correct:

suppose fit is an object obtained from stan(...). then the posterior for any percentile is obtained from:

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

where the number in square brackets is the chain i guess. in case this hasn't been clear: i use rstan

Answer 2

As from rstan v1.0.3 (not released yet), you will be able to utilize the workhorse apply() function directly on an object of stanfit class that is produced by the stan() function. If fit is an object obtained from stan(), then for example,

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

or

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

The former applies FUN to each parameter in each chain, while the latter combines the chains before applying FUN to each parameter. If you were only interested in one parameter, then something like

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

is an option.