# 2D Gaussian Mixture Model

**URL:** <https://discourse.pymc.io/t/2d-gaussian-mixture-model/286>\
**Category:** Questions\
**Created:** [August 24, 2017, 9:41pm UTC](https://discourse.pymc.io/t/2d-gaussian-mixture-model/286 "2017-08-24T21:41:20Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![cmlakhan](https://avatars.discourse-cdn.com/v4/letter/c/278dde/32.png) [@cmlakhan](https://discourse.pymc.io/u/cmlakhan)\
**Post date:** [August 24, 2017, 9:41pm UTC](https://discourse.pymc.io/t/2d-gaussian-mixture-model/286/1 "2017-08-24T21:41:20Z")

</div>

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](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](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.

---

<div class="post-metadata">

**Author:** ![AustinRochford](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/austinrochford/32/92_2.png) [@AustinRochford](https://discourse.pymc.io/u/AustinRochford)\
**Post date:** [August 25, 2017, 12:08pm UTC](https://discourse.pymc.io/t/2d-gaussian-mixture-model/286/2 "2017-08-25T12:08:13Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![cmlakhan](https://avatars.discourse-cdn.com/v4/letter/c/278dde/32.png) [@cmlakhan](https://discourse.pymc.io/u/cmlakhan)\
**Post date:** [August 25, 2017, 6:17pm UTC](https://discourse.pymc.io/t/2d-gaussian-mixture-model/286/3 "2017-08-25T18:17:40Z")

</div>

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)
```

---

<div class="post-metadata">

**Author:** ![cmlakhan](https://avatars.discourse-cdn.com/v4/letter/c/278dde/32.png) [@cmlakhan](https://discourse.pymc.io/u/cmlakhan)\
**Post date:** [August 25, 2017, 6:24pm UTC](https://discourse.pymc.io/t/2d-gaussian-mixture-model/286/4 "2017-08-25T18:24:14Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![AustinRochford](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/austinrochford/32/92_2.png) [@AustinRochford](https://discourse.pymc.io/u/AustinRochford)\
**Post date:** [August 25, 2017, 6:29pm UTC](https://discourse.pymc.io/t/2d-gaussian-mixture-model/286/5 "2017-08-25T18:29:13Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![cmlakhan](https://avatars.discourse-cdn.com/v4/letter/c/278dde/32.png) [@cmlakhan](https://discourse.pymc.io/u/cmlakhan)\
**Post date:** [August 25, 2017, 6:40pm UTC](https://discourse.pymc.io/t/2d-gaussian-mixture-model/286/6 "2017-08-25T18:40:43Z")

</div>

Yes here you go

> <https://github.com/cmlakhan/notebooks/blob/master/multivariate_normal.ipynb>

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](http://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:

---

<div class="post-metadata">

**Author:** ![cmlakhan](https://avatars.discourse-cdn.com/v4/letter/c/278dde/32.png) [@cmlakhan](https://discourse.pymc.io/u/cmlakhan)\
**Post date:** [August 25, 2017, 6:42pm UTC](https://discourse.pymc.io/t/2d-gaussian-mixture-model/286/7 "2017-08-25T18:42:37Z")

</div>

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)
```

---

<div class="post-metadata">

**Author:** ![cmlakhan](https://avatars.discourse-cdn.com/v4/letter/c/278dde/32.png) [@cmlakhan](https://discourse.pymc.io/u/cmlakhan)\
**Post date:** [August 25, 2017, 9:40pm UTC](https://discourse.pymc.io/t/2d-gaussian-mixture-model/286/8 "2017-08-25T21:40:04Z")

</div>

I also tried it using pm.Mixture

> <https://github.com/cmlakhan/notebooks/blob/master/multivariate2.ipynb>

essentially just copying what you have done here

> <https://gist.github.com/AustinRochford/41109579c8be23a10e2cd167209dbe25>

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\>)