How to build a multi-covariance model in pymc3?

Questions
1

Hey community,

I am trying to study the covariance between several variables coming from two(or more) different engines. Variables are temperature, vibration, humidity, etc. I built a simple covariance model using 2 variables of one engine using this tutorial ([LKJ](https://docs.pymc.io/notebooks/LKJ.html))

Currently, I want to use 2 engines and study the independent covariance of the variables. But, I’m stuck with not understanding the shape parameter. The error is given at the bottom. The idea is to eventually build a hierarchical model. But, I am not sure if my strategy is correct.

with pm.Model() as model:    
      p = pm.Dirichlet('p', a = np.ones(2), shape = (2, 2))
      category = pm.Categorical('category', p=p, shape= 2) 

     # Covariance matrix model for engine 1
     # LKJ prior for covariance matrix 
     packed_L1 = pm.LKJCholeskyCov('packed_L1', n=2, eta=2., sd_dist=pm.HalfCauchy.dist(6))    
     L1 = pm.expand_packed_triangular(2, packed_L1)    
     cov1 = pm.Deterministic('cov1', L1.dot(L1.T))


     # Covariance matrix model for engine 2
     # LKJ prior for covariance matrix 
     packed_L2 = pm.LKJCholeskyCov('packed_L2', n=2, eta=2., sd_dist=pm.HalfCauchy.dist(6))        
     L2 = pm.expand_packed_triangular(2, packed_L2)    
     cov2 = pm.Deterministic('cov2', L2.dot(L2.T))

     # stack covariance of model 1 and model 2
     cov_both = t.stack([cov1,cov2])

     # Let's assume each component has an independent mean
     mu = pm.Normal('mu', mu=0, sigma=6, shape=(2, 2)) 

     # f(x|mu, cov)
     components = pm.MvNormal('components', mu[category], cov_both[category], shape = (2, 2))      

     obs = pm.Mixture('obs', w=p,
                  comp_dists = components, 
                  observed = [df_bearing.iloc[:, k:k+2] for k in [0, 2]], # to get 2 variables as a df at a time
                   shape=(2, 2))
    
     n_samples = 1000

     step = pm.Metropolis()
     trace_ = pm.sample(n_samples, step, cores=1)

I get this error:
—> 53 components = pm.MvNormal(‘components’, np.array([0, 0]), np.array([[[1, 0], [0,1]], [[1,2], [2,1]]]), shape = 2)
54
55

3 frames
/usr/local/lib/python3.6/dist-packages/pymc3/distributions/multivariate.py in init(self, mu, cov, chol, tau, lower, *args, **kwargs)
59 cov = tt.as_tensor_variable(cov)
60 if cov.ndim != 2:
—> 61 raise ValueError(‘cov must be two dimensional.’)
62 self.chol_cov = cholesky(cov)
63 self.cov = cov

ValueError: cov must be two dimensional.

1 Like
2

Hi rahul.ai and community, a colleague and I am struggling with a similar issue. If anybody had idea that would be super helpful. All the best, Martin

3

Unfortunately, currently the multivariate distributions in PyMC3 are a bit incomplete as they does not support “batch” of multivariate random variable - so you will need to build the MvNormal in a for loop and then use pm.Mixture.

You can find some example here: Mixed multivariate Gauss distribution

1 Like