Bookmark and Share

Notice: On March 31, it was announced that Statalist is moving from an email list to a forum. The old list will shut down at the end of May, and its replacement, statalist.org is already up and running.


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

st: Random draws from a negative binomial distribution


From   Dirk Enzmann <dirk.enzmann@uni-hamburg.de>
To   statalist@hsphsun2.harvard.edu
Subject   st: Random draws from a negative binomial distribution
Date   Wed, 19 Jun 2013 00:59:50 +0200

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/


© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   Site index