Help the Stat Consulting Group by giving a gift

Simulating Discrete (Geometric, Poisson and Zero-Inflated Poisson, Negative Binomial and Zero-Inflated Negative Binomial) Random Variables

Geometric Random Variable:

It can be shown that a Geometric random variable can be
simulated using the following argument (int(ln(u)/ln(1-p)) + 1)
where **u** is a uniform(0,1) random variable and **p** is the probability
of observing a success (Simulation by Ross, 2003). In this example we are going to generate
a Geometric
random variable with 1000 observations with probability of success **p** = 0.25.

clear set obs 1000 local p = .25 gen u = uniform( ) gen g = int(ln(u)/ln(1-`p')) + 1

Negative Binomial Random Variable as a sum of independent Geometric Random
Variables:

To generate a Negative Binomial random variable we make use
of the fact that a Negative Binomial random variable is sum of **r**
independent Geometric random variables, where **r** is the of trials required
to observe the r^{th} success and **p** is the probability of a success. In this
example, we are going to simulate **p** = 0.25 and let **r** = 4. Note
that unlike in the geometric case, the uniform(0,1) is defined within the **
foreach** argument so that the r trials are independent.

clear set obs 1000 local r = 4 local p = .25 foreach var of newlist nb1-nb`r' { gen `var' = int(ln(uniform( ))/ln(1-`p')) + 1 } egen nb = rsum(nb1-nb`r') drop nb1-nb`r'

Poisson and Zero-Inflated Poisson Random Variable:

We can simulate a Poisson and Zero-Inflated Poisson random variable
by recoding a uniform(0,1) random variable in terms of the cumulative
distribution. What allows us to simulate both types of variables in the same
code fragment is when we specify **p0**--the percent of zeros not explained
by a Poisson distribution--to be 0, the Zero-Inflated Poisson reduces to the
traditional Poisson distribution. The **while** loop continues until the
cumulative probability exceeds 0.99999 (a problem may exist if a value from the
uniform random variable is greater than 0.99999). What is required for this code
fragment is the specification of the Poisson parameter **lambda** and **p0.
**We need to note that the program was stable for values of
lambda less than 100. We recommend checking that

clear set obs 1000 gen x = uniform( ) gen u = x local lambda = 10 local p0 = .40 local cp=0 local n =0 while `cp'<=.99999 { local p = `p0'*(`n'==0) + (1-`p0')*exp(-1*`lambda')*(`lambda'^`n')/exp(lnfact(`n')) local cp= `cp'+`p' local pcp = `cp'-`p' quietly recode x (`pcp'/`cp' = `n') local n = `n' + 1 }

Negative Binomial and Zero-Inflated Negative Binomial Random Variable that allows for over dispersion.

The hallmark of the Poisson distribution is that the mean is
equal to the variance. Many times that assumption is not satisfied and the
variance is greater than the mean. To account for the over-dispersion, we can
use a Negative
Binomial as a Gamma mixture of Poisson random variable that accounts for
over-dispersion by adding a parameter alpha. Like the prior code fragment,
we can simulate both types of variables within the same code fragment (the prior
code was more stable when we specified extreme values of **lambda** and **p0**,
this one is not. We recommend using conservative values of the parameter and
checking that **x** was recoded for all values). When alpha is equal to zero, we end up with a
Poisson random variable (this
program blows up if you set alpha = 0, instead we recommend setting alpha =0.01).
The code fragment requires the specification of the **mean, p0** and **alpha**. The density
formula of the negative binomial distribution was taken from Regression Analysis
of Count Data by Cameron and Trivedi, formula (4.14).

clear set obs 1000 gen x = uniform( ) gen u = x local mean = 20 local alpha = .1 local p0 = .20 local a = 1/`alpha' local cp=0 local n =0 while `cp'<.99999 { local p = `p0'*(`n'==0) + (1-`p0')*(exp(lngamma(`n'+`a')) / (exp(lngamma(`a'))*(exp(lnfact(`n')))))*((`a'/(`a'+`mean'))^`a') /// *((`mean'/(`mean'+`a'))^`n') local cp= `cp'+`p' local pcp = `cp'-`p' quietly recode x (`pcp'/`cp' = `n') local n = `n' + 1 }

The content of this web site should not be construed as an endorsement of any particular web site, book, or software product by the University of California.