Multivariate normal with changing covariance over observations


I would like to construct a model where the observed data consists of a set of samples - assumed to be multivariate Normal, but with a different covariance matrix for each sample, because these cov matrices are functions of other variables in the model.

As far as I can tell, if I use MvNormal - then mu values can be specified for each observed point, but cov values cannot - is that correct ? (if not, how do you set up cov for this scenario?)

As an alternative, I have tried writing a ‘logp’ function which is used in a pm.Potential - this function calls the logp() of pm.MvNormal.dist in a loop, passing in the appropriate mu and cov for every point - but this seems to be incredibly slow and memory intensive.

Is there a better or recommended way of doing something like this ?

Since the current target problem is 2-dimensional, if the issue is down to the matrix operations in the density - then I could code the analytic solutions directly in this case - but I don’t want to end up ‘reinventing the wheel’ if there is already a good way of doing this somewhere in the package.

I happen to be using PyMC3 at the moment, but I’m guessing that isn’t very relevant.