/** * Convert a variance-covariance matrix to a correlation matrix. * * @param cov A variance-covariance matrix. * @return The corresponding correlation matrix. */ functions { matrix cov2cor(matrix cov) { // Compute inverse of standard deviations of the covariance matrix. matrix[rows(cov), cols(cov)] stdinv = diag_matrix(1.0 ./ sqrt(diagonal(cov))); // Compute the correlation matrix. return stdinv * cov * stdinv; } } data { int N; // number of observations int J; // dimension of observations // Observations of either one bivariate random normal variable (the base) or // the some of the base and a second random normal variable (the condition). vector[J] y[N]; // Binary condition. // Zero when the observation is of the base only. // One when we are observing the sum of the base + cond variables. int cond[N]; // binary condition } parameters { // Error scales cholesky_factor_corr[J] L_corr_base; cholesky_factor_corr[J] L_corr_sum; vector[J] sigma_base; vector[J] sigma_sum; // Means vector[J] beta_base; vector[J] beta_sum; } model { // Cholesky factor variance-covariance matrices. matrix[J, J] L_Sigma_base; matrix[J, J] L_Sigma_sum; L_Sigma_base = diag_pre_multiply(sigma_base, L_corr_base); L_Sigma_sum = diag_pre_multiply(sigma_sum, L_corr_sum); // Non-informative priors. sigma_base ~ cauchy(0, 5); sigma_sum ~ cauchy(0, 5); L_corr_base ~ lkj_corr_cholesky(1); L_corr_sum ~ lkj_corr_cholesky(1); beta_base ~ normal(0,10); beta_sum ~ normal(0,10); // Sampling distribution of the observations. for (i in 1:N) { if (cond[i]) y[i] ~ multi_normal_cholesky(beta_sum, L_Sigma_sum); else y[i] ~ multi_normal_cholesky(beta_base, L_Sigma_base); } } generated quantities { // Convert cholesky factors back to correlation and variance-covariance. matrix[J, J] Omega_base = multiply_lower_tri_self_transpose(L_corr_base); matrix[J, J] Omega_sum = multiply_lower_tri_self_transpose(L_corr_sum); matrix[J, J] Sigma_base = quad_form_diag(Omega_base, sigma_base); matrix[J, J] Sigma_sum = quad_form_diag(Omega_sum, sigma_sum); // Infer correlation matrix of unobserved, condition-only variable. matrix[J, J] Sigma_cond = Sigma_sum - Sigma_base; matrix[J, J] Omega_cond = cov2cor(Sigma_cond); }