data { int N; // # of samples vector [N] Y; // weight of sample int NLOTS; // # of lots int lot [N]; // lot of sample real meanProcessMean; real tauProcessMean; real dfsigmaA; real meansigmaA; } transformed data { real sigmaProcessMean; real alpha; real beta; sigmaProcessMean = 1 / sqrt (tauProcessMean); alpha = dfsigmaA / 2; beta = alpha * meansigmaA ^ 2; } parameters { vector [NLOTS] mu_raw; real processMean; real tauA; real sigmaB_unif; real sigmaW_unif; vector [N] Y_true_raw; vector [NLOTS] delta; vector [NLOTS] nu; vector [N] om; vector [N] gamma_raw; } transformed parameters { vector [N] Y_true; vector [NLOTS] mu; real sigmaA; real sigmaB; real sigmaW; vector [N] mu_b; vector [N] tau; vector [N] abs_gamma_raw; for (n in 1 : N) abs_gamma_raw [n] = fabs (gamma_raw [n]); for (n in 1 : N) tau [n] = 1 / sqrt (om [n]); sigmaA = 1 / sqrt (tauA); sigmaB = 2.5 * tan (sigmaB_unif); // sigmaB ~ cauchy(0, 2.5) sigmaW = 2.5 * tan (sigmaW_unif); // sigmaW ~ cauchy(0, 2.5) // assume normal inter-lot variation mu = processMean + sigmaB * mu_raw; // assume skew-t intra-lot variation mu_b = mu [lot] + delta [lot] .* abs_gamma_raw .* tau; Y_true = mu_b + sigmaW * tau .* Y_true_raw; } model { tauA ~ gamma (alpha, beta); processMean ~ normal (meanProcessMean, sigmaProcessMean); mu_raw ~ normal (0, 1); Y_true_raw ~ normal (0, 1); Y ~ normal (Y_true, sigmaA); om ~ gamma (nu [lot] / 2, nu [lot] / 2); gamma_raw ~ normal (0, 1); delta ~ normal (0, 5); nu ~ exponential (0.1); } generated quantities { vector[N] log_lik; vector[N] Y_rep; for (n in 1 : N) { real tmp; log_lik [n] = normal_lpdf (Y [n] | Y_true [n], sigmaA); tmp = normal_rng (mu_b [n], sigmaW); Y_rep [n] = normal_rng (tmp, sigmaA); } }