/**********************************************************************
 * hierarchical normal model 
 *
 * math[i, j] ~ normal(theta[i], sigma)
 * theta[i] ~ normal(mu, tau)
 * sigma^2 ~ scaled_inv_chi_square(1, 10)
 * mu ~ normal(50, 5)
 * tau^2 ~ scaled_inv_chi_square(1, 10)
 **********************************************************************/

data {
    int<lower=0> n;               // number of observations 
    int<lower=0> m;               // number of groups 
    real math[n];                 // y 
    int<lower=0> school[n];       // school ID, 1, ..., m 
}
parameters {
    real theta[m];                // school mean 
    real mu; 
    real<lower=0> tau; 
    real<lower=0> sigma;
}
model {
    sigma^2 ~ scaled_inv_chi_square(1, 10);
    tau^2 ~ scaled_inv_chi_square(1, 10); 
    mu ~ normal(50, 5); 

    theta ~ normal(mu, tau);

    for(i in 1:n) {
        math[i] ~ normal(theta[school[i]], sigma); 
    }
}

/**********************************************************************
 * THE END 
 **********************************************************************/