Motivação

Em Inferência, o processo de estimação muitas vezes envolve a necessidade do cálculo de integrais não triviais; por exemplo, no cálculo de distribuições preditivas e momentos. Para estes casos, é conveniente usar métodos para aproximar integrais complexas (p.ex., de alta dimensão) e/ou que não possuem soluções analíticas.

Importante: Sempre que possível devemos utilizar soluções exatas, isto é, não aproximadas, se elas existem.

Método de Monte Carlo Simples

Suponha que queremos calcular a integral \[\displaystyle I = \int_0^1 g(x)dx.\] De fato, se \(U\) é uma v.a. uniformemente distribuída no intervalo \((0,1)\), podemos escrever \[\displaystyle I = \int_0^1 g(x)dx = E[g(U)].\]

Se \(U_1,...,U_n\) formam uma amostra i.i.d. de \(U\), \(g(U_1),...,g(U_n)\) também são i.i.d., com média \(E[g(U)]\); portanto, a Lei Forte dos Grandes Números garante que, com probabilidade 1,
\[\displaystyle \hat{I} = \frac{1}{n} \sum_{i=1}^n g(U_i) \stackrel{n \rightarrow \infty}{\longrightarrow} E[g(U)] = I.\] Portanto, podemos aproximar a integral \(I\) gerando um grande número de valores \(u_i\) de \(U \sim \text{U}(0,1)\) e calculando a média dos valores \(g(u_i)\) calculados ao longo da amostra.

Esta aproximação é chamada de Aproximação de Monte Carlo e possui usos um pouco mais gerais, conforme veremos a seguir.

A ideia do método é escrever a integral que se deseja calcular como um valor esperado, ao longo de valores dos quais sabemos simular.

Suponha que queremos calcular \[I = \int_a^b g(x)dx.\]

Note que \[\displaystyle I = \int_a^b g(x)dx = \int_a^b (b-a)g(x)\frac{1}{b-a}dx = (b-a)E[g(X)], \] onde \(X \sim \text{U}(a,b)\) (ou seja, transformamos o problema matemático de avaliar \(I\) em um problema estatístico de avaliar uma esperança).

Se temos uma amostra aleatória de tamanho \(n\) da \(U(a,b)\), digamos \(x_1,...,x_n\), podemos estimar \(I\) como \[\displaystyle \hat{I} = (b-a) \frac{1}{n}\sum_{i=1}^ng(X_i).\]

De fato, note que \(\hat{I}\) é uma estimativa não-viesada de \(I\), pois \[\displaystyle E[\hat{I}] = (b-a)\frac{1}{n}\sum_{i=1}^n E[g(X_i)] = (b-a)\frac{1}{n}n E[g(X)] = (b-a) E[g(X)] = I.\]

Algoritmo

\[I = \int_a^b g(x)dx\]

  1. Gere uma amostra \(x_1,...,x_n\) de \(U(a,b)\).

  2. calcule \(g(x_i)\), \(i=1,...,n\).

  3. Aproxime \(I\) por \(\displaystyle \hat{I} =\frac{1}{n} \sum_{i=1}^n g(x_i)\).

Exemplo 1

Calcule \(\displaystyle \int_1^3 e^{-x}dx\).

Solução:

\[I=\displaystyle \int_1^3 e^{-x}dx = (3-1)\displaystyle \int_1^3 e^{-x} \frac{1}{3-1}dx = 2\; E[e^{-X}], \qquad \mbox{onde } X \sim \text{U}(1,3).\]

Algoritmo 1

  1. Gere \(n\) valores da distribuição \(U(1,3)\): \(x_1,...,x_n\).

  2. Calcule \(g(x_i)=e^{-x_i}\), \(i=1,...,n\).

  3. Aproxime \(I\) por \(\displaystyle \hat{I} = (3-1) \frac{1}{n} \sum_{i=1}^n g(x_i) = \frac{2}{n} \sum_{i=1}^n e^{-x_i}\).

g = function(x){
  exp(-x)
}
  
I.chapeu = function(n,a,b,g){
  x = runif(n,a,b)
  I = (b-a)*mean(g(x))
  return(I)
}
  
n = 100
a = 1
b = 3
I.est = I.chapeu(n,a,b,g)
I.est
## [1] 0.3072662

Exercício 1

Observe que \(\displaystyle \int_1^3 e^{-x}dx = P(X\leq3)-P(X\leq1) = 0,3181\), pois \(e^{-x}\) (para \(x>0\)) é a função de densidade de probabilidade da distribuição exponencial com média 1.

Conhecendo o valor verdadeiro da integral, faça os seguintes exercícios:

  1. Compare o valor obtido de forma aproximada com o valor exato da integral.

  2. Aumente o tamanho da amostra (\(n\)) e repita a comparação.

  3. Repita o cálculo 100 vezes com \(n\) fixo e faça um histograma dos resultados encontrados.

  4. Repita o exercício anterior com 1000 cálculos e observe o que ocorre.

Observação 1

A Lei Forte dos Grandes Números garante que \[\displaystyle \hat{I} = \frac{1}{n} \sum_{i=1}^n g(X_i) \stackrel{n \rightarrow \infty}{\longrightarrow} E[g(X)] = I \qquad \mbox{q.c.};\] além disso, a variância do estimador \(\hat{I}\) é \[\displaystyle Var(\hat{I}) = Var\left(\frac{1}{n} \sum_{i=1}^n g({X}_i)\right) = \frac{1}{n^2} \sum_{i=1}^n Var\left( g({X}_i)\right);\] e, portanto, podemos melhorar a aproximação tanto quanto se queira, bastando aumentar \(n\).

Generalizando o Algoritmo

Note que a estratégia de olhar a integral como a esperança de uma função avaliada em valores de uma a.a.s. uniforme não se aplica diretamente para integrais definidas em intervalos ilimitados.

Uma possível solução é generalizar o método, assumindo a geração de pontos através de outros modelos probabilísticos, cuja imagem seja mais adequada ao problema.

Suponha que \(g(x)\) pode ser decomposta da forma \(g(x) = h(x)p(x)\), onde \(p(x)\) representa uma modelo probabilístico do qual sabemos simular. Neste caso, podemos utilizar o algoritmo a seguir para aproximar a integral

\[I = \int g({x})d{x} = \int h({x})p(x)d{x}.\]

  1. Gere uma amostra \({x}_1,...,{x}_n\) de \(p({x})\).

  2. Calcule \(h({x}_i)\), \(i=1,...,n\).

  3. Aproxime \(I\) por \(\displaystyle \hat{I} =\frac{1}{n} \sum_{i=1}^n h({x}_i)\).

Exemplo 2

Calcule a função gama, denotada por \(\Gamma(\alpha)\) para diferentes valores de \(\alpha\).

\[I=\displaystyle \Gamma(\alpha) = \int_0^\infty x^{\alpha-1}e^{-x}dx = E[X^{\alpha-1}], \qquad \mbox{onde } X \sim \text{Exp}(1).\]

Solução:

Algoritmo 2

  1. Gere \(n\) valores da distribuição \(\text{Exp}(1)\): \(x_1,...,x_n\).

  2. Calcule \(h(x_i)=x_i^{\alpha-1}\), \(i=1,...,n\).

  3. Aproxime \(I\) por \(\displaystyle \hat{I} = \frac{1}{n} \sum_{i=1}^n h(x_i) = \frac{1}{n} \sum_{i=1}^n x_i^{\alpha-1}\).

h = function(x,alpha){
  x^(alpha-1)
}

I.chapeu = function(n,alpha,h){
  x = rexp(n)
  I = mean(h(x,alpha))
  return(I)
}
  
n = 100
alpha = 1.7
I.est = I.chapeu(n,alpha,h)
I.est  
## [1] 0.8536099
gamma(alpha)
## [1] 0.9086387

Observação 2

Note que a definição vale tanto no caso univariado quanto no caso multivariado no qual \(\boldsymbol{x}=(x_1,...,x_k)\) e \(p(\boldsymbol{x})\) representa a distribuição conjunta de todas as componentes de \(\boldsymbol{x}\).

Exemplo 3

(Ross, 2022) Estimação de \(\pi\)

Considere o quadrado \(D=[-1,1]\times[-1,1]\), centrado na origem e, dentro desse quadrado, o círculo unitário

\[ C=\{(x,y)\in D:x^2+y^2\leq 1\}. \]

Deseja-se utilizar o Método de Monte Carlo simples para estimar a integral dupla \[ I= \int_{-1}^{1} \int_{-1}^{1} \mathbf{1}_{{x^2+y^2\leq1}} dxdy, \] em que \(\mathbf{1}_{A}\) é a função indicadora do evento \(A\): \[ \mathbf{1}_{A}= \begin{cases} 1,&\text{se }A\text{ ocorre},\\ 0,&\text{caso contrário}. \end{cases} \] Determine uma estimativa de \(I\) utilizando amostras aleatórias uniformes e compare o resultado obtido com o valor exato da integral.

De fato, note que \(I\) representa exatamente a área do círculo unitário; portanto, \(I=\pi\). Assim, o problema consiste em utilizar simulação aleatória para obter uma aproximação de \(\pi\).

Solução:

Note que \[ I = \int_{-1}^{1} \int_{-1}^{1} \mathbf{1}_{{x^2+y^2\leq1}} dxdy = 4 \int_{-1}^{1} \int_{-1}^{1} \mathbf{1}_{{x^2+y^2\leq1}} \frac{1}{4} dxdy = E[\mathbf{1}_{{X^2+Y^2\leq1}}], \qquad (X,Y) \sim \text{U}(D). \]

Algoritmo 3

  1. Gere \(n\) valores da distribuição \(\text{U}(D)\): \((x_1,y_1),...,(x_n,y_n)\). Isto é equivalente a simular \(X \sim \text{U}(-1,1)\) e \(Y \sim \text{U}(-1,1)\) independentes.

  2. Calcule \(g(x_i,y_i)=\mathbf{1}_{{x_i^2+y_i^2\leq1}}\), \(i=1,...,n\).

  3. Aproxime \(I\) por \(\displaystyle \hat{I} = \frac{4}{n} \sum_{i=1}^n g(x_i,y_i) = \frac{4}{n} \sum_{i=1}^n \mathbf{1}_{{x_i^2+y_i^2\leq1}}\).

n = 10000

x = runif(n,-1,1)
y = runif(n,-1,1)

indicadora = (x^2 + y^2 <= 1)

pi.mc = 4*sum(indicadora)/n
pi.mc
## [1] 3.1416

O resultado não será exatamente igual a \(\pi\), pois estamos utilizando uma aproximação numérica baseada em uma amostra aleatória, porém o valor obtido deve estar próximo de \[ \pi\approx3.14159265. \]

Exercício 2

Construa uma função que aproxime integrais duplas do tipo \[ I = \int_a^b \int_c^d g(x,y)dydx. \] através do Método de Monte Carlo.

Método de Monte Carlo via Função de Importância

Suponha que queremos calcular

\[\displaystyle I=\int h({x})p({x})d{x},\] mas não sabemos gerar de \(p({x})\).

Neste caso, se sabemos gerar de uma distribuição \(q({x})\) definida no mesmo espaço de variação de \(p({x})\), podemos fazer \[I = \displaystyle \int h({x})p({x})d{x}=\int \frac{h({x})p({x})}{q({x})}q({x})d{x},\] que é a esperança de \(\displaystyle \frac{h({X})p({X})}{q({X})}\) com relação a distribuição \(q({x})\).

A função \(q({x})\) é chamada de função de importância.

Algoritmo

  1. Gere uma amostra \({x}_1,...,{x}_n\) de \(q({x})\).

  2. Calcule \(\displaystyle \frac{h({x}_i)p({x}_i)}{q({x}_i)}\), \(i=1,...,n\), \(i=1,...,n\).

  3. Aproxime \(I\) por \(\displaystyle \hat{I} =\frac{1}{n} \sum_{i=1}^n \frac{h({x}_i)p({x}_i)}{q({x}_i)}\).

Pode-se mostrar que a escolha ótima no sentido de minimizar a variância do estimador consiste em tomar \(q({x}) \propto h({x})p({x})\).

Exemplo 4

Estimar a probabilidade

\[ I=P(X>5),\qquad \mbox{onde }X\sim N(0,1). \]

Solução:

Podemos escrever essa probabilidade como uma integral:

\[ I=\int_5^\infty \frac{1}{\sqrt{2\pi}}e^{-x^2/2} dx = \int_{-\infty}^\infty \mathbb{I}(x>5)\frac{1}{\sqrt{2\pi}}e^{-x^2/2} dx = E[\mathbb{I}(X>5)], \qquad X \sim \text{N}(0,1). \] Como o evento \(X>5\) é bastante raro, o Método de Monte Carlo simples pode apresentar uma estimativa pouco eficiente. De fato, se gerarmos valores diretamente de \(N(0,1)\), é muito provável que uma parcela muito pequena, ou até nenhuma, das observações esteja na região \(x>5\).

Para melhorar a eficiência, utilizaremos uma Função de Importância. Escolhemos uma distribuição normal deslocada para a região de interesse:

\[ X\sim N(\mu,1), \qquad \text{com} \qquad \mu=5. \]

A densidade da distribuição original é

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

enquanto a densidade da função de importância é

\[ q(x)=\frac{1}{\sqrt{2\pi}}e^{-(x-5)^2/2}. \]

Podemos multiplicar e dividir o integrando por \(q(x)\):

\[\begin{align*} I&=\int_{-\infty}^\infty \mathbb{I}(x>5)\frac{p(x)}{q(x)}q(x) dx = \int_{-\infty}^\infty \mathbb{I}(x>5)\frac{\frac{1}{\sqrt{2\pi}}e^{-x^2/2}}{\frac{1}{\sqrt{2\pi}}e^{-(x-\mu)^2/2}}q(x) dx \\ &= \int_{-\infty}^\infty \mathbb{I}(x>5)e^{(-2\mu x+\mu^2)/2} \frac{1}{\sqrt{2\pi}}e^{-(x-\mu)^2/2} dx = \int_{-\infty}^\infty \mathbb{I}(x>5)e^{-5 x+25/2} \frac{1}{\sqrt{2\pi}}e^{-(x-5)^2/2} dx \\ &= E[\mathbb{I}(X>5)e^{-5 X+25/2}], \qquad X \sim \text{N}(5,1). \end{align*}\]

Algoritmo 4.1 (Método de Monte Carlo Simples)

  1. Gere \(n\) valores da distribuição \(\text{N}(0,1)\): \(x_1,...,x_n\).

  2. Calcule \(h_1(x_i)=\mathbb{I}(x_i>5)\), \(i=1,...,n\).

  3. Aproxime \(I\) por \(\displaystyle \hat{I} = \frac{1}{n} \sum_{i=1}^n h_1(x_i) = \frac{1}{n} \sum_{i=1}^n \mathbb{I}(x_i>5)\).

Algoritmo 4.2 (Método de Monte Carlo via Função de Importância)

  1. Gere \(n\) valores da distribuição \(\text{N}(5,1)\): \(x_1,...,x_n\).

  2. Calcule \(h_2(x_i)=\mathbb{I}(x_i>5)e^{-5 x_i+25/2}\), \(i=1,...,n\).

  3. Aproxime \(I\) por \(\displaystyle \hat{I} = \frac{1}{n} \sum_{i=1}^n h_2(x_i) = \frac{1}{n} \sum_{i=1}^n \mathbb{I}(x_i>5)e^{-5 x_i+25/2}\).

n = 10000

h1 = function(x,mu){
  I(x > 5)
}

I.chapeu.1 = function(n,h){
  x = rnorm(n,mean=0,sd=1)
  I = mean(h(x))
  return(I)
}

I.MCS = I.chapeu.1(n,h1)
I.MCS
## [1] 0
h2 = function(x){ 
  (x > 5)*exp(-5*x + 25/2) 
} 

I.chapeu.2 = function(n,h){
  x = rnorm(n,mean=5,sd=1) 
  I = mean(h(x)) 
  return(I) 
}

I.MCFI = I.chapeu.2(n,h2)
I.MCFI
## [1] 3.015429e-07
# valor exato
pnorm(5,lower.tail=FALSE)
## [1] 2.866516e-07

Referências

ROSS, Sheldon M. Simulation. 6. ed. London: Academic Press, 2022.