Bookmark and Share

Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.


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

st: Random draws from a negative binomial distribution


From   Dirk Enzmann <[email protected]>
To   [email protected]
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: [email protected]
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–2018 StataCorp LLC   |   Terms of use   |   Privacy   |   Contact us   |   Site index