5
$\begingroup$

See iPython notebook for full example

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:

enter image description here

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.

See ipython notebook with issue reproduced.

enter image description here

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – twiecki
    Aug 29 '15 at 10:33
  • $\begingroup$ Sure thing @twiecki, I've created a sample NB with error reproduced: nbviewer.ipython.org/gist/markdregan/e8376040266a7eb42c34. Thanks. $\endgroup$
    – Mark Regan
    Aug 29 '15 at 17:15
  • $\begingroup$ 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 $\endgroup$
    – inversion
    Aug 30 '15 at 20:24
3
$\begingroup$

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:

enter image description here

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

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

enter image description here

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.

enter image description here

$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – inversion
    Aug 30 '15 at 20:26
  • $\begingroup$ 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 $\endgroup$
    – Mark Regan
    Aug 31 '15 at 15:33
  • $\begingroup$ 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. $\endgroup$
    – inversion
    Aug 31 '15 at 16:35
  • $\begingroup$ 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. $\endgroup$
    – Mark Regan
    Sep 1 '15 at 9:29
3
$\begingroup$

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/

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.