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?