Multivariate normal with diagonal covariance matrix

In my Stan time series model, the s-dimensional vector-valued latent variable x[k] has a multivariate normal distribution with mean vector returned by the function getEXdiag(), and a covariance matrix that is a scalar multiple of the identity matrix. The function getVXdiag() returns the required scalar. I’ve skipped the details of the arguments to these functions. Both the following sampling statements compile. I had expected them to give the same answer. On simulated data, the multi_normal version leads to sensible parameter estimates, while the normal version leads to estimates that are completely wrong.

x[k] ~ multi_normal(getEXdiag(eBT, x[k-1], b1, a[j]), diag_matrix(rep_vector(1.0, s)) * getVXdiag(eBT, b1, sigma2));

x[k] ~ normal(getEXdiag(eBT, x[k - 1], b1, a[j]), getVXdiag(eBT, b1, sigma2));

Did I just do something stupid with the syntax in the normal version? I’d like to use it if I could, because it’s much faster (although I expect that I could speed up the multi_normal version a bit by not first making an identity matrix and then scaling it).

Seems the same to me. Can you post the rest of your model here? I’m curious about the exact types floating around. Might be something else going on (or maybe one of these types is funky).

Just do

real sigma = sqrt(sigma2); // I suspect you forgot that normal accepts a standard deviation
for (k in 1:K)
  x[k] ~ normal(getEXdiag(eBT, x[k - 1], b1, a[j]), sigma);
1 Like

Thank you! I had forgotten that normal wants the standard deviation, not the variance. Fixed that and now the univariate version is working nicely, and much faster than the multivariate version.