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

Question

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

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

泊松逆变换算法有两种:常规版本和改进版本。

普通版步骤如下:

1. 模拟观察 U来自U(0,1)​ .
2. 设置i=0​和​F=F(0)=p0=e^{-λ}​ .
3. 如果U<F​ ，选择​X=​i并终止。
4. 如果U >= F ​，获取i=i+1, F=F+pi​并返回上一步。

我编写并测试上述步骤如下:

### 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) .

• 现在生成泊松随机变量 X平均值 λ通过生成随机数 U ，注意是否 X <= I​通过查看是否​ U <= F(I) ​.然后从​I开始向下搜索​ 在 X <= 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。然后“在中间”不需要 if 语句。

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