1

Estou tentando sobrepor dois gráficos, um da distribuição binomial e outro de poisson. Eles deveriam ter valores próximos para N grande. N aqui foi chamado de x e guardados os valores gerados por runif dentro de um vetor x. Da forma que estou plotando, a função lines fica sempre limitando Y a cerca de 170 e não sei o porquê. Qual outra função posso usar?

Minha intenção é simular até n=x=1000.inserir a descrição da imagem aqui

k<-2
p<-0.0001
i<-1
x<-0
while(i<=1000){
  x[i]<-runif(1,min=5,max=200)
  i<-i+1
}
fun<-function(x){
  (exp(-x*p))*((x*p)**k)/factorial(k)
}

plot(x, fun(x), main="b(k=2;n=60:180,p=0.2)", col="red", type="l")

func<-function(x){
  (factorial(x)/(factorial(k)*factorial(x-k)))*(p**k)*(1-p)**(x-k)
}
lines(x, func(x), col="blue")

2 Respostas 2

3

O R sempre vai plotar os pontos na ordem em que eles aparecem no vetor. Veja o exemplo abaixo, no qual ploto a função x^2:

x <- c(1, 5, 3, 6, 2, 4)
plot(x, x^2, type = "l")

inserir a descrição da imagem aqui

Veja como o gráfico parece não fazer sentido. A melhor maneira de organizar esses resultados e visualizar o gráfico desejado é ordenando os valores de x usando o comando sort:

x.ord <- sort(x)
plot(x.ord, x.ord^2, type = "l")

inserir a descrição da imagem aqui

Veja que agora o gráfico faz muito mais sentido.

Os teus dados sofrem do mesmo problema. Como os valores de x são gerados aleatoriamente, eles não estão ordenados. Portando, basta ordená-los para ter um gráfico que faça sentido:

k<-2
p<-0.0001
i<-1
x<-0
while(i<=1000){
    x[i]<-runif(1,min=5,max=200)
    i<-i+1
}

fun<-function(x){
    (exp(-x*p))*((x*p)**k)/factorial(k)
}

x.ord <- sort(x)

plot(x.ord, fun(x.ord), main="b(k=2;n=60:180,p=0.2)", col="red", type="l")

func<-function(x){
    (factorial(x)/(factorial(k)*factorial(x-k)))*(p**k)*(1-p)**(x-k)
}
lines(x.ord, func(x.ord), col="blue")

inserir a descrição da imagem aqui

O gráfico da binomial é cortado em torno de 170 porque o R tem um problema numérico para o fatorial acima de 170. Ele não consegue calcular (dá overflow, risos), o que implica que o valor da densidade da binomial não pode ser calculada.

factorial(171)/(factorial(171-2)*factorial(2))
## Inf

Uma maneira de resolver isso é usando a função choose, que calcula a combinação entre dois números:

choose(171, 2)
## [1] 14535

Outra forma de resolver isso é usar a função lfactorial, que calcula o logaritmo do fatorial. Assim, o overflow é evitado. Basta aplicar algumas propriedades de logaritmo e o resultado desejado sai:

exp(lfactorial(171) - (lfactorial(171-2)+lfactorial(2)))
## [1] 14535
0

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

inserir a descrição da imagem aqui

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

inserir a descrição da imagem aqui

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .