Hey all,

I understand how to **sample** from multivariate densities, using something like: M + SLZ, where Z is some standard random variate from the distribution(s), L is a cholesky correlation matrix, and S is a diagonal scale matrix. One can use this to sample parameters from various combinations of multivariate densities, or multivariate densities not yet included in Stan.

However, I’m curious whether one could use a similar approach to modeling multivariate likelihoods. Assume for the moment that the multivariate normal doesn’t exist in Stan, but we’re predicting multiple normally distributed outcomes. We want to assume that residual covariance exists, and is accounted for. Is it possible to do something like:

M = XB (M is an NxJ matrix of predicted values for N observations and J variables; B is a coefficient matrix; X is a design matrix).

Mhat = M + SLZ; where Z is a standard normal variate matrix; L is a cholesky correlation matrix; S is a scale matrix.

y[,1] ~normal(Mhat[,1], S[1,1]); ?

In other words, can you use the same trick for simulating from multivariate densities to create “multivariate likelihoods” that are really univariate likelihoods that account for covariance?

My future usecase would be to use multivariate likelihoods for things like the skew-normal, skew-t, etc; sampling from those multivariate forms is easy, but the multivariate density is terrible, so if I could model the covariance of the residuals AND still use the defined univariate skew likelihoods, that would be fantastic; I just have no idea if that’s possible, or if I should suck it up and code up the multivariate density.