Sum of cholesky factors vs sum of covariance

Modeling
3

Thanks @RJTK, that was exactly my problem. I didn’t even stop to think that the Cholesky factored covariance matrix is basically the square root of the variance-covariance matrix! Parameterizing on covariance fixes the problem, although you then lose the efficiency of the Cholesky parameterization:

parameters {
  corr_matrix[J] Omega;
  corr_matrix[J] Omega_cond;
  vector<lower=0>[J] sigma;
  vector<lower=0>[J] sigma_cond;
  ...
}
model {
  sigma ~ cauchy(0, 5);
  sigma_cond ~ cauchy(0, 5);
  Omega ~ lkj_corr(1);
  Omega_cond ~ lkj_corr(1);
  ...
  for (i in 1:N) {
      y[i] ~ multi_normal(beta + cond[i] * beta_cond, Sigma + cond[i] * Sigma_cond);
  }
}

It is possible to continue using the Cholesky parameterization with the code below, but you lose most of the efficiency gain. I am not sure if there is a more efficient way to sum the Cholesky factors…

functions {
  matrix sum_chol_factors(matrix chol1, matrix chol2) {
    // Re-factoring the covariance matrix doesn't seem efficient...
    return cholesky_decompose(chol1 * chol1' + chol2 * chol2');
  }
}
...
model {
  ...
  for (i in 1:N) {
    if (cond[i])
      y[i] ~ multi_normal_cholesky(
        beta + beta_cond,
        sum_chol_factors(L_Sigma, L_Sigma_cond));
    else
      y[i] ~ multi_normal_cholesky(beta, L_Sigma);
  } 
}