# Advice with hierarchical model

**URL:** <https://discourse.pymc.io/t/advice-with-hierarchical-model/1925>\
**Category:** Questions\
**Created:** [September 16, 2018, 9:46pm UTC](https://discourse.pymc.io/t/advice-with-hierarchical-model/1925 "2018-09-16T21:46:12Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![rpgoldman](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/rpgoldman/32/1800_2.png) [@rpgoldman](https://discourse.pymc.io/u/rpgoldman)\
**Post date:** [September 16, 2018, 9:46pm UTC](https://discourse.pymc.io/t/advice-with-hierarchical-model/1925/1 "2018-09-16T21:46:12Z")

</div>

I have a three level hierarchical model that has:

1. hyper-hyper parameters (for populations)
2. hyper parameters (for sub-populations)
3. measurements

The measurements are N(mu, sigma) where the mu and sigma are sampled from the hyper parameters.  
The hyper-parameters for each sub-population are mu ~ N(hyper-mu, hyper-sigma) and sigma ~ Gamma(mu=hyper-sd-mean, sigma=hyper-sd-sd).  
I have been having troubles with divergence, so first I reformulated the mus as recommended using

```auto
offset=pm.Normal('offset', mu=0, sd=1)
mu = pm.Deterministic('mu', hyper_mu + offset * hyper_sigma)
pm.Normal('measurement', mu=mu, sd=sd)

```

That worked, but didn’t solve my divergence problem. So I moved to the next step, and reformulated the `sd`:

```auto
sigma_offset = pm.Normal('sigma offset', mu=0, sd=1)
sigma = pm.Deterministic('sigma', largest(hyper_sigma_mu + sigma_offset * hyper_sigma_sd, 0))

```

But when I added this my chains started failing, with this error:

```auto
ValueError: Bad initial energy: inf. The model might be misspecified.

```

It’s quite possible either (a) there’s a deep reason behind this which I don’t know or (b) I have done something stupid (like a typo).  
Can anyone suggest either what is wrong with reformulating the standard deviation hyper-parameters like this or how to find whatever bug there is? I suppose the latter would be to somehow look at the probability function so I can see what is wrong with my code?  
I’m lost here and would be very grateful for any guidance.

---

<div class="post-metadata">

**Author:** ![\_eigenfoo](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/_eigenfoo/32/3059_2.png) [@\_eigenfoo](https://discourse.pymc.io/u/_eigenfoo)\
**Post date:** [September 16, 2018, 11:11pm UTC](https://discourse.pymc.io/t/advice-with-hierarchical-model/1925/2 "2018-09-16T23:11:31Z")

</div>

Perhaps you could try sampling from the prior? `pm.sample_prior_predictive` will generate samples from your specified priors, and you can manually inspect them to see if they make sense. There could be something wrong with how your model specifies the priors.

---

<div class="post-metadata">

**Author:** ![\_eigenfoo](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/_eigenfoo/32/3059_2.png) [@\_eigenfoo](https://discourse.pymc.io/u/_eigenfoo)\
**Post date:** [September 16, 2018, 11:14pm UTC](https://discourse.pymc.io/t/advice-with-hierarchical-model/1925/3 "2018-09-16T23:14:34Z")

</div>

In fact, this is actually one of our FAQs, perhaps you could look here for more info: [Frequently Asked Questions](https://discourse.pymc.io/t/frequently-asked-questions/74/5)

---

<div class="post-metadata">

**Author:** ![rpgoldman](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/rpgoldman/32/1800_2.png) [@rpgoldman](https://discourse.pymc.io/u/rpgoldman)\
**Post date:** [September 17, 2018, 3:55pm UTC](https://discourse.pymc.io/t/advice-with-hierarchical-model/1925/4 "2018-09-17T15:55:13Z")

</div>

Thanks. I have been doing that. I have seen that my variance parameters seemed like they could hit zero which could cause some of my divergences and chain failures. I have eliminated those, but my divergences remain. In `traceplot()` output, in the frequency charts, I see some of the variables, for some of the chains have ‘spikes’ – very tight points with high frequencies. Am I right in thinking that these are suggestive of divergences?

---

<div class="post-metadata">

**Author:** ![rpgoldman](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/rpgoldman/32/1800_2.png) [@rpgoldman](https://discourse.pymc.io/u/rpgoldman)\
**Post date:** [September 17, 2018, 3:59pm UTC](https://discourse.pymc.io/t/advice-with-hierarchical-model/1925/5 "2018-09-17T15:59:19Z")

</div>

Thanks. I looked at the FAQ (might be a good idea to link from the PyMC3 docs to it), but am not sure I see how to proceed further, now that I have removed the chain failure (by ensuring that the variance parameters never go to zero), but my divergences remain.

Is it correct to think that this could be because my variance parameters are still too low, and that’s why the model gets to a point where the gradient is enormous?

Thanks

---

<div class="post-metadata">

**Author:** ![\_eigenfoo](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/_eigenfoo/32/3059_2.png) [@\_eigenfoo](https://discourse.pymc.io/u/_eigenfoo)\
**Post date:** [September 18, 2018, 4:17am UTC](https://discourse.pymc.io/t/advice-with-hierarchical-model/1925/6 "2018-09-18T04:17:49Z")

</div>

+1 to making a FAQ to the docs! I’ll keep that in mind for a future PR 🙂

Here are some tips on divergences that seem applicable to your situation taken from [my blog post](https://eigenfoo.xyz/bayesian-modelling-cookbook/). I might’ve already linked you to the post, but I can’t exactly recall, so better safe than sorry.

- Inspect the `pairplot` of your variables one at a time (if you have a plate of variables, it’s fine to pick a couple at random). It’ll tell you if the two random variables are correlated, and help identify any troublesome neighborhoods in the parameter space (divergent samples will be colored differently, and will cluster near such neighborhoods).
- Maybe trying increasing `target_accept`? Usually 0.9 is a good number (currently the default in PyMC3 is 0.8). This will help get rid of false positives from the test for divergences. However, divergences that _don’t_ go away are cause for alarm.
- Increasing `tune` can sometimes help as well: this gives the sampler more time to 1) find the typical set and 2) find good values for step sizes, scaling factors, etc. If you’re running into divergences, it’s always possible that the sampler just hasn’t started the mixing phase and is still trying to find the typical set.

---

<div class="post-metadata">

**Author:** ![rpgoldman](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/rpgoldman/32/1800_2.png) [@rpgoldman](https://discourse.pymc.io/u/rpgoldman)\
**Post date:** [September 18, 2018, 9:17pm UTC](https://discourse.pymc.io/t/advice-with-hierarchical-model/1925/7 "2018-09-18T21:17:53Z")

</div>

Thanks for the advice. It seems like nudging my variance parameters (and hyper-parameters) up was necessary to fix the problem. Once I ruled out sigma of 0, the divergences went away.  
One remaining problem that I haven’t read about – I am getting complaints about there being not many effective samples for the variables that I added in the reparameterization (i.e., the normals that are used as the offset). Any intuitions about whether this is a real problem or not?

---

<div class="post-metadata">

**Author:** ![rpgoldman](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/rpgoldman/32/1800_2.png) [@rpgoldman](https://discourse.pymc.io/u/rpgoldman)\
**Post date:** [September 18, 2018, 9:23pm UTC](https://discourse.pymc.io/t/advice-with-hierarchical-model/1925/8 "2018-09-18T21:23:25Z")

</div>

By the way, sampling from the priors was very helpful. But it raised an issue with PyMC3 (at least I _think_ it’s an issue). I don’t believe that there’s a plotting function for the output of `sample_prior_predictive` (if there is one, I didn’t find it in the docs). The sampler was neat and its output helpful, but I had to spend a few hours getting the plotting set up. In particular, auto-detecting and dealing with vector-valued random variables was a pain – good things I didn’t have any matrix-valued ones!  
This would be a great feature to add. If it would help to see mine, I’d be happy to share it, but I’m not a good enough pyplot user to contribute it as a PR.

---

<div class="post-metadata">

**Author:** ![\_eigenfoo](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/_eigenfoo/32/3059_2.png) [@\_eigenfoo](https://discourse.pymc.io/u/_eigenfoo)\
**Post date:** [September 19, 2018, 4:07am UTC](https://discourse.pymc.io/t/advice-with-hierarchical-model/1925/9 "2018-09-19T04:07:40Z")

</div>

Again quoting straight from the blog post (we really should port all these tips into an FAQ, shouldn’t we?):

`The number of effective samples is smaller than XYZ for some parameters.`

- Quoting [Junpeng Lao on discourse.pymc3.io](https://discourse.pymc.io/t/the-number-of-effective-samples-is-smaller-than-25-for-some-parameters/1050/3): “A low number of effective samples is usually an indication of strong autocorrelation in the chain.”
- Make sure you’re using an efficient sampler like NUTS. (And not, for instance, Metropolis-Hastings. (I mean seriously, it’s the 21st century, why would you ever want Metropolis-Hastings?))
- Tweak the acceptance probability ( `target_accept` ) - it should be large enough to ensure good exploration, but small enough to not reject all proposals and get stuck.

Also, I think that prior visualization will be moved to ArviZ, see [this issue](https://github.com/pymc-devs/pymc3/issues/3104). But we appreciate any code you think would be good to contribute!

---

<div class="post-metadata">

**Author:** ![rpgoldman](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/rpgoldman/32/1800_2.png) [@rpgoldman](https://discourse.pymc.io/u/rpgoldman)\
**Post date:** [October 8, 2018, 6:08pm UTC](https://discourse.pymc.io/t/advice-with-hierarchical-model/1925/10 "2018-10-08T18:08:14Z")

</div>

Thanks, again. Sorry to have taken so long to respond. Your advice about divergence was very helpful, and between your [site](https://eigenfoo.xyz/bayesian-modelling-cookbook/) and the two articles and [the lecture video](https://www.youtube.com/watch?v=DJ0c7Bm5Djk&feature=youtu.be&t=4h40m9s) it pointed me to, I was able to solve my divergence problems.

A quick follow-up – do you know of resources like that for dealing with failure to get sufficient effective samples?

Here’s my case – I did the divergence trick of replacing

```auto
sigmaRVs = pm.Gamma(self.sigma[i],
                   mu=hyper_sigmaRVs_mu[i],
                   sd=hyper_sigmaRVs_sigma[i],
                   shape=self.num_replicates)

```

with this:

```auto
sigmaRVs = [pm.Deterministic('%s replicate %d' % (self.sigma[i], j),
                            largest(hyper_sigmaRVs_mu[i] +
                                    (sigma_offset[j] * hyper_sigmaRVs_sigma[i]),
                                    0.3))
                    for j in range(self.num_replicates)] 

```

The _good news_ is that this solved my divergence problem (well, that and raising my variance priors, which were too low).

The _bad news_ is that the `sigmaRVs` in the new models are getting hardly any effective samples. The `sigma_offset` RVs get between 75-150 effective samples in chains of 3000, and the `sigmaRVs` approximately 2!

These `sigmas`, and corresponding `mus`, fit as hyperparameters into a normal RV that is observed.

The `mu` RVs show more effective samples, but the same relative relationship: their offset variables have two orders of magnitude more effective samples. The offset RVs have about 2200 effective samples in 4 chains of 1000 samples, each, and the mu RVs only 20-25.

Can anyone suggest something I might do to this model to improve things, or do I just need more iterations?

Thanks