Exercício 1

Comparação de Métodos de Monte Carlo

Considere a integral

\[ I= \int_0^1 \frac{e^{-x}}{1+x^2}\,dx. \]

O objetivo deste exercício é comparar diferentes estratégias de Monte Carlo para aproximar a integral acima.

Estime a integral pelo Método de Monte Carlo simples.
Escolha uma distribuição de importância e construa um estimador para \(I\) via Função de Importância.
Estime a integral utilizando variáveis antitéticas.
Estime a integral utilizando uma variável de controle.
Compare todos os resultados com uma aproximação numérica obtida por meio da função integrate().

Para cada método citado

Utilize \(N=10.000\) observações;
repita cada procedimento \(500\) vezes;
compare as distribuições das estimativas obtidas utilizando gráficos e medidas resumo.

Solução:

1. Método de Monte Carlo Simples

\[I=\displaystyle \int_0^1 \frac{e^{-x}}{1+x^2}\,dx = E\left[\frac{e^{-x}}{1+x^2}\right], \qquad \mbox{onde } X \sim \text{U}(0,1).\]

Neste caso, vamos assumir \[\hat{I}_{MC}=\frac{1}{N}\sum_{i=1}^N g(X_i) = \frac{1}{N} \sum_{i=1}^N \frac{e^{-X_i}}{1+X_i^2}.\]

Algoritmo 1

Gere \(N\) valores da distribuição \(U(0,1)\): \(x_1,...,x_N\).
Calcule \(g(x_i)=\frac{e^{-x}}{1+x^2}\), \(i=1,...,N\).
Aproxime \(I\) por \(\displaystyle \hat{I} = \frac{1}{N} \sum_{i=1}^N g(x_i) = \frac{1}{N} \sum_{i=1}^N \frac{e^{-x_i}}{1+x_i^2}\).

n = 10000

g = function(x){
  exp(-x)/(1+x^2)
}
  
I.chapeu.MC = function(n,g){
  x = runif(n)
  I = mean(g(x))
  return(I)
}

I.est.MC = I.chapeu.MC(n,g)
I.est.MC

## [1] 0.5261481

integrate(g,0,1)

## 0.5247971 with absolute error < 5.8e-15

2. Variáveis antitéticas

Como \(1-X\sim U(0,1)\), utilizamos

\[ \hat{I}_A = \frac1N \sum_{i=1}^N \frac{g(X_i)+g(1-X_i)}{2}. \]

I.chapeu.A = function(n,g){
  x = runif(n)
  y = 1-x
  I = mean((g(x)+g(y))/2)
  return(I)
}
  
I.est.A = I.chapeu.A(n,g)
I.est.A

## [1] 0.5249907

integrate(g,0,1)

## 0.5247971 with absolute error < 5.8e-15

3. Variável de controle

Podemos utilizar \(Y \sim U(0,1)\) como variável de controle, pois conhecemos \(E[Y]=\frac12\).

Definimos \[ \hat I_{C} = \frac1N \sum_{i=1}^N [g(X_i)+c(Y_i-1/2)], \qquad c = -\frac{\widehat{Cov}(g(X),Y)}{\widehat{Var}(Y)} \]

I.chapeu.C = function(n,g){
  x = runif(n)
  y = runif(n)
  c = -cov(g(x),y)/var(y)
  I = mean(g(x)+c*(y-0.5))
  return(I)
}
  
I.est.C = I.chapeu.C(n,g)
I.est.C

## [1] 0.5261769

integrate(g,0,1)

## 0.5247971 with absolute error < 5.8e-15

4. Função de Importância

Note que \[g(x)=\frac{e^{-x}}{1+x^2} = \frac{1}{1+x^2} \; e^{-x};\] portanto, \[I=\displaystyle \int_0^1 \frac{e^{-x}}{1+x^2}\,dx = E\left[\frac{1}{1+x^2}\mathbb{I}(X\leq1)\right], \qquad \mbox{onde } X \sim \text{Exp}(1).\]

Note que a integral original está definida no intervalo \([0,1]\); por isso introduzimos a função indicadora \(\mathbb{I}(X\leq1)\).

Neste caso, o estimador será \[ \hat{I}_{F} = \frac{1}{N} \sum_{i=1}^N \frac{1}{1+X_i^2}\mathbb{I}(X_i\leq1) . \]

h = function(x){
  (1/(1+x^2)) * I(x<=1)
}
  
I.chapeu.F = function(n,h){
  x = rexp(n)
  I = mean(h(x))
  return(I)
}
  
I.est.F = I.chapeu.F(n,h)
I.est.F

## [1] 0.5355342

integrate(g,0,1)

## 0.5247971 with absolute error < 5.8e-15

Fazendo o estudo com base nas r = 500 replicações:

# Número de observações em cada replicação
N = 10000

# Número de replicações
nrep = 500

# Matriz para armazenar os resultados
resultados = matrix(
  NA,
  nrow = nrep,
  ncol = 4
)

colnames(resultados) = c(
  "MC.Simples",
  "V.Antitética",
  "V.Controle",
  "F.Importância"
)

for(r in 1:nrep){
  
  # ---------------------------
  # Monte Carlo simples
  # ---------------------------
  
  X = runif(N)
  resultados[r, "MC.Simples"] = mean(g(X))
  
  # ---------------------------
  # Variáveis antitéticas
  # ---------------------------
  
  valores_antiteticos = (g(X) + g(1 - X)) / 2
  resultados[r, "V.Antitética"] = mean(valores_antiteticos)
  
  # ---------------------------
  # Variável de controle
  # ---------------------------
  
  c = -cov(g(X),X)/var(X)
  valores_controle = g(X)+c*(X-0.5)
  resultados[r, "V.Controle"] = mean(valores_controle)
  
  # ---------------------------
  # Função de Importância
  # ---------------------------
  
  Y = rexp(N)
  valores_importancia = h(Y)
  resultados[r, "F.Importância"] = mean(valores_importancia)
}

# Valor de referência
I_exato = integrate(g, 0, 1)$value

# Resumo dos resultados
medias = colMeans(resultados)

variancias <- apply(
  resultados,
  2,
  var
)

resumo = data.frame(
  Metodo = colnames(resultados),
  Media = medias,
  Variancia = variancias,
  Erro = abs(medias - I_exato)
)

cat("Valor exato:", I_exato, "\n\n")

## Valor exato: 0.5247971

resumo

##                      Metodo     Media    Variancia         Erro
## MC.Simples       MC.Simples 0.5248728 6.077728e-06 7.563834e-05
## V.Antitética   V.Antitética 0.5248110 1.113645e-07 1.381126e-05
## V.Controle       V.Controle 0.5248091 1.129885e-07 1.192847e-05
## F.Importância F.Importância 0.5248357 1.882543e-05 3.857835e-05

Comparando graficamente as distribuições das estimativas:

boxplot(
  resultados,
  col = c(
    "aliceblue",
    "lightgreen",
    "lightyellow",
    "lightpink"
  ),
  main = "Comparação dos Métodos de Monte Carlo",
  ylab = "Estimativa da Integral"
)

abline(
  h = I_exato,
  col = "red",
  lwd = 2,
  lty = 2
)

A comparação das variâncias permite avaliar empiricamente o ganho obtido pelas técnicas de redução de variância. Quanto menor a variância das estimativas, maior é a eficiência do método para um mesmo tamanho amostral.

Integração de Monte Carlo - Exercícios