How to use the posterior predictive distribution for checking a model from PyMC

If you grab the posterior predictive samples like this:

ppc = pm.sample_posterior_predictive(idata)

Then you should be able to get the max of each simulated data set like this:

az.extract(ppc.posterior_predictive.max(dim=["likelihood_dim_2"]))["likelihood"]

ppc.posterior_predictive.max(dim=["likelihood_dim_2"]) takes the max over simulated data sets. az.extract() flattens this down into a 1 dimensional array. The ["likelihood"] grabs the variable of interest (i.e., the observed data in your original model). If you prefer to work with a numpy array, you can take the values:

az.extract(ppc.posterior_predictive.max(dim=["likelihood_dim_2"]))["likelihood"].values

If you pass extend_inferencedata=True to pm.sample_posterior_predictive(), then you should be able to do the same with the idata object:

az.extract(idata.posterior_predictive.max(dim=["likelihood_dim_2"]))["likelihood"]