Variance-covariance matrix multivariate normal

General
1

I’m trying to model using a covariance matrix but I’m not sure whether i should be using the variance-covariance matrix or the square root of the var-cov matrix( specially, the square root of the Covariance will be used to define the matrix). In the univariate case for a normal distribution, I know that we use standard deviations, but just want to be absolutely clear that it is the standard deviation of the variance covariance matrix that is being used for the multivariate case. Does this also mean when decomposing the covarinace matrix into a scale and a correlation matrix, the specification of the correalation matrix would be different, than if we had specified the variance-covariance matrix ?

2

There are multiple ways to define the square root of a variance-covariance matrix. The easiest one to use in Stan is the Cholesky factor, which can (and often should) be decomposed further into a vector of standard deviations and a Cholesky factor of a correlation matrix, which can be multiplied together with diag_pre_multiply. In the usual case that it is being used with the multivariate normal distribution, you can pass that to the last argument of multi_normal_cholesky_lpdf if you are modeling data. If the multivariate normal is used as a prior, it is often best to just use the stochastic representation and pre-multiply a vector of standard normal variates by the Cholesky factor and shift by the mean vector (if it is not all zeros).

3

I’ve actually used Omega = quad_form_diag(rho,omega), where rho is the correlation matrix and omega is a vector of standard deviations. But this will return Omega to be a matrix with standard deviations on the diagonal and covariance on the off-diagonals. Then passed this in as a multivariate normal prior ~ multi_normal(beta,Omega). Can it be done this way? Or is the Omega matrix essentially wrong- should the diagonals be variances?

4

If you do Omega = quad_form_diag(rho,omega) and omega is a correlation matrix, then the diagonal elements of Omega will be rho squared, i.e. variances. That is valid, but then multi_normal has to do a Cholesky decomposition of Omega internally, so you can avoid that by calling multi_normal_cholesky with diag_pre_multiply(rho, L) as the last argument. But better still do

parameters {
  vector[K] z;
  cholesky_factor_corr[K] L;
  ...
}
transformed parameters {
  vector[K] prior = beta + diag_pre_multiply(rho, L) * z;
}
model {
  target += normal_lpdf(z | 0, 1); // implies prior is MVN
  ...
}
5

Brilliant, thanks for the response!

1 Like