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.16067318897452204
0.033349619046247295
0.060702521234135595
0.20527470654846877
0.4827442694605983
0.4212568112173069
0.22506353725848197
0.31423611447458877
0.09785290475075122
0.5862559913111378
⋮
0.3577854206128256
0.1419468470244178
0.37435776256621733
0.16280599822437028
0.6562827314351217
0.24576904134618077
0.1076270956133472
0.16586266836271255
0.06137027711633528
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.20768136762713474
0.6804525084507673
0.7894702684775973
0.6672636685448472
0.25011572331795295
0.3092144185444903
0.1561267924552292
0.9260514845801957
0.25368488104584397
0.8601731998444724
⋮
0.3372750472126288
0.8828745586754533
0.34917597743737205
0.3914565129712684
0.4557472318693915
0.19860397090206872
0.3936468931973939
0.7314791360468379
0.18168776722936253
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.