[Math] Algorithm to produce random number with a gamma distribution

Question

I'd like to produce pseudo-random numbers with different distributions for a Monte Carlo simulation.

I've got the poisson distribution working nicely with an algorithm from Knuth. I'm having trouble getting a nice easy and fast algorithm for a power distribution. The gamma distribution should do, but the article in wikipaedia gives an algorithm, but remarks that it's not a good one, without providing a link to a better one.

http://en.wikipedia.org/wiki/Gamma_distribution#Generating_gamma-distributed_random_variables

Is there a good, fast algorithm for a gamma distribution?

Answer 1

The difficulty mentioned in Wikipedia refers to gamma distributions with small shape parameter; this has been addressed in arXiv:1302.1884:

The gamma distribution with small shape parameter can be difficult to characterize. For this reason, standard algorithms for sampling from such a distribution are often unsatisfactory. In this paper, we first obtain a limiting distribution for a suitably normalized gamma distribution when the shape parameter tends to zero. Then this limiting distribution provides insight to the construction of a new, simple, and highly efficient acceptance--rejection algorithm. Pseudo-code and an R implementation for this new sampling algorithm are provided.