O problema
Tal como o Marcus Nunes diz na sua resposta, o problema vem do cálculo da função mass de probabilidade da distribuição binomial. Os fatoriais crescem muito depressa e dá-se uma condição de overflow.
Ora, o problema de calcular a função dbinom
é um problema muito difícil. O código R é uma implementação de
Catherine Loader (2000). Fast and Accurate Computation of Binomial
Probabilities.
Veja um thread recente do r-devel sobre o cálculo de dbinom
:
Mais à frente vou calcular os valores de func
com dbinom
. Essa nova função será chamada de func2
.
A pergunta
O código para simular as variáveis uniformes pode ser muito simplificado. Basta ver que em R muitas funções são vetorizadas, e uma delas é runif
. Qualquer dos três modos de gerar o vetor x
(ou x2
ou y
) que se segue dá resultados idênticos.
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
Além disso, os conjuntos suporte das distribuições binomial e de Poisson são conjuntos de números inteiros. Não faz sentido usar o 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")
Agora com 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")