我正在阅读4.2节本文由Sheldon M. Ross在Simulation (2006, 4ed., Elsevier)中介绍了用逆变换方法生成泊松随机变量。
令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，选择X=i并终止。 - 如果是
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)。现在通过生成随机数
U生成一个泊松随机变量X，其平均值为λ，通过观察U <= F(I)是否存在来判断X <= 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))的模拟结果来看，方差部分是无意义的。我应该如何解决这个问题?
主要问题是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