为什么我在模拟泊松随机变量的“for”循环中得到许多 NA答案

【问题标题】：Why do I get many NA's in a "for" loop that simulates Poisson random variables为什么我在模拟泊松随机变量的“for”循环中得到许多 NA
【发布时间】：2019-03-08 04:15:33
【问题描述】：

我不断收到错误消息，说我在for 循环中创建了数百个NA。那些NA 来自哪里？任何帮助将不胜感激！

drip <- function(rate = 1, minutes = 120) {
  count <- 0
  for(i in 1:(minutes)) {
    count <- count + rpois(1, rate)
    rate <- rate * runif(1, 0, 5)
  }
  count
}
drip()

【问题讨论】：

标签： r loops for-loop na poisson

【解决方案1】：

你得到整数溢出。试试

set.seed(0)
rpois(1, 1e+8)
#[1] 100012629
rpois(1, 1e+9)
#[1] 999989683
rpois(1, 1e+10)
#[1] NA
#Warning message:
#In rpois(1, 1e+10) : NAs produced

只要lambda 太大，整数的32 位表示就不够了，就会返回NA。（回想一下，泊松随机变量是整数）。

您的循环在rate (lambda) 上有动态增长，最终可能变得太大。使用较小的minutes（例如10）运行您的函数很好。

相比之下，产生双精度浮点数的ppois 和dpois 适合较大的lambda。

dpois(1e+8, 1e+8)
#[1] 3.989423e-05
dpois(1e+9, 1e+9)
#[1] 1.261566e-05
dpois(1e+10, 1e+10)
#[1] 3.989423e-06
dpois(1e+11, 1e+11)
#[1] 1.261566e-06

ppois(1e+8, 1e+8)
#[1] 0.5000266
ppois(1e+9, 1e+9)
#[1] 0.5000084
ppois(1e+10, 1e+10)
#[1] 0.5000027
ppois(1e+11, 1e+11)
#[1] 0.5000008

每经过一分钟，速率参数会增加x%，其中x 是区间[0, 5] 上均匀分布的随机值。

rate 增加了x% 而不是x。所以你应该使用

rate <- rate * (1 + runif(1, 0, 5) / 100)

【讨论】：