Add sample_generative function

[ColCarroll]

This replaces #2876, since @springcoil's recent improvements to draw_values (#2979) makes this somewhat easier. Includes a test, plus two more I had written for draw_values. Both of them actually failed without changes here.

This method is quite fast now.

[springcoil]

LGTM - needs an addition to the RELEASE NOTES though.

[springcoil]

Not for this PR - but a future PR would include some benchmarks I think. Other than that awesome work @ColCarroll

[junpenglao]

forward_val is the numpy version of the forward function and should be faster.

[junpenglao]

Should use .transformation.forward here.

[junpenglao]

Oops i meant .transformation.forward_val should be used here... Should be slightly faster:
model[untransformed].transformation.forward(data[untransformed]).eval() --> model[untransformed].transformation.forward_val(data[untransformed])

[junpenglao]

the size args is not working yet:

from pymc3.distributions.distribution import draw_values

with pm.Model() as m:
    p = pm.Beta('p', 1., 1.)
draw_values([m['p']], size=10)
[array(0.72159403)]

[ColCarroll]

your previous example is still a little funny too:

X = theano.shared(np.arange(3))
with pm.Model() as m:
    ind = pm.Categorical('i', np.ones(3)/3)
    x = pm.Deterministic('X', X[ind])
    prior=pm.sample_generative()

[junpenglao]

I am really good at discovering edge cases.

[twiecki]

should we rename this sample_prior_predictive?

[twiecki]

and sample_ppc should then probably be sample_posterior_predictive (it's not really doing a check).

[ColCarroll]

I think this is now actually the generative function (see the test case). The prior predictive would be nice to also implement.

[twiecki]

I looked at the test but I don't understand the difference between generative and prior predictive here.

[ColCarroll]

Ugh... I think I still don't. Will work out a simple model by hand later to try to help my intuition.

[springcoil]

I agree with the renaming

…

On Sat, 19 May 2018, 2:43 pm Thomas Wiecki, ***@***.***> wrote: ***@***.**** commented on this pull request. ------------------------------ In pymc3/sampling.py <#2983 (comment)>: > @@ -1208,6 +1209,49 @@ def sample_ppc_w(traces, samples=None, models=None, weights=None, return {k: np.asarray(v) for k, v in ppc.items()} +def sample_generative(samples=500, model=None, vars=None, random_seed=None): and sample_ppc should then probably be sample_posterior_predictive (it's not really doing a check). — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#2983 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AA8DiLV0bYDwo4DMtt3WzNFjxvzTuEHOks5t0BNygaJpZM4UFdN9> .

[ferrine]

I think that isinstance check should me replaced with checker for name

[twiecki]

That seems a bit crude, couldn't we just test for the mean and sd?

[AllenDowney]

I am coming to this discussion late, but I was talking with Chris Fonnesbeck about this issue and he pointed me to this PR.

This functionality would be very helpful to me. I am working on a new volume of Think Bayes that will use PyMC, and I plan to demonstrate a development process that uses a pm.Model to generate a prior predictive distribution before running an inference. I think this is useful for both validation (model makes sense) and verification (model specification is correct).

Especially for verification, it would be ideal if generation and inference are as identical as possible. For example, I would love to be able to run

with pm.Model() as model:
    mu = pm.Gamma('mu', alpha, beta)
    goals = pm.Poisson('goals', mu)
    trace = pm.sample(1000)

To generate the prior predictives, and then

with pm.Model() as model:
    mu = pm.Gamma('mu', alpha, beta)
    goals = pm.Poisson('goals', mu, observed=6)
    trace = pm.sample(1000)

If the first one validates, it would be really hard to mess up the second one.

This API would require pm.sample to analyze the model, see that there is no data, and run sample_prior_predictive rather than any of the MCMC samplers. I assume that would not be hard, but I don't really know.

This notebook contains more examples of what I have in mind.

https://github.com/AllenDowney/BayesMadeSimple/blob/master/zigzag.ipynb

In those examples I use pm.sample to generate, but (you will not be surprised to hear) it doesn't work very well. @ColCarroll , since you saw the talk where I presented this notebook, you might have thoughts to add.

Thank you all!

[springcoil]

Fascinating idea!!!

…

On Fri, 1 Jun 2018, 6:40 pm Allen Downey, ***@***.***> wrote: I am coming to this discussion late, but I was talking with Chris Fonnesbeck about this issue and he pointed me to this PR. This functionality would be very helpful to me. I am working on a new volume of Think Bayes that will use PyMC, and I plan to demonstrate a development process that uses a pm.Model to generate a prior predictive distribution before running an inference. I think this is useful for both validation (model makes sense) and verification (model specification is correct). Especially for verification, it would be ideal if generation and inference are as identical as possible. For example, I would love to be able to run with pm.Model() as model: mu = pm.Gamma('mu', alpha, beta) goals = pm.Poisson('goals', mu) trace = pm.sample(1000) To generate the prior predictives, and then with pm.Model() as model: mu = pm.Gamma('mu', alpha, beta) goals = pm.Poisson('goals', mu, observed=6) trace = pm.sample(1000) If the first one validates, it would be really hard to mess up the second one. This API would require pm.sample to analyze the model, see that there is no data, and run sample_prior_predictive rather than any of the MCMC samplers. I assume that would not be hard, but I don't really know. This notebook contains more examples of what I have in mind. https://github.com/AllenDowney/BayesMadeSimple/blob/master/zigzag.ipynb In those examples I use pm.sample to generate, but (you will not be surprised to hear) it doesn't work very well. @ColCarroll <https://github.com/ColCarroll> , since you saw the talk where I presented this notebook, you might have thoughts to add. Thank you all! — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#2983 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AA8DiOoNyaRuR_1nXGUanaVFfORluf86ks5t4W52gaJpZM4UFdN9> .

[fonnesbeck]

So I wonder if sample could infer that there are no observed RVs and just do forward sampling. Can anyone think of a reason for doing MCMC without any likelihoods?

[kyleabeauchamp]

I think there are lots of situations where we want to sample from a joint (unobserved) distribution that we describe using the machinery of pymc---pymc allows one to describe models in terms of random variables and functions of them, which is the right level of abstraction for that kind of work.

One example: Monte Carlo based probabilistic sensitivity analysis of an economic model.

[kyleabeauchamp]

Another example: for performance reasons, one may sometimes apply Bayes' rule analytically with a known conjugate prior. E.g., perhaps it's computationally faster to model the posterior RVs directly, rather than include the prior and observations. In such cases, the model won't have observed / likelihood nodes.

[fonnesbeck]

Right, but those situations don’t require MCMC, do they? Just forward simulation using independent sampling.

[kyleabeauchamp]

Agreed, but my point is that such folks ought to use pymc rather than write their own model description / independent sampling code.

[fonnesbeck]

Yeah I think we are talking about the same thing. Following on Allen’s point, we could engineer sample such that it uses simple forward sampling when there are no likelihoods in the model, without the user having to do anything. My question was related to whether you would ever not want to do this.

[ColCarroll]

@AllenDowney I actually started using your notebook as a test case on the train home last night!

I think the issues brought up here with size and shape may be too hard to fix in a single go, and I think making the sort of analysis you did "just work" would be a good start (#2984 is something of a blocker), even if some edge-cases remain.

My goal for the api is more like

with pm.Model() as model:
    mu = pm.Gamma('mu', alpha, beta)
    goals = pm.Poisson('goals', mu, observed=6)

with model:
    prior = pm.sample_prior()

with model:
    trace = pm.sample()

with model:
    posterior_pred = pm.sample_ppc(trace)

which allows for model reuse. There's no reason your approach could not also be supported, though it makes me a little nervous to overload a function like that: in particular, silently ignoring arguments like tune or step.

[AllenDowney]

@ColCarroll

That makes sense: the parameters of sample_prior, sample, and sample_ppc are different, so it makes sense for them to be different functions. And I like the idea of separating the specification of the model from the sampling functions.

To be completely systematic, you might want four functions:

1. sample_prior
2. sample_prior_predictive
3. sample_posterior (synonym for sample)
4. sample_posterior_predictive (synonym for sample_ppc)

But maybe (1) is not necessary; if you're going to run the model forward, you might as well generate both the priors and the prior predictives.

If you like shorter names, maybe sample_prior_pred, sample_post, and sample_post_pred.

[twiecki]

I like this naming scheme. For completeness we should indeed add pm.sample_posterior and then make pm.sample() an alias for that function for backwards compatibility.

[springcoil]

I'm at a talk by Mike Betancourt tonight. Is this implemented yet.

[ColCarroll]

Have been making surprisingly steady progress. I'll push a more complete description of the changes when it is ready (tomorrow?), but will push something partial tonight.

Has turned into a partial rewrite of all shape handling, which runs through generate_samples. I think what I have is much cleaner now, and ironically most of what is left is cleaning up all the code that deals with bad behavior in the old/current implementation. See, for example, Multinomial._random.

I do have it so that @AllenDowney's notebook runs cleanly without "abusing" PyMC3 😄 !

[AllenDowney]

That's awesome. Thank you, Colin!

…

On Thu, Jun 14, 2018 at 3:26 PM, Colin ***@***.***> wrote: Have been making surprisingly steady progress. I'll push a more complete description of the changes when it is ready (tomorrow?), but will push something partial tonight. Has turned into a partial rewrite of all shape handling, which runs through generate_samples. I think what I have is much cleaner now, and ironically most of what is left is cleaning up all the code that deals with bad behavior in the old/current implementation. See, for example, Multinomial._random. I do have it so that @AllenDowney <https://github.com/AllenDowney>'s notebook runs cleanly without "abusing" PyMC3 😄 ! — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#2983 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABy37cTXRb_cYoyamShDFpCK5_hzX02rks5t8rjfgaJpZM4UFdN9> .

[springcoil]

Wicked!!!

…

On Thu, 14 Jun 2018, 8:26 pm Colin, ***@***.***> wrote: Have been making surprisingly steady progress. I'll push a more complete description of the changes when it is ready (tomorrow?), but will push something partial tonight. Has turned into a partial rewrite of all shape handling, which runs through generate_samples. I think what I have is much cleaner now, and ironically most of what is left is cleaning up all the code that deals with bad behavior in the old/current implementation. See, for example, Multinomial._random. I do have it so that @AllenDowney <https://github.com/AllenDowney>'s notebook runs cleanly without "abusing" PyMC3 😄 ! — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#2983 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AA8DiMQsRwcyib1Fi62pUCQK8ykK89Gdks5t8rjagaJpZM4UFdN9> .

[ColCarroll]

Oddly, I think this is nearly ready to merge. There were a few unexpected tests failing, and I think a few distributions (notably, Dirchlet) are not tested as well as they should be. In any case, I expect some edge cases to pop up after merging.

Everything is quite fast for .random and .draw_values (see the screen shot below). Indeed, this may now be the most flexible Python library for generating random numbers from a network of variables. I am a little worried that I accidentally slowed down a hot path for actual MCMC sampling, though, so checking the http://pandas.pydata.org/speed/pymc3/ after merging may also be good.

I do have a modestly embarrassing question about nomenclature, and what the prior vs prior predictive is. In the example below, would I be able to define

prior = {k: v for k, v in generated.items() if k != 'obs'}
prior_predictive = {'obs': generated['obs']}

[junpenglao]

My understanding is that prior predictive will return:

prior_predictive = {k: v for k, v in generated.items()}

The current output from the function is what I would expect.

[twiecki]

Found some small things. I'm very excited about having this functionality. It definitely needs a release note as well as an example, though.

[ColCarroll]

Thanks for those! I'd prefer to add an example in a separate PR, if that is ok. This touches so much code in subtle ways...

[AllenDowney]

Thank you all for your work on this, especially @ColCarroll! Do you know when this will make it into a release (or has it already)?

…

On Mon, Jun 18, 2018 at 1:01 AM, Junpeng Lao ***@***.***> wrote: Merged #2983 <#2983>. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#2983 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABy37S5y87VqK85Pt936ijxQwJPs98KJks5t9zQkgaJpZM4UFdN9> .

[junpenglao]

It will be in the next release. @fonnesbeck @twiecki maybe we should prepare a release actually, we got quite some new features.

[ColCarroll]

I've been using this function and the graph function daily-ish and ironing out the bugs. Which are shape related, obviously.

[twiecki]

Are there any non-shape bugs?

The answer to should we do a new release is always yes.

[fonnesbeck]

OK, I will get a release ready over the weekend.