Drawing a RV from a multi-variate normal mixture

Questions
1

Hi again,

I need a bit of help with drawing a random variable from a mixture model. I have to sets of latent parameters, A and B, which are correlated. There are outliers in A, but not in B, so I need to draw the latent variables from a mixture of one population drawn from p(A, B | \mu, \Sigma), and another drawn from p(A, B | \mu, \Sigma_o), where \Sigma_o is a modified covariance matrix.

I’ve got code set up that I’ve adapted from a similar example for univariate normal distributions, but I’m consistently getting the same error. Does anybody know what direction I should be looking for the solution? Will I need to write a pm.Densitydist of the mixture-model likelihood to make this work?

import pymc3 as pm
model = pm.Model()
with model:
    # Draw the covariance matrix and mean values
    packed_L = pm.LKJCholeskyCov('packed_L', n = 2,
                                eta = 2., sd_dist=pm.HalfCauchy.dist(beta=2.5)) 
    L = pm.expand_packed_triangular(2, packed_L)
    Sigma = pm.Deterministic('Sigma', L.dot(L.T))
    mu_ = pm.Normal('mu', mu = 0., sigma = 1., shape=2, testval=np.zeros(2))
                    
    # Outlier treatment
    ## Create the outlier covariance matrix
    s = pm.Lognormal('s', pm.math.log(5), sigma=5.)
    so = tt.eye(2,2,0)
    so1 = tt.set_subtensor(so[1][1], (s+1))
    so = tt.set_subtensor(so[1], so1)    
    Sigmao = (so.dot(Sigma).dot(so))
    
    ## Create the mixture distributions and weighting factor
    inB = pm.MvNormal.dist(mu=mu_, chol=L, shape=[nstars ,2])
    outB = pm.MvNormal.dist(mu=mu_, cov=Sigmao, shape=[nstars ,2])
    Q = pm.Beta('Q', 8,2)
    
    # Draw the latent variables from this mixture
    latents = pm.Mixture('latents', w=[Q, 1-Q], comp_dists = [inB, outB])
    
    obsA = pm.Normal('obsA', mu=latents[:,0], sigma=sigma_obs[0,0], observed=dfout['A'].values)
    obsB = pm.Normal('obsB', mu=latents[:,1], sigma=sigma_obs[0,0], observed=dfout['B'].values)
    
    outtrace = pm.sample()

and the error message I get is:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-140-f08ef490e6e1> in <module>
     22     Q = pm.Beta('Q', 8,2)
     23 
---> 24     latents = pm.Mixture('latents', w=[Q, 1-Q], comp_dists = [inB, outB])
     25 
     26     obsA = pm.Normal('obsA', mu=latents[:,0], sigma=sigma_obs[0,0], observed=dfout['A'].values)

~/Library/Python/3.7/lib/python/site-packages/pymc3/distributions/distribution.py in __new__(cls, name, *args, **kwargs)
     44             total_size = kwargs.pop('total_size', None)
     45             dist = cls.dist(*args, **kwargs)
---> 46             return model.Var(name, dist, data, total_size)
     47         else:
     48             raise TypeError("Name needs to be a string but got: {}".format(name))

~/Library/Python/3.7/lib/python/site-packages/pymc3/model.py in Var(self, name, dist, data, total_size)
    824                 with self:
    825                     var = FreeRV(name=name, distribution=dist,
--> 826                                  total_size=total_size, model=self)
    827                 self.free_RVs.append(var)
    828             else:

~/Library/Python/3.7/lib/python/site-packages/pymc3/model.py in __init__(self, type, owner, index, name, distribution, total_size, model)
   1272             self.tag.test_value = np.ones(
   1273                 distribution.shape, distribution.dtype) * distribution.default()
-> 1274             self.logp_elemwiset = distribution.logp(self)
   1275             # The logp might need scaling in minibatches.
   1276             # This is done in `Factor`.

~/Library/Python/3.7/lib/python/site-packages/pymc3/distributions/mixture.py in logp(self, value)
    413         w = self.w
    414 
--> 415         return bound(logsumexp(tt.log(w) + self._comp_logp(value), axis=-1),
    416                      w >= 0, w <= 1, tt.allclose(w.sum(axis=-1), 1),
    417                      broadcast_conditions=False)

~/Library/Python/3.7/lib/python/site-packages/pymc3/distributions/mixture.py in _comp_logp(self, value)
    275         else:
    276             return tt.squeeze(tt.stack([comp_dist.logp(value)
--> 277                                         for comp_dist in comp_dists],
    278                                        axis=-1))
    279 

~/Library/Python/3.7/lib/python/site-packages/pymc3/distributions/mixture.py in <listcomp>(.0)
    275         else:
    276             return tt.squeeze(tt.stack([comp_dist.logp(value)
--> 277                                         for comp_dist in comp_dists],
    278                                        axis=-1))
    279 

~/Library/Python/3.7/lib/python/site-packages/pymc3/distributions/multivariate.py in logp(self, value)
    316         TensorVariable
    317         """
--> 318         quaddist, logdet, ok = self._quaddist(value)
    319         k = value.shape[-1].astype(theano.config.floatX)
    320         norm = - 0.5 * k * pm.floatX(np.log(2 * np.pi))

~/Library/Python/3.7/lib/python/site-packages/pymc3/distributions/multivariate.py in _quaddist(self, value)
     84         mu = self.mu
     85         if value.ndim > 2 or value.ndim == 0:
---> 86             raise ValueError('Invalid dimension for value: %s' % value.ndim)
     87         if value.ndim == 1:
     88             onedim = True

ValueError: Invalid dimension for value: 0

Cheers!
Oli

2

Hi! I managed to fix the problem but passing a shape kwarg to the pm.Mixture call. The functional code is:

model = pm.Model()
with model:
    # Draw the covariance matrix and mean values
    packed_L = pm.LKJCholeskyCov('packed_L', n = 2,
                                eta = 2., sd_dist=pm.HalfCauchy.dist(beta=2.5)) 
    L = pm.expand_packed_triangular(2, packed_L)
    Sigma = pm.Deterministic('Sigma', L.dot(L.T))
    mu_ = pm.Normal('mu', mu = 0., sigma = 1., shape=2, testval=np.zeros(2))
                    
    # Outlier treatment
    ## Create the outlier covariance matrix
    s = pm.Lognormal('s', tt.log(5), sigma=5.)
    so = tt.eye(2,2,0)
    so1 = tt.set_subtensor(so[1][1], (s+1))
    so = tt.set_subtensor(so[1], so1)    
    Sigmao = pm.Deterministic('Sigmao', (so.dot(Sigma).dot(so)))
    
    ## Create the mixture distributions and weighting factor
    inB = pm.MvNormal.dist(mu = mu_, chol = L, shape = [nstars,2])
    outB = pm.MvNormal.dist(mu = mu_, cov = Sigmao, shape = [nstars,2])
    Q = pm.Beta('Q', 8,2)
    
    latents = pm.Mixture('latents', w=[Q, 1-Q], comp_dists = [inB, outB], shape=(nstars,2))
    
    obsA = pm.Normal('obsA', mu=latents[:,0], sigma=sigma_obs[0,0], observed=dfout['A'].values)
    obsB = pm.Normal('obsB', mu=latents[:,1], sigma=sigma_obs[0,0], observed=dfout['B'].values)
    
    outtrace = pm.sample()
1 Like
3

@ojhall94, can you elaborate on the shape definition, what parts of the component distributions were used.
Can you elaborate on where nstars is defined and what it is supposed to represent, ie is it a number proportional to the number of MvNormals, if so what is the 2 in the shape of the Mvnormal snd mixture.