13  Amostragem por Importância

O Método de Monte Carlo, embora funcione, nem sempre é eficiente. Veremos a seguir um exemplo disso.

13.1 Exemplo 1: Estimativa de \(\mathbb{P}(Z > 4.5)\)

Neste exemplo, vamos estimar a probabilidade \(\mathbb{P}(Z > 4.5)\) onde \(Z \sim N(0,1)\). Este é um evento raro, com probabilidade pequena. O valor exato dessa probabilidade é: \[ \log(\mathbb{P}(Z > 4.5)) = -12.59242 \] Aplicaremos o método de Monte Carlo simulando \(n = 100.000\) valores de \(Z \sim N(0,1)\) e calculando a proporção dos valores que são maiores que 4.5.

Mostrar código 
set.seed(0)
# Valor exato de P(Z > 4.5)
valor_exato <- log(pnorm(4.5, 0, 1, lower.tail = FALSE))
cat("Valor exato (log(P(Z > 4.5))):", valor_exato, "\n")
Valor exato (log(P(Z > 4.5))): -12.59242
Mostrar código 
# Simulação Monte Carlo
n <- 100000
z <- rnorm(n, 0, 1)

# Estimativa Monte Carlo de P(Z > 4.5)
p_estimada <- mean(1 * (z > 4.5))
cat("Estimativa Monte Carlo de P(Z > 4.5):", p_estimada, "\n")
Estimativa Monte Carlo de P(Z > 4.5): 0
Mostrar código 
# Logaritmo da estimativa Monte Carlo
log_p_estimada <- log(p_estimada)
cat("Logaritmo da estimativa Monte Carlo:", log_p_estimada, "\n")
Logaritmo da estimativa Monte Carlo: -Inf
Mostrar código 
import numpy as np
from scipy.stats import norm

# Valor exato de P(Z > 4.5)
valor_exato = np.log(norm.sf(4.5))  # sf is the survival function, equivalent to (1 - cdf)
print(f"Valor exato (log(P(Z > 4.5))): {valor_exato}")
Valor exato (log(P(Z > 4.5))): -12.592419735713081
Mostrar código 
# Simulação Monte Carlo
n = 100000
z = np.random.normal(0, 1, n)

# Estimativa Monte Carlo de P(Z > 4.5)
p_estimada = np.mean(z > 4.5)
print(f"Estimativa Monte Carlo de P(Z > 4.5): {p_estimada}")
Estimativa Monte Carlo de P(Z > 4.5): 0.0
Mostrar código 
# Logaritmo da estimativa Monte Carlo
log_p_estimada = np.log(p_estimada)
print(f"Logaritmo da estimativa Monte Carlo: {log_p_estimada}")
Logaritmo da estimativa Monte Carlo: -inf

O método de Monte Carlo exige uma grande quantidade de amostras para fornecer uma estimativa precisa devido à raridade do evento \(Z > 4.5\). A abordagem pode se tornar ineficiente para eventos raros.

13.2 Amostragem por Importância

A Amostragem por Importância é uma técnica que nos permite simular valores \(Y_1, \dots, Y_n\) a partir de uma densidade auxiliar \(g(y)\).

A ideia chave é que, se \(X\) tem densidade \(f(x)\), a esperança \(\mathbb{E}[h(X)]\) pode ser reescrita como:

\[ \mathbb{E}[h(X)] = \int h(x) f(x) dx = \int h(x) \frac{f(x)}{g(x)} g(x) dx = \mathbb{E}\left[ h(Y) \frac{f(Y)}{g(Y)} \right], \] em que \(Y\) tem densidade \(g\). Isso vale desde que \(g\) possua suporte tão ou mais amplo que a região de integração. Assim, podemos estimar \(\mathbb{E}[h(X)]\) usando \(n\) valores simulados de \(Y\) independentemente pela fórmula: \[ \hat\theta_n = \frac{1}{n}\sum_{i=1}^n \frac{f(y_i)}{g(y_i)} h(y_i)= \frac{1}{n}\sum_{i=1}^n w(y_i) h(y_i), \] onde os pesos \(w(y_i):=\frac{f(y_i)}{g(y_i)}\) ajustam a importância de cada amostra, podendo tornar a simulação mais eficiente. De fato, a variância de \(\hat\theta_n\) é \[ \mathrm{Var}(\hat\theta_n) =\frac{1}{n}\mathrm{Var}_g\left[h(Y)\frac{f(Y)}{g(Y)}\right] =\frac{1}{n}\left\{\int h(x)^2\frac{f(x)^2}{g(x)}dx-\theta^2\right\}= \frac{1}{n}\left\{\int h(x)^2\frac{f(x)^2}{g(x)}dx-\theta^2\right\}, \] que pode ser menor que a variância do estimador Monte Carlo usual: \[ \widetilde\theta_n=\frac1n\sum_{i=1}^n h(X_i)\quad (X_i\stackrel{iid}{\sim} f), \qquad \mathrm{Var}(\widetilde\theta_n)=\frac{1}{n}\mathrm{Var}_f[h(X)]=\frac{1}{n}\left\{\int h(x)^2f(x)dx-\theta^2\right\}, \] especialmente quando \(g\) concentra probabilidade nas regiões onde \(|h(x)|f(x)\) é grande. Assim, quando possível, devemos escolher \(g\) com maior massa nas regiões que mais contribuem para a integral—isto é, onde \(|h(x)|f(x)\) é grande. Assim, reduzimos \(\int h(x)^2\frac{f(x)^2}{g(x)}dx\).

13.3 Revistando o Exemplo 1 via Amostragem por Importância

Para melhorar a eficiência no Exemplo 1, vamos agora usar Amostragem por Importância. Em vez de simular diretamente de \(Z \sim N(0,1)\), simulamos de uma distribuição exponencial truncada \(Exp(1)\) a partir de 4.5. A densidade de importância é: \[ g(y) = e^{-(y - 4.5)}, \quad y > 4.5 \]

O estimador de Monte Carlo ponderado é, portanto, \[ \hat{\theta}_n = \frac{1}{n} \sum_{i=1}^n \frac{f(y_i)}{g(y_i)}I(y_i>4.5)=\frac{1}{n} \sum_{i=1}^n \frac{f(y_i)}{g(y_i)}, \] onde \(f(y_i)\) é a densidade da normal \(N(0,1)\).

Mostrar código 
n <- 100
# Gerar n valores da distribuição exponencial truncada
y <- rexp(n) + 4.5

# Estimativa por importância
estimativa_importancia <- mean(dnorm(y) / exp(-(y - 4.5)))
cat("Estimativa por amostragem por importância:", estimativa_importancia, "\n")
Estimativa por amostragem por importância: 3.723979e-06
Mostrar código 
# Logaritmo da estimativa por importância
log_estimativa_importancia <- log(estimativa_importancia)
cat("Logaritmo da estimativa por amostragem por importância:", log_estimativa_importancia, "\n")
Logaritmo da estimativa por amostragem por importância: -12.50072
Mostrar código 
import numpy as np
from scipy.stats import expon, norm

n = 100
# Gerar n valores da distribuição exponencial truncada
y = expon.rvs(scale=1, size=n) + 4.5

# Estimativa por importância
estimativa_importancia = np.mean(norm.pdf(y) / np.exp(-(y - 4.5)))
print(f"Estimativa por amostragem por importância: {estimativa_importancia}")
Estimativa por amostragem por importância: 3.4849152097977753e-06
Mostrar código 
# Logaritmo da estimativa por importância
log_estimativa_importancia = np.log(estimativa_importancia)
print(f"Logaritmo da estimativa por amostragem por importância: {log_estimativa_importancia}")
Logaritmo da estimativa por amostragem por importância: -12.567066844091448

Com essa técnica, concentramos as simulações na região relevante, melhorando a eficiência da estimativa para eventos raros.

13.4 Exemplo 2: Estimativa de \(\mathbb{E}[\sqrt{X}]\) com Amostragem por Importância

Queremos estimar \[ \theta=\mathbb{E}[\sqrt{X}],\quad X\sim \mathrm{Exp}(1). \] Está equivale à seguinte integral: \[ \theta=\int_0^\infty \sqrt{x}e^{-x}dx =\Gamma \Big(\tfrac32\Big)=\frac{\sqrt{\pi}}{2}. \] Para comparação, a variância do MC direto (amostrar \(X\sim\mathrm{Exp}(1)\) e tirar a média de \(\sqrt{X}\)) é \[ \mathrm{Var}(\sqrt{X})=\mathbb{E}[X]-\theta^2=1-\frac{\pi}{4}, \qquad \Rightarrow\quad \mathrm{Var}(\hat\theta_{\text{MC}})=\frac{1}{n}\left( 1-\frac{\pi}{4} \right). \]

Para implementar a amostragel por importância, escolhemos a proposta \(g_\lambda(x)=\lambda e^{-\lambda x}\mathbf{1}{x\ge 0}\) para algum \(0<\lambda<2\) fixo. Com essa proposta, a amostragem por importância fara a média de termos do tipo \[ A_\lambda(x)=\frac{h(x)f(x)}{g_\lambda(x)} =\frac{\sqrt{x}e^{-x}}{\lambda e^{-\lambda x}} =\frac{1}{\lambda}\sqrt{x}e^{-(1-\lambda)x}. \] O estimador por importância é, portanto, \[ \hat\theta_{\text{IS}}=\frac1n\sum_{i=1}^n A_\lambda(X_i),\qquad X_i\stackrel{iid}{\sim}\mathrm{Exp}(\lambda). \]

Iremos agora otimizar \(\lambda\), isto é, escolher \(\lambda\) que retorna a menor variância. Para isso, vamos calcular o segundo momento sob \(g_\lambda\) de \(A_\lambda(X)\): \[ \mathbb{E}_{g_\lambda}[A_\lambda(X)^2] =\frac{1}{\lambda}\int_0^\infty x e^{-(2-\lambda)x}dx =\frac{1}{\lambda(2-\lambda)^2}. \] Logo, \[ \mathrm{Var}(\hat\theta_{\text{IS}}) =\frac{1}{n}\Big(\frac{1}{\lambda(2-\lambda)^2}-\theta^2\Big). \] Minimizando \(1/[\lambda(2-\lambda)^2]\) em \(0<\lambda<2\), obtemos a melhor taxa na família: \[ \lambda^\star=\tfrac{2}{3}. \] Para esta escolha, \[ \mathrm{Var}(\hat\theta_{\text{IS},,\lambda^\star}) =\frac{1}{n}\Big(0{,}84375-\frac{\pi}{4}\Big). \] um ganho de aproximadamente \(0{,}2146/0{,}05835\approx 3{,}68\) vezes em relação ao MC direto.

Finalmente, o estimador por amostragem por importância consiste em:

  1. Escolha \(\lambda\) (use \(\lambda=2/3\) para ficar próximo do ótimo).
  2. Gere \(X_1,\dots,X_n\stackrel{iid}{\sim}\mathrm{Exp}(\lambda)\).
  3. Calcule \(A_i=\dfrac{\sqrt{X_i}e^{-X_i}}{\lambda e^{-\lambda X_i}}\).
  4. Estime \(\hat\theta=\dfrac1n\sum_{i=1}^n A_i\).
  5. Erro-padrão: \(\widehat{\mathrm{SE}}=\sqrt{\widehat{\mathrm{Var}}(A)/n}\), onde \(\widehat{\mathrm{Var}}(A)\) é a variância amostral dos \(A_i\).
Mostrar código 
set.seed(1)

theta_true <- sqrt(pi)/2
n <- 1000

is_exp <- function(n, lambda) {
  x  <- rexp(n, rate = lambda)
  A  <- sqrt(x) * exp(-x) / (lambda * exp(-lambda * x))
  th <- mean(A)
  se <- sqrt(var(A) / n)
  list(theta = th, se = se)
}

# Comparando lambdas: 1 (referência) vs 2/3 (ótimo)
out_1   <- is_exp(n, lambda = 1)
out_opt <- is_exp(n, lambda = 2/3)

cat(sprintf("Valor exato           : %.6f\n", theta_true))
Valor exato           : 0.886227
Mostrar código 
cat(sprintf("IS Exp(1):    theta=%.6f  SE≈%.5f\n", out_1$theta,   out_1$se))
IS Exp(1):    theta=0.906481  SE≈0.01448
Mostrar código 
cat(sprintf("IS Exp(2/3):  theta=%.6f  SE≈%.5f\n", out_opt$theta, out_opt$se))
IS Exp(2/3):  theta=0.882012  SE≈0.00763
Mostrar código 
import numpy as np

rng = np.random.default_rng(1)
theta_true = np.sqrt(np.pi)/2
n = 1000

def is_exp(n, lam, rng):
    x = rng.exponential(scale=1/lam, size=n)
    A = (np.sqrt(x) * np.exp(-x)) / (lam * np.exp(-lam*x))
    theta = np.mean(A)
    se = np.sqrt(np.var(A, ddof=1)/n)
    return {"theta": theta, "se": se}

out_1   = is_exp(n, 1.0, rng)
out_opt = is_exp(n, 2/3, rng)

print(f"Valor exato          : {theta_true:.6f}")
Valor exato          : 0.886227
Mostrar código 
print(f"IS Exp(1):   theta={out_1['theta']:.6f}  SE≈{out_1['se']:.5f}")
IS Exp(1):   theta=0.888850  SE≈0.01478
Mostrar código 
print(f"IS Exp(2/3): theta={out_opt['theta']:.6f}  SE≈{out_opt['se']:.5f}")
IS Exp(2/3): theta=0.891175  SE≈0.00773

13.5 Exercícios

Exercício 1. Usando a Amostragem por Importância (AI), estime a probabilidade de que uma variável aleatória \(Z \sim N(0,1)\) seja maior que 5.0. Compare os resultados obtidos com a amostragem direta usando Monte Carlo simples e discuta a eficiência da Amostragem por Importância em relação ao método direto.

Exercício 2. Estime o valor da integral \(\int_0^2 \sqrt{x} e^{-x} dx\) utilizando a AIcom duas distribuições de importância diferentes. Justifique a escolha das distribuições e compare a eficiência de cada uma.

Exercício 3. Utilize a AI para estimar \(\mathbb{E}[\log(1 + X)]\), onde \(X \sim Exp(1)\). Escolha uma densidade auxiliar apropriada para melhorar a eficiência da estimativa e justifique sua escolha. Compare os resultados com a aproximação via Monte Carlo simples.

Exercício 4 Estime \(\mathbb{P}(Z>6)\) com \(Z\sim\mathcal N(0,1)\) usando AI com proposta \(g_\mu=\mathcal N(\mu,1)\).

  1. Mostre que o estimador é \(\hat p=\frac1n\sum I(Y_i>6)\frac{\phi(Y_i)}{\phi(Y_i-\mu)}\), \(Y_i\sim\mathcal N(\mu,1)\).

  2. Encontre o melhor \(\mu\).

  3. Compare MC ingênuo, AI com alguns valores diferentes de \(\mu\) (do seu exemplo) e AI com \(\mu^\star\) em termos de erro-padrão