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.