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).