grow with the size of the dataset. Also, if the data represent outcomes of multiple interdepen- dent processes (such as the timing and magnitude of purchases), both the posterior parameters and successive MCMC samples tend to be correlated, requiring a larger, yet unknown, number of iterations. We believe that the most important reason Bayesian methods have not been embraced “in the field” nearly as much as classical approaches is that they are difficult and expensive to implement routinely using MCMC, even with semi-automated software procedures. Practitioners simply do not have an academician’s luxury of letting an MCMC chain run for days or weeks with no guarantee that the chain has converged to produce “correct” answers at the end of the process. 

With recent developments in multiple core processing and distributed computing systems, it is reasonable to look to parallel computing technology as a solution to the convergence problem. However, each MCMC cycle depends on the outcome of the previous one, so we cannot collect posterior samples in parallel by allocating the work across distributed processing cores. Using parallel processors to generate one draw from a target distribution, or running several MCMC chains in parallel, is not the same as generating all the required independent samples in parallel. Hence, extant parallel MCMC methods are also subject to the same pesky question of convergence; indeed, now one has to ensure that all of the parallel chains have converged. On the other hand, non-MCMC methods like rejection sampling have the advantage of being able to generate samples from the correct target posterior in parallel, but these methods are beset with their own set of implementation issues. For instance, the inability to find efficient “envelope” distributions renders standard rejection sampling almost impractical for all but the smallest problems. 

In this paper, we propose a solution to sample from Bayesian parametric, hierarchical models that is inspired by two pre-MCMC approaches: rejection sampling, and sampling from a multi- variate normal (MVN) approximation around the posterior mode. Our contribution is an algo- rithm that recasts traditional rejection sampling in a way that circumvents the difficulties associ- ated with these two approaches. The algorithm requires that one be able to compute the unnor- malized log posterior of the parameters (or a good approximation of it); that the posterior distri- bution is bounded from above over the parameter space; and that available computing resources can locate any local maxima of the log posterior. There is no need to derive conditional posterior distributions (as with blockwise Gibbs sampling), and there are no conjugacy requirements. 

We present the details of our method in Section 2, and in Section 3, we share some examples 

