Efficient and Accurate Parallel Inversion of the Gamma Distribution

Abstract

A method for parallel inversion of the gamma distribution is described. This is very desirable for random number generation in Monte Carlo simulations where gamma variates are required. Let $\alpha$ be a fixed but arbitrary positive real number. Explicitly, given a list of uniformly distributed random numbers our algorithm applies the quantile function (inverse CDF) of the gamma distribution with shape parameter $\alpha$ to each element. The result is, therefore, a list of random numbers distributed according to the said distribution. The output of our algorithm has accuracy close to a choice of single- or double-precision machine epsilon. Inversion of the gamma distribution is traditionally accomplished using some form of root finding. This is known to be computationally expensive. Our algorithm departs from this paradigm by using an initialization phase to construct, on the fly, a piecewise Chebyshev polynomial approximation to a transformation function, which can be evaluated very quickly during variate generation. The Chebyshev polynomials are high order, for good accuracy, and generated via recurrence relations derived from nonlinear second order ODEs. A novelty of our approach is that the same change of variable is applied to each uniform random number before evaluating the transformation function. This is particularly amenable to implementation on SIMD architectures, whose performance is sensitive to frequently diverging execution flows due to conditional statements (branch divergence). We show the performance of a CUDA GPU implementation of our algorithm (called Quantus) is within an order of magnitude of the time to compute the normal quantile function.