Integração Monte Carlo
Prof. Lupercio F. Bessegato
27/07/2021
Introdução
As duas maiores classes de problemas numéricos em inferência estatística relacionam-se com otimização e integração. Em geral, os problemas associados com a solução de equações implícitas podem ser reformulados como un problema de otimização. Embora não haja classificações estritas dessas classes de problemas, os problemas de otimização estão geralemnte associados com a abordagem de verossimilhança e, os de integração, com a abordagem bayesiana.
Entretanto, nem sempre é possível calcular analiticamente os estimadores associados a um dado problema (máxima verossimilhança, método dos momentos, etc.). Assim, seja qual for o tipo de inferência estatística, freqüentemente é necessário considerar soluções numéricas. Uma solução geral é usar a simulação da distribuição verdadeira ou de alguma distribuição substituta para calcular as quantidades de interesse. A possibilidade de produzir um número quase infinito de variáveis aleatórias distribuídas de acordo com uma dada distribuição, dá acesso ao uso de resultados frequentistas e assintóticos muito mais facilmente do que nos cenários inferenciais usuais. Pode-se assim aplicar resultados probabilísticos como a Lei dos Grandes Números ou o Teorema do Limite Central, uma vez que é possível avaliar a convergência dos métodos de simulação (a qual é equivalente aos limites determinísticos usados pelas abordagens numéricas).
Os métodos numéricos, como o de Simpson e as quadraturas gaussianas, são uma alternativa óbvia ao uso de métodos de simulação para aproximar integrais da forma \(\int_ {R_X} h (x) f(x) \operatorname{d}\! x\), em que f é uma função de densidade de probabilidade. O R oferece duas funções que executam integração simples, area() (pacote MASS) e integrate(). No entanto, a função area() não trabalhar com limites infinitos na integral, exigindo, portanto, um conhecimento prévio da região de integração. Por outro lado, a função integrate() aceita limites infinitos, mas infelizmente não é robusta e pode produzir saídas sem confiabilidade.
Exemplo
Compare o uso da função integrate() para aproximar o valor de \(\Gamma(\alpha)\), definido como: ## _0x{- 1}^{-x} x, ##
com o cálculo efetuado pela função gamma().
# integração numérica da função gama - integrate()
f.gama <- function(a){integrate(function(x){x^(a - 1)*exp(-x)},0,Inf)$val}
# valores do parâmetro alfa da função gama
alfa <- seq(0.01, 10, le = 100)
# gráfico do log da função gama - log(integrate()) vs. lgamma()
plot(lgamma(alfa), log(sapply(alfa, f.gama)),
xlab = "log(integrate(f))", ylab = expression(log(Gamma(alpha))),
pch = 19, cex = 0.6)
Não se percebe diferenças nas quantidades aproximadas por integrate(), mesmo para pequenos valores de \(\alpha\) (todos os pontos situam-se praticamente sobre a identidade).
A principal dificuldade dos métodos de integração numérica, como integrate(), é que muitas vezes eles não conseguem identificar as regiões de importantes para a integração dafunção. Por outro lado, os métodos de simulação focam estas regiões naturalmente, ao explorar as informações fornecidas pela função de densidade de probabilidade associada às integrais.
Exemplo
Considere uma amostra \(x_1, x_2, \dots, x_n\), de variáveis aleatórias independentes e com distribuição de Cauchy, com parâmtero de locação \(\theta = 350\). Deseja-se obter a quantidade
\[ m(\mathbf{x}) = \int_{-\infty}^\infty \prod_{i=1}^{10} \frac{1}{\pi} \frac{1}{1 + (x_i - \theta)^2} \operatorname{d}\! \theta. \]
# amostra Cauchy
amostra <- rcauchy(10)+350
# função de verossimilhança
veros <- function(teta){
u <- dcauchy(amostra[1]- teta)
for (i in 2:10)
u <- u*dcauchy(amostra[i]- teta)
return(u)
}
# integral em toda a reta
integrate(veros,-Inf,Inf)## 5.899823e-44 with absolute error < 1.2e-43 # integral região da cauda superior
integrate(veros, 200, 400)## 3.729636e-11 with absolute error < 6.9e-11A função integrate() fornece valor errado, falhando em sinalizar a dificuldade de aproximar o valor exato, já que o erro de aproximação é muito pequeno. O resultado obtido pela função area() é muito melhor
# ----- Comparação integrate() e area()
# amostra - locação zero
amostra <- rcauchy(10)
# integração em região simétrica - integrate()
nin <- function(a){integrate(veros, -a, a)$val}
# integração em região simétrica - area()
nan <- function(a){MASS::area(veros,-a, a)}
# limites de integração
x <- seq(1, 1E3, le = 1E4)
# integrais por integrate()
y <- log(apply(as.matrix(x),1,nin))
# integrais por area()
z <- log(apply(as.matrix(x),1,nan))
# comparação entre integrate() e area() - logescala
plot(x, y, type = "l", ylim = range(cbind(y,z)), lwd = 2)
lines(x, z, lty=2, col = "sienna", lwd = 2)
Na figura acima, observa-se a comparação entre integrate() e area() na aproximação da integral de \(m(\mathbf{x})\), em escala logarítmica. (a curva tracejada corresponde ao resultado obtido com a função area()). Percebe-se, nesse caso, que o uso de area() oferece uma avaliação mais confiável, já que a area(veros, -a, a) aumenta à medida com o aumento de a. Isso obviamente exige algum conhecimento prévio sobre a localização da moda do integrando. Observe abaixo o gráfico do integrando e perceba as situações relacionadas com as aproximações obtidas anteriormente.
# plot do integrando
curve(veros, from = -0.2E2, to = 0.2E2)
Por último, os procedimentos de integração numérica não resolvem facilmente as integrais multidimensionais (com quantidade moderada ou alta de variáveis), as quais são a regra em problemas estatísticos. O desenvolvimento de ferramentas de integração específicas para esses problemas seria muito caro computacionalmente, até porque podemos tirar proveito da natureza probabilística dessas integrais.
Eficiência na Programação em R:
É importante usar as funções em R de maneira completa. Por exemplo, desejamos usar simulação Monte Carlo para estimar a integral:
\[ \int_0^1 x^2 \operatorname{d}\! x. \]
Basicamente, seria como se atirássemos em direção à região do plano que contém a curva e contássemos a quantidade de acertos abaixo da curva. Observe o gráfico abaixo:
# Exemplo de simulação Monte Carlo para estimar área
# função a ser integrada no intervalo (0, 1)
curve(x^2, from = 0, to = 1, lwd = 2, col = "blue", ylab = expression(x^2))
# geração dos pontos a serem plotados
pontos <- 30
u1 <- runif(30)
u2 <- runif(30)
# aceitação/rejeição dos pontos
aceita <- u1^2 > u2
# aceita = T; rejeita = F
codigos <- as.factor(aceita)
# cores dos pontos conforme região
cores <- c("deepskyblue", "darkblue")
# anotação dos pontos gerados
points(u1, u2, col = cores[codigos], pch = 19, cex = 1.3)
# identificação das regiões
text(0.9, 0, pos = 3, font = 2, label = "Aceita")
text(0.1, 1, pos = 1, font = 2, label = "Rejeita")
Para estimar a área abaixo da curva são plotados no gráfico pontos aleatórios, contando-se a quantidade de pontos que estão abaixo da curva.
Um algoritmo clássico para efetuar esse tipo de estimativa seria:
1. Inicialização: rodadas = 0
2. for i in 1:n`
3. Gerar U1 e U2. dois números ao acaso entre 0 e 1.
4. Se U2 < U1^2, então rodadas = rodadas +1
5. Final do for.
6. Área estimada: rodadas/NO seguinte código seria uma possível implementação para esse algoritmo de simulação Monte Carlo:
# função não vetorizada
MonteCarlo <- function(N) {
rodadas = 0
for (i in seq_len(N)) {
U1 <- runif(1)
U2 <- runif(1)
if (U1^2 > U2)
rodadas = rodadas + 1
}
return(rodadas / N)
}A função criada pelo código acima não é vetorizada e verifica-se seu tempo de execução.
# tempo de execução
N <- 5E5
system.time(MonteCarlo(N))## user system elapsed
## 2.51 0.00 3.12Uma abordagem mais voltada ao R seria:
# função vetorizada
MonteCarlo.vet <- function(N) sum(runif(N)^2 > runif(N)) / NA função MonteCarlo.vet() contém pelo menos quatro aspectos de vetorização:
- A função
runif()é integralmente vetorizada. - Elevamos vetores inteiros ao quadrado com o operador
^. - As comparações com o operador
>são vetorizadas. - A função
sum()é mais rápida que seu equivalente em um loop
# tempo de execução
N <- 5E5
system.time(MonteCarlo.vet(N))## user system elapsed
## 0.04 0.00 0.05A função MonteCarlo.vet é aproximadamente 30 vezes mais rápida que a função MonteCarlo.
Uma outra questão é sobre a precisão da estimativa. Uma pergunta pertinente é se a aproximação da integral fica mais precisa se forem usados mais pontos aleatórios.
Exemplo
vamos obter uma estimativa Monte Carlo para a integral:
\[ \int_0^1 sen(10x^2)^2sen(x) + 0.1 \operatorname{d}\! x \]
# função do integrando
f <- function(x) ((sin(10*x^2))^2*sin(x))*x + 0.1
# gráfico da função
curve(f, from = 0, to = 1, lwd = 2, col = "blue", ylab = "", cex.axis = 0.8)
title(ylab = expression(sin(10*x^2)^2+sin(x)), line = 2, cex.lab = 1.1)
# área abaixo da curva - exata
(int.ex <- integrate (f, lower = 0, upper = 1)$value)## [1] 0.2426328Embora a integral possa ser calculada analiticamente, a função integrate() fornece o valor de 0.2426328, que consideraremos como o valor exato da área sob a curva no intervalo \((0, 1)\). Vamos proceder à simulação Monte Carlo de utilizando 1000 pontos escolhidos inteiramente ao acaso na região \((0,1)^2\),
# Geração dos pontos
N <- 1E3
x1 <- runif(N, min = 0 , max =1 )
y1 <- runif(N, min = 0 , max =1 )
# aceitação/rejeição dos pontos
aceita <- f(x1) > y1
# aceita = T; rejeita = F
codigos <- as.factor(aceita)
# cores dos pontos conforme região
cores <- c("deepskyblue", "darkblue")
# gráfico da função
curve(f, from = 0, to = 1, lwd = 2, col = "blue", ylab = "", cex.axis = 0.8)
title(ylab = expression(sin(10*x^2)^2+sin(x)), line = 2, cex.lab = 1.1)
# anotação dos pontos gerados
points(x1, y1, col = cores[codigos], pch = 20, cex = 0.75)
# aproximação área abaixo da curva
(int.est <- sum(aceita)/N)## [1] 0.265 # erro relativo da estimativa
(int.err <- (int.est - int.ex)/int.ex)## [1] 0.09218561est.1e3 <- c(int.est, int.err)A estimativa obtida foi 0.265, com um erro relativo de 0.0921856. Vamos verificar a estimativa com \(N=10000\).
# -------- 10000 pontos
# Geração dos pontos
N <- 1E4
x1 <- runif(N, min = 0 , max =1 )
y1 <- runif(N, min = 0 , max =1 )
# aceitação/rejeição dos pontos
aceita <- f(x1) > y1
# aceita = T; rejeita = F
codigos <- as.factor(aceita)
# cores dos pontos conforme região
cores <- c("deepskyblue", "darkblue")
# gráfico da função
curve(f, from = 0, to = 1, lwd = 2, col = "blue", ylab = "", cex.axis = 0.8)
title(ylab = expression(sin(10*x^2)^2+sin(x)), line = 2, cex.lab = 1.1)
# anotação dos pontos gerados
points(x1, y1, col = cores[codigos], pch = 20, cex = 0.75)
# aproximação área abaixo da curva
(int.est <- sum(aceita)/N)## [1] 0.2466 # erro relativo da estimativa
(int.err <- (int.est - int.ex)/int.ex)## [1] 0.01635083est.10e3 <- c(int.est, int.err)A estimativa obtida foi 0.2466, com um erro relativo de 0.0163508. Vamos verificar a estimativa com \(N=15000\).
# -------- 15000 pontos
# Geração dos pontos
N <- 1.5E4
x1 <- runif(N, min = 0 , max =1 )
y1 <- runif(N, min = 0 , max =1 )
# aceitação/rejeição dos pontos
aceita <- f(x1) > y1
# aceita = T; rejeita = F
codigos <- as.factor(aceita)
# cores dos pontos conforme região
cores <- c("deepskyblue", "darkblue")
# gráfico da função
curve(f, from = 0, to = 1, lwd = 2, col = "blue", ylab = "", cex.axis = 0.8)
title(ylab = expression(sin(10*x^2)^2+sin(x)), line = 2, cex.lab = 1.1)
# anotação dos pontos gerados
points(x1, y1, col = cores[codigos], pch = 20, cex = 0.75)
# aproximação área abaixo da curva
(int.est <- sum(aceita)/N)## [1] 0.2377333 # erro relativo da estimativa
(int.err <- (int.est - int.ex)/int.ex)## [1] -0.02019274est.15e3 <- c(int.est, int.err)As estimativas apresentaram os seguintes erros relativos ao estimar a quantidade exata (0.2426328). São significativas as modificações no erro relativo à medida que o tamanho de pontos gerados aumenta.
## estimativa Erro relativo
## 1000 0.2650000 0.09218561
## 10000 0.2466000 0.01635083
## 15000 0.2377333 -0.02019274Temos entretanto algumas questões com relação à esta abordagem: o esforço computacional e a precisão, já que apenas uma parte dos pontos gerados foi utilizado na estimação. Além disso, surge uma pergunta: como estimar integrais impróprias?
Integração Monte Carlo para integrais definidas no intervalo \((0,1)\):
Há um outro tipo de abordagem que leva a resultados mais precisos, com um menr esforço computacional. Seja uma função \(g(x)\) e suponha que desejamos calcular \(\theta\), em que: \[ \theta = \int_0^1g(x) \,\operatorname{d}\! x %\mathrm{d}x. \]
Perceba que, se \(U \sim \mathcal{U}(0,1)\), podemos expressar o valor de \(\theta\) como:
\[ \theta = \operatorname{E}[g(U)]. \]
Se \(U_1, U_2, \dots, U_n\) são variáveis aleatórias independentes e uniformente distribuídas sobre o intervalo \((0, 1)\), então, \(g(U_1), g(U_2), \dots, g(U_n)\) são variáveis aleatórias independentes e identicamente distribuídas, com média \(\theta\). Portanto, pela Lei Forte dos Grandes Números, segue que, com probabilidade 1: \[ \bar{\theta}_n = \frac{1}{n}\sum_{i=1}^n g(U_i) \to \operatorname{E}[g(U)]=\theta, \, \text{para } n \to \infty. \]
Ou seja, podemos obter uma aproximação de \(\theta\) gerando um grande número de números aleatórios \(u_i\) e usar a média dos valores \(g(u_i)\) como uma aproximação de \(\theta\). Essa abordagem para estimação de integrais definidas é denominada Integração Monte Carlo.
Exemplo
Estimar a integral definida: \[ B(x, y) = \int_0^1 (1-t)^{y-1} t^{x-1} \, \operatorname{d}\! t, \, x>0, \, y>0 \]
Integração Monte Carlo para integrais definidas com suporte compacto:
Integrais do tipo:
\[ \theta = \int_a^b g(x) \operatorname{d}\! x. \]
Para se usar o mesmo procedimento (aproximação pela média dos valores da função \(g\) avaliada em valores aletaórios de uma uniforme), pode-se usar uma substituição conveniente para transfomar os limites de integração do intervalo \((a, b)\) para o intervalo \((0, 1)\). Se substituímos a variável \(x\) por: \[ u = \frac{(x-a)}{b-a}, \]
então, \(\operatorname{d}\! u = \frac{\operatorname{d}\! x}{(b-a)}\) e os extremos de integração tornam-se \(0\) e \(1\). Logo: \[ \theta = \int_0^1 g(a + [b-a]y)\, (b-a) \, \operatorname{d}\! y, \]
ou, de maneira mais simplificada
\[\begin{equation} \theta = \int_a^b g(x) \operatorname{d}\! x = \int_0^1 h(u) \, \operatorname{d}\! u, \end{equation}\]
com \(h(u) = (b - a)g(a + [b-a]u)\), com \(0 < u < 1\).
Dessa maneira \(\theta\) pode ser aproximada pela geração de números aleatórios uniformente distribuídos no intervalo \((0, 1)\), estimando-se \(\theta\) pelo valor médio da função \(h\), calculadas nesses números aleatórios.
Exemplo
Usar integração Monte Carlo para estimar a probabilidade \(\operatorname{P}(1 < Z < 2)\), sendo \(Z\) a variável aleatória normal padrão. Deseja-se então estimar a quantidade:
\[ \theta = \operatorname{P}(1 < Z < 2) = \int_1^2 \frac{1}{\sqrt(2\pi)} \operatorname{e}^{-\frac{z^2}{2}} \, \operatorname{d}\! z. \]
A transformação adequada é para a função \(h(u) = (2 - 1)g(1 + [2-1]u\), ou seja:
\[ \theta = \int_0^1 \frac{1}{\sqrt(2\pi)} \operatorname{e}^{-\frac{(u + 1)^2}{2}} \, \operatorname{d}\! u. \]
Integração Monte Carlo para integrais definidas com suporte na semi-reta:
Integrais do tipo:
\[ \theta = \int_0^\infty g(x) \operatorname{d}\! x. \]
Uma substituição conveniente para transfomar os limites de integração do intervalo \((0, \infty)\) para o intervalo \((0, 1)\) é substituímos a variável \(x\) por: \[ u = \frac{1}{x + 1}, \]
então, \(\operatorname{d}\! u = \frac{-\operatorname{d}\! x}{(x+1)^2}\) e os extremos de integração tornam-se \(0\) e \(1\). Logo: \[ \theta = \int_0^1 \frac{g(\frac{1}{u}-1)}{u^2}\, \operatorname{d}\! u, \]
ou seja, \[ \theta = \int_0^\infty g(x) \operatorname{d}\! x = \int_0^1 h(u) \, \operatorname{d}\! u, \] em que \(h(u) = \frac{g\left(\frac{1}{u}-1\right)}{u^2}\), com \(0 < u < 1\).
Exemplo
Usar integração Monte Carlo para estimar a quantidade: \[ \theta = \int_0^\infty \operatorname{e}^{-x^2} \,\operatorname{d}\! x \]
Uma transformação conveniente é a função \(h(u) = \frac{\operatorname{e^{-\left(\frac{1}{u}-1 \right)^2}}}{u^2}\), ou seja:
\[ \theta = \int_0^1 \frac{\operatorname{e^{-\left(\frac{1}{u}-1 \right)^2}}}{u^2} \, \operatorname{d}\! u. \]
Há outras parametrizações possíveis, como por exemplo \(u = \frac{x}{1 + x}\), em que \(\operatorname{d}\! u = \frac{\operatorname{d}\! x}{(1 - x)^2}\). Dessa maneira: \[ \theta = \int_0^\infty f(x) \, \operatorname{d}\! x = \int_0^1 \frac{f\left(\frac{u}{1 - u} \right)}{(1 - u)^2} \, \operatorname{d}\! u \]
A quantidade \(\theta\) pode também ser expressa como: \[ \theta = \int_0^\infty f(x) \operatorname{d}\! x = \int_0^1 f(x) \operatorname{d}\! x + \int_1^\infty f(x) \operatorname{d}\! x \]
e, dessa maneira, a substituição \(u = \frac{1}{x}\), com \(\operatorname{d}\! u = - \frac{\operatorname{d}\! x}{x^2}\) permite obtermos duas integrais definidas no intervalo \((0, 1)\).
Integral Clássica de Monte Carlo:
Denominamos o procedimento em estudo como Método Clássico de Integração de Monte Carlo. Seja o problema genérico de avaliar a integral
\[ \tag{*} \operatorname{E}_f[h(X)] = \int_{R_X} h(x) f(x) \, \operatorname{d}\! x, \]
em que, \(f(x)\) é a função de densidade de probabilidade da variável aleatória \(X\) e \(h\) uma função.
Seja, por exemplo, \(\operatorname{E}(X) = \mu\) e \(h(x) = (x - \mu)^2\), então \(\operatorname{E}[h(X)]= \operatorname{E}(X-\mu)^2 = \operatorname{Var}(X)\). Uma outra situação é para \(h(x) = \mathbb{I}_A(x)\), em que:
\[ \mathbb{I}_A(x) = \begin{cases} 1 &, \text{ se } x \in A \\ 0 &, \text{ se } x \not\in A, \end{cases} \]
logo, \(\operatorname{E}[h(X)] = \operatorname{P}(A)\).
Nesse contexto, para estimar \(\operatorname{E}[h(x)]\) \((*)\) , é natural propor o uso de uma sequência de números aleatórios \((X_1, X_2, \dots, X_n)\), uma amostra da variável aleatória \(X\), com função de densidade de probabilidade \(f\), para aproximar pela média empírica:
\[ \bar{h}_n = \frac{1}{n}\sum_{j=1}^n h(x_j) \]
Pela Lei Forte dos Grandes Números, a média amostral \(\bar{h}_n\) converge quase certamente para \(\operatorname{E}_f[h(x)]\). Quando \(h^2\) tem esperança finita sob \(f\), a velocidade de convergência de \(\bar{h}_n\) pode ser calculada, já que:
\[ \operatorname{Var}(\bar{h}_n) = \frac{\operatorname{Var}[h(X)]}{n} = \frac{1}{n} \int_{R_X} (h(x) - \operatorname{E}_f[h(x)])^2 f(x) \operatorname{d}\! x, \]
e pode ser estimada pela amostra \((X_1, X_2, \dots, X_n)\) por meio da estatística:
\[ v_n = \frac{1}{n} \sum_{j=1}^n (h(x_j) - \bar{h}_n)^2, \]
já que \[ \hat{\operatorname{Var}[h(X)]} = \frac{\sum_{j=1}^n(h(x_j) - \bar{h}_n)^2}{m-1} \approx \frac{\sum_{j=1}^n(h(x_j) - \bar{h}_n)^2}{m}. \]
Assim, o erro padrão associado aos erros padrão da estimativa Monte Carlo são:
\[ \operatorname{ep}_{\bar{h}} = \sqrt{\operatorname{Var}(\bar{h}_n)} \approx \sqrt{\frac{\hat{\operatorname{Var}[h(X)]}}{m}}. \]
Por outro lado, para \(n\) suficientemente grande
\[ \frac{\bar{h}_n - \operatorname{E}_f[h(X)]}{\sqrt{v_n}} \sim \mathcal{N}(0,1). \]
Esse fato permite a construção de teste de convergência e de limites de confiança para a aproximação de \(\operatorname{E}_f[h(X)]\).
Exemplo
A função de distribuição acumulada de uma normal não pode ser expressa de forma explícita. É possível utilizarmos simulações na construção de tabelas com probabilidades acumuladas dessa distribuição.
Considere uma amostra aleatória de tamanho n, \((x_1, x_2, \dots, x_n)\), gerada pelo algoritmo de Box-Muller. A o valor da função de probabilidade acumulada em \(t\):
\[ \Phi(t) = \int_{-\infty}^t \frac{1}{\sqrt{2 \pi}}\operatorname{e}^{- \frac{z^2}{2}} \operatorname{d}\! z, \]
pode ser aproximada pelo método de Monte Carlo por meio de: \[ \hat{\Phi}(t) = \frac{1}{n}\sum_{j=1}^n \mathbb{I}_{x_j \leq t}. \]
A variância exata de \(\hat{\Phi}(t)\) é:
\[ \frac{\Phi(t)(1- \Phi(t))}{n}, \]
já que as variáveis aleatórias \(\mathbb{I_{x_j \leq t}}\) são independentes e identicamente distribuídas de acordo a uma distribuição de Bernoulli, com probabilidade de sucesso \(\Phi(t)\).
Percebe-se que, para valores de \(t\) em torno de zero, a variância é aproximadamente \(\frac{1}{4n}\). Assim, para obter-se uma precisão de 4 casa decimais, a aproximação de Monte Carlo exige uma média \(n=(\sqrt{2}\, 10^4)^2\) iterações, ou seja, aproximadamente 200 milhões de iterações.
Robert e Casella (1999) avaliaram a estimação por integração de Monte Carlo de probabilidades de uma normal padrão para vários valores de t, apresentando os resultado da tabela abaixo:
A avaliação é precisa (considerada exata) para 100 milhões de iterações, percebendo-se maior precisão nas caudas. Há métodos mais eficientes para serem usados na estimação dos valores dessas probabilidade, como, por exemplo a importance sampling.
Exemplo - Toy function
Seja a função:
\[ h(x) = [cos(50x) + sin(20x)]^2, \]
cujo gráfico encontra-se no painel superior da figura abaixo. Embora essa função possa ser integrada analiticamente, pode-se usar o método Monte carlo para estimadar a quantidade \(\int_0^1h(x) \operatorname{d}\! x\). Essa integral pode ser vista como \(\operatorname{E}_U[h(x)]\), como \(U \sim \mathcal{U}(0,1)\). Dessa maneira, podemos obter uma aproximação do valor da integral por meio da média amostral \(\sum_{j=1}^n h(U_i)\), em que \(U_1, U_2, \dots, U_n\) é uma sequência de números aleatórios independentes e identicamente distribuídas, gerados ao acaso no intervalo \((0, 1)\).
# ------------------------------------------
# Integração da Toy Function
# ------------------------------------------
# função h(x)
h <- function(x){(cos(50*x)+sin(20*x))^2}
# parâmetros do gráfico
par(mar = c(2, 2, 2, 1), mfrow = c(2, 1))
# Gráfico da função
curve(h, xlab = "Função", ylab = "", lwd = 2)
# valor exato da integral
(int.ex <- integrate(h, 0, 1)$value)## [1] 0.9652009 # ----- Aproximação Monte Carlos
n <- 1E4
U <- runif(n)
x <- h(U)
# estimativa móvel da integral
int.est <- cumsum(x)/(1:n)
# erro padrão móvel da estimativa
err.est <- sqrt(cumsum((x - int.est)^2))/(1:n)
# gráfico da estimativa móvel da integral
plot(int.est, xlab = "Média e intervalo de erro", type = "l", lwd= 2,
ylim = mean(x) + 20*c(-err.est[n], err.est[n]), ylab="")
# limites do intervalo de errro
lines(int.est + 2*err.est, col = "gold", lwd = 2)
lines(int.est - 2*err.est, col = "gold", lwd = 2)
# objetivo - valor exato da integral
abline(h = int.ex, lty = 2, col = "blue")
O painel inferior apresenta a média móvel e os limites obtidos a partir dos erros padrão estimados versus a quantidade de simulações (n). A banda de confiança obtida não é a clássica banda com 95% de confiança.
Estimação Monte Carlo de Integrais Múltiplas
A utilidade de usar números aleatórios para aproximar integrais se torna mais aparente no caso de integrais multidimensionais. Seja \(g\) uma função em \(\mathbb{R}^p\) e estamos interessados em calcular:
\[ \theta = \int_0^1\int_0^1 \dots \int_0^1 g(x_1, x_2, \dots, x_p) \operatorname{d}\! x_1 \operatorname{d}\! x_2 \dots \operatorname{d}\! x_p. \]
pode ser expresso como a seguinte esperança:
\[ \operatorname{E}[g(U_1, U_2, \dots, U_p)], \]
em que \(U_1, U_2, \dots, U_p\) são variáveis aleatórias independentes e uniformemente distribuídas no intervalo \((0,1)\).
Seja a geração de \(n\) conjuntos independentes, com \(p\) variáveis aleatórias uniformes independentes:
\[\begin{align} (U_1^{(1)}, U_2^{(1)}, & \dots, U_p^{(1)}) \\ (U_1^{(2)}, U_2^{(2)}, & \dots, U_p^{(2)}) \\ & \vdots \\ (U_1^{(n)}, U_2^{(n)}, & \dots, U_p^{(n)}) . \end{align}\]
Por outro lado, as variáveis aleatórias \(g(U_1^{(j)}, U_2^{(j)}, \dots, U_p^{(j)})\), \(j=1, 2, \dots, p\), são variáveis aleatórias independentes e identicamente distribuídas, com média \(\theta\). Então podemos estimar \(\theta\) por:
\[ \bar{\theta} = \frac{\sum_{j=1}^n g(U_1^{(j)}, U_2^{(j)}, \dots, U_p^{(j)})}{n}. \]
Exemplo
Duas pessoas têm encontro no Poleiro do Galo no Valentine’s Day. P2 chega ao acaso entre 22:00 e 23:30 e P1, entre 22:30 e 24:00.
- Qual a probabilidade de que
P2chegue antes deP1?
Sejam \(T_1\) e \(T_2\) as avriáveis aleatórias associadas ao tempo de chegada (medidos em horas após o meio-dia) de P1 e P2, respectivamente. Assume-se que \(T_1\) e \(T_2\) são independentes, com \(T_1\), distribuída uniformemente no intervalo \((10.5, 12.0)\) e \(T_2\), no intervalo \((10.0, 11.5)\). Deseja-se obter a probabilidade:
\[ \operatorname{P}\{T_1 < T_2\} = \int_{t_1 < t_2} f_{T_1}(t_1) f_{T_2}(t_2) \operatorname{d}\! t_1 \operatorname{d}\! t_2, \]
em que \(f_{T_1}\) e \(f_{T_2}\) denotam as funções de densidade de probabilidade uniformes para \(T_1\) e \(T_2\), respectivamente. Dada a suposição de independência, a função de densidade de probabilidade conjunta de \((T_1, T_2)\) é \(f_{T_1, T_2}(t_1 t_2) = f_{T_1}(t_1)f_{T_2}(t_2)\). A probabilidade de interesse pode ser representada pela esperança:
\[ \operatorname{P}\{T_1 < T_2\} = \operatorname{E}(\mathbb{I}(T_1<T_2)) = \int \mathbb{I}(T_1<T_2) f_{T_1}(t_1) f_{T_2}(t_2) \operatorname{d}\! t_1 \operatorname{d}\! t_2, \]
em que \(\mathbb{I}(T_1<T_2)\) é a função indicadora de \(T_1<T_2\), ou seja, é iguala a um quando \(t_1 < t_2\) e zero, caso contrário. Para aproximar essa integral pelo método de Monte Carlo, simula-se uma grande número, digamos 1000, de valores da distribuição de \((T_1, T_2)\). Como \(T_1\) e \(T_2\) são independentes, simulamos 1000 valores de \(T_1\), da distribuição uniforme no intervalo \((10.5, 12)\) e, independentemente, outros 1000 valores da distribuição uniforme no intervalo \((10, 11.5)\).
# geração dos tempos de chegada
n <- 1E3
t1 <- runif(n, 10.5, 12)
t2 <- runif(n, 10, 11.5)
# estimativa da probabilidade
prob <- sum(t1 < t2)/n
prob## [1] 0.225A estimativa da probabilidade \(\operatorname{P}(T_1 < T_2)\) é dada pela proporção dos pares ordenados simulados \((t_1, t_2)\), nos quais \(t_1 < t_2\). Estima-se que a probabilidade de P1 chegar antes de P2 é 0.225.
# Figura da simulação
cores <- c("deepskyblue", "darkblue")
aceita <- as.factor(t1 < t2)
plot(t2, t1, xlab = "P2", ylab = "P1", type = "n")
# área de integração
polygon(c(10.5, 11.5, 11.5, 10.5), c(10.5, 10.5, 11.5, 10.5), density = 10,
angle = 135)
# área descartada
polygon(c(10, 10.5, 11.5, 11.5, 10), c(10.5, 10.5, 11.5, 12, 12), col = "cyan")
# pontos simulados
points(t2, t1, pch = 19, col=cores[aceita])
A figura ilusta o processo de estimação da probabilidade. A região hachurada representa a região em que \((T_1 < T_2)\) e a estimativa Monte Caro é a proporção de pontos simulados que pertence à área hachurada. A função polygon() adiciona a área hachurada (região de aceitação) e a área colorida (região de rejeição).
Como essa estimativa Monte Carlo da probabilidade é uma proporção, podemos usar a fórmula usual do erro padrão de uma proporção para calcular o erro padrão dessa estimativa:
\[ \operatorname{ep}_\hat{p} = \sqrt(\frac{\hat{p}(1-\hat{p})}{n}), \]
em que \(n\) é a quantidade de pares ordenados simulados. Para nossa estimativa, obtemos:
# erro padrão
(ep <- sqrt(prob * (1 - prob)/n))## [1] 0.01320511Aplicando a aproximação normal para a distribuição amostral de \(\hat{p}\), com 95% de confiança P1 chegará mais cedo que P2 entre 23.17% e 23.6% das vezes.
- Qual a diferença esperada entre os dois tempos de chegada
Como é mais provável que P1 chegue mais tarde que P2, desejamos calcular a esperança \(\operatorname{E}(T_1 - T_2)\), que pode ser calculada por:
\[ \operatorname{E}(T_1 - T_2) = \int(t_1 - t_2) f_{T_1}(t_1) f_{T_2}(t_2) \operatorname{d}\! t_1 \operatorname{d}\! t_2 \]
A estimativa Monte Carlo dessa esperança é dada por:
\[ \widehat{\operatorname{E}(T_1 - T_2)} = \frac{\sum_{j=1}^n(T_1^{(j)} - T_2^{(j)})}{n}, \]
em que $(T_1^{(j)}, T_2^{(j)}) é um ponto simulado da densidade conjunta de \((T_1, T_2)\), ou seja, a estimativa Monte Carlo da diferença é dada por:
# diferença de chegadas em cada encontro
diferenca <- t1 - t2
# estimativa Monte Carlo da esperança
(E.mc <- mean(diferenca))## [1] 0.4861744 # Erro padrão da estimativa Monte Carlo da esperança
(ep.mc <- sd(diferenca)/sqrt(n))## [1] 0.01919389Estimamos que P1 vai chegar, em média, 0.49 horas mais tarde que P2. Como o erro padrão é 0.02, estamos 95% confiantes que a verdadeira diferença está num intervalo de 0.04 dessa estimativa. Ressalta-se que o valor exato dessa diferença é 0.5.
Exemplo - Estimação de \(\pi\)
Suponha que o vetor aleatório \((X,Y)\) esteja uniformemente distribuído em um quadrado de área 4, centrado na origem. Considere agora a probabilidade de um ponto aleatório \((X,Y)\) estar contido em um círculo de raio 1, inscrito nesse quadrado.
par(mar = c(4.1, 4.1, 4.1, 4.1))
plot(NA, xlim = c(-1, 1), ylim = c(-1, 1), type = "n", xaxs="i",yaxs="i",
ylab = "", xlab = "")
# círculo
plotrix::draw.circle(0, 0, 1, col = "cyan")
Como o vetor aleatório \((X,Y)\) está uniformemente distribuído no quadrado, tem-se que:
\[ \operatorname{P}\{(X, Y) \in \text{ círculo} \} = \operatorname{P}\{X^2 + Y^2 \leq 1 \} = \frac{\text{área círculo}}{\text{área quadrado}} = \frac{\pi}{4} \]
Portanto, se gerarmos um grande número de pontos aleatórios no quadrado, a proporção de pontos que pertencerá ao círculo será aproximadamente \(\pi/4\). Assim, se \(X\) e \(Y\) forem independentes e uniformemente distribuídas entre \((-1, 1)\), sua função de densidade de probabilidade conjunta será:
\[\begin{align} f_{XY}(x,y) & =f_X(x) f_Y(y) \\ & = \frac{1}{2}\cdot \frac{1}{2} \\ & = \frac{1}{4} \text{, }-1 \leq x \leq 1 \text{, }-1 \leq y \leq 1. \end{align}\]
Como a densidade de \((X, Y)\) é constante no quadrado, por definição, \((X,Y)\) é uniformemente distribuída no quadrado. Por outro lado, se \(U\) é uniformemente distribuída no intervalo \((0, 1)\), então \(2U\) é uniforme em \((0, 2)\) e \((2U -1)\) é uniforme em \((-1,1)\). Então, podemos implementar o seguinte algoritmo para efetuar essa estimativa:
1. Gerar ao acaso os números U1 e U2
2. Calcular X = 2U1 -1 e Y = 2U2 - 1
3. Calcular a proporção de pares em que X^2 + Y^2 <= 1.Perceba que se \[ \mathbb{I}(x, y) = \begin{cases} 1 & \text{, se } X^2 + Y^2 \leq 1 \\ 0 & \text{, caso contrário,} \end{cases} \]
então
\[\begin{align} \operatorname{E}[\mathbb{I}(x,y)] & = \int_\Omega \mathbb{I}(x,y) \operatorname{d}\! x \operatorname{d}\! y \\ & = \operatorname{P}\{ X^2 + Y^2 \leq 1 \} \\ & = \frac{\pi}{4}, \end{align}\]
em que \(\Omega = [-1, 1] \times [-1, 1]\)
O seguinte código seria uma possível implementação para esse algoritmo de simulação Monte Carlo:
# Estimação de pi
# geração da uniforme (0,1)
n <- 1E3
u1 <- runif(n)
u2 <- runif(n)
# geração dos pares ordenados
x <- 2*u1 - 1
y <- 2*u2 - 1
circulo <- x^2 + y^2 <= 1
# Proporção de pontos no círculo
(mean(circulo))## [1] 0.779 # Estimativa de pi
mean(circulo)*4## [1] 3.116Essa estimativa apresenta os seguintes erro amostral e erro relativo:
# erro amostral
mean(circulo)*4 - pi## [1] -0.02559265 # erro amostral
(mean(circulo)*4 - pi)/pi## [1] -0.008146395A simulação apresentou o seguinte comportamento, sendo que 77.9% dos pontos situaram-se dentro do cícrulo.
# Pontos plotados
par(mar = c(4.1, 4.1, 4.1, 4.1))
cores <- c("deepskyblue", "darkblue")
aceita <- as.factor(circulo)
plot(NA, xlim = c(-1, 1), ylim = c(-1, 1), type = "n", xaxs = "i", yaxs = "i",
ylab = "", xlab = "")
# círculo
plotrix::draw.circle(0, 0, 1, col = "cyan")
# pontos gerados
points(x, y, pch = 19, col=cores[aceita], cex = 0.75)