A simple method for generating gamma variables

Abstract

We offer a procedure for generating a gamma variate as the cube of a suitably scaled normal variate. It is fast and simple, assuming one has a fast way to generate normal variables. In brief: generate a normal variate x and a uniform variate U until In(U) < 0.5x2 + d - dv + dln(u), then return dv. Here, the gamma parameter is α ≥ 1, and v = (1 + x/√9d)3, with d = α - 1/3. The efficiency is high, exceeding 0.951, 0.981, 0.992, 0.996 at α = 1, 2, 4, 8. The procedure can be made to run faster by means of a simple squeeze that avoids the two logarithms most of the time: return dv if U < 1 - 0.0331x4. We give a short C program for any α ≥ 1, and show how to boost an α < 1 into an α > 1. The gamma procedure is particularly fast for C implementation if the normal variate is generated in-line, via the #define feature. We include such an inline version, based on our ziggurat method. With it, and an inline uniform generator, gamma variates can be produced in 400MHz CPUs at better than 1.3 million per second, with the parameter a changing from call to call. Categories and Subject Descriptors: G.4 [Mathematics of Computing]: Mathematical Software; 1.6 [Computing Methodologies]: Simulation and Modeling.