# Mixture of hierarchical model

**URL:** <https://discourse.pymc.io/t/mixture-of-hierarchical-model/1298>\
**Category:** Questions\
**Created:** [May 30, 2018, 3:02am UTC](https://discourse.pymc.io/t/mixture-of-hierarchical-model/1298 "2018-05-30T03:02:34Z")\
**Posts on this page:** 1\
**Showing post:** 3

<div class="post-metadata">

**Author:** ![peppe](https://avatars.discourse-cdn.com/v4/letter/p/d26b3c/32.png) [@peppe](https://discourse.pymc.io/u/peppe)\
**Post date:** [May 30, 2018, 9:04am UTC](https://discourse.pymc.io/t/mixture-of-hierarchical-model/1298/3 "2018-05-30T09:04:31Z")

</div>

Hi thanks for the reply.  
I am trying to rewrite the model using pm.Mixtures, by following other discussions:

> [@Gaussian Mixture of regression](https://discourse.pymc.io/t/gaussian-mixture-of-regression/537/2):
>
> Did you check the trace? I am getting the same result from Model 1 (also much faster with NUTS): This is the code I am using (small rewrite so it is more compact): df = pd.read\_csv('mixture.csv') X = df['X'].values[:, None] Y = df['Y'].values k = 2 with pm.Model() as Mixture\_regression1: # Priors for weight parameter π = pm.Dirichlet('π', np.array([1]\*k), testval=np.ones(k)/k) # Priors for unknown model parameters α = pm.Normal('α', mu=0, sd=100, shape=(1, k))…

However, I find it difficult to handle the fact that observations are “nested” in subjects. I tried several models (attached only two, i.e. Model 1 and Model2), but none works. Could you be a bit more specific about the way of rewriting the model?

Thanks.

**Model 1**

```
k=2
#x_t = time
#x_time = shared(x_t)
x_t = time[:,np.newaxis]
x_time = shared(x_t,broadcastable=(False,True))
with Model() as varying_intercept:
    # Priors
    mu_a = Normal('mu_a', mu=0., sd=0.1,shape=k)
    sigma_a = HalfCauchy('sigma_a', 5)
    mu_b = Normal('mu_b', mu=0., sd=0.1,shape=k)
    sigma_b = HalfCauchy('sigma_b', 5)

    # Random intercepts
    a = Normal('a', mu=mu_a, sd=sigma_a, shape=(n_sbj,k))
    b = pm.Normal('b',mu=mu_b,sd=sigma_b,shape=(n_sbj,k))

    # Model error
    sd_y = HalfCauchy('sd_y', 5)

    # Expected value
    y_hat = a[subjects,:] + b[subjects,:] * x_time

    # Data likelihood   
    p = pm.Dirichlet('p', np.ones(k))
    likelihood = pm.NormalMixture('likelihood', p, y_hat, sd=sd_y, observed=observed)

```

> Auto-assigning NUTS sampler…  
> Initializing NUTS using jitter+adapt\_diag…  
> Multiprocess sampling (4 chains in 4 jobs)  
> NUTS: [p\_stickbreaking\_\_, sd\_y\_log\_\_, b, a, sigma\_b\_log\_\_, mu\_b, sigma\_a\_log\_\_, mu\_a]  
> 100%|██████████| 6000/6000 [01:26\<00:00, 69.28it/s]  
> The gelman-rubin statistic is larger than 1.4 for some parameters. The sampler did not converge.  
> The estimated number of effective samples is smaller than 200 for some parameters.

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/1X/d81751c4b062d2ae85fe5e25c66ba2793621c9ad.jpg)

**Model 2**

```
k=2
x_t = time
x_time = shared(x_t)
with Model() as varying_intercept:

    # cluster sizes
    p = pm.Dirichlet('p', np.ones(k))

    mu_a = pm.Normal('mu_a', mu=0, sd=0.1, shape=k) # intercept
    mu_b = pm.Normal('mu_b', mu=0, sd=0.1, shape=k) # slope

    sd_a = HalfCauchy('sd_a', 5,shape=k)
    sd_b = HalfCauchy('sd_b', 5,shape=k)

    a = pm.NormalMixture('mix_mu_a', p, mu_a, sd=sd_a,shape=n_sbj)
    b = pm.NormalMixture('mix_mu_b', p, mu_b, sd=sd_b,shape=n_sbj)

    # Expected value
    y_hat = a[subjects] + b[subjects] * x_time

    #Model error
    sd_y = HalfCauchy('sd_y', 5)

    # Data likelihood
    y_like = Normal('y_like', mu=y_hat, sd=sd_y, observed=observed)

```

> Auto-assigning NUTS sampler…  
> Initializing NUTS using jitter+adapt\_diag…  
> Multiprocess sampling (4 chains in 4 jobs)  
> NUTS: [sd\_y\_log\_\_, mix\_mu\_b, mix\_mu\_a, sd\_b\_log\_\_, sd\_a\_log\_\_, mu\_b, mu\_a, p\_stickbreaking\_\_]  
> 100%|██████████| 6000/6000 [00:37\<00:00, 160.68it/s]  
> There were 1 divergences after tuning. Increase `target_accept` or reparameterize.  
> The gelman-rubin statistic is larger than 1.4 for some parameters. The sampler did not converge.  
> The estimated number of effective samples is smaller than 200 for some parameters.

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/1X/2cf225e12e82879ee89329ee432fec2bcff9bc32.jpg)

---

_[View the full topic](https://discourse.pymc.io/t/mixture-of-hierarchical-model/1298)._