5

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!

4
  • 3
    Não sou especialista, mas uma busca rápida parece indicar que o problema é mais matemático/estatístico do que de programação. Aparentemente a solução envolve manipular a sua função de forma que ela possa ser alimentada com valores aleatórios (uniformes), e daí gerar os valores com tal distribuição.
    – Molx
    9/12/2015 às 23:59
  • Segundo este link, você pode obter uma aproximação disso simplesmente ao inverter a função e resolver para valores com distribuição uniforme (se a função for discreta). O ideal seria você postar a função para que possamos tentar resolver o problema.
    – Molx
    10/12/2015 às 0:05
  • 2
    Cada distribuição possui uma função específica no R para obtermos amostras. Por exemplo, para obter valores de uma distribuição binomial, a função é 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. 10/12/2015 às 2:42
  • Se a função de densidade é conhecida, vale o comentário do @AiltonAndradedeOliveira, pois há uma função para cada distribuição. Caso seja desconhecida, mas você tenha uma amostra, a questão se complica (eg. qual o domínio da função? Reais? Reais não-negativos? Intervalo definido?).
    – Bernardo
    11/12/2015 às 18:56

1 Resposta 1

5

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 inserir a descrição da imagem aqui segue a distribuição definida por F.

A prova desse resultado é simples:

inserir a descrição da imagem aqui

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:

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 .