st: Random draws from a negative binomial distribution

[Dirk Enzmann <dirk.enzmann@uni-hamburg.de>]

Unfortunately, I am not able to solve the following problem in Stata 
which I can solve easily using R:

As far as I can see Stata does not allow to draw random values from a 
negative binomial distribution if "size" (= 1/alpha) is less than 0.1 
(see -h rnbinomial-). I tried to circumvent this problem by (1) creating 
random draws from a gamma distribution with shape parameter = size and 
scale parameter = (1-prob)/prob, with prob = size/(size+mu), and 
subsequently creating random draws from a poisson distribution with 
parameter m = the result of the previous random draws from the gamma 
distribution. However, if size is small, this does not help either.

Here an example which works, followed by an example which does not:

* ---- begin Stata example -------------
* a) The following works because size > 0.1:

clear
input x freq
0 9316
1  601
2   61
3   15
4    5
5    1
7    1
end

expand freq
nbreg x, irr
local mu = exp(_b[_cons])
local size = 1/e(alpha)
local prob = `size'/(`size'+`mu')
local scale = (1-`prob')/`prob'

* directly via -rnbinomial-:
gen xrnbinom = rnbinomial(`size',`prob')
nbreg xrnbinom, irr
di "size = " 1/e(alpha) ", prob = " ///
   1/e(alpha)/(1/e(alpha)+exp(_b[_cons]))

* indirectly via -rgamma- and -rpoisson-:
gen xg = rgamma(`size',`scale')
gen xnb = poisson(xg)
nbreg xrnbinom, irr
di "size = " 1/e(alpha) ", prob = " ///
   1/e(alpha)/(1/e(alpha)+exp(_b[_cons]))

* ---------------------------------------
* b) The following does not work because size < 0.1:

clear
input x freq
0 2041
1   79
2   22
3   13
4    5
6    1
7    1
8    1
10   1
13   1
end

expand freq
nbreg x, irr
local mu = exp(_b[_cons])
local size = 1/e(alpha)
local prob = `size'/(`size'+`mu')
local scale = (1-`prob')/`prob'

* directly via -rnbinomial-:
gen xrnbinom = rnbinomial(`size',`prob')
nbreg xrnbinom, irr
di "size = " 1/e(alpha) ", prob = " ///
   1/e(alpha)/(1/e(alpha)+exp(_b[_cons]))

* indirectly via -rgamma- and -rpoisson-:
gen xg = rgamma(`size',`scale')
gen xnb = poisson(xg)
nbreg xrnbinom, irr
di "size = " 1/e(alpha) ", prob = " ///
   1/e(alpha)/(1/e(alpha)+exp(_b[_cons]))

* --- End Stata example. --------------------

If this were possible I could use Stata for analyses of count data, if 
not I have to switch to R which I am trying to avoid for consistency 
reasons.

# --- Begin R example: ---------------------
# b)

library(MASS)

x = rep(c(0,1,2,3,4,6,7,8,10,13),c(2041,79,22,13,5,1,1,1,1,1))
table(x)
fit = fitdistr(x,densfun="negative binomial")
fit

xnb = rnbinom(length(x),size=fit$estimate[1],mu=fit$estimate[2])
table(xnb)
fitdistr(xnb,densfun="negative binomial")
fit

# --- End R example. ---------------------


Dirk

========================================
Dr. Dirk Enzmann
Institute of Criminal Sciences
Dept. of Criminology
Rothenbaumchaussee 33
D-20148 Hamburg
Germany

phone: +49-(0)40-42838.7498 (office)
       +49-(0)40-42838.4591 (Mrs Billon)
fax:   +49-(0)40-42838.2344
email: dirk.enzmann@uni-hamburg.de
http://www2.jura.uni-hamburg.de/instkrim/kriminologie/Mitarbeiter/Enzmann/Enzmann.html
========================================
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/