Just to give more context, this is for calculating the Bayesian posterior predictive p-value (here’s one of the first papers about it, but its use has been advocated by Gelman among others). The idea is to use it as a measure of how well a posterior fits the observed data, irrespective of the other models, rather than using it as a decision criterion. Basically, we choose some statistic of the data, get a trace, get simulated data for each trace sample, and then check for each trace sample whether the statistic of the simulated data is at least as extreme as that of the actual data (under that trace sample). That is why I would like to get for each posterior sample the loglikelihood of data simulated from the sample. However, the method above with pm.set_data is so slow!