Hi,

Let’s say that I am trying to fit a multivariate normal distribution of a vector of parameters named `Vec`

with all means being zero and the covariance matrix being `Sigma`

(other parts of the code have been left out of the example). Usually I would do something like this:

```
data {
int<lower = 1> K; // number of variables
vector[K] Zero; // vector of zeros for the parameter means
}
parameters {
row_vector[K] Vec;
corr_matrix(K) corr; // correlations of parameters
vector<lower=0> sigma; // variances of parameters
}
transformed parameters{
cov_matrix[K] Sigma;
Sigma = quad_form_diag(corr, sigma);
}
model{
Vec ~ multi_normal(Zero, Sigma);
}
```

However, I would like to switch to using Cholesky factors because it makes the estimation more stable and I am currently having problems with model convergence and the culprit seems to be this covariance matrix. But, I don’t understand how exactly should the code look like. The documentation says to put ` cholesky_factor_corr[K] L;`

in the parameters block and `L ~`

`lkj_corr_cholesky`

`(eta)`

in the model block, but I don’t understand if this will give me a full matrix, with the diagonals included, or should I also use `quad_form_diag`

to join the variances with the covariances? Furthermore, should the resulting matrix be used in the `multi_normal`

to specify the covariance given that it’s a correlation matrix and not the covariance matrix? I would appreciate any help on how my code should look like.