Trying to understand the inverse transform method

2019-08-27 13:42发布

问题:

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?