# Need code for sampling negative binomial with non-integer $r$

I'm trying to write code (in C) to sample from a negative binomial distribution parameterized by $r,p$, where $r$ is not necessarily integer (also called Polya distribution).

I've found a number of explanations, and it looks like I could sample $y$ from a Gamma(r,1) distribution, transform $y$ somehow, and then use it to sample from a Poisson distribution, but that means I need to find a Gamma distribution with non-integer $r$.

All of this is a bit scary, and I was wondering if anyone could point me to some source code to do it. I don't much care what language. I have access to gamma functions, which I'm pretty sure will be needed.

• This is a general comment (and perhaps not helpful), but did you consider sampling from uniform(0,1) and using an inverse probability integral transform (like here)? – Richard Hardy Apr 11 '15 at 6:19
• @RichardHardy: I was hoping for something faster. To use the inverse-transformation, you need to construct the CDF for the specific $r,p$ arguments, and then use some sort of simple search function to invert it. Nevermind, I managed to find info for the Gamma-Poisson method and get it working. Perhaps I will answer my own question. The tricky part is sampling the Gamma for non-integer shape. In Wikipedia I found something called Ahrens-Dieter for that. – Mike Dunlavey Apr 11 '15 at 13:12

The code for sampling a negative binomial is in rnbinom.c. It does it in two steps. First it samples from the gamma distribution in rgamma.c. That number is used as the mean of a poisson, defined in rpois.c.