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

Yi

I added a stan tag for you. – Celebrity

@DirkEddelbuettel Thanks a lot! I proposed a wiki :) probably Rstan could be a choice as well - for the r-implementation – Yi

I think this question is better suited for crossvalidated.com, as it doesn't concern a programming technical question. – Fiedling

@SachaEpskamp well actually i thought it does since the problem how to tell the programm to include more percentiles in the output ;) – Yi

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

Yi

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.

Bridie

The thing about apply(fit, MARGIN = "parameters", FUN = quantile, probs = (1:100) / 100)) is that it will include the stuff in the generated quantities{} block as well. So extract might be better. – Kerplunk