GPU-accelerated hierarchical Bayesian estimation of luminosity functions using flux-limited observations with photometric noise

We describe a C++ framework that uses graphical processing units (GPUs) to accelerate basic hierarchical Bayesian computation, and use it to compute the posterior distribution for the parameters of a galaxy luminosity function based on data with photometric noise and selection effects arising from a flux-limited detection threshold. In addition to estimating the population-level luminosity function parameters, the framework also provides estimates of the absolute magnitude of each object, which automatically correct for Eddington bias. To sample the posterior, we implement a Metropolis–Hastings-within-Gibbs algorithm that alternates between exploring the population-level and member-level parameters. The algorithm exploits conditional independence in hierarchical Bayesian models, which makes member-level exploration readily parallelizable on a GPU. The framework uses adaptive MCMC, automatically tuning Metropolis–Hastings proposal distributions on-the-fly. We present a simulation study demonstrating the accuracy and computational scaling of the algorithm. In addition, we compare hierarchical Bayesian estimation with maximum likelihood estimation (known to provide inconsistent estimates in this setting), providing a concrete demonstration of the benefits of using hierarchical models for luminosity function inference.