themselves in parallel. The sampling steps of our method can be run in their entirety in parallel. Another attractive feature of our method is that the model is fully specified by the log poste- rior, and possibly its gradient and Hessian. The tools that we discuss in Section 4 are components of a reusable infrastructure. Only the function that returns the log posterior changes from model to model. This feature is unlike a blockwise Gibbs sampler, for which we need to derive and imple- ment a set of conditional densities for each model. A small change in the model specification can require a substantial change in the sampling strategy. For example, a change to a prior might mean that the sampler is no longer conditionally conjugate. So while it might be possible to construct a highly efficient Gibbs sampler for a particular model (e.g., through blocking, data augmentation, or parameter transformation), there can be considerable upfront investment in doing so. We now provide three examples of our method in action. In the first example, we simulate data from a basic, conditionally conjugate model, and compare the marginal posterior distributions that were generated by our method with those from a Gibbs sampler. The second example is a non-conjugate hierarchical model of store-level retail sales. In that example, we compare esti- mates from our method with those from Hamiltonian Monte Carlo (HMC), using the Stan software package (Stan Development Team 2014). For these first two examples, the MCMC methods are ef- ficient, so it is likely that they generate good estimates of the true posterior densities. Thus, we can use those estimates as benchmarks against which we can assess the accuracy of our method. The third example is a more complicated model for which MCMC performs poorly. We use this example to not only assess the accuracy of our method (for those parameters for which we think the MCMC estimates are reasonable), but to illustrate some of the computational problems that are inherent in MCMC methods. In all of our examples, we implemented our method using the bayesGDS package (Braun 2013a) for R . In our first example, we simulated T = 10 observations for each of N = 1500 heterogeneous units. Each observation for unit i = 1 . . . N is a sample from a normal distribution, with mean q i 12 3.1 Simulated data, conditionally conjugate model 3 Examples