CODEWITHSUNDEEP
rloopsfor-loopnapoisson
Tanner
Tanner

Reputation: 25

Why do I get many NA's in a "for" loop that simulates Poisson random variables

I keep getting error messages saying that I am creating hundreds of NA's in my for loop. Where do those NA come from? Any help would be greatly appreciated!

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()

Upvotes: 1

Views: 390

Answers (1)

Zheyuan Li
Zheyuan Li

Reputation: 73295

You get integer overflow. Try

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

As soon as lambda is too large, 32-bit representation of integer is insufficient and NA is returned. (Recall that Poisson random variables are integers).

Your loop has a dynamic growth on rate (lambda), which can eventually become too big. Running your function with a smaller minutes, say 10, is fine.

By contrast, ppois and dpois which produce double-precision floating point numbers are fine with large 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

With each passing minute, the rate parameter increases by x%, where x is a random value from the uniform distribution on the interval [0, 5].

The rate increases by x% not x. So you should use

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

Upvotes: 2

Related Questions