How to generate the posterior predictive distribution for hierarchal model in PYMC3

The below stochastic node y_pred enables me to generate the posterior predictive distribution:

with pm.Model() as model:
alpha = pm.Gamma('alpha', alpha=.1, beta=.1)
mu = pm.Gamma('mu', alpha=.1, beta=.1)
y_pred = pm.NegativeBinomial('y_pred', mu=mu, alpha=alpha)
y_est = pm.NegativeBinomial('y_est', mu=mu, alpha=alpha, observed=messages['time_delay_seconds'].values)

start = pm.find_MAP()
step = pm.Metropolis(start=start)
trace = pm.sample(20000, step, start=start, progressbar=True)


And the posterior predictive distribution plot:

Solution #1: @inversion points out that the shape parameter must be provided.

Problem #2: The samples from the posterior predictive are not correct. They get stuck at 0 or 1. Any ideas how to resolve this? I have only been able to reproduce this issue with hierarchal models.

• I think the traceback is cut off. Can you paste the full one? Also, if you upload a NB that can run stand-alone it would go a long way in enabling others to help you. Aug 29 '15 at 10:33
• Sure thing @twiecki, I've created a sample NB with error reproduced: nbviewer.ipython.org/gist/markdregan/e8376040266a7eb42c34. Thanks. Aug 29 '15 at 17:15
• I've added an answer below, but I'm interested to know why Metropolis(start=start) works in the above code, and fails when I try to use it (unexpected kwd error). pm.__version__ = 3.0 Aug 30 '15 at 20:24

The hierarchical model works if you specify the shape of y_pred.

y_pred = pm.NegativeBinomial('y_pred',
mu=mu[people_idx],
alpha=alpha[people_idx],
shape=people.shape)


Here's the trace:

And the posterior predictive plot (after the appropriate flattening):

y_pred = trace[burn::thin].get_values('y_pred').ravel()


You can, of course, use people_idx to compare posterior predictive plots of specific individuals.

EDIT: Noting that y_pred takes a long time to snap out of a value of 1.

• Note: I worked off of the linked IPython notebook, which is why the distributions appear to be different and the variable names are slightly different. Aug 30 '15 at 20:26
• Thanks @inversion, But the traceplots show issues. Have a look at this ipython notebook. All values for y_est are 1. See: gist.github.com/markdregan/61581d05c14a30e3a6bc Aug 31 '15 at 15:33
• I don't understand exactly why, but y_pred takes a very long time to break out from a value of 1. See my full trace in the edit. Aug 31 '15 at 16:35
• I believe the issue is caused by the priors. Gamma(.1, .1) suggests that alpha and mu are relatively large. When I constrict priors to Uniform(0,100) y_pred becomes unstuck. Sep 1 '15 at 9:29

Second @inversion's answer. Note that there is (as of 5 minutes ago) proper support for drawing random samples directly from random variables. Here is an example of how to do this: http://pymc-devs.github.io/pymc3/posterior_predictive/