data { int nmetro; int nx; vector[nmetro] index; vector[nmetro] index_se; matrix[nmetro, nx] x; matrix[nmetro, nx] x_se; } transformed data { matrix[nmetro, nx] x_var; vector[nmetro] index_var; x_var = x_se .* x_se; index_var = index_se .* index_se; } parameters { real alpha; vector[nx] beta; vector[nx] mu; vector[nx] sigma2; real tau2; } model { for(i in 1:nx){ x[,i] ~ normal(mu[i], sqrt(sigma2[i] + x_var[,i])); } for(i in 1:nmetro){ real mubar; real sig2bar; vector[nx] lambda; lambda = sigma2 ./ (sigma2 + to_vector(x_var[i]) ); mubar = alpha + dot_product(mu + lambda .* (to_vector(x[i]) - mu), beta); sig2bar = index_var[i] + tau2 + dot_product(sigma2 .* (1 - lambda), beta .* beta); index ~ normal(mubar, sqrt(sig2bar)); } }