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

# st: RE: RE: Re: Random draws from a negative binomial distribution

 From Timothy Mak <[email protected]> To "[email protected]" <[email protected]> Subject st: RE: RE: Re: Random draws from a negative binomial distribution Date Wed, 19 Jun 2013 14:44:04 +0800

```<>

This is very interesting. Usually there is a good reason why Stata doesn't allow you to do something, but I haven't been able to find why in this case. I've looked at the help file for -rpoisson- and also -rpois- in R. The Stata algorithm appears to be more advanced than that of R (Stata 12 vs R2.15).

It would be good if StataCorp can comment on this.

I've also found a blog that gives fairly good details on the intricacies of simulating Poisson variables.

http://evolvedmicrobe.com/blogs/?p=6http://evolvedmicrobe.com/blogs/?p=6

There's no mention of difficulties with small lambda's...

Tim

-----Original Message-----
From: [email protected] [mailto:[email protected]] On Behalf Of Kieran McCaul
Sent: 19 June 2013 12:28
To: [email protected]
Subject: st: RE: Re: Random draws from a negative binomial distribution

...
The problem is with rpoisson().
the  smallest mu that it will accept is 1e-6

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'

* indirectly via -rgamma- and -rpoisson-:
gen double xg = rgamma(`size',`scale')

replace xg = 1e-6 if xg< 1e-6
gen double xnb = rpoisson(xg)

nbreg xnb, irr
di "size = " 1/e(alpha) ", prob = " ///
1/e(alpha)/(1/e(alpha)+exp(_b[_cons]))

-----Original Message-----
From: [email protected] [mailto:[email protected]] On Behalf Of Dirk Enzmann
Sent: Wednesday, 19 June 2013 7:06 AM
To: [email protected]
Subject: st: Re: Random draws from a negative binomial distribution

Because there was a mistake in the Stata example syntax of my previous
mail, here again my question (please ignore my previous mail):

=================================================================

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'

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

* 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]))

* ---------------------------------------
* 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'

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

* 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]))

* --- 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)
fitx = fitdistr(x,densfun="negative binomial")
fitx

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

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

Dirk

=================================================================

Am 19.06.2013 00:59, schrieb Dirk Enzmann:
> Unfortunately, I am not able to solve the following problem in Stata
> which I can solve easily using R: ...

========================================
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/

*
*   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/

*
*   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/
```