using Distributions
using StatsPlotsSimulação de variáveis aleatórias com os métodos de Inversão e Rejeição
Veja a implementação dos métodos de Rejeição e Inversão em Julia, ferramentas poderosas para a simulação de variáveis aleatórias de diversas distribuições de probabilidade..
Introdução
Em muitos problemas de Estatística e Simulação, surge a necessidade de gerar amostras de distribuições de probabilidade que não estão diretamente implementadas em pacotes prontos. Nesse contexto, os métodos da Inversão e da Rejeição são ferramentas fundamentais, pois permitem construir algoritmos simples e eficientes para esse fim.
Neste post apresentaremos ambos os métodos e daremos exemplos de implementação ideias para pegar e usar.
Método da Inversão
Em primeiro lugar, vamos carregar os pacotes que serão utilizados:
O Método da Inversão consiste, basicamente, em aplicar a inversa da função de distribuição acumulada sobre uma amostra de variáveis aleatórias uniformes. Mais detalhes podem ser encontrados neste link.
Podemos definir uma função amostra_inversao() que automatize esse processo. Essa função recebe como argumentos:
n: Tamanho da amostra desejada;inv_acumulada: Inversa da função de distribuição acumulada da distribuição de interesse.
function amostra_inversao(n, inv_acumulada)
return amostras = inv_acumulada.(rand(Uniform(0,1), n))
endamostra_inversao (generic function with 1 method)
Para exemplificar, vamos obter uma amostra de tamanho 1000 da distribuição Exponencial(2). Basta definirmos a inversa da sua acumulada e utilizar a função criada:
inv_exp(x) = log(1 - x)/(-2);
amostras_exponencial = amostra_inversao(1000, inv_exp)1000-element Vector{Float64}:
0.7162636586004516
0.23362229189684883
0.7769275712962066
0.6534392047804466
0.4367840693268625
0.5733200782273268
0.11026774022591004
0.04536942839111509
0.08568063825850655
0.0230905198412386
⋮
0.8496533595421486
0.088497549870055
0.33127922386395997
0.0062714504013609215
0.29479722722589996
0.18434677620863832
0.3517994291741801
0.3816870026331362
0.06904407400308499
Para verificar que de fato obtivemos uma amostra da distribuição correta, podemos comparar o histograma das amostras com a densidade teórica:
densidade_exponencial(x) = 2*exp(-2*x);
x = range(0, maximum(amostras_exponencial), length = 1000);
histogram(amostras_exponencial, normalize=:pdf, label= "Amostra")
plot!(x, densidade_exponencial.(x), color=:red, label= "Distribuição")
Método da Rejeição
O Método da Rejeição consiste em escolher uma distribuição auxiliar cuja densidade (\(g(x)\)), multiplicada por uma constante C, seja sempre maior ou igual à densidade da distribuição de interesse(\(f(x)\)), ou seja \(f(x) \leq C g(x)\). Uma explicação detalhada pode ser encontrada neste link.
Podemos definir uma função amostra_rejeicao() para esse procedimento. Ela recebe os seguintes argumentos:
n: Tamanho da amostra desejada;f: Densidade da distribuição de interesse;g: Densidade da distribuição auxiliar;amostra_g: Função que gera uma observação da distribuição auxiliar;C: Constante a ser utilizada.
function amostra_rejeicao(n, f, g, amostra_g, C)
amostras = Vector{Float64}(undef, n)
n_amostras = 1
while n_amostras <= n
Y = amostra_g()
U = rand()
if U <= f(Y)/(C*g(Y))
amostras[n_amostras] = Y
n_amostras = n_amostras + 1
end
end
return amostras
endamostra_rejeicao (generic function with 1 method)
Para exemplificar, vamos obter uma amostra de tamanho 1000 da distribuição Beta(2,2) utilizando como distribuição auxiliar a Uniforme(0,1). Após fazer as contas, a constante \(C\) necessária é 3/2.
\(C = \displaystyle max_{0\leq x \leq 1} \dfrac{f(x)}{g(x)} = max_{0\leq x \leq 1} 6 \times (x - x^2) / 1\). O ponto de máximo acontece em \(x = 1/2\) e f(1/2)/g(1/2) = 3/2.
densidade_beta(x) = 6*(x - x^2); #Beta(2, 2)densidade_beta (generic function with 1 method)
densidade_uniforme(x) = 1;
amostras_beta = amostra_rejeicao(1000, densidade_beta, densidade_uniforme, rand, 3/2)1000-element Vector{Float64}:
0.9383794162282901
0.8080231892548565
0.3559616842217317
0.7120191990173921
0.8555349671169985
0.37322122662036317
0.7752462960292195
0.7472584632958286
0.14177554541002357
0.25336391057040286
⋮
0.2449935273635483
0.500684502505503
0.6632273199686185
0.5481639258176497
0.8334555174235013
0.7088198652843193
0.7070527922216966
0.33066485034624726
0.26231317228320195
Assim como no caso anterior, para verificar que de fato obtivemos uma amostra da distribuição correta, podemos comparar o histograma das amostras com a densidade teórica:
x = range(0, 1, length = 1000)0.0:0.001001001001001001:1.0
histogram(amostras_beta, normalize=:pdf, label= "Amostra")
plot!(x, densidade_beta.(x), color=:red, label= "Distribuição")
Conclusão
Os métodos da Inversão e da Rejeição são estratégias fundamentais para gerar amostras de distribuições arbitrárias quando não dispomos de funções prontas em pacotes estatísticos. Portanto, com essas implementações, passamos a ter ferramentas práticas e flexíveis para simulação, reforçando o poder do Julia como linguagem para Estatística Computacional.
Ferramentas de IA foram utilizadas para correção ortográfica e aprimoramento do texto.