试读产生泊松随机变量的逆变换法

问题描述

Wikipedia给出的使用逆变换方法生成泊松分布随机变量的算法是：

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.

我在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)
}

但我不明白这是如何生成这个随机变量的值的，我也认为代码是不完整的，因为输出总是1。

有人能解释一下吗？

推荐答案

对于泊松分布的逆变换抽样，我建议(例如)阅读谢尔顿·M·罗斯Simulation (2006, 4ed., Elsevier)中的第四章。然后按照前面提到的wiki步骤，您可以尝试进行R编码

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

这篇关于试读产生泊松随机变量的逆变换法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！