Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.
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 |
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/