I sample my model with two chains and 600 draws each. Then I sample the
posterior predictive, using pm.sample_posterior_predictive() (or its speedy cousin pm.fast_sample_posterior_predictive()). The posterior predictive dataset has one chain and 1200 draws.
Why does sampling the posterior predictive use only one chain instead of the two in the posterior? Maybe I shouldn’t care about this: 1x1200 is maybe equivalent to 2x600. But a subsequent call to pm.summary() on the posterior predictive group makes numerous warnings about the change in shape:
Should these warnings be ignored?
1 Like
Hi David,
I don’t think that’s a problem: that’s probably because sample_posterior_predictive spits out a Multitrace object, which stacks draws and chains, whereas summary, which is actually the ArviZ function, expects an InferenceData object with seperate dimensions for chains and draws.
I’m guessing that if you inject your posterior predictive samples in the idata object where you already have the trace, the warnings should go away – speaking under the control of @OriolAbril here
2 Likes
To have sample_posterior_predictive return a dict whose arrays keep the chain, draw shape, you need to use the keep_size argument, otherwise the chains and draws get flattened as you have already seen. It’s not because only one chain is used however, they are all used but the samples are stored as nchain*ndraw, ... arrays
2 Likes
When I use keep_size=True on pm.fast_sample_size_predictive(), I still
get a posterior dataset with only a single chain rather than two. And now the dataset has 2 draws (instead of 600 it should have, or 1200 it has with keep_size=False), and an extra (unwanted) dimension (of 600) on all my result variables.
When I try pm.summary() on the posterior predictive, it gives me thousands of error messages, instead of the 19 message before:
and so on, for longer than I wanted to wait.
Have I discovered a bug in the keep_size logic of pm.fast_sample_size_predictive()? Or perhaps I am making a different mistake now?
Mmmh, as if there were some mismatch in the dimension attributions in fast_sample_posterior_predictive 
Pinging @rpgoldman here, as he might have some insights
1 Like
Is there a return_inferencedata argument to sample_posterior_predictive? We have had some issues retrieving the original shape due to this internal reshaping when sampling the posterior predictive. Here it looks like it is ignoring the info in the provided trace (or posterior dataset) and independently of the arguments it adds a dummy dimension at the beggining to play the role of chain dimension.
Try getting a dict from sample_posterior_predictive using keep_size and then use from_dict instead of from_pymc3.
2 Likes
Not yet: it seems that there is a PR open to do just that (by @rpgoldman BTW).
Exactly, I have the same impression
1 Like
from_dict(). I tried reading the documentation to understand the difference between it and from_pymc3(), but the site 404ed.
from_dict
You can manually add a dict of ndarrays (assumes dim order as chain,draw,*shape) for each group you want to populate.
This will return a new InferenceData object, which can then be merged with the original InferenceData (one from from_pymc).
You might want to try a few times how it works in Jupyter Notebook (html_repr makes it easy to inspect). <-- this is especially true when trying dims and coords.
2 Likes
See this link (I think we just upgraded the site?)
1 Like
I’m not sure what is going wrong here without an example to run.
What happens in fast_sample_posterior_predictive is that we create a single chain which is a dict of variable names to arrays of shape n, ... where n = cs for c the number of chains and s the length of each chain in the posterior multi-trace (and heaven help you if the chains aren’t all of equal length).
If it is desired to keep size, then for each ary in the dictionary we do
ary.reshape((nchains, ndraws, *ary.shape[1:]))
… which looks like it should do the right thing, assuming that nchains and ndraws are being computed correctly. Again, given an example, I could check this. It would be interesting to know whether there is any difference when you pass an InferenceData or a MultiTrace as the posterior trace to this function.
This is likely a problem of pm.summary() expecting something from the posterior predictive that the latter is not prepared to give. TBH, I didn’t even know that pm.summary() could be used with a posterior predictive – I thought it could only be used on a posterior trace, which has sampling statistics that are not present in the posterior predictive (which is sampled with an entirely different procedure).
@OriolAbril – Should pm.summary() – which is really just az.summary() even work on a posterior predictive distribution? In PyMC3 these are very different things from posterior distributions.
I don’t think we can make more progress without a runnable example. @bridgeland would you see if you can provide one?
Returning inference_data from posterior predictive sampling could be fraught with issues, because PyMC3 uses the same functions (spp and fspp) to predict both in and out of sample.
Or, more accurately, PyMC3 uses these functions to predict in sample, and if you torture the input trace, you can use them to predict out of sample, as well.
I can create an example to show the problem, but not until Wednesday, as I am out of town for a few days, visiting my elderly father in Michigan.
That would be great, thanks. Right now, I’m not able to tell for certain if this is an error in fast_sample_posterior_predictive (and you say you get a related error in sample_posterior_predictive, right?), in pm.summary() (which is really arviz.summary()), or in the way the two interact.
I’ll look forward to seeing your test case.
Best,
r
1 Like
Example:
import pymc3 as pm
import arviz as az
with pm.Model() as model:
foo = pm.Beta('foo', 1.3, 1.4)
bar = pm.Beta('bar', 1 + foo, 2 - foo, shape=50)
trace = pm.sample(draws=600, chains=2, tune=500, return_inferencedata=True)
trace_ppc = pm.fast_sample_posterior_predictive(
trace, var_names=['bar'])
trace_aux = az.from_pymc3(posterior_predictive=trace_ppc)
trace.extend(trace_aux)
with model:
summary = pm.summary(trace, group='posterior_predictive', var_names=['bar'])
summary
Results:
Another example, with keep_size=True:
with pm.Model() as model2:
foo = pm.Beta('foo', 1.3, 1.4)
bar = pm.Beta('bar', 1 + foo, 2 - foo, shape=50)
trace = pm.sample(draws=600, chains=2, tune=500, return_inferencedata=True)
trace_ppc = pm.fast_sample_posterior_predictive(
trace, var_names=['bar'], keep_size=True)
trace_aux = az.from_pymc3(posterior_predictive=trace_ppc)
trace.extend(trace_aux)
bar in the posterior predictive has a bogus shape:
There appear to be 30,000 error messages:
with model2:
summary = pm.summary(trace, group='posterior_predictive', var_names=['bar'])
summary
etc.
I can replicate this error. I will see if I can track it down.
2 Likes
Meanwhile I am suppressing the thousands of warnings via:
import logging
class DisableLogger():
def __enter__(self):
logger = logging.getLogger('arviz.stats.stats_utils')
logger.disabled = True
def __exit__(self, exit_type, exit_value, exit_traceback):
logger = logging.getLogger('arviz.stats.stats_utils')
logger.disabled = False
with DisableLogger():
function_that_produces_warnings()
