Tenho conhecimento que, por exemplo runif(1000,0,3)
gera 1000 valores aleatórios para uma distribuição uniforme para x compreendido entre 0 e 3. Mas como fazer isso para uma função densidade de probabilidade qualquer? Qualquer pista é agradecida!
1 Resposta
Toda distribuição estatística pode ser definida por uma função de distribuição acumulada F(x).
Um resultado bastante conhecido afirma que, se você possui uma variável aleatória U com distribuição uniforme no intervalo (0,1). Então segue a distribuição definida por F.
A prova desse resultado é simples:
Portanto se você consegue gerar números aleatórios com a distribuição uniforme e conhece a função de distribuição acumulada, você também consegue gerar números aleatórios de acordo com essa distribuição.
No R
, isso já está programado para diversas distribuições: ver a lista que o Ailton postou. Mas, não é muito complicado programar para uma outra distribuição se você conseguir inverter a sua F(x).
Se você conhecer apenas a função de distribuição acumulada, você pode escrever como um problema de otimização. Defina a distribuição acumulada ( Aqui é a distribuição acumulada da exponencial):
dF <- function(x, lambda = 0.5){
1 - exp(-lambda*x)
}
Sorteie um número aleatório entre 0 e 1 com distribuição uniforme. NO meu caso obtive:
n <- 0.917915
Em seguida, você tem que encontrar x
da distribuição dF
que mais se aproxima de n. Isso pode ser feito da seguinte forma:
op <- optim(runif(1), function(x){
abs(dF(x, 0.5) - n)
}, method = "BFGS")
Veja que:
> op$par
[1] 5
E que :
> dF(5, 0.5)
[1] 0.917915
Esse artigo da wikipedia explica de forma mais completa o que falei aqui: https://en.wikipedia.org/wiki/Inverse_transform_sampling
Esse procedimento pode ser repetido para, a partir de uma amostra aleatória de uma variável aleatória com distribuição uniforme obtermos uma amostra com distribuição acumulada definida no objeto dF
. No caso usamos a distribuição exponencial:
dF <- function(x, lambda = 0.5){
1 - exp(-lambda*x)
}
amostra <- runif(1000)
Para obter a amostra com distribuição exponencial é necessário encontrar o valor da inversa de dF
para cada número aleatório gerado. Isso pode ser obtido da seguinte maneira:
inverter_distribuicao <- function(x){
m <- nlm(function(y){
abs(dF(y, 0.5) - x)
}, p = 1
)
return(m$estimate)
}
amostra_exp <- sapply(amostra, inverter_distribuicao)
Veja agora o histograma da amostra gerada:
rbinom
, para a distribuição de Poisson,rpois
. Ambas no pacote stats que já vem instalado com o R. No site do projeto (em inglês) há uma lista com várias distribuições e os pacotes que dão suporte às mesmas. O que está listado em Base functionality está incluído em qualquer versão atualizada do R.