# Using a mixture of multivariate normal distributions as a prior

**URL:** <https://discourse.pymc.io/t/using-a-mixture-of-multivariate-normal-distributions-as-a-prior/2858>\
**Category:** Questions\
**Created:** [March 5, 2019, 11:35am UTC](https://discourse.pymc.io/t/using-a-mixture-of-multivariate-normal-distributions-as-a-prior/2858 "2019-03-05T11:35:10Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![mattpitkin](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mattpitkin/32/1031_2.png) [@mattpitkin](https://discourse.pymc.io/u/mattpitkin)\
**Post date:** [March 5, 2019, 11:35am UTC](https://discourse.pymc.io/t/using-a-mixture-of-multivariate-normal-distributions-as-a-prior/2858/1 "2019-03-05T11:35:10Z")

</div>

Following on from [Using a Multivariate Normal as a prior](https://discourse.pymc.io/t/using-a-multivariate-normal-as-a-prior/2851) is it possible to use a mixture of multivariate normals as a prior? Can the `pm.Mixture` even be used in defining a prior distribution, as it often seems to require `observed` values?

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [March 5, 2019, 12:09pm UTC](https://discourse.pymc.io/t/using-a-mixture-of-multivariate-normal-distributions-as-a-prior/2858/2 "2019-03-05T12:09:49Z")

</div>

The current mixture Class is not optimal for multivariate distributions (although it has gotten so much better lately thanks to @lucianopaz). If you are using is a prior, make sure you specify the shape carefully.

---

<div class="post-metadata">

**Author:** ![lucianopaz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/lucianopaz/32/2186_2.png) [@lucianopaz](https://discourse.pymc.io/u/lucianopaz)\
**Post date:** [March 5, 2019, 2:37pm UTC](https://discourse.pymc.io/t/using-a-mixture-of-multivariate-normal-distributions-as-a-prior/2858/3 "2019-03-05T14:37:49Z")

</div>

It is possible in the current master branch on GitHub. You can do something like the following

```python
with pm.Model():
    component0 = pm.MvNormal('c0', mu=[0, 0], cov=[[1, 1], [1, 1]], shape=2)
    component1 = pm.MvNormal('c1', mu=[0, 0], cov=[[1, -1], [-1, 1]], shape=2)
    w = pm.Beta('w', 1, 1)
    mix = pm.Mixture('mix', w=[w, 1 - w], comp_dist=[component0.distribution, component1.distribution], shape=2)
    obs = pm.Normal('obs', mu=mix[0] * x + mix[1], sigma= 1, observed=y)

```

---

<div class="post-metadata">

**Author:** ![mattpitkin](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mattpitkin/32/1031_2.png) [@mattpitkin](https://discourse.pymc.io/u/mattpitkin)\
**Post date:** [March 6, 2019, 11:01am UTC](https://discourse.pymc.io/t/using-a-mixture-of-multivariate-normal-distributions-as-a-prior/2858/4 "2019-03-06T11:01:20Z")

</div>

Thanks very much. That works very well. Any idea when the next release will be that includes this?

---

<div class="post-metadata">

**Author:** ![lucianopaz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/lucianopaz/32/2186_2.png) [@lucianopaz](https://discourse.pymc.io/u/lucianopaz)\
**Post date:** [March 6, 2019, 2:18pm UTC](https://discourse.pymc.io/t/using-a-mixture-of-multivariate-normal-distributions-as-a-prior/2858/5 "2019-03-06T14:18:16Z")

</div>

We had a brief discussion last week on whether we could bundle together a nightly-release. I’m not sure when that would be, nor when 3.7 will be released. @colcarroll, @twiecki, @ferrine, do you have any updates on nightly and 3.7 releases?

---

<div class="post-metadata">

**Author:** ![mattpitkin](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mattpitkin/32/1031_2.png) [@mattpitkin](https://discourse.pymc.io/u/mattpitkin)\
**Post date:** [March 6, 2019, 2:21pm UTC](https://discourse.pymc.io/t/using-a-mixture-of-multivariate-normal-distributions-as-a-prior/2858/6 "2019-03-06T14:21:35Z")

</div>

Thanks.

Sorry to add another question, but do you know if it’s also possible to use the mixture within a `pm.Bounded` object, to set lower and upper bounds on the distribution?

---

<div class="post-metadata">

**Author:** ![twiecki](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/twiecki/32/6930_2.png) [@twiecki](https://discourse.pymc.io/u/twiecki)\
**Post date:** [March 6, 2019, 2:36pm UTC](https://discourse.pymc.io/t/using-a-mixture-of-multivariate-normal-distributions-as-a-prior/2858/7 "2019-03-06T14:36:29Z")

</div>

I think we should push out 3.7 asap.

---

<div class="post-metadata">

**Author:** ![lucianopaz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/lucianopaz/32/2186_2.png) [@lucianopaz](https://discourse.pymc.io/u/lucianopaz)\
**Post date:** [March 7, 2019, 6:34am UTC](https://discourse.pymc.io/t/using-a-mixture-of-multivariate-normal-distributions-as-a-prior/2858/8 "2019-03-07T06:34:33Z")

</div>

> [@mattpitkin](#):
>
> Sorry to add another question, but do you know if it’s also possible to use the mixture within a `pm.Bounded` object, to set lower and upper bounds on the distribution?

I haven’t tried, but in principle it should be possible. Let us know if you run in to any trouble.

---

<div class="post-metadata">

**Author:** ![mattpitkin](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mattpitkin/32/1031_2.png) [@mattpitkin](https://discourse.pymc.io/u/mattpitkin)\
**Post date:** [March 8, 2019, 11:18am UTC](https://discourse.pymc.io/t/using-a-mixture-of-multivariate-normal-distributions-as-a-prior/2858/9 "2019-03-08T11:18:35Z")

</div>

I’ve tried adding bounds and it seems to lead to something a bit odd. Below is some example code:

```python
%matplotlib inline

from matplotlib import pyplot as pl
import numpy as np
import pymc3 as pm

# set some example data (which is uniformative)
x = np.linspace(0, 9, 10)
sigma = 100.
data = sigma*np.random.randn(10) # broad (uninformative) Gaussian data

mus = [[0., 0.], [5., 5.]] # means of Gaussians
covs = [np.array([[1., 0.], [0., 1.]]), np.array([[1., 0.9], [0.9, 1.]])] # covariances

weights = [0.7, 0.3] # weights of each mode

lowerbounds = [0., 0.]
upperbounds = [5., 5.]

shape = len(mus[0]) # number of dimensions

# set test values in the centre of the bounds
testvals = np.add(lowerbounds, np.subtract(upperbounds, lowerbounds)/2.)

with pm.Model():
    # define bounded MvNormal
    BoundedMvN = pm.Bound(pm.MvNormal, lower=lowerbounds, upper=upperbounds)

    comp_dists = []
    for i, mu, cov in zip(range(shape), mus, covs):
        comp_dists.append(BoundedMvN('c{}'.format(i), mu=mu, cov=cov, shape=shape).distribution)
    mix = pm.Mixture('mix', w=weights, comp_dists=comp_dists, shape=shape, testval=testvals)

    obs = pm.Normal('obs', mu=mix[0] * x + mix[1], sd=sigma, observed=data)

    trace = pm.sample(2000)

```

If I plot the output:

```auto
pl.plot(trace['mix'][:,0], trace['mix'][:,1], '.')

```

It gives:

 ![example](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/b/b9c76f3c42eeee56f8f40b72c9497f8cf1021e16.png)

which shows it only seems to have sampled from one of the modes of the mixture.

If I remove the bounds, e.g., using:

```auto
with pm.Model():
    comp_dists = []
    for i, mu, cov in zip(range(shape), mus, covs):
        comp_dists.append(pm.MvNormal('c{}'.format(i), mu=mu, cov=cov, shape=shape).distribution)
    mix = pm.Mixture('mix', w=weights, comp_dists=comp_dists, shape=shape, testval=testvals)

    obs = pm.Normal('obs', mu=mix[0] * x + mix[1], sd=sigma, observed=data)

    trace = pm.sample(2000)

```

It samples from both modes as expected:

 ![example](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/5/5a9a0e4214c21b31af2b2301449919a0fb8b247e.png)

Do you have any idea what might be the issue?

---

<div class="post-metadata">

**Author:** ![mattpitkin](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mattpitkin/32/1031_2.png) [@mattpitkin](https://discourse.pymc.io/u/mattpitkin)\
**Post date:** [March 8, 2019, 11:27am UTC](https://discourse.pymc.io/t/using-a-mixture-of-multivariate-normal-distributions-as-a-prior/2858/10 "2019-03-08T11:27:44Z")

</div>

I should add that if I plot the the individual components `c0` and `c1` they do seem to be sampled correctly, its just the mixture `mix` that only shows one component.

---

<div class="post-metadata">

**Author:** ![lucianopaz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/lucianopaz/32/2186_2.png) [@lucianopaz](https://discourse.pymc.io/u/lucianopaz)\
**Post date:** [March 8, 2019, 3:18pm UTC](https://discourse.pymc.io/t/using-a-mixture-of-multivariate-normal-distributions-as-a-prior/2858/11 "2019-03-08T15:18:27Z")

</div>

And if you use `sample_prior_predictive` instead, do you still get just a single component?

---

<div class="post-metadata">

**Author:** ![mattpitkin](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mattpitkin/32/1031_2.png) [@mattpitkin](https://discourse.pymc.io/u/mattpitkin)\
**Post date:** [March 8, 2019, 3:20pm UTC](https://discourse.pymc.io/t/using-a-mixture-of-multivariate-normal-distributions-as-a-prior/2858/12 "2019-03-08T15:20:49Z")

</div>

Using `sample_prior_predictive` gives the following error:

```auto
ValueError: Drawing samples from distributions with array-valued bounds is not supported.

```

---

<div class="post-metadata">

**Author:** ![lucianopaz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/lucianopaz/32/2186_2.png) [@lucianopaz](https://discourse.pymc.io/u/lucianopaz)\
**Post date:** [March 8, 2019, 5:14pm UTC](https://discourse.pymc.io/t/using-a-mixture-of-multivariate-normal-distributions-as-a-prior/2858/13 "2019-03-08T17:14:20Z")

</div>

Hmm, I’ll try to look into the problem when I have some time. The `sample_prior_predictive` part will take some time so I’ll start with the inference part with no data being stuck in a single component.

---

<div class="post-metadata">

**Author:** ![colcarroll](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/colcarroll/32/9_2.png) [@colcarroll](https://discourse.pymc.io/u/colcarroll)\
**Post date:** [March 8, 2019, 8:33pm UTC](https://discourse.pymc.io/t/using-a-mixture-of-multivariate-normal-distributions-as-a-prior/2858/14 "2019-03-08T20:33:46Z")

</div>

I’m actually running the snippet you supplied and it works fine for me with a conda-installed pymc3-3.6. It feels weird to say that the supplied example runs _too_ well 😃 :

```python
import arviz as az

data = az.from_pymc3(trace)
az.plot_joint(data.posterior, var_names=['mix'], coords={'mix_dim_0': np.array([0, 1])});

```

![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/5/53781f2f58e0a26a469c33aba00991ec2fa632da.png)

---

<div class="post-metadata">

**Author:** ![colcarroll](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/colcarroll/32/9_2.png) [@colcarroll](https://discourse.pymc.io/u/colcarroll)\
**Post date:** [March 8, 2019, 9:06pm UTC](https://discourse.pymc.io/t/using-a-mixture-of-multivariate-normal-distributions-as-a-prior/2858/15 "2019-03-08T21:06:03Z")

</div>

Playing around a bit more, it is possible NUTS was just having trouble mixing on your first run. I reran with `target_accept=0.99` since there were a lot of divergences, and get a better mix of samples. Still plenty of divergences, I assume from the boundary conditions.

Here’s the new marginal plot. The scatter is hard to distinguish in the first plot, but you can see differences in the marginal distributions. I’m using a custom stylesheet here!

![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/5/5a847a5e41ee609c2054ce92e017faa1328dced4.png)

---

<div class="post-metadata">

**Author:** ![mattpitkin](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mattpitkin/32/1031_2.png) [@mattpitkin](https://discourse.pymc.io/u/mattpitkin)\
**Post date:** [March 11, 2019, 11:13am UTC](https://discourse.pymc.io/t/using-a-mixture-of-multivariate-normal-distributions-as-a-prior/2858/16 "2019-03-11T11:13:39Z")

</div>

Cool! Trying again it looks like I get decent sampling of both modes now. It must have just been that my initial attempt got a bit stuck!

---

<div class="post-metadata">

**Author:** ![mattpitkin](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mattpitkin/32/1031_2.png) [@mattpitkin](https://discourse.pymc.io/u/mattpitkin)\
**Post date:** [March 12, 2019, 9:42am UTC](https://discourse.pymc.io/t/using-a-mixture-of-multivariate-normal-distributions-as-a-prior/2858/17 "2019-03-12T09:42:57Z")

</div>

Just a quick note - I think part (or maybe all!) of the reason I originally only got one mode was because the likelihood wasn’t entirely uninformative and therefore led to one mode being favoured. If I make the data “noise free” then _I think_ it more consistently returns both modes.