# Multinomial hierarchical regression with multiple observations per group ("Bad energy issue")

**URL:** <https://discourse.pymc.io/t/multinomial-hierarchical-regression-with-multiple-observations-per-group-bad-energy-issue/2755>\
**Category:** Questions\
**Created:** [February 19, 2019, 2:46pm UTC](https://discourse.pymc.io/t/multinomial-hierarchical-regression-with-multiple-observations-per-group-bad-energy-issue/2755 "2019-02-19T14:46:25Z")\
**Posts on this page:** 1\
**Showing post:** 1

<div class="post-metadata">

**Author:** ![AlexAndorra](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/alexandorra/32/9142_2.png) [@AlexAndorra](https://discourse.pymc.io/u/AlexAndorra)\
**Post date:** [February 19, 2019, 2:46pm UTC](https://discourse.pymc.io/t/multinomial-hierarchical-regression-with-multiple-observations-per-group-bad-energy-issue/2755/1 "2019-02-19T14:46:25Z")

</div>

Hi everyone! To try my hand on hierarchical models with PyMC, I’m trying to adapt @twiecki’s [model](https://docs.pymc.io/notebooks/GLM-hierarchical.html#Hierarchical-Model) to a multinomial case with 7 continuous regressors encoded in the variable `X`. The multinomial observations are in `y`, with 7 possible categories.  
There are ~2000 observations in total, broke down by counties (~85), which constitute the groups in the model. Thus, there are multiple observations _per_ group.  
Here is my work in progress:

```
with pm.Model() as multinom_multivariate:
   
    mu_a = pm.Normal('mu_a', mu=0., sd=100.)
    sigma_a = pm.HalfCauchy('sigma_a', 5.)
    mu_b = pm.Normal('mu_b', mu=0., sd=100.)
    sigma_b = pm.HalfCauchy('sigma_b', 5.)   
   
    a = pm.Normal('a', mu=mu_a, sd=sigma_a, shape=(n_counties, n_categs))
    b = pm.Normal('b', mu=mu_b, sd=sigma_b, shape=(n_counties, n_categs, n_regressors))

    results_est = a[counties_idx, :] + pm.math.sum(b[counties_idx, :, :] * np.expand_dims(X, axis=1), axis=2)
    p = pm.Deterministic('p', tt.nnet.softmax(results_est))
   
    likelihood = pm.Multinomial('likelihood', n=y.sum(axis=1), p=p, observed=y)
   
    trace = pm.sample(cores=4)

```

I can’t find what’s misspecified, although it obviously is, since I can’t even sample:

 ![error-model](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/e/e4a5fc7884329e14eaeaec481461787d524c8f12.png)  
I went back and checked: my X and y don’t have any NaNs or 0s, so I can’t see why there would be an infinite value in the log-probability…  
Anyone has any thought? That would be much appreciated! A huge thanks in advance 🙏

---

_[View the full topic](https://discourse.pymc.io/t/multinomial-hierarchical-regression-with-multiple-observations-per-group-bad-energy-issue/2755)._