R中泊松随机变量生成的改进逆变换方法

Question

我正在阅读 Sheldon M. Ross Simulation (2006, 4ed., Elsevier) 的 Section 4.2，其中介绍了通过逆变换方法生成泊松随机变量。

表示pi =P(X=xi)=e^{-λ} λ^i/i!, i=0,1,...和F(i)=P(X<=i)=Σ_{k=0}^i pi分别为泊松的PDF和CDF，可以通过R中的dpois(x,lambda)和ppois(x,lambda)计算。

泊松逆变换算法有两种：普通版和改进版。

普通版步骤如下：

模拟来自 U(0,1) 的观察 U。
设置i=0和F=F(0)=p0=e^{-λ}。
如果U<F，select X=i并终止。
如果U >= F ，得到i=i+1, F=F+pi 和return到上一步。

我将以上步骤编写测试如下：

### write the regular R code
pois_inv_trans_regular = function(n, lambda){
  X = rep(0, n) # generate n samples
  for(m in 1:n){
    U = runif(1)
    i = 0; F = exp(-lambda) # initialize
    while(U >= F){
      i = i+1; F = F + dpois(i,lambda) # F=F+pi
    }
  X[m] = i
  }
X
}
### test the code (for small λ, e.g. λ=3)
set.seed(0); X = pois_inv_trans_regular(n=10000,lambda=3); c(mean(X),var(X))
# [1] 3.005000 3.044079

请注意，Poisson(λ)的均值和方差都是λ，因此常规代码的编写和测试是有意义的！

接下来我尝试了改进的，它是为大型λ设计的，根据书上的描述如下：

正则算法需要进行1+λ次搜索，即O(λ)计算复杂度，λ小的时候还好，大的时候可以当 λ 较大时得到改进。
事实上，由于具有均值 λ 的泊松随机变量最有可能取最接近 λ 的两个整数值之一，因此更有效的算法将首先检查这些值之一，而不是从 0 开始并向上工作。例如，令I=Int(λ) 并递归确定F(I)。
现在通过生成一个随机数U生成一个均值λ的泊松随机变量X，注意是否X <= I通过查看是否 U <= F(I) 。然后在X <= I 的情况下从I 开始向下搜索，否则从 I+1开始向上搜索。
据说改进后的算法只需要1+0.798√λ次搜索，即O(√λ)复杂度

我尝试编写改进后的R代码如下：

### write the improved R code
pois_inv_trans_improved = function(n, lambda){
  X = rep(0, n) # generate n samples
  p = function(x) {dpois(x,lambda)} # PDF: p(x) = P(X=x) = λ^x exp(-λ)/x!
  F = function(x) {ppois(x,lambda)} # CDF: F(x) = P(X ≤ x)
  I = floor(lambda) # I=Int(λ)
  F1 = F(I); F2 = F(I+1) # two close values
  for(k in 1:n){
    U = runif(1)
    i = I
    if ( F1 < U  &  U <= F2 ) { 
      i = I+1 
    } 
    while (U <= F1){ # search downward
      i = i-1; F1 = F1 - p(i)
    }
    while (U > F2){ #  search upward
      i = i+1; F2 = F2 + p(i)
    }
    X[k] = i
  }
  X
}
### test the code (for large λ, e.g. λ=100)
set.seed(0); X = pois_inv_trans_improved(n=10000,lambda=100); c(mean(X),var(X))
# [1] 100.99900000   0.02180118

从 [1] 100.99900000 0.02180118 的模拟结果来看 c(mean(X),var(X))，这表明方差部分毫无意义。我应该如何解决这个问题？

Answer 1

主要问题是F1和F2是在循环中修改的，没有重置，所以最终很大范围的U被认为在中间。
第二个问题是向下搜索时使用的p(i)应该是原来的i，因为F(x) = P(X <= x)。如果没有这个，代码会因低 U 而挂起。最简单的解决方法是开始 i = I + 1。如果不需要语句，则“在中间”。

pois_inv_trans_improved = function(n, lambda){
  X = rep(0, n) # generate n samples
  p = function(x) {dpois(x,lambda)} # PDF: p(x) = P(X=x) = λ^x exp(-λ)/x!
  `F` = function(x) {ppois(x,lambda)} # CDF: F(x) = P(X ≤ x)
  I = floor(lambda) # I=Int(λ)
  F1 = F(I); F2 = F(I+1) # two close values
  for(k in 1:n){
    U = runif(1)
    i = I + 1
    # if ( F1 < U  &  U <= F2 ) { 
    #   i = I + 1
    # }
    F1tmp = F1
    while (U <= F1tmp){ # search downward
      i = i-1; F1tmp = F1tmp - p(i);  
    }
    F2tmp = F2
    while (U > F2tmp){ #  search upward
      i = i+1; F2tmp = F2tmp + p(i)
    }
    X[k] = i
  }
  X
}

这给出：

[1] 100.0056 102.2380

R中泊松随机变量生成的改进逆变换方法

Improved inverse transform method for Poisson random variable generation in R

r

poisson