Generate pseudo-random overdispersed Poisson numbers I have multiple sets of data which conform to overdispersed Poisson distributions which I can model with the alternative parameterization of a negative binomial distribution ($\mu$ and $D$ instead of $p$ and $r$).
I am interested in being able to generate pseudo-random numbers from these distributions but cannot find the ICDF to do so.  I also tried to look for how R did it but it's buried in the Cpp deep enough that I get lost (not a C++ programmer).
How can I convert a psuedo-random uniform number into a pseudo-random negative binomial number with parameters $\mu$=mean and $D$=dispersion?
clarification and thanks
Thank you for the answer and comments pointing me to rnbinom in R.  That in combination with fitdistr from the MASS package do exactly what I want.  However my interest is in generating these random numbers in SQL Server without leveraging a .NET language.
That being said, I found the stata code interesting, but it's parameterized for $p$ and $r$ and perhaps I'm just very dense, but I'm not sure how to convert my $\mu$ and $D$ back to $p$ and $r$.
All that being said, I have found a JSTOR article that about negative binomial distributions with the specific application as an overdispersed Poisson.  I'll see if I can conjure a CDF from there.
 A: R package MASS provides the function rnegbin 
From the documentation:

Description Function to generate random outcomes from a Negative
  Binomial distribution, with mean mu and variance mu + mu^2/theta.
  Usage rnegbin(n, mu = n, theta = stop("’theta’ must be specified"))
  road 131 Arguments n If a scalar, the number of sample values
  required. If a vector, length(n) is the number required and n is used
  as the mean vector if mu is not specified. mu The vector of means.
  Short vectors are recycled. theta Vector of values of the theta
  parameter. Short vectors are recycled. Details The function uses the
  representation of the Negative Binomial distribution as a continuous
  mixture of Poisson distributions with Gamma distributed means. Unlike
  rnbinom the index can be arbitrary. Value Vector of random Negative
  Binomial variate values. Side Effects Changes .Random.seed in the
  usual way. 

Sample code:
 # Negative Binomials with means fitted(fm) and theta = 4.5
 fm <- glm.nb(Days ~ .,data = quine) dummy <- rnegbin(fitted(fm), theta = 4.5)

