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]

From |
Timothy Mak <tshmak@hku.hk> |

To |
"statalist@hsphsun2.harvard.edu" <statalist@hsphsun2.harvard.edu> |

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

**References**:**st: Random draws from a negative binomial distribution***From:*Dirk Enzmann <dirk.enzmann@uni-hamburg.de>

**st: Re: Random draws from a negative binomial distribution***From:*Dirk Enzmann <dirk.enzmann@uni-hamburg.de>

**st: RE: Re: Random draws from a negative binomial distribution***From:*Kieran McCaul <kieran.mccaul@uwa.edu.au>

- Prev by Date:
**st: Variable name length limit** - Next by Date:
**st: RE: Variable name length limit** - Previous by thread:
**st: RE: Re: Random draws from a negative binomial distribution** - Next by thread:
**st: Re: Random draws from a negative binomial distribution** - Index(es):