The problem
As Marcus Nunes says in his reply, the problem comes from calculating the mass probability function of the binomial distribution. Factorials grow very fast and an overflow condition occurs.
The problem of calculating the function dbinom
is a very difficult problem. R code is an implementation of
Catherine Loader (2000). Fast and Accurate Computation of Binomial
Probabilities.
See a recent thread of r-devel on the calculation of dbinom
:
Later I will calculate the values of func
with dbinom
. This new function will be called func2
.
The question
The code to simulate uniform variables can be greatly simplified. Just see that in R many functions are vectorized, and one of them is runif
. Any of the three modes of generating the vector x
(or x2
or y
) which gives identical results.
set.seed(2020)
i <- 1
x <- 0
while(i <= 1000){
x[i] <- runif(1, min = 5, max = 200)
i <- i + 1
}
set.seed(2020)
x2 <- numeric(1000)
for(i in seq.int(1000)){
x2[i] <- runif(1, min = 5, max = 200)
}
set.seed(2020)
y <- runif(1000, min = 5, max = 200)
identical(x, x2) #[1] TRUE
identical(x, y) #[1] TRUE
In addition, the support sets of binomial and Poisson distributions are sets of integer numbers. It makes no sense to use the runif
.
fun <- function(x, p){
(exp(-x*p))*((x*p)**k)/factorial(k)
}
func <- function(x, p){
choose(x, k)*(p**k)*(1 - p)**(x-k)
}
n <- 1000
k <- 2
p <- 0.0001
set.seed(2020)
x <- sample(n, n, TRUE)
x.ord <- sort(x)
plot(x.ord, fun(x.ord, p), main="b(k=2;n=60:180,p=0.2)", col="red", type="l")
lines(x.ord, func(x.ord, p), col="blue", lty = "dotted")
Now with dbinom
.
func2 <- function(x, p) dbinom(k, size = x, prob = p)
plot(x.ord, fun(x.ord, p), main="b(k=2;n=60:180,p=0.2)", col="red", type="l")
lines(x.ord, func2(x.ord, p), col="blue", lty = "dotted")