pystan_eight_schools_bayesian_hierarchical_model_sampling.py

1import stan
2
3schools_code = """
4data {
5  int<lower=0> J;         // number of schools
6  real y[J];              // estimated treatment effects
7  real<lower=0> sigma[J]; // standard error of effect estimates
8}
9parameters {
10  real mu;                // population treatment effect
11  real<lower=0> tau;      // sd of treatment effects
12  vector[J] eta;          // unscaled deviation from mu by school
13}
14transformed parameters {
15  vector[J] theta = mu + tau * eta;        // school treatment effects
16}
17model {
18  eta ~ std_normal();       // prior log-density
19  y ~ normal(theta, sigma); // log-likelihood
20}
21"""
22
23schools_data = {"J": 8,
24                "y": [28,  8, -3,  7, -1,  1, 18, 12],
25                "sigma": [15, 10, 16, 11,  9, 11, 10, 18]}
26
27posterior = stan.build(schools_code, data=schools_data)
28fit = posterior.sample(num_chains=4, num_samples=1000)
29eta = fit["eta"]  # array with shape (8, 4000)
30df = fit.to_frame()  # pandas `DataFrame`
31print(df)