试图理解生成泊松随机变量的逆变换方法答案

【问题标题】：Trying to understand the inverse transform method for generating a Poisson random variable试图理解生成泊松随机变量的逆变换方法
【发布时间】：2019-08-20 22:40:02
【问题描述】：

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。

谁能解释一下？

【问题讨论】：

您是否熟悉并熟悉连续随机变量的逆变换方法？除了泊松之外的离散随机变量怎么样？
@pjs 我还不熟悉连续随机变量的逆变换方法。并且以某种方式熟悉离散随机变量:)。为什么？
我试图找出解决您问题的合理起点。您可以查看this paper 的第 4.3.1 节和附录 A。如果这对您有意义，我将探讨离散反演如何适用于泊松。
@pjs 我已经检查过上面提到的论文，是的，它有某种意义。

标签： r algorithm statistics simulation distribution

【解决方案1】：

对于泊松分布的逆变换采样，我建议（例如）阅读 Sheldon M. Ross 在Simulation (2006, 4ed., Elsevier) 中的第 4 章。然后按照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

【讨论】：