In [2]:

    

library(rstan)


    

    




Loading required package: Rcpp
Loading required package: inline

Attaching package: ‘inline’

The following object is masked from ‘package:Rcpp’:

    registerPlugin

rstan (Version 2.7.0-1, packaged: 2015-07-17 18:12:01 UTC, GitRev: 05c3d0058b6a)
For execution on a local, multicore CPU with excess RAM we recommend calling
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

In [3]:

    

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 mu; 
  real<lower=0> tau;
  real eta[J];
}
transformed parameters {
  real theta[J];
  for (j in 1:J)
    theta[j] <- mu + tau * eta[j];
}
model {
  eta ~ normal(0, 1);
  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, chains = 4)


    

    




COMPILING THE C++ CODE FOR MODEL 'schools_code' NOW.

SAMPLING FOR MODEL 'schools_code' NOW (CHAIN 1).

Chain 1, Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 1, Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 1, Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 1, Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 1, Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 1, Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 1, Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 1, Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 1, Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 1, Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 1, Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 1, Iteration: 1000 / 1000 [100%]  (Sampling)
#  Elapsed Time: 0.015087 seconds (Warm-up)
#                0.015038 seconds (Sampling)
#                0.030125 seconds (Total)


SAMPLING FOR MODEL 'schools_code' NOW (CHAIN 2).

Chain 2, Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 2, Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 2, Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 2, Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 2, Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 2, Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 2, Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 2, Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 2, Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 2, Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 2, Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 2, Iteration: 1000 / 1000 [100%]  (Sampling)
#  Elapsed Time: 0.012598 seconds (Warm-up)
#                0.01365 seconds (Sampling)
#                0.026248 seconds (Total)


SAMPLING FOR MODEL 'schools_code' NOW (CHAIN 3).

Chain 3, Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 3, Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 3, Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 3, Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 3, Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 3, Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 3, Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 3, Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 3, Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 3, Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 3, Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 3, Iteration: 1000 / 1000 [100%]  (Sampling)
#  Elapsed Time: 0.012619 seconds (Warm-up)
#                0.014506 seconds (Sampling)
#                0.027125 seconds (Total)


SAMPLING FOR MODEL 'schools_code' NOW (CHAIN 4).

Chain 4, Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 4, Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 4, Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 4, Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 4, Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 4, Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 4, Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 4, Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 4, Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 4, Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 4, Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 4, Iteration: 1000 / 1000 [100%]  (Sampling)
#  Elapsed Time: 0.014785 seconds (Warm-up)
#                0.014077 seconds (Sampling)
#                0.028862 seconds (Total)

In [4]:

    

print(fit)


    

    




Inference for Stan model: schools_code.
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
mu        7.89    0.16 4.97  -2.88  4.85  7.99 11.23 16.98   957    1
tau       6.75    0.21 5.63   0.22  2.49  5.39  9.58 22.53   752    1
eta[1]    0.39    0.02 0.91  -1.46 -0.19  0.40  0.99  2.20  1539    1
eta[2]    0.01    0.02 0.85  -1.69 -0.51  0.02  0.55  1.76  1652    1
eta[3]   -0.21    0.03 0.97  -2.07 -0.87 -0.22  0.43  1.80  1328    1
eta[4]   -0.02    0.02 0.89  -1.79 -0.60 -0.02  0.54  1.70  1314    1
eta[5]   -0.37    0.02 0.86  -1.97 -0.94 -0.40  0.18  1.45  1309    1
eta[6]   -0.20    0.02 0.89  -1.86 -0.83 -0.21  0.36  1.58  1439    1
eta[7]    0.32    0.02 0.89  -1.43 -0.26  0.34  0.93  2.02  1618    1
eta[8]    0.03    0.02 0.92  -1.78 -0.56  0.03  0.62  1.80  1661    1
theta[1] 11.54    0.25 8.49  -1.86  6.14 10.46 15.26 32.43  1129    1
theta[2]  7.89    0.14 6.23  -4.37  4.02  7.86 11.68 20.20  2000    1
theta[3]  5.97    0.22 8.04 -12.42  1.93  6.72 10.74 20.74  1309    1
theta[4]  7.75    0.16 6.39  -5.08  3.85  7.79 11.60 20.54  1530    1
theta[5]  5.05    0.15 6.16  -8.29  1.24  5.41  9.27 16.42  1601    1
theta[6]  6.17    0.17 6.86  -8.84  2.04  6.44 10.65 18.93  1567    1
theta[7] 10.68    0.18 6.74  -1.26  6.36 10.06 14.56 25.31  1419    1
theta[8]  8.22    0.21 8.08  -7.49  3.52  8.13 12.36 25.72  1485    1
lp__     -4.79    0.12 2.69 -10.77 -6.40 -4.52 -3.01 -0.16   520    1

Samples were drawn using NUTS(diag_e) at Sat Aug  1 18:27:17 2015.
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).

In [6]:

    

plot(fit)


    

In [ ]: