Cheng Li, Sanvesh Srivastava, and David Dunson write:

We propose a new scalable algorithm for posterior interval estimation. Our algorithm first runs Markov chain Monte Carlo or any alternative posterior sampling algorithm in parallel for each subset posterior, with the subset posteriors proportional to the prior multiplied by the subset likelihood raised to the full data sample size divided by the subset sample size. To obtain an accurate estimate of a posterior quantile for any one-dimensional functional of interest, we simply calculate the quantile estimates in parallel for each subset posterior and then average these estimates.

Wow—does this really work? This could be awesome. Definitely worth trying out in Stan.

Looks like they’re using Stan (actually, “STAN” :Op) already; see pages 11 & 35/36.

I’d be happy to evaluate this for a variety of models in Stan and report back the results if someone could help me discern what precisely needs to be done in the stan code for each subset. So far as I can tell one needs to use increment_log_prob(p) (or the newere “target += p” syntax), but I believe that p needs to be multiplied by a scaling factor that is presumably a function of the proportion of the total data that is present in the subset?

I think I’ve done this by accident before? I.e. when learning how to parallelise MCMCglmm runs, I accidentally drew new a new subset for each parallel chain. I did not immediately notice my mistake, because the chains had converged fairly well.