# Hierarchical Mixture Model

**URL:** <https://discourse.pymc.io/t/hierarchical-mixture-model/2898>\
**Category:** Questions\
**Created:** [March 11, 2019, 4:52pm UTC](https://discourse.pymc.io/t/hierarchical-mixture-model/2898 "2019-03-11T16:52:06Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![rz5](https://avatars.discourse-cdn.com/v4/letter/r/f0a364/32.png) [@rz5](https://discourse.pymc.io/u/rz5)\
**Post date:** [March 11, 2019, 4:52pm UTC](https://discourse.pymc.io/t/hierarchical-mixture-model/2898/1 "2019-03-11T16:52:06Z")

</div>

Hi,

I’m new to working with PyMC3 and I’m trying to specify a hierarchical mixture model to cluster house types (ie one story vs two story or something of that nature) based on their sale prices, using the house’s county and town as nested covariates. In my model, each town has mean house price `beta_j`, but I would like to add mixture priors to the model, where each town’s house prices are drawn from one of `K` normal distributions. I am able to include the hierarchical portion of the model without issues, but haven’t been successful trying to code the mixture priors for each town.

Below is an example of one of the models I have tried that hasn’t worked. I am getting an error because the number of observations being fed in is different than the number of distributions:

ValueError: Input dimension mis-match. (input[0].shape[0] = 10000, input[1].shape[0] = 3)

```auto
model = pm.Model()

with model:
  mu = pm.Normal('mu', mu=GM, sd=SD) # overall mean of house prices
  sigma = pm.HalfCauchy('sigma', beta=1) # overall standard

  alpha = pm.Normal('alpha', mu=mu, sd=sigma, shape=I) # county level mean
  gamma = pm.Uniform('gamma', lower=0, upper=sigma, shape=I) # county level standard deviation

  beta = pm.Normal('beta', # town level mean
                   mu=alpha[county_indices],
                   sd=gamma[county_indices],
                   shape=county_count)
  omega = pm.Uniform('omega', # town level standard deviation
                     lower=0,
                     upper=gamma[county_indices],
                     shape=county_count)

  p = pm.Dirichlet('p', np.ones(K), shape=(town_count, K)) # component mixing proportions
  theta = pm.Normal('theta', # component level means for each town
                    mu=beta[town_indices],
                    sd=omega[town_indices],
                    shape=(K, town_count))
  eta = pm.Uniform('eta', # component level standard deviations for each town
                   lower=0,
                   upper=omega[town_indices],
                   shape=(K, town_count))
  y = pm.NormalMixture('y',
                 p[indexed_transactions.town_idx.values, :],
                 mu=theta[:, indexed_transactions.town_idx.values],
                 sd=eta[:, indexed_transactions.town_idx.values],
                 observed=indexed_transactions.SalePrice)

  # approx = pm.fit(n=50000)
  trace = pm.sample(10000, chains=1)

```

I suspect that my issues may stem from either how I am specifying the Dirichlet priors, or how I am slicing the component-level mean and standard deviation matrices, but at this point I’m not too sure.

The code I am using the generate the data follows [this post](http://www.databozo.com/deep-in-the-weeds-complex-hierarchical-models-in-pymc3) and is included below.

```auto
import pandas as pd
import pymc3 as pm
import numpy as np
import matplotlib.pyplot as plt

I = 4 # number of counties
town_mult = 5 # number of towns per county
J = I*town_mult # total number of towns
K = 3 # number of mixture components
GM = 250000 # grand mean of all house prices
SD = 50000 # standard deviation of all house prices
N = 10000 # number of observations
spread = 200000 # difference between component means

comp_centers = np.array([GM-spread, GM, GM+spread])
comp_vals = np.random.randint(0, K, N)
alpha = np.random.normal(loc=comp_centers, scale=SD/2, size=(I, K)) # component mean for each county
county_vals = np.random.randint(0, I, N)
town_vals = np.ones(N)
for i in range(len(county_vals)):
  town_vals[i] = np.random.randint(county_vals[i]*town_mult,
                                   county_vals[i]*town_mult+I)

alpha = np.reshape(alpha, I*K)
beta = np.random.normal(loc=alpha, scale=SD/2, size=(J, I*K)) # component mean for each town
y = np.ones(N)
for i in range(N): # simulate individual house prices
  y[i] = np.random.normal(loc=beta[int(town_vals[i]),
                                   int((county_vals[i]*K)+comp_vals[i])],
                          scale=SD/2)

transactions = pd.DataFrame({'County': county_vals,
                             'Town': town_vals,
                             'SalePrice': y})

county_idx = transactions.groupby(['County']).all().reset_index().reset_index()[['index', 'County']]
county_idx.columns = ['county_idx', 'County']
county_town_idx = transactions.groupby(['County', 'Town']).all().reset_index().reset_index()[['index', 'County', 'Town']]
county_town_idx.columns = ['town_idx', 'County', 'Town']

county_idx_df = pd.merge(county_idx, county_town_idx, how='inner', on='County')

indexed_transactions = pd.merge(transactions,
                                county_idx_df,
                                how='inner',
                                on=['County', 'Town']).reset_index()

county_indices = county_idx_df.county_idx.values
county_count = len(county_indices)

town_indices = county_idx_df.town_idx.values
town_count = len(town_indices)

```

I appreciate any help you can offer!

---

<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 12, 2019, 11:36am UTC](https://discourse.pymc.io/t/hierarchical-mixture-model/2898/2 "2019-03-12T11:36:31Z")

</div>

You need to pass the correct `shape` and `comp_shape` into `pm.NormalMixture`, because it is a mixture of multivariate normals.

I’m a bit lost by the many index values you have, so I’ll tell you the idea of how to structure the `shape` and `comp_shape` values. `shape` should be the resulting random variable’s shape, this means that it should not contain any information on the number of mixture components. `comp_shape` should be the shape of the multidimensional normal random variables that go into the mixture, plus the number of mixture components as the **last** axis.

Without the `indexed_transactions.town_idx.values` indexing, your mixture should be defined somewhat like this:

```python
  theta = pm.Normal('theta', # component level means for each town
                    mu=beta[town_indices, None], # Maybe you should tt.stack due to #3147
                    sd=omega[town_indices, None], # Maybe you should tt.stack due to #3147
                    shape=(town_count, K)) # The mixture components must be the last index
  eta = pm.Uniform('eta', # component level standard deviations for each town
                   lower=0,
                   upper=omega[town_indices, None], # Maybe you should tt.stack due to #3147
                   shape=(town_count, K)) # The mixture components must be the last index
  y = pm.NormalMixture('y',
                 p,
                 mu=theta,
                 sd=eta,
                 comp_shape=(town_count, K),
                 shape=(town_count,),
                 observed=...)

```

I hope this helps you figure out how to write up the correct shapes for your use case with the `indexed_transactions` indexing.

---

<div class="post-metadata">

**Author:** ![rz5](https://avatars.discourse-cdn.com/v4/letter/r/f0a364/32.png) [@rz5](https://discourse.pymc.io/u/rz5)\
**Post date:** [March 13, 2019, 8:12pm UTC](https://discourse.pymc.io/t/hierarchical-mixture-model/2898/3 "2019-03-13T20:12:29Z")

</div>

It’s running now. Thanks for your help!