library(rstan) 
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)

print(fit)
plot(fit)