# Mixture Model Metropolis vs. NUTS

**URL:** <https://discourse.pymc.io/t/mixture-model-metropolis-vs-nuts/4783>\
**Category:** Questions\
**Created:** [April 1, 2020, 3:53pm UTC](https://discourse.pymc.io/t/mixture-model-metropolis-vs-nuts/4783 "2020-04-01T15:53:41Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![tvwenger](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/tvwenger/32/7369_2.png) [@tvwenger](https://discourse.pymc.io/u/tvwenger)\
**Post date:** [April 1, 2020, 3:53pm UTC](https://discourse.pymc.io/t/mixture-model-metropolis-vs-nuts/4783/1 "2020-04-01T15:53:41Z")

</div>

New user of PyMC here! I’m attempting to sample a mixture model that is marginalized over the categorical membership of each data point. For example, given some independent data x, observed values y, models M\_i, and model parameters \theta\_i, the likelihood is

P(y | M, x, \theta, \sigma) = \sum\_i q\_i \mathcal{N}(y - M\_i(x, \theta\_i), \sigma)

where q\_i is the Dirichlet membership prior. I’ve set up a minimum working example, which is trying to generate the slopes and intercepts of two lines (M\_i(x, m\_i, b\_i) = m\_i x + b\_i):

```
import numpy as np
np.random.seed(1234)
import pymc3 as pm

num_data = 1000
xdata = np.random.uniform(0.0, 10.0, num_data)
obs = np.zeros(num_data)
model = np.random.choice([0, 1], p=[0.75, 0.25], size=num_data)
m_trues = np.array([6.0, 2.0])
b_trues = np.array([-12.0, 8.0])
sigma = 1.0
obs[model == 0] = m_trues[0]*xdata[model == 0] + b_trues[0]
obs[model == 1] = m_trues[1]*xdata[model == 1] + b_trues[1]
obs += np.random.randn(num_data)*sigma

```

![data](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/e/e734b650c1efb5c65fec6eee714d6545a8ad1084.png)

Using Metropolis steps to sample the posterior, I get reasonable results (i.e., the chains converge to the expected values):

```
with pm.Model() as model:
    p = pm.Dirichlet('q', a=np.ones(2))
    slopes = pm.Normal('slopes', mu=0.0, sigma=10.0, shape=2, testval=[0.0, 0.0])
    intcps = pm.Normal('intcps', mu=0.0, sigma=25.0, shape=2, testval=[-10.0, 10.0], transform=pm.transforms.Ordered())
    mus = slopes * xdata[:, np.newaxis] + intcps
    like = pm.NormalMixture('like', w=p, mu=mus, sigma=sigma, observed=obs)

with model:
    step = pm.Metropolis()
    trace = pm.sample(5000, step=step, tune=1000, cores=8)

```

 ![metropolis](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/0/0b978aa4856bdd1deb84582e644d1968dbc22818.png)

But, with NUTS, I get unreasonable results (i.e., the chains converge, but _not_ to the expected values):

```
with model:
    step = pm.NUTS()
    trace = pm.sample(5000, step=step, tune=1000, cores=8)

```

 ![nuts](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/9/9a6521be3e51174d6d2dcf6f8b5d6e2297b863a5.png)

Is there something simple that I’m missing? Thanks in advance for your help.

---

<div class="post-metadata">

**Author:** ![tvwenger](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/tvwenger/32/7369_2.png) [@tvwenger](https://discourse.pymc.io/u/tvwenger)\
**Post date:** [April 1, 2020, 11:02pm UTC](https://discourse.pymc.io/t/mixture-model-metropolis-vs-nuts/4783/2 "2020-04-01T23:02:11Z")

</div>

Solved the problem! The issue was simply poor prior constraints. The NUTS chains were getting stuck at a local maximum.