// student_grade_model_ppc.stan

data {
    int<lower=0> n; // number of students
    int<lower=0> p; // number of tests
    array[n, p] int<lower=0, upper=100> X; // student test grades
    real<lower=0> tau;
    real<lower=0> a;
    real<lower=0> b;
}
parameters {
    array[n] real z;
    real mu;
    real<lower=0> sigma_sq;
}
transformed parameters {
    array[n] real<lower=0, upper=1> theta;
    real sigma;
    theta = inv_logit(z);
    sigma = sqrt(sigma_sq);
}
model {
    sigma_sq ~ inv_gamma(a,b);
    mu ~ normal(0, sigma * tau);
    z ~ normal(mu, sigma);
    for (i in 1:n)
      X[i] ~ binomial(100,theta[i]);
}
generated quantities{
    array[n, p] int<lower=0, upper=100> X_rep;
    real log_lhd = 0;
    real log_lhd_rep = 0;
    real ppp;
    for (i in 1:n){
      for (j in 1:p){
	log_lhd += binomial_lpmf(X[i][j] | 100,theta[i]);
        X_rep[i][j] = binomial_rng(100,theta[i]);
	log_lhd_rep += binomial_lpmf(X_rep[i][j] | 100,theta[i]);
      }
    }
    ppp = log_lhd >= log_lhd_rep ? 1 : 0;
}