CODEWITHSUNDEEP
ralgorithmstatisticssimulationdistribution
user9802913
user9802913

Reputation: 245

Trying to understand the inverse transform method for generating a Poisson random variable

The algorithm Wikipedia gives for generating Poisson-distributed random variables using the inverse transform method is:

init:
     Let x ← 0, p ← e^−λ, s ← p.
     Generate uniform random number u in [0,1].
while u > s do:
     x ← x + 1.
     p ← p * λ / x.
     s ← s + p.
return x.

I implemented it in R:

f<-function(lambda)
{
x<-0
p<-exp(-lambda)
s<-p    
u<-runif(1)    
while (u>s )
{
   x<-x+1    
   p<-p*lambda/x    
   s<-s+p
}    
return(x)
}

but I don't understand how this generates values of this random variable, and I also think the code is incomplete because the output is always 1.

Could someone please explain?

Upvotes: 0

Views: 3458

Answers (1)

John Stone
John Stone

Reputation: 715

For Inverse transform sampling of poisson distribution, I suggest (e.g.) reading Chapter 4 in Simulation (2006, 4ed., Elsevier) by Sheldon M. Ross. Then according to the wiki steps mentioned, you can try following R coding

pois_wiki = function(n, lambda){
  ### Inverse transform sampling for poisson distribution
  ### ref: https://en.wikipedia.org/wiki/Poisson_distribution
  ### note: work efficiently for small λ
  X = rep(0, n) # generate positions for n samples
  for(j in 1:n){ # for each j-sample we do following wiki steps:
    x = 0; p = exp(-lambda); s = p
    u = runif(1) # Generate uniform random number u in [0,1]
    while(u > s){
      x = x + 1
      p = p*lambda/x
      s = s + p
    }
    X[j] = x # put generated sample for X in j-th element position
  }
  X # return n samples
}

### test in R
### note: for X~possion(λ), E(X)=Var(X)=λ
### set.seed(0) for reproducibility
n = 10000; lambda = 3
set.seed(0); mean(pois_wiki(n,lambda)); var(pois_wiki(n,lambda))
# [1] 3.005
# [1] 2.983695

Upvotes: 1

Related Questions