data{ int NvarsX; // num independent variables int NvarsY; // num dependent variables int N; // Total number of rows vector[NvarsX] x [N]; // data for independent vars vector[NvarsY] y [N]; // data for dependent vars } parameters{ cholesky_factor_corr[NvarsY] L_Omega; // L form of correlation matrix for Yvars vector[NvarsY] L_sigma; // Vector of scales for covariance matrix for Yvars vector[NvarsY] alpha; // Intercepts for each Y matrix[NvarsY, NvarsX] beta ; // Coefficients for each X vs each Y } transformed parameters{ vector[NvarsY] mu[N]; matrix[NvarsY, NvarsY] L_Sigma; L_Sigma = diag_pre_multiply(L_sigma, L_Omega); // L form of Sigma covariance matrix for Yvars for (n in 1:N){ mu[n] = alpha + beta * x[n]; } } model{ //Priors L_Omega ~ lkj_corr_cholesky(1); L_sigma ~ normal(0, 5); alpha ~ normal(0, 3); to_vector(beta) ~ normal(0, 2); //likelihood y ~ multi_normal_cholesky(mu, L_Sigma); } generated quantities{ matrix[NvarsY, NvarsY] Omega; // Square form of correlation matrix matrix[NvarsY, NvarsY] Sigma; // Square form of covariance matrix vector[NvarsY] mu_ppc[N]; // model predictions for posterior predictive check vector[NvarsY] y_ppc [N]; // model predictions for posterior predictive check // Generate correlation matrix Omega = multiply_lower_tri_self_transpose(L_Omega); Sigma = quad_form_diag(Omega, L_sigma); // Generate y_ppc for (n in 1:N){ mu_ppc[n] = alpha + beta * x[n]; } y_ppc = multi_normal_cholesky_rng(mu_ppc, L_Sigma); }