8  Método da Rejeição para Variáveis Contínuas

No capítulo anterior usamos o método da rejeição para gerar v.a. discretas: sorteávamos candidatos a partir de uma distribuição fácil de simular e aceitávamos cada candidato com probabilidade proporcional à razão entre o alvo e a proposta. A versão contínua é a mesma ideia, palavra por palavra, trocando as funções de probabilidade \(p_j\) e \(q_j\) por densidades \(f\) e \(g\).

Aqui o método é ainda mais valioso. No Capítulo 4 vimos que a inversão exige uma fórmula explícita para \(F^{-1}\), e avisamos que isso é a exceção: para a normal, por exemplo, nem mesmo \(F\) tem forma fechada. O método da rejeição não precisa de \(F\) nem de \(F^{-1}\) — basta saber calcular a densidade alvo em um ponto e saber simular de alguma outra densidade parecida com ela. Como veremos no Exemplo 2, isso já é suficiente para gerar valores da normal.

8.1 Algoritmo

Seja \(X\) uma v.a. contínua com densidade \(f\); essa é a distribuição alvo. O método supõe que sabemos simular uma segunda v.a. \(Y\), com densidade \(g\), chamada de distribuição proposta.

Atenção: a proposta precisa cobrir o suporte do alvo

É indispensável que \(g(x) > 0\) sempre que \(f(x) > 0\). Se a proposta nunca sugere valores de uma região onde o alvo tem densidade positiva, essa região jamais aparecerá na amostra — e o algoritmo devolve, silenciosamente, uma distribuição errada.

Além disso, supomos conhecida uma constante \(c\) tal que

\[ \frac{f(x)}{g(x)} \leq c \quad \text{para todo } x \text{ com } f(x) > 0. \]

Exatamente como no caso discreto, necessariamente \(c \geq 1\): integrando a desigualdade \(f(x) \leq c\, g(x)\) sobre o conjunto onde \(f\) é positiva, obtemos

\[ 1 = \int_{\{f > 0\}} f(x)\, dx \;\leq\; c \int_{\{f > 0\}} g(x)\, dx \;\leq\; c, \]

pois a integral de \(g\) sobre um subconjunto da reta é, no máximo, \(1\).

Pseudo-algoritmo: rejeição para v.a. contínuas

  1. Gere \(Y\) a partir da densidade proposta \(g\).

  2. Gere \(U \sim \text{Unif}(0,1)\), independente de \(Y\).

  3. Se \(U \leq \dfrac{f(Y)}{c\, g(Y)}\), devolva \(X = Y\). Caso contrário, descarte \(Y\) e volte ao passo 1.

Como no caso discreto, o passo 3 é apenas um sorteio de Bernoulli: aceitamos o candidato com probabilidade \(f(Y) / (c\, g(Y))\), e é a escolha de \(c\) que garante que esse número esteja entre \(0\) e \(1\).

8.2 Interpretação geométrica

A figura a seguir ilustra o passo 3 no caso em que a proposta é uma \(\text{Unif}(0,1)\), de modo que \(c\,g(x)\) é uma reta horizontal na altura \(c\). Sorteado um candidato \(Y\), olhamos a haste vertical que vai de \(0\) até \(c\, g(Y)\): a parte verde (abaixo de \(f(Y)\)) corresponde a aceitar, e a parte vermelha a rejeitar. O papel de \(U\) é sortear um ponto uniformemente ao longo dessa haste.

Mostrar código 
library(ggplot2)

# Densidade alvo: f(x) = 20 x (1-x)^3, para 0 < x < 1
f <- function(x) 20 * x * (1 - x)^3

# Densidade proposta: Unif(0,1). O rep() faz a função servir tanto para um
# número quanto para um vetor de valores
g <- function(x) rep(1, length(x))

# Menor constante com f(x) <= cte * g(x). Chamamos de `cte` (e não de `c`)
# para não esconder a função c() do R
cte <- 135 / 64

# Candidato usado apenas para ilustrar o sorteio de aceitação
Y <- 0.6

grade <- seq(0, 1, length.out = 400)
df_curvas <- data.frame(x = grade, f = f(grade), teto = cte * g(grade))

ggplot(df_curvas, aes(x = x)) +
  geom_line(aes(y = f, color = "f(x)"), linewidth = 1) +
  geom_line(aes(y = teto, color = "c g(x)"), linewidth = 1) +
  # A haste em Y vai de 0 até o teto c*g(Y); abaixo de f(Y) aceitamos
  annotate("segment", x = Y, xend = Y, y = 0, yend = f(Y),
           color = "darkgreen", linewidth = 1.5) +
  annotate("segment", x = Y, xend = Y, y = f(Y), yend = cte * g(Y),
           color = "red", linewidth = 1.5) +
  annotate("text", x = Y + 0.03, y = f(Y) / 2, label = "aceita",
           color = "darkgreen", hjust = 0) +
  annotate("text", x = Y + 0.03, y = (f(Y) + cte) / 2, label = "rejeita",
           color = "red", hjust = 0) +
  annotate("text", x = Y, y = -0.09, label = "Y") +
  scale_color_manual(values = c("f(x)" = "black", "c g(x)" = "blue"), name = "") +
  coord_cartesian(ylim = c(-0.12, 2.6)) +
  labs(x = "x", y = "Densidade") +
  theme_minimal() +
  theme(legend.position = "top")

Mostrar código 
import numpy as np
import matplotlib.pyplot as plt

# Densidade alvo: f(x) = 20 x (1-x)^3, para 0 < x < 1
def f(x):
    return 20 * x * (1 - x)**3

# Densidade proposta: Unif(0,1). O ones_like faz a função servir tanto para um
# número quanto para um vetor de valores
def g(x):
    return np.ones_like(x, dtype=float)

# Menor constante com f(x) <= cte * g(x)
cte = 135 / 64

# Candidato usado apenas para ilustrar o sorteio de aceitação
Y = 0.6

grade = np.linspace(0, 1, 400)

plt.figure(figsize=(8, 5))
plt.plot(grade, f(grade), color="black", label="f(x)")
plt.plot(grade, cte * g(grade), color="blue", label="c g(x)")

# A haste em Y vai de 0 até o teto c*g(Y); abaixo de f(Y) aceitamos
plt.vlines(Y, 0, f(Y), color="darkgreen", linewidth=3)
plt.vlines(Y, f(Y), cte * g(Y), color="red", linewidth=3)

plt.text(Y + 0.03, f(Y) / 2, "aceita", color="darkgreen")
plt.text(Y + 0.03, (f(Y) + cte) / 2, "rejeita", color="red")
plt.text(Y, -0.09, "Y", horizontalalignment="center")

plt.ylim(-0.12, 2.6)
(-0.12, 2.6)
Mostrar código 
plt.xlabel("x")
plt.ylabel("Densidade")
plt.legend(loc="upper center", ncol=2)
plt.show()

Essa leitura sugere uma segunda maneira de enxergar o método, que costuma ser mais esclarecedora do que a conta. Chame de \(V = U \cdot c\, g(Y)\) a altura sorteada na haste. O par \((Y, V)\) é um ponto sorteado uniformemente na região abaixo da curva \(c\,g\), e o passo 3 mantém apenas os pontos que caem abaixo de \(f\). O gráfico abaixo mostra 500 candidatos e o que acontece com cada um deles.

Mostrar código 
set.seed(42)

n_candidatos <- 500

# Passo 1: candidatos da proposta Unif(0,1)
Y <- runif(n_candidatos)

# Passo 2: os uniformes que decidem a aceitação
U <- runif(n_candidatos)

# Altura sorteada na haste: um ponto uniforme entre 0 e c*g(Y)
V <- U * cte * g(Y)

# Passo 3: ficamos com os pontos que caíram abaixo da curva f
aceito <- V <= f(Y)

df_pontos <- data.frame(
  x = Y,
  altura = V,
  status = ifelse(aceito, "aceito", "rejeitado")
)

ggplot(df_pontos, aes(x = x, y = altura, color = status)) +
  geom_point(size = 1.5, alpha = 0.7) +
  geom_line(data = df_curvas, aes(x = x, y = f), inherit.aes = FALSE,
            linewidth = 1) +
  scale_color_manual(values = c("aceito" = "darkgreen", "rejeitado" = "red"),
                     name = "") +
  coord_cartesian(ylim = c(0, 2.6)) +
  labs(x = "x", y = "Altura sorteada") +
  theme_minimal() +
  theme(legend.position = "top")

Mostrar código 
np.random.seed(42)

n_candidatos = 500

# Passo 1: candidatos da proposta Unif(0,1)
Y = np.random.uniform(0, 1, n_candidatos)

# Passo 2: os uniformes que decidem a aceitação
U = np.random.uniform(0, 1, n_candidatos)

# Altura sorteada na haste: um ponto uniforme entre 0 e c*g(Y)
V = U * cte * g(Y)

# Passo 3: ficamos com os pontos que caíram abaixo da curva f
aceito = V <= f(Y)

plt.figure(figsize=(8, 5))
plt.scatter(Y[aceito], V[aceito], s=12, alpha=0.7, color="darkgreen",
            label="aceito")
plt.scatter(Y[~aceito], V[~aceito], s=12, alpha=0.7, color="red",
            label="rejeitado")
plt.plot(grade, f(grade), color="black", linewidth=1)

plt.ylim(0, 2.6)
(0.0, 2.6)
Mostrar código 
plt.xlabel("x")
plt.ylabel("Altura sorteada")
plt.legend(loc="upper center", ncol=2)
plt.show()

Outra maneira de ver o método

Os pontos verdes estão espalhados uniformemente na região sob a curva \(f\). E um ponto uniforme sob o gráfico de uma densidade tem uma propriedade notável: sua abscissa tem exatamente essa densidade. A razão é que a região é mais alta onde \(f\) é maior, então mais pontos caem sobre esses valores de \(x\) — na proporção certa.

Ou seja, o método da rejeição pode ser lido assim: para simular de \(f\), sorteie um ponto uniforme sob o gráfico de \(f\) e devolva sua coordenada horizontal. A proposta \(g\) e a constante \(c\) servem só para construir uma região maior, e fácil de sortear, que contenha a região de interesse. O Exercício 8 pede a demonstração desse fato.

8.3 Exemplo 1: uma densidade em \((0,1)\) com proposta uniforme

Queremos gerar valores da densidade

\[ f(x) = 20x(1 - x)^3, \quad 0 < x < 1, \]

que é a densidade de uma \(\text{Beta}(2,4)\). Como o suporte é o intervalo \((0,1)\), a escolha mais simples de proposta é \(g(x) = 1\) para \(0 < x < 1\), isto é, \(Y \sim \text{Unif}(0,1)\).

Nesse caso a razão \(f(x)/g(x)\) é a própria \(f(x)\), e a menor constante possível é o valor máximo da densidade. Derivando,

\[ f'(x) = 20\left[(1-x)^3 - 3x(1-x)^2\right] = 20(1-x)^2(1 - 4x), \]

que se anula em \(x = 1/4\) (e em \(x = 1\), onde \(f\) vale zero). Logo

\[ c = f(1/4) = 20 \cdot \frac{1}{4} \cdot \left(\frac{3}{4}\right)^3 = \frac{135}{64} \approx 2{,}11. \]

Pseudo-algoritmo: Beta(2,4) com proposta uniforme

  1. Gere \(Y \sim \text{Unif}(0,1)\).

  2. Gere \(U \sim \text{Unif}(0,1)\).

  3. Se \(U \leq \dfrac{20\,Y(1-Y)^3}{135/64}\), devolva \(X = Y\). Caso contrário, volte ao passo 1.

O código a seguir gera \(B = 1000\) valores. Além das amostras aceitas, contamos quantos candidatos foram necessários, para comparar a taxa de aceitação observada com o valor teórico \(1/c\).

Mostrar código 
library(ggplot2)

set.seed(42)

f <- function(x) 20 * x * (1 - x)^3
g <- function(x) 1          # densidade da Unif(0,1)
cte <- 135 / 64

B <- 1000            # quantos valores queremos gerar
amostras <- c()      # guarda os valores aceitos
n_propostas <- 0     # conta quantas vezes o passo 1 foi executado

while (length(amostras) < B) {
  # Passo 1: gerar Y da proposta Unif(0,1)
  Y <- runif(1)
  n_propostas <- n_propostas + 1

  # Passo 2: gerar U ~ Unif(0,1)
  U <- runif(1)

  # Passo 3: aceitar Y com probabilidade f(Y) / (cte * g(Y))
  if (U <= f(Y) / (cte * g(Y))) {
    amostras <- c(amostras, Y)
  }
}

cat("Candidatos gerados:", n_propostas, "\n")
Candidatos gerados: 2081
Mostrar código 
cat("Taxa de aceitação observada:", round(B / n_propostas, 3), "\n")
Taxa de aceitação observada: 0.481
Mostrar código 
cat("Taxa de aceitação teórica (1/c):", round(1 / cte, 3), "\n")
Taxa de aceitação teórica (1/c): 0.474
Mostrar código 
df <- data.frame(x = amostras)

ggplot(df, aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30,
                 fill = "skyblue", color = "black") +
  stat_function(fun = f, color = "red", linewidth = 1) +
  labs(title = "Rejeição para f(x) = 20x(1-x)^3",
       x = "Valor de X", y = "Densidade") +
  theme_minimal()

Mostrar código 
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

def f(x):
    return 20 * x * (1 - x)**3

def g(x):
    return 1.0              # densidade da Unif(0,1)

cte = 135 / 64

B = 1000             # quantos valores queremos gerar
amostras = []        # guarda os valores aceitos
n_propostas = 0      # conta quantas vezes o passo 1 foi executado

while len(amostras) < B:
    # Passo 1: gerar Y da proposta Unif(0,1)
    Y = np.random.uniform(0, 1)
    n_propostas += 1

    # Passo 2: gerar U ~ Unif(0,1)
    U = np.random.uniform(0, 1)

    # Passo 3: aceitar Y com probabilidade f(Y) / (cte * g(Y))
    if U <= f(Y) / (cte * g(Y)):
        amostras.append(Y)

amostras = np.array(amostras)

print("Candidatos gerados:", n_propostas)
Candidatos gerados: 2084
Mostrar código 
print("Taxa de aceitação observada:", round(B / n_propostas, 3))
Taxa de aceitação observada: 0.48
Mostrar código 
print("Taxa de aceitação teórica (1/c):", round(1 / cte, 3))
Taxa de aceitação teórica (1/c): 0.474
Mostrar código 
# Malha usada só para desenhar a densidade teórica
grade = np.linspace(0, 1, 400)

plt.figure(figsize=(8, 5))
plt.hist(amostras, bins=30, density=True, color="skyblue", edgecolor="black")
plt.plot(grade, f(grade), color="red", linewidth=2)
plt.title("Rejeição para f(x) = 20x(1-x)^3")
plt.xlabel("Valor de X")
plt.ylabel("Densidade")
plt.show()

8.4 Exemplo 2: gerando uma normal a partir de exponenciais

Este é o exemplo que a inversão não conseguia resolver. Queremos gerar \(X \sim N(0,1)\), cuja f.d.a. não tem forma fechada — muito menos sua inversa.

A ideia é gerar primeiro o módulo \(|X|\) e depois sortear o sinal. Se \(X \sim N(0,1)\), a densidade de \(|X|\) é o dobro da densidade da normal restrita aos valores positivos:

\[ f(x) = \sqrt{\frac{2}{\pi}}\, e^{-x^2/2}, \quad x > 0. \]

Como proposta, usamos \(Y \sim \text{Exp}(1)\), com densidade \(g(x) = e^{-x}\) para \(x > 0\) — uma distribuição que já sabemos simular por inversão (Capítulo 4) e que, como o alvo, tem suporte em \((0, \infty)\) e cai a zero na cauda direita. A razão entre as duas densidades é

\[ \frac{f(x)}{g(x)} = \sqrt{\frac{2}{\pi}}\, e^{x - x^2/2}. \]

O expoente \(x - x^2/2\) é máximo em \(x = 1\), onde vale \(1/2\). Portanto

\[ c = \sqrt{\frac{2}{\pi}}\, e^{1/2} = \sqrt{\frac{2e}{\pi}} \approx 1{,}3155, \]

e a taxa de aceitação é \(1/c \approx 0{,}76\): descartamos menos de um quarto dos candidatos.

O critério de aceitação fica particularmente simples. Substituindo \(c\),

\[ \frac{f(Y)}{c\, g(Y)} = \frac{\sqrt{2/\pi}\, e^{Y - Y^2/2}}{\sqrt{2/\pi}\, e^{1/2}} = e^{-(Y-1)^2/2}, \]

uma expressão que não envolve nem \(\pi\) nem raízes.

Pseudo-algoritmo: normal padrão por rejeição

  1. Gere \(Y \sim \text{Exp}(1)\).

  2. Gere \(U \sim \text{Unif}(0,1)\).

  3. Se \(U > e^{-(Y-1)^2/2}\), volte ao passo 1. Caso contrário, siga: \(Y\) é um valor de \(|X|\).

  4. Sorteie um sinal \(S\), igual a \(+1\) ou \(-1\) com probabilidade \(1/2\) cada, e devolva \(X = S \cdot Y\).

O passo 4 é o que transforma o módulo em uma normal: como a densidade da \(N(0,1)\) é simétrica em torno de zero, metade da massa está de cada lado.

Mostrar código 
library(ggplot2)

set.seed(42)

# A constante c só é usada no relatório final: o critério do passo 3 já está
# escrito na forma simplificada
cte <- sqrt(2 * exp(1) / pi)

B <- 5000            # quantos valores queremos gerar
modulos <- c()       # guarda os valores aceitos de |X|
n_propostas <- 0

while (length(modulos) < B) {
  # Passo 1: Y ~ Exp(1), gerada por inversão (Capítulo 4)
  U1 <- runif(1)
  Y <- -log(1 - U1)
  n_propostas <- n_propostas + 1

  # Passo 2: o uniforme que decide a aceitação
  U <- runif(1)

  # Passo 3: a razão f(Y) / (c * g(Y)) se simplifica para exp(-(Y-1)^2 / 2)
  if (U <= exp(-(Y - 1)^2 / 2)) {
    modulos <- c(modulos, Y)
  }
}

# Passo 4: cada módulo recebe um sinal + ou - com probabilidade 1/2
sinais <- sample(c(-1, 1), size = B, replace = TRUE)
X <- sinais * modulos

cat("Candidatos gerados:", n_propostas, "\n")
Candidatos gerados: 6610
Mostrar código 
cat("Taxa de aceitação observada:", round(B / n_propostas, 3), "\n")
Taxa de aceitação observada: 0.756
Mostrar código 
cat("Taxa de aceitação teórica (1/c):", round(1 / cte, 3), "\n")
Taxa de aceitação teórica (1/c): 0.76
Mostrar código 
cat("Média e desvio padrão amostrais:", round(mean(X), 3), round(sd(X), 3), "\n")
Média e desvio padrão amostrais: 0 1.01
Mostrar código 
df <- data.frame(x = X)

ggplot(df, aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)), bins = 40,
                 fill = "lightcoral", color = "black") +
  stat_function(fun = dnorm, color = "red", linewidth = 1) +
  labs(title = "Normal padrão gerada por rejeição",
       x = "Valor de X", y = "Densidade") +
  theme_minimal()

Mostrar código 
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

np.random.seed(42)

# A constante c só é usada no relatório final: o critério do passo 3 já está
# escrito na forma simplificada
cte = np.sqrt(2 * np.exp(1) / np.pi)

B = 5000             # quantos valores queremos gerar
modulos = []         # guarda os valores aceitos de |X|
n_propostas = 0

while len(modulos) < B:
    # Passo 1: Y ~ Exp(1), gerada por inversão (Capítulo 4)
    U1 = np.random.uniform(0, 1)
    Y = -np.log(1 - U1)
    n_propostas += 1

    # Passo 2: o uniforme que decide a aceitação
    U = np.random.uniform(0, 1)

    # Passo 3: a razão f(Y) / (c * g(Y)) se simplifica para exp(-(Y-1)^2 / 2)
    if U <= np.exp(-(Y - 1)**2 / 2):
        modulos.append(Y)

modulos = np.array(modulos)

# Passo 4: cada módulo recebe um sinal + ou - com probabilidade 1/2
sinais = np.random.choice([-1, 1], size=B)
X = sinais * modulos

print("Candidatos gerados:", n_propostas)
Candidatos gerados: 6558
Mostrar código 
print("Taxa de aceitação observada:", round(B / n_propostas, 3))
Taxa de aceitação observada: 0.762
Mostrar código 
print("Taxa de aceitação teórica (1/c):", round(1 / cte, 3))
Taxa de aceitação teórica (1/c): 0.76
Mostrar código 
print("Média e desvio padrão amostrais:",
      round(X.mean(), 3), round(X.std(ddof=1), 3))
Média e desvio padrão amostrais: 0.008 1.002
Mostrar código 
# Malha usada só para desenhar a densidade teórica
grade = np.linspace(-4, 4, 400)

plt.figure(figsize=(8, 5))
plt.hist(X, bins=40, density=True, color="lightcoral", edgecolor="black")
plt.plot(grade, norm.pdf(grade), color="red", linewidth=2)
plt.title("Normal padrão gerada por rejeição")
plt.xlabel("Valor de X")
plt.ylabel("Densidade")
plt.show()

8.5 Por que o método funciona

Proposição

Suponha que \(g(x) > 0\) sempre que \(f(x) > 0\) e que \(f(x)/g(x) \leq c\) para todo \(x\) com \(f(x) > 0\). Então:

  1. o valor \(X\) devolvido pelo método da rejeição tem densidade \(f\);

  2. o número \(N\) de candidatos gerados até o primeiro aceite tem distribuição geométrica com probabilidade de sucesso \(1/c\). Em particular, \(\mathbb{E}[N] = c\).

Seja \(Y\) o candidato gerado em um passo qualquer do algoritmo e \(U\) o uniforme correspondente. Condicionando no valor de \(Y\), a probabilidade de aceitar é

\[ \mathbb{P}\left(U \leq \frac{f(Y)}{c\, g(Y)} \,\Big|\, Y = y\right) = \frac{f(y)}{c\, g(y)}, \]

porque \(\mathbb{P}(U \leq a) = a\) para \(a \in [0,1]\). Logo, para qualquer \(x\),

\[ \begin{aligned} \mathbb{P}\left(Y \leq x,\ \text{aceitar } Y\right) &= \int_{-\infty}^{x} \frac{f(y)}{c\, g(y)}\, g(y)\, dy \\ &= \frac{1}{c} \int_{-\infty}^{x} f(y)\, dy \;=\; \frac{F(x)}{c}, \end{aligned} \]

em que \(F\) é a f.d.a. do alvo. Fazendo \(x \to \infty\) e usando que \(F(\infty) = 1\),

\[ \mathbb{P}(\text{aceitar}) = \frac{1}{c}. \]

Como os passos do algoritmo são independentes e cada um aceita com probabilidade \(1/c\), o número \(N\) de candidatos até o primeiro aceite tem distribuição geométrica com probabilidade de sucesso \(1/c\), e portanto \(\mathbb{E}[N] = c\). Isso prova (ii).

Para (i), o valor devolvido é o candidato do passo em que houve aceite, ou seja, \(Y\) condicionado ao aceite. Pela definição de probabilidade condicional,

\[ \mathbb{P}(X \leq x) = \mathbb{P}(Y \leq x \mid \text{aceitar}) = \frac{\mathbb{P}(Y \leq x,\ \text{aceitar})}{\mathbb{P}(\text{aceitar})} = \frac{F(x)/c}{1/c} = F(x). \]

Portanto \(X\) tem f.d.a. \(F\), isto é, densidade \(f\). \(\square\)

Repare que a demonstração usa \(f\) e \(g\) apenas através da razão \(f/g\). Isso tem uma consequência prática grande, explorada no Exercício 3: se conhecemos o alvo apenas a menos de uma constante de normalização — sabemos calcular \(f^*\), mas não a constante \(k\) com \(f = k f^*\) — o método continua funcionando. Basta achar \(c'\) com \(f^*(x) \leq c' g(x)\) e usar o critério \(U \leq f^*(Y)/(c'\,g(Y))\), que é idêntico ao original com \(c = k\,c'\).

8.6 Eficiência e escolha da distribuição proposta

O item (ii) diz que, para gerar um valor do alvo, o algoritmo gasta em média \(c\) candidatos; equivalentemente, a fração aceita é \(1/c\). Toda a eficiência do método está nessa única constante, e \(c\) depende só de \(f\) e \(g\) — não do tamanho da amostra que queremos gerar.

Na leitura geométrica da seção anterior, \(1/c\) é uma razão de áreas: a região sob \(f\) tem área \(1\), a região sob \(c\,g\) tem área \(c\), e aceitamos os pontos que caem na primeira. Quanto mais justo o “envelope” \(c\,g\) ficar em volta de \(f\), menos área desperdiçada e mais eficiente o algoritmo. No Exemplo 1, \(c \approx 2{,}11\): o retângulo de altura \(2{,}11\) tem mais que o dobro da área sob \(f\), e por isso rejeitamos mais da metade dos candidatos. No Exemplo 2 a proposta exponencial acompanha bem melhor o formato do alvo, e \(c \approx 1{,}32\).

Há duas armadilhas frequentes na escolha de \(g\).

Atenção: a proposta precisa ter caudas mais pesadas que o alvo

Se \(g\) vai a zero mais rápido que \(f\) quando \(|x| \to \infty\), a razão \(f/g\) explode e não existe constante \(c\) válida. É o que acontece, por exemplo, ao tentar usar uma normal como proposta para uma Cauchy: a densidade da normal decai como \(e^{-x^2/2}\), e a da Cauchy apenas como \(1/x^2\). Na direção contrária — Cauchy como proposta para gerar uma normal — a razão é limitada e o método funciona bem (Exercício 6).

Atenção: uma proposta espalhada demais custa caro

Mesmo quando \(c\) existe, ele pode ser grande. No Exemplo 1, se usássemos \(Y \sim \text{Unif}(0,3)\) — que também cobre o suporte de \(f\) — teríamos \(g(x) = 1/3\) e \(c = 3 \times 135/64 \approx 6{,}33\): o mesmo resultado a um custo três vezes maior, porque dois terços dos candidatos cairiam fora de \((0,1)\) e seriam rejeitados na certa.

Por fim, nem sempre o máximo de \(f/g\) sai no papel. Quando não sair, calcule a razão em uma grade fina de pontos, ou use uma rotina de otimização numérica (optimize em R, scipy.optimize.minimize_scalar em Python), como pede o Exercício 7.

8.7 Exercícios

Exercício 1. Seja \(X \sim \text{Gama}(3/2, 1)\), com densidade

\[ f(x) = \frac{1}{\Gamma(3/2)} x^{1/2} e^{-x}, \quad x > 0. \]

Use como proposta a distribuição \(\text{Exp}(2/3)\).

Observação: o parâmetro \(2/3\) foi escolhido para que a proposta tenha a mesma média do alvo: \(\mathbb{E}[X] = 3/2 = \mathbb{E}[Y]\).

  1. Mostre que a razão \(f(x)/g(x)\) é máxima em \(x = 3/2\) e calcule a constante \(c\).

  2. Escreva um pseudo-algoritmo para simular um valor de \(X\).

  3. Implemente uma função que gere um valor de \(X\) e devolva também quantos candidatos foram necessários. Compare esse número com \(\mathbb{E}[N] = c\).

  4. Implemente uma função que gere uma amostra de tamanho \(B\) de \(X\), e use-a com \(B = 5000\).

  5. Qual foi o número médio de candidatos por valor gerado? Compare com \(c\).

  6. Compare o histograma dos valores simulados com a densidade de \(X\) (dgamma em R, scipy.stats.gamma.pdf em Python).

Exercício 2. Queremos gerar valores da densidade

\[ f(x) = \frac{1}{2}\sin(x), \quad 0 \leq x \leq \pi. \]

  1. Usando como proposta a \(\text{Unif}(0, \pi)\), isto é, \(g(x) = 1/\pi\), encontre a menor constante \(c\) tal que \(f(x) \leq c\, g(x)\) para todo \(x \in [0, \pi]\).

  2. Implemente o método e gere \(B = 5000\) valores. Compare o histograma com a densidade teórica e a taxa de aceitação observada com \(1/c\).

  3. Considere agora a proposta

\[ g_2(x) = \frac{6x(\pi - x)}{\pi^3}, \quad 0 \leq x \leq \pi, \]

que é uma \(\text{Beta}(2,2)\) reescalada para o intervalo \([0,\pi]\). Verifique que \(g_2\) integra \(1\) e mostre que a razão \(f(x)/g_2(x)\) é máxima em \(x = \pi/2\), com \(c_2 = \pi/3\).

  1. Compare \(1/c\) com \(1/c_2\). Qual das duas propostas é mais eficiente, e por quê? Olhe o formato das duas propostas ao lado do de \(f\).

  2. Para usar \(g_2\) seria preciso saber simular dela. Explique por que isso não é imediato, e sugira um caminho (dica: qual método deste livro gera uma \(\text{Beta}(2,2)\)?).

Exercício 3. Muitas vezes conhecemos a densidade alvo apenas a menos de uma constante de normalização: sabemos calcular \(f^*\), mas não a constante \(k\) tal que \(f = k f^*\). Considere

\[ f^*(x) = e^{-x} \sin^2(x), \quad x > 0, \]

e use como proposta a \(\text{Exp}(1)\), ou seja, \(g(x) = e^{-x}\) para \(x > 0\).

  1. Mostre que \(f^*(x)/g(x) = \sin^2(x)\) e conclua que \(c' = 1\) serve.

  2. Escreva o pseudo-algoritmo usando apenas \(f^*\), \(g\) e \(c'\), e explique por que ele não depende de \(k\).

  3. Implemente o algoritmo e gere \(B = 5000\) valores. Construa o histograma de densidade.

  4. Calcule \(k\) exatamente, usando que \(\sin^2 x = (1 - \cos 2x)/2\) e que \(\int_0^\infty e^{-x}\cos(2x)\, dx = 1/5\). Sobreponha \(f = k f^*\) ao histograma do item (c).

  5. (Desafio) Pela proposição, a taxa de aceitação estima \(1/c = 1/(k c')\). Use a taxa observada no item (c) para construir uma estimativa \(\hat{k}\) e compare com o valor exato. Repare que acabamos de estimar uma integral sem nunca tê-la calculado.

Exercício 4. Seja \(Z \sim N(0,1)\) e considere como alvo a distribuição de \(Z\) condicionada a \(Z \geq 2\), cuja densidade é

\[ f(x) = \frac{\varphi(x)}{1 - \Phi(2)}, \quad x \geq 2, \]

onde \(\varphi\) e \(\Phi\) são a densidade e a f.d.a. da \(N(0,1)\).

  1. Usando a própria \(N(0,1)\) como proposta, mostre que \(f(x)/\varphi(x) = 1/(1 - \Phi(2))\) para \(x \geq 2\) e \(0\) caso contrário. Conclua que \(c = 1/(1 - \Phi(2))\) e calcule seu valor (pnorm em R, scipy.stats.norm.cdf em Python).

  2. Mostre que, com esse \(c\), a probabilidade de aceitação vale \(1\) quando \(y \geq 2\) e \(0\) caso contrário — ou seja, o algoritmo se reduz a “gere normais até obter uma maior ou igual a 2”. Implemente e gere \(B = 2000\) valores, registrando o número médio de candidatos.

  3. Considere agora a proposta exponencial deslocada, com densidade \(g_\lambda(x) = \lambda e^{-\lambda(x - 2)}\) para \(x \geq 2\), que se simula fazendo \(Y = 2 + \text{Exp}(\lambda)\). Mostre que \(f(x)/g_\lambda(x) \propto e^{-x^2/2 + \lambda x}\) e que essa razão é limitada para todo \(\lambda > 0\). Onde está o máximo quando \(\lambda \geq 2\)?

  4. Calcule \(c(\lambda)\) para \(\lambda\) em uma grade de \(2\) a \(5\) e faça um gráfico de \(c(\lambda)\) contra \(\lambda\). Implemente o método com o \(\lambda\) que minimiza \(c\) e compare a taxa de aceitação com a do item (b).

  5. (Desafio) Mostre analiticamente que o \(\lambda\) ótimo é \(\lambda^* = 1 + \sqrt{2}\). O que esse exercício sugere sobre usar rejeição para amostrar de eventos raros?

Exercício 5. Este exercício explora diretamente a leitura geométrica do método: em vez de uma densidade na reta, vamos sortear um ponto uniformemente em uma região do plano.

Queremos gerar um ponto \((X_1, X_2)\) uniformemente distribuído no disco unitário \(D = \{(x_1,x_2): x_1^2 + x_2^2 \leq 1\}\), isto é, com densidade conjunta \(f(x_1,x_2) = 1/\pi\) em \(D\). Como proposta, sorteamos um ponto uniforme no quadrado \([-1,1]^2\), cuja densidade conjunta é \(g(x_1,x_2) = 1/4\).

  1. Mostre que \(f/g = 4/\pi\) dentro do disco e \(0\) fora dele, e conclua que \(c = 4/\pi\). Verifique que o critério do passo 3 se reduz a “aceite se o ponto caiu dentro do disco”.

  2. Implemente o algoritmo, gere \(B = 5000\) pontos e faça um gráfico de dispersão.

  3. A taxa de aceitação é \(1/c = \pi/4\). Use isso ao contrário: a partir da taxa observada, construa um estimador \(\hat{\pi}\) e compare com o valor de \(\pi\).

  4. (Desafio) O mesmo raciocínio vale em \(d\) dimensões: sorteie um ponto no cubo \([-1,1]^d\) e aceite se ele cair na bola unitária. O volume da bola é \(V_d = \pi^{d/2}/\Gamma(d/2 + 1)\), de modo que a taxa de aceitação é \(V_d/2^d\). Calcule essa taxa para \(d = 2, 5, 10\) e \(20\). O que isso diz sobre usar rejeição em dimensão alta?

Exercício 6. Vamos verificar a afirmação sobre caudas feita no texto. Sejam \(f\) a densidade da \(N(0,1)\) e \(g\) a densidade da Cauchy padrão,

\[ g(x) = \frac{1}{\pi(1 + x^2)}, \quad x \in \mathbb{R}. \]

  1. Mostre que

\[ \frac{f(x)}{g(x)} = \sqrt{\frac{\pi}{2}}\,(1 + x^2)\, e^{-x^2/2} \]

e que essa razão é máxima em \(x = \pm 1\), com \(c = \sqrt{2\pi/e} \approx 1{,}52\).

  1. Sabendo que a f.d.a. da Cauchy é \(G(x) = \frac{1}{2} + \frac{1}{\pi}\arctan(x)\), mostre que \(G^{-1}(u) = \tan\!\left(\pi(u - 1/2)\right)\) e escreva o pseudo-algoritmo completo (inversão para a proposta, rejeição para o alvo).

  2. Implemente e gere \(B = 5000\) valores. Compare o histograma com a densidade da normal e a taxa de aceitação com \(1/c\).

  3. Troque os papéis: alvo Cauchy, proposta \(N(0,1)\). Mostre que a razão \(f(x)/g(x)\) é ilimitada e que, portanto, não existe \(c\) válido. Faça um gráfico dessa razão para \(x \in [0, 6]\) para visualizar o problema.

  4. Compare a eficiência do item (c) com a do Exemplo 2 do capítulo. Qual das duas propostas se ajusta melhor à normal?

Exercício 7. Nos exemplos do capítulo conseguimos maximizar \(f/g\) no papel. Aqui não. Considere a densidade conhecida a menos de constante

\[ f^*(x) = e^{-x^2/2}\left(\sin^2(6x) + 3x^2\cos^2(x) + 1\right), \quad x \in \mathbb{R}, \]

e use como proposta uma \(N(0, 2^2)\), isto é, uma normal de média \(0\) e desvio padrão \(2\).

  1. Faça um gráfico de \(f^*\) no intervalo \([-6, 6]\). Quantos modos ela tem? Convença-se de que inverter a f.d.a. aqui está fora de questão.

  2. Escreva a razão \(r(x) = f^*(x)/g(x)\) e mostre que ela é limitada. (Dica: compare o expoente \(-x^2/2\) da alvo com o \(-x^2/8\) da proposta.)

  3. Calcule \(r(x)\) em uma grade fina de \([-10, 10]\) e tome \(c'\) como o máximo observado. Confirme o resultado com optimize (R) ou scipy.optimize.minimize_scalar (Python).

  4. Implemente o método e gere \(B = 5000\) valores. Sobreponha ao histograma a curva \(f^*\) reescalada de modo a integrar \(1\) (você pode obter a constante por integração numérica, com integrate em R ou scipy.integrate.quad em Python).

  5. Qual a taxa de aceitação observada? Explique por que a \(N(0,1)\) não serviria como proposta aqui.

  6. (Desafio) Se a grade do item (c) fosse grosseira, \(c'\) poderia ficar abaixo do máximo verdadeiro. O que aconteceria com a distribuição gerada? (Compare com o Exercício 7(d) do capítulo anterior.)

Exercício 8. (Desafio) Este exercício formaliza a leitura geométrica do método. Seja \(f\) uma densidade e

\[ A = \{(x, v) \in \mathbb{R}^2 : 0 \leq v \leq f(x)\} \]

a região sob o seu gráfico, que tem área \(\int f(x)\, dx = 1\).

  1. Suponha que \((X, V)\) seja um ponto sorteado uniformemente em \(A\). Mostre que a densidade marginal de \(X\) é \(f\). (Dica: a densidade conjunta de \((X,V)\) é constante e igual a \(1\) em \(A\); integre em \(v\).)

  2. No algoritmo do capítulo, seja \(Y\) o candidato e \(V = U \cdot c\, g(Y)\). Mostre que \((Y, V)\) é uniforme na região sob \(c\,g\).

  3. Conclua que, condicionado ao aceite, o par \((Y, V)\) é uniforme em \(A\) e que, portanto, o valor devolvido tem densidade \(f\) — uma demonstração alternativa da parte (i) da proposição.

  4. Interprete \(1/c\) como razão entre as áreas de \(A\) e da região sob \(c\,g\).

  5. Volte ao Exemplo 2 e considere a proposta \(\text{Exp}(\lambda)\), com \(\lambda > 0\) qualquer. Mostre que

\[ c(\lambda) = \sqrt{\frac{2}{\pi}}\, \frac{e^{\lambda^2/2}}{\lambda}, \]

e que essa função é minimizada em \(\lambda = 1\) — ou seja, a escolha feita no capítulo é a melhor possível dentro da família exponencial. Faça o gráfico de \(c(\lambda)\) para confirmar.