PyMC3 passing stochastic covariance matrix to pm.MvNormal()

Question

I've tried to fit a simple 2D gaussian model to observed data by using PyMC3.

import numpy as np
import pymc3 as pm

n = 10000;
np.random.seed(0)
X = np.random.multivariate_normal([0,0], [[1,0],[0,1]], n);

with pm.Model() as model:
    # PRIORS
    mu = [pm.Uniform('mux', lower=-1, upper=1), 
          pm.Uniform('muy', lower=-1, upper=1)]
    cov = np.array([[pm.Uniform('a11', lower=0.1, upper=2), 0],
                    [0, pm.Uniform('a22', lower=0.1, upper=2)]])

    # LIKELIHOOD
    likelihood = pm.MvNormal('likelihood', mu=mu, cov=cov, observed=X)

with model:
    trace = pm.sample(draws=1000, chains=2, tune=1000)

while I can do this in 1D by passing the sd to pm.Normal, I have some trouble in passing the covariance matrix to pm.MvNormal.

Where am I going wrong?

Answer 1

PyMC3 distribution objects are not simple numeric objects or numpy arrays. Instead, they are nodes in a theano computation graph and often require operations from either pymc3.math or theano.tensor to manipulate them. Moreover, placing PyMC3 objects in numpy arrays is unnecessary since they already are multidimensional.

Original Model

Keeping with the intent of your code, a working version would go something like

import numpy as np
import pymc3 as pm
import theano.tensor as tt

N = 10000
np.random.seed(0)
X = np.random.multivariate_normal(np.zeros(2), np.eye(2), size=N)

with pm.Model() as model:
    # use `shape` argument to define tensor dimensions
    mu = pm.Uniform('mu', lower=-1, upper=1, shape=2)

    # diagonal values on covariance matrix
    a = pm.Uniform('a', lower=0.1, upper=2, shape=2)

    # convert vector to a 2x2 matrix with `a` on the diagonal
    cov = tt.diag(a)

    likelihood = pm.MvNormal('likelihood', mu=mu, cov=cov, observed=X)

Alternative Model

I assume the example you provided is just a toy to communicate the problem. But just in case, I'll mention that the main advantage of using a multivariate normal (modeling covariance between parameters) is lost when restricting the covariance matrix to be diagonal. Furthermore, the theory of priors for covariance matrices is well-developed, so it's worth one's time considering existing solutions. In particular, there is a PyMC3 example using the LKJ prior for covariance matrices.

Here's a simple application of that example in this context:

with pm.Model() as model_lkj:
    # use `shape` argument to define tensor dimensions
    mu = pm.Uniform('mu', lower=-1, upper=1, shape=2)

    # LKJ prior for covariance matrix (see example)
    packed_L = pm.LKJCholeskyCov('packed_L', n=2,
                                 eta=2., sd_dist=pm.HalfCauchy.dist(2.5))
    # convert to (2,2)
    L = pm.expand_packed_triangular(2, packed_L)

    likelihood = pm.MvNormal('likelihood', mu=mu, chol=L, observed=X)