Método da Rejeição

Há ocasiões em que o método da transformação inversa ou mesmo as transformações gerais não poderão ser usado para construir geradores de números aleatórios de muitas distribuições. Nesse casos, recorremos a métodos indiretos; ou seja, métodos nos quais geramos uma variável aleatória candidata e testamo a solução candidata para acaeitá-la ou não. Essa classe de métodos é bastante poderosa e nos permitirá simular praticamente qualquer distribuição. Dentre esses, os métodos de Aceitação/Rejeição exigem apenas que conheçamos a forma funcional da densidade \(f\) de interesse (chamada densidade alvo). Usamos uma densidade \(g\), cuja simulação é mais simples e que é denominada densidade candidata. Ela será usada para gerar a variável aleatória de interesse. As únicas restrições impostas a esta densidade candidata \(g\) são:

  1. \(f\) e \(g\) devem ter suportes compatíveis, ou seja, \(g(x) > 0\), quando \(f(x) > 0\).

  2. Há uma constante \(M\) com \(\frac{f(x)}{g(x)} \leq M\) para \(\forall x\)

Neste caso, podemos contruir um gerador de números aleatórios da variável de \(f_X\) usando o procedimento a seguir. Primeiramente geramos \(Y \sim g_Y\) e, independentemente, geramos \(U \sim \mathcal{Unif.}(0, 1)\). Se \[ U \leq \frac{1}{M}\frac{f_X(Y)}{g_Y(Y)}, \]

então \(X = Y\) Se a desigualdade não for satisfeita, descarta-se \(Y\) e \(U\) e tenta-se novamente, com o mesmo procedimento.

Resumidamente, o método de Aceitação/Rejeição pode ser representado no pseudocódigo a seguir:

Algoritmo - Método de Aceitação/Rejeição

  1. Gerar \(Y \sim g_Y\) e \(U \sim \mathcal{Unif.}(0,1)\).

  2. Aceitar \(X = Y\), se \(U f_X(Y) \leq M g_X(Y)\).

  3. Caso contrário, retornar ao passo 1.

Por exemplo, queremos gerar números da variável aleatória com função de densidade de probabilidade alvo dada por:

\[ f_X(x) = \begin{cases} \frac{2}{\sqrt{2\pi}} \operatorname{e}^{- \frac{x^2}{2}} & \text{, } \quad x \geq 0 \\ 0 &\text{, }\quad x<0 \end{cases} \]

Uma função candidata \(g\) é função de densidade exponencial padrão que tem o mesmo suporte de \(f\) e é fácil de simular (sua inversa é bastante simples). Sua densidade é: \[ g_Y(y) = \begin{cases} \operatorname{e}^{- y} & \text{, } \quad x \geq 0 \\ 0 &\text{, }\quad x<0 \end{cases} \] Temos de determinar o valor da constante \(M\). Na realidade queremos que um múltiplo de \(g\) esteja sempre envolvendo a função alvo \(f\). A razão entre as funções alvo e candidata são: \[ \frac{f_X(x)}{g_Y(x)} = \frac{\frac{2}{\sqrt{2\pi}} \operatorname{e}^{- \frac{x^2}{2}}}{\operatorname{e}^{-x}} = \frac{2}{\sqrt{2\pi}} \operatorname{e}^{x - \frac{x^2}{2}} \]

O máximo é atingido quando a derivada desta razão é zero, ou seja: \[ \frac{\operatorname{d}}{\operatorname{d} \! x} \left( \frac{f_X(x)}{g_Y(x)} \right) = \frac{2}{\sqrt{2\pi}}(1 - x) \operatorname{e}^{x - \frac{x^2}{2}} , \]

que é zero, para \(x = 1\). O ponto crítico então será: \[ \frac{f_x(1)}{g_Y(1)}= \frac{2}{\sqrt{2\pi}} \operatorname{e}^{1 - \frac{1}{2}} = \frac{2}{\sqrt{2\pi}} \operatorname{e}^{\frac{1}{2}}, \]

logo,
\[ \frac{f_x(x)}{g_Y(x)} \leq \frac{2}{\sqrt{2\pi}} \operatorname{e}^{\frac{1}{2}}=M \text{ , } \quad \forall x \] Observe o gráfico da densidade alvo e da candidata ajustada pela constante \(M\)

Neste caso a estrutura do algoritmo será

while(f.X <= Y){
    U = runif(1)
    V = runif(1)
    X = - log(U)
    Y = V * exp(-X) * M
}

Nessa estrutura exp(-X) * M é \(M \, g_Y(X)\) e se \(V = 1\), \(Y\) está situado na envoltória. Tentem implementar.