I’m fitting for the parameters of a multi-variate (3D) Gaussian: the 3 means, the 3 standard deviations, and only one correlation coefficient for the covariance matrix (I’m allowing a covariance between only two variables, the other off-diagonal terms of the covariance matrix are assumed to be zero). I am new to PyMC. I have a working version based on scipy+numpy where I’ve used the pytensor @as_op decorator on a forward_model() function that takes in the vector of means and the covariance matrix, and then spits out the probability for a bunch of points according to the current Gaussian. This works but it’s very slow to sample and PyMC chooses the Slice sampler by default – I think because I’m using things like scipy.stats.multivariate_normal and np.tensordot inside my forward_model() function, instead of using pymc.MvNormal and pytensor.tensor operations (which I think would allow for NUTS sampling).

I have a very basic question – how do I convert the following 2 lines of scipy code to instead use pm.MvNormal and pytensor.tensor.dotproduct so that I can put it all directly into my “with pm.Model() as model” context without needing a separate decorated forward_model function?

`mvn = scipy.stats.multivariate_normal(mean=vector_means,cov=matrix_covariances)`

`probs = mvn.pdf(points)`

Here points is a list of size-3 [x,y,z] pairs at which we will evaluate the probability given the prescribed 3D Gaussian parameters vector_means and matrix_covariances (the latter are concatenated from pm.Uniform Stochastic Random Variable objects).

I tried to do

`mvn = pm.MvNormal('mvn',mu=vector_means,cov=matrix_covariances)`

`probs = pm.Deterministic('probs', pt.exp(pm.logp(mvn,points)))`

but while pm.sample does work, it is treating mvn as another random variable to sample (whereas it should only be sampling the 3 mean, 3 standard deviation and 1 correlation coefficient parameters).

So I guess I want mvn to be some kind of Deterministic (or Potential?) object which is fully determined by the 7 stochastic RVs in its inputs vector_means and matrix_covariances. Finally I tried

`mvn = pm.Deterministic('mvn',pm.MvNormal('mvn_inside',mu=vector_means,cov=matrix_covariances))`

`probs = pm.Deterministic('probs', pt.exp(pm.logp(mvn,points)))`

but then I get an error for the pm.logp call:

`TypeError: When RV is not a pure distribution, value variable must have the same type`

Any suggestions? Is there a standard way of evaluating MvNormal as part of a PyMC model?