# Can I avoid sampling irrelevant bits of the distribution?

**URL:** <https://discourse.pymc.io/t/can-i-avoid-sampling-irrelevant-bits-of-the-distribution/4443>\
**Category:** Questions\
**Created:** [February 2, 2020, 2:46pm UTC](https://discourse.pymc.io/t/can-i-avoid-sampling-irrelevant-bits-of-the-distribution/4443 "2020-02-02T14:46:30Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![TheDerivator](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/thederivator/32/2453_2.png) [@TheDerivator](https://discourse.pymc.io/u/TheDerivator)\
**Post date:** [February 2, 2020, 2:46pm UTC](https://discourse.pymc.io/t/can-i-avoid-sampling-irrelevant-bits-of-the-distribution/4443/1 "2020-02-02T14:46:30Z")

</div>

I discovered pymc3 last month and started playing with it to fit some simple models. For instance, imagine I have:

\alpha, \beta \sim \mathrm{Exp}(-1/10)  
p\_i \sim \mathrm{Beta}(\alpha, \beta)  
N\_i \sim \mathrm{Binom}(ss\_i, p\_i) _(observed data)_  
M\_i \sim \mathrm{Binom}(ps\_i, p\_i)

where i runs from 1\dots N. There could be a story associated with it, such as there are N hospitals with 100 patients each and in each hospital we ask 10 patients if they are happy. We would be interested in statistics of N\_1 + N\_2 + M\_1 + M\_2 because we want to know how many people are happy in all hospitals, or we want to see if we can be more than 90% certain that the overall satisfaction rate is higher than some number.

I found it really easy to build this model with PyMc3 and sample from the posterior. I will also generate some data first to have a working example

```auto
np.random.seed(17610407)
sample_sizes = [10, 20, 30, 40]
n = len(sample_sizes)
population_sizes = [100] * n

alpha_, beta_ = pm.Exponential.dist(0.1).random(size=2)
p_ = pm.Beta.dist(alpha_, beta_).random(size=n)
N_ = pm.Binomial.dist(sample_sizes, p_).random()

with pm.Model() as demo:
    alpha = pm.Exponential('alpha', 0.1)
    beta = pm.Exponential('beta', 0.1)
    p = pm.Beta('p', alpha, beta, shape=n)
    N = pm.Binomial('N', sample_sizes, p, observed=N_)
    M = pm.Binomial('M', population_sizes, p, shape=n) # ?
    
pm.model_to_graphviz(demo)

```

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

But unfortunately, this will sometimes run into problems. I think it is because we are mixing discrete M\_i distributions with continuous p\_i distributions. Namely if I sample from the posterior …

```auto
with demo:
    trace = pm.sample(10000, chains=5, tune=2000, cores=5)
pm.summary(trace)

```

… this gives me very few effective samples for the p’s: around 200. It also complains about many divergences, although I am not sure if these are a problem. However if I remove the M node from the diagram and do everything again, I get about 40k-50k effective samples and far fewer divergences.

So it seems to me that I can get better convergences by first removing M, and sampling from the joint distribution over \alpha, \beta and p. Afterwards, I can compute M for each sample by drawing from the binomial distribution.

- is this a valid approach?
- is this ‘trick’ something pymc3 can do for me if I specify some options to the sampler?

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [February 2, 2020, 3:45pm UTC](https://discourse.pymc.io/t/can-i-avoid-sampling-irrelevant-bits-of-the-distribution/4443/2 "2020-02-02T15:45:27Z")

</div>

> [@TheDerivator](#):
>
> is this a valid approach?

Yes, and it is usually a better approach - NUTS would interact weirdly with discrete samplers sometimes.

One trick is to run samples, and then create a new node for posterior predictive sampling:

```python
with pm.Model() as demo:
    alpha = pm.Exponential('alpha', 0.1)
    beta = pm.Exponential('beta', 0.1)
    p = pm.Beta('p', alpha, beta, shape=n)
    N = pm.Binomial('N', sample_sizes, p, observed=N_)
    trace = pm.sample()
    
with demo:
    M = pm.Binomial('M', population_sizes, p, shape=n)
    posterior_predictive = pm.sample_posterior_predictive(trace, vars=[M])

```