  Gibbs sampling is a widely used Markov chain Monte Carlo (MCMC) method for
numerically approximating integrals of interest in Bayesian statistics and
other mathematical sciences. Many implementations of MCMC methods do not extend
easily to parallel computing environments, as their inherently sequential
nature incurs a large synchronization cost. In the case study illustrated by
this paper, we show how to do Gibbs sampling in a fully data-parallel manner on
a graphics processing unit, for a large class of exchangeable models that admit
latent variable representations. Our approach takes a systems perspective, with
emphasis placed on efficient use of compute hardware. We demonstrate our method
on a Horseshoe Probit regression model and find that our implementation scales
effectively to thousands of predictors and millions of data points
simultaneously.
