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

[Timothy Mak <tshmak@hku.hk>]

<>

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: owner-statalist@hsphsun2.harvard.edu [mailto:owner-statalist@hsphsun2.harvard.edu] On Behalf Of Kieran McCaul
Sent: 19 June 2013 12:28
To: statalist@hsphsun2.harvard.edu
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: owner-statalist@hsphsun2.harvard.edu [mailto:owner-statalist@hsphsun2.harvard.edu] On Behalf Of Dirk Enzmann
Sent: Wednesday, 19 June 2013 7:06 AM
To: statalist@hsphsun2.harvard.edu
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: 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/

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