data { int C; int R; int S[R]; int RS; int V; int Y[RS]; matrix[C, V] X[RS]; vector[V] prior_mean; vector[V] prior_sd; } parameters { vector[V] sigma; row_vector[V] theta; cholesky_factor_corr[V] L_omega; matrix[V, R] standard_normal; } transformed parameters { matrix[V, V] L_sigma; matrix[R, V] beta; L_sigma = diag_pre_multiply(sigma, L_omega); beta = rep_matrix(theta, R) + (L_sigma * standard_normal)'; } model { int rs = 1; vector[V * R] vector_normal; sigma ~ gamma(1, 1); theta ~ normal(prior_mean, prior_sd); L_omega ~ lkj_corr_cholesky(4); vector_normal = to_vector(standard_normal); vector_normal ~ normal(0, 1); for (r in 1:R) { for (s in 1:S[r]) { Y[rs] ~ categorical_logit(X[rs] * to_vector(beta[r])); rs += 1; } } }