# Sampling from nested models

**URL:** <https://discourse.pymc.io/t/sampling-from-nested-models/7511>\
**Category:** Questions\
**Created:** [May 26, 2021, 8:33pm UTC](https://discourse.pymc.io/t/sampling-from-nested-models/7511 "2021-05-26T20:33:36Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![non87](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/non87/32/3859_2.png) [@non87](https://discourse.pymc.io/u/non87)\
**Post date:** [May 26, 2021, 8:33pm UTC](https://discourse.pymc.io/t/sampling-from-nested-models/7511/1 "2021-05-26T20:33:36Z")

</div>

Say we generate data by nesting a latent variable from a logistic model within a Poisson model (perhaps “nesting” is not the best word here, but I do not know a better one). For example:

```auto
import pymc3 as pm
import numpy as np

# generate data
rng = np.random.default_rng(12345678)

# Parameters of interest
a = 0.9
b = -0.7
c = -0.3
d = 0.5
x_1 = rng.random(1000)
# Latent is the logistic reg. latent variable
latent = a + x_1*b
p = 1/(1+np.exp(-latent))
# observed trials
trial_success = np.int_( rng.random(x_1.shape[0]) < p)
# Now the number of events
lam = np.exp(c + d*latent)
num_event = rng.poisson(lam=lam)

```

Intuitively, the way I would model this is to nest a logistic regression inside a Poisson model:

```auto
# Now both at the same time
with pm.Model() as simple_log:
    a_ = pm.Normal('a', 0, sd=100)
    b_ = pm.Normal('b', 0, sd=100)
    c_ = pm.Normal('c', 0, sd=100)
    d_ = pm.Normal('d', 0, sd=100)
    # Latent specification
    # Latent for the bernoulli
    latent_ = pm.Deterministic(name='latent_bin', var= a_ + b_ * x_1)
    # Latent for the poisson
    latent_2 = pm.Deterministic(name='latent_pois', var=c_ + latent_*d_ )
    # Bernoulli random vector with probability of success
    pm.Bernoulli(name='logit', p=pm.invlogit(latent_), observed=trial_success)
    # Poisson random vector
    pm.Poisson(name='pois', mu=pm.math.exp(latent_2), observed=numbers)

```

This is the DAG of the model, which corresponds to the generating process above (as far as I can understand):

 ![Digraph.gv](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/c/c096cecb3c8652819eddd21a6ed9e32459e3c486.png)

However, I am having a really hard time to sample from this model, as I anticipated. Is there a good way to re-parametrize the model – or any miraculous receipt for sampling?

To stabilize sampling, I can think about ways to independently sample the logistic and Poisson part of the model. First, I can fit a standard logistic regression. Then, if I fix the a and b parameters to their posterior expected value in the Poisson part, the Poisson part becomes a standard GLM again. Yet, if I fix the a and b parameters, I am losing all information about the variability in my posteriors. So I am wondering if there is any way to incorporate that information.

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [May 27, 2021, 5:45am UTC](https://discourse.pymc.io/t/sampling-from-nested-models/7511/2 "2021-05-27T05:45:47Z")

</div>

What problems are you seeing?

---

<div class="post-metadata">

**Author:** ![non87](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/non87/32/3859_2.png) [@non87](https://discourse.pymc.io/u/non87)\
**Post date:** [May 27, 2021, 1:14pm UTC](https://discourse.pymc.io/t/sampling-from-nested-models/7511/3 "2021-05-27T13:14:41Z")

</div>

All the signs of non-convergence you can get, starting from a _very_ slow sampling:

```auto
with simple_log:
    log_trace = pm.sample(4000, chains=3, tune=1000, random_seed=101010, return_inferencedata=True, target_accept=0.9)

```

Output:

```auto
Sampling 3 chains for 1_000 tune and 4_000 draw iterations (3_000 + 12_000 draws total) took 750 seconds.
There was 1 divergence after tuning. Increase `target_accept` or reparameterize.
There were 9 divergences after tuning. Increase `target_accept` or reparameterize.
The acceptance probability does not match the target. It is 0.8152433470462593, but should be close to 0.9. Try to increase the number of tuning steps.
There were 4 divergences after tuning. Increase `target_accept` or reparameterize.
The chain reached the maximum tree depth. Increase max_treedepth, increase target_accept or reparameterize.
The rhat statistic is larger than 1.4 for some parameters. The sampler did not converge

```

The sampling is clearly biased – or very slow to mix. For example the d parameter has an expected posterior of -13 (should be 0.5).

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [May 27, 2021, 1:31pm UTC](https://discourse.pymc.io/t/sampling-from-nested-models/7511/4 "2021-05-27T13:31:41Z")

</div>

Your priors also look a bit too broad.

Other than that you might have some redundancy going on in your linear models, I would try to set some of the intercepts / coefficients to the true values and see if it samples better, to identify the culprits.

Can you show some trace plots and pairwise correlation plots of the posterior parameters?

---

<div class="post-metadata">

**Author:** ![non87](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/non87/32/3859_2.png) [@non87](https://discourse.pymc.io/u/non87)\
**Post date:** [May 27, 2021, 1:52pm UTC](https://discourse.pymc.io/t/sampling-from-nested-models/7511/5 "2021-05-27T13:52:25Z")

</div>

I have never seen such bad plots 🙂

 ![trace](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/b/baeaaa8652da377e30580eab128fd76d30a44e04.png)  
 ![pairwise](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/8/8839300c7402463ccfb758148d6c2feba9602650.png)

I would conclude that yes, the sample does not converge in the right region. The pairwise plot shows that the parameters are really correlated. In particular d, c and b (well, 3 out of 4).

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [May 27, 2021, 1:59pm UTC](https://discourse.pymc.io/t/sampling-from-nested-models/7511/6 "2021-05-27T13:59:40Z")

</div>

It seems that the `latent` does not vary enough to inform the `lam` parameter. In particular you can reverse `c` and `d` and get basically the same `lam` parameter

```python
sns.kdeplot(c + latent * d)
sns.kdeplot(-c + latent * (-d))

```

![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/7/7de1a2648284f8baa73535b61fcc887d58467226.png)

---

<div class="post-metadata">

**Author:** ![non87](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/non87/32/3859_2.png) [@non87](https://discourse.pymc.io/u/non87)\
**Post date:** [May 27, 2021, 2:04pm UTC](https://discourse.pymc.io/t/sampling-from-nested-models/7511/7 "2021-05-27T14:04:24Z")

</div>

Agreed, but I do not understand why. latent\_ shouldn’t be that uninformative, it has variation…  
In any case, 1. this is screaming for re-parametrization (but I would not know how, exactly) and 2. I will follow your advice to fix some parameters to see if it gets better.

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [May 27, 2021, 2:07pm UTC](https://discourse.pymc.io/t/sampling-from-nested-models/7511/8 "2021-05-27T14:07:39Z")

</div>

Not enough variation for this model, it seems. It spans the range `0.2 - 0.9` which is not much for a Poisson GLM with uninformative priors (and / or this dataset size / coefficient ranges).

---

<div class="post-metadata">

**Author:** ![non87](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/non87/32/3859_2.png) [@non87](https://discourse.pymc.io/u/non87)\
**Post date:** [May 27, 2021, 2:10pm UTC](https://discourse.pymc.io/t/sampling-from-nested-models/7511/9 "2021-05-27T14:10:35Z")

</div>

So you think that if I had more variation in the `latent_` variables the model should be fine?

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [May 27, 2021, 2:22pm UTC](https://discourse.pymc.io/t/sampling-from-nested-models/7511/10 "2021-05-27T14:22:26Z")

</div>

I think so. Here I increased the variance of x\_1 by a factor of 4 and it seems to sample just fine: [Google Colaboratory](https://colab.research.google.com/drive/1HlWicirOOBW2R_NbgEK5bhwnu8DL8cBA?usp=sharing)

Now you might be unlucky in that your data actually looks like your original parameters, in which case you might need a simpler model (or much more regularization)

---

<div class="post-metadata">

**Author:** ![non87](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/non87/32/3859_2.png) [@non87](https://discourse.pymc.io/u/non87)\
**Post date:** [May 27, 2021, 2:24pm UTC](https://discourse.pymc.io/t/sampling-from-nested-models/7511/11 "2021-05-27T14:24:19Z")

</div>

You are right that sampling looks much better with more variability. Thanks for all your help with this.