Simulating draws from a negative binomial I am looking for a way to simulate draws from a negative binomial distribution for a computational experiment on biological sequencing data. I am using a high performance package which only has certain distributions however, and though I know that gamma+poisson draws would give me the required simulation, the package lacks the latter. Here is a list of the distributions I have: [uniform, triangular, beta, exponential, gaussian, lognormal, vonmises, weibull, pareto, gamma]. Is there a simple approximation I can use from any of these functions that would require the least computational complexity (i.e. dont loop over creating many random draws for a single sample).
 A: One possible way is to generate a large number of value with frequencies their probabilities and then sample from a discrete uniform choice without replacement. By example, suppose your probability of success is $p=1/2$ and the number of failure until the experience has stop is $r=5$.
Then you start by computing the probability of each value using the negative binomial distribution. For the example, you will have $$P(X=0) = 0.03125, P(X=1) = 0.07813, P(X=2)=0.11719, \dots$$
Now, you create a large list of frequency corresponding to those probability. If you want to create a vector of length 100000, then this vector will have 3125 times $0$, 7813 times $1$; 11719 times $2$ and so one.
After, you just choose without replacement from that list the number $n$ of trials you want.
A: Here is a screenshot of a reference by NL Johnson et al. (2005) summarising the possible approaches based on a Geometric sampling:

Devroye (1985) indeed does not cover any other algorithm than the Gamma-Poisson version. Note that generating a Poisson variate is feasible from the Exponential generator, due to the connection between the Poisson process and the Exponential distribution. It is however not the recommended algorithm.
A: Why don't you use the inverse CDF method? The CDF of the negative binomial is the regularized incomplete beta function, which is just a ratio of beta functions so it should not be a problem to obtain it with your setup. You can see the detailshere and here. If you want something more fancy tailored specifically to the negative binomial distribution, the standard reference is Devroye's book Non-uniform random variate generation.
