2D Gaussian Mixture Model

Questions
1

I am a newbie to PyMC3 but I was wondering if there is an example of how one would model a full 2D mixture of Gaussian model. In a sense I would like to combine the example from here

https://pymc-devs.github.io/pymc3/notebooks/LKJ.html

with the marginalized mixture of Gaussian example

https://pymc-devs.github.io/pymc3/notebooks/marginalized_gaussian_mixture_model.html

my sense is that pm.NormalMixture only fits 1D data, is that correct?

If someone could give me pointers for how to start that would be wonderful.

2

You are correct; unfortunately NormalMixture only supports 1D observations at this time. This should be possible in an unmarginalized way by sampling the mixture components from a categorical distribution. This approach, however, will take much longer to mix.

3

I am a little confused about how to group covariance matrices here is my model, where i create two covariance matrices cov_1 and cov_2. Is this not correct? It should be a 2D mixture of gaussians with two components.

import pandas as pd
import numpy as np
import pymc3 as pm
import matplotlib.pyplot as plt
import seaborn as sns
import theano.tensor as tt

nclus=2

Count_Row=x.shape[0]

with pm.Model() as model:

p = pm.Dirichlet('p', a=np.array([1., 1.]), shape=nclus)    
category = pm.Categorical('category', p=p, shape=Count_Row)

with model:

sigma_1 = pm.Lognormal('sigma_1', np.zeros(2), np.ones(2), shape=2)
nu_1 = pm.Uniform('nu_1', 0, 5)
C_triu_1 = pm.LKJCorr('C_triu_1', nu_1, 2)
C_1 = pm.Deterministic('C_1', tt.fill_diagonal(C_triu_1[np.zeros((2, 2), dtype=np.int64)], 1.))
sigma_diag_1 = pm.Deterministic('sigma_mat_1', tt.nlinalg.diag(sigma_1))
cov_1 = pm.Deterministic('cov_1', tt.nlinalg.matrix_dot(sigma_diag_1, C_1, sigma_diag_1))

with model:

sigma_2 = pm.Lognormal('sigma_2', np.zeros(2), np.ones(2), shape=2)
nu_2 = pm.Uniform('nu_2', 0, 5)
C_triu_2 = pm.LKJCorr('C_triu_2', nu_2, 2)
C_2 = pm.Deterministic('C_2', tt.fill_diagonal(C_triu_2[np.zeros((2, 2), dtype=np.int64)], 1.))
sigma_diag_2 = pm.Deterministic('sigma_mat_2', tt.nlinalg.diag(sigma_2))
cov_2 = pm.Deterministic('cov_2', tt.nlinalg.matrix_dot(sigma_diag_2, C_2, sigma_diag_2))    
cov_both = tt.stack([cov_1,cov_2])

with model:

mu_both = pm.MvNormal('mu', mu=np.zeros(2), shape=(2,2))
x_ = pm.MvNormal('x', mu_both[category], cov_both[category], observed=x)
4

If I use

cov_both = np.stack([cov_1,cov_2], axis=0)

then I get the following error

IndexError Traceback (most recent call last)
in ()
34 with model:
35 mu_both = pm.MvNormal(‘mu’, mu=np.zeros(2), shape=(2,2))
—> 36 x_ = pm.MvNormal(‘x’, mu_both[category], cov_both[category], observed=x)
37

IndexError: only integers, slices (:), ellipsis (...), numpy.newaxis (None) and integer or boolean arrays are valid indices

5

What goes wrong with the first code you posted? Can you link to a notebook in a Gist to make reading it easyer?

6

Yes here you go

I get the following error

AssertionError Traceback (most recent call last)
in ()
32 with model:
33 mu_both = pm.MvNormal(‘mu’, mu=np.zeros(2), shape=(2,2))
—> 34 x_ = pm.MvNormal(‘x’, mu_both[category], cov_both[category], observed=x)

/Users/chirag/anaconda/lib/python3.6/site-packages/pymc3/distributions/distribution.py in new(cls, name, *args, **kwargs)
28 if isinstance(name, string_types):
29 data = kwargs.pop(‘observed’, None)
—> 30 dist = cls.dist(*args, **kwargs)
31 return model.Var(name, dist, data)
32 else:

/Users/chirag/anaconda/lib/python3.6/site-packages/pymc3/distributions/distribution.py in dist(cls, *args, **kwargs)
39 def dist(cls, *args, **kwargs):
40 dist = object.new(cls)
—> 41 dist.init(*args, **kwargs)
42 return dist
43

/Users/chirag/anaconda/lib/python3.6/site-packages/pymc3/distributions/multivariate.py in init(self, mu, cov, tau, *args, **kwargs)
92 DeprecationWarning)
93 self.mean = self.median = self.mode = self.mu = mu
—> 94 self.tau, self.cov = get_tau_cov(mu, tau=tau, cov=cov)
95
96 def random(self, point=None, size=None):

/Users/chirag/anaconda/lib/python3.6/site-packages/pymc3/distributions/multivariate.py in get_tau_cov(mu, tau, cov)
49 tau = np.eye(len(mu))
50 else:
—> 51 tau = tt.nlinalg.matrix_inverse(cov)
52
53 else:

/Users/chirag/anaconda/lib/python3.6/site-packages/theano/gof/op.py in call(self, *inputs, **kwargs)
613 “”"
614 return_list = kwargs.pop(‘return_list’, False)
–> 615 node = self.make_node(*inputs, **kwargs)
616
617 if config.compute_test_value != ‘off’:

/Users/chirag/anaconda/lib/python3.6/site-packages/theano/tensor/nlinalg.py in make_node(self, x)
71 def make_node(self, x):
72 x = as_tensor_variable(x)
—> 73 assert x.ndim == 2
74 return Apply(self, [x], [x.type()])
75

AssertionError:

7

I have a feeling is is a problem with how the mean and covariance matrix is being indexed

problem seems to be this line

x_ = pm.MvNormal('x', mu_both[category], cov_both[category], observed=x)
8

I also tried it using pm.Mixture

essentially just copying what you have done here

but get the following error:

AsTensorError Traceback (most recent call last)
in ()
1 with model:
----> 2 trace = pm.sample(1000, pm.Metropolis())[500::5]

/Users/chirag/anaconda/lib/python3.6/site-packages/pymc3/step_methods/arraystep.py in new(cls, *args, **kwargs)
58 # If we don’t return the instance we have to manually
59 # call init
—> 60 step.init([var], *args, **kwargs)
61 # Hack for creating the class correctly when unpickling.
62 step.__newargs = ([var], ) + args, kwargs

/Users/chirag/anaconda/lib/python3.6/site-packages/pymc3/step_methods/metropolis.py in init(self, vars, S, proposal_dist, scaling, tune, tune_interval, model, **kwargs)
96
97 shared = pm.make_shared_replacements(vars, model)
—> 98 self.delta_logp = delta_logp(model.logpt, vars, shared)
99 super(Metropolis, self).init(vars, shared)
100

/Users/chirag/anaconda/lib/python3.6/site-packages/pymc3/step_methods/metropolis.py in delta_logp(logp, vars, shared)
424
425 def delta_logp(logp, vars, shared):
–> 426 [logp0], inarray0 = pm.join_nonshared_inputs([logp], vars, shared)
427
428 tensor_type = inarray0.type

/Users/chirag/anaconda/lib/python3.6/site-packages/pymc3/theanof.py in join_nonshared_inputs(xs, vars, shared, make_shared)
198 replace.update(shared)
199
–> 200 xs_special = [theano.clone(x, replace, strict=False) for x in xs]
201 return xs_special, inarray
202

/Users/chirag/anaconda/lib/python3.6/site-packages/pymc3/theanof.py in (.0)
198 replace.update(shared)
199
–> 200 xs_special = [theano.clone(x, replace, strict=False) for x in xs]
201 return xs_special, inarray
202

/Users/chirag/anaconda/lib/python3.6/site-packages/theano/scan_module/scan_utils.py in clone(output, replace, strict, share_inputs, copy_inputs)
256 [],
257 strict,
–> 258 share_inputs)
259
260 return outs

/Users/chirag/anaconda/lib/python3.6/site-packages/theano/compile/pfunc.py in rebuild_collect_shared(outputs, inputs, replace, updates, rebuild_strict, copy_inputs_over, no_default_updates)
230 else:
231 if isinstance(outputs, Variable):
–> 232 cloned_v = clone_v_get_shared_updates(outputs, copy_inputs_over)
233 cloned_outputs = cloned_v
234 # computed_list.append(cloned_v)

/Users/chirag/anaconda/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
91 if owner not in clone_d:
92 for i in owner.inputs:
—> 93 clone_v_get_shared_updates(i, copy_inputs_over)
94
95 clone_d[owner] = owner.clone_with_new_inputs(

/Users/chirag/anaconda/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
91 if owner not in clone_d:
92 for i in owner.inputs:
—> 93 clone_v_get_shared_updates(i, copy_inputs_over)
94
95 clone_d[owner] = owner.clone_with_new_inputs(

/Users/chirag/anaconda/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
91 if owner not in clone_d:
92 for i in owner.inputs:
—> 93 clone_v_get_shared_updates(i, copy_inputs_over)
94
95 clone_d[owner] = owner.clone_with_new_inputs(

/Users/chirag/anaconda/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
91 if owner not in clone_d:
92 for i in owner.inputs:
—> 93 clone_v_get_shared_updates(i, copy_inputs_over)
94
95 clone_d[owner] = owner.clone_with_new_inputs(

/Users/chirag/anaconda/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
91 if owner not in clone_d:
92 for i in owner.inputs:
—> 93 clone_v_get_shared_updates(i, copy_inputs_over)
94
95 clone_d[owner] = owner.clone_with_new_inputs(

/Users/chirag/anaconda/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
91 if owner not in clone_d:
92 for i in owner.inputs:
—> 93 clone_v_get_shared_updates(i, copy_inputs_over)
94
95 clone_d[owner] = owner.clone_with_new_inputs(

/Users/chirag/anaconda/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
91 if owner not in clone_d:
92 for i in owner.inputs:
—> 93 clone_v_get_shared_updates(i, copy_inputs_over)
94
95 clone_d[owner] = owner.clone_with_new_inputs(

/Users/chirag/anaconda/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
91 if owner not in clone_d:
92 for i in owner.inputs:
—> 93 clone_v_get_shared_updates(i, copy_inputs_over)
94
95 clone_d[owner] = owner.clone_with_new_inputs(

/Users/chirag/anaconda/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
94
95 clone_d[owner] = owner.clone_with_new_inputs(
—> 96 [clone_d[i] for i in owner.inputs], strict=rebuild_strict)
97 for old_o, new_o in zip(owner.outputs, clone_d[owner].outputs):
98 clone_d.setdefault(old_o, new_o)

/Users/chirag/anaconda/lib/python3.6/site-packages/theano/gof/graph.py in clone_with_new_inputs(self, inputs, strict)
240 remake_node = True
241 if remake_node:
–> 242 new_node = self.op.make_node(*new_inputs)
243 new_node.tag = copy(self.tag).update(new_node.tag)
244 else:

/Users/chirag/anaconda/lib/python3.6/site-packages/theano/tensor/elemwise.py in make_node(self, *inputs)
574 using DimShuffle.
575 “”"
–> 576 inputs = list(map(as_tensor_variable, inputs))
577 out_dtypes, out_broadcastables, inputs = self.get_output_info(
578 DimShuffle, *inputs)

/Users/chirag/anaconda/lib/python3.6/site-packages/theano/tensor/basic.py in as_tensor_variable(x, name, ndim)
169 if not isinstance(x.type, TensorType):
170 raise AsTensorError(
–> 171 “Variable type field must be a TensorType.”, x, x.type)
172
173 if ndim is None:

AsTensorError: (‘Variable type field must be a TensorType.’, mu0_shared, <theano.gof.type.Generic object at 0x111f763c8>)