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

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

 From Dirk Enzmann To statalist@hsphsun2.harvard.edu Subject st: Re: Random draws from a negative binomial distribution Date Wed, 19 Jun 2013 01:06:01 +0200

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