O método de integração de Monte Carlo

Desejamos calcular uma integral definida, como \(I = \int_{a}^{b}f(x)dx\). Onde \(x\in[a,b]\) e \(f(x)\in[c,d]\). Sendo \(A\) a área abaixo da curva limitada pela caixa \([a,b]\times[c,d]\), então podemos estimar o valor da integral da seguinte forma: \[I = A + c(b - a)\] Portanto, basta estimar A que teremos um valor estimado da integral.

Para estimar \(A\) basta usar o mesmo princípio utilizado no relatório 1. Geramos \(n\) pontos aleatórios dentro da caixa \([a,b]\times[c,d]\), contando quantos desses pontos estão abaixo da curva da função. A proporção de pontos abaixo da curva será em média equivalente a razão entre a área \(A\) e a área da caixa, isto é \[\frac{A}{(b-a)(d-c)}\].

Sendo \(X\sim U[a,b]\) e \(Y\sim U[c,d]\) então \((X,Y)\) segue uma distribuição uniforme na caixa \([a,b]\times[c,d]\). Sendo \(Z\) uma variável que assume 1 caso o ponto esteja abaixo da curva, ou seja \(Z = 1\) se \(Y \leq f(x)\), e 0 caso contrário. A esperança de \(Z\) será igual a proporção de pontos abaixo da curva, isto é: \[E(Z) = \frac{A}{(b-a)(d-c)}\] Portanto: \[A = \frac{E(Z)}{(b-a)(d-c)}\]

Então podemos estimar a integral da seguinte forma: \[I = E(Z)(b-a)(d-c) + c(b-a)\]

Implementação

Utilizando como exemplo um quarto de uma circunferência (ou seja, apenas o primeiro quadrante) de raio 1, temos que \(x^2 + y^2 = 1\), então \(y = \sqrt{1-x^2}\) onde \(x\in[0,1]\) e consequentemente \(f(x)\in[0,1]\). O valor real da integral é: \[\int_{0}^{1}\sqrt{1-x^2}dx = \frac{\pi}{4} \approx 0,7854\]

Agora, vamos estimar o valor e comparar com o valor real da integral.

#número de pontos gerados
n <- 1000 

#limites
a <- 0
b <- 1
c <- 0
d <- 1

soma <- 0
for (i in 1:n) {
  X <- runif(1, a, b)
  Y <- runif(1, c, d)
  Z <- (sqrt(1-X^2) >= Y)
  soma <- soma + Z
}
I <- (b - a)*c + (soma/n)*(b - a)*(d - c)
print(I)
## [1] 0.792
#número de pontos gerados
n <- 10000 

#limites
a <- 0
b <- 1
c <- 0
d <- 1

soma <- 0
for (i in 1:n) {
  X <- runif(1, a, b)
  Y <- runif(1, c, d)
  Z <- (sqrt(1-X^2) >= Y)
  soma <- soma + Z
}
I <- (b - a)*c + (soma/n)*(b - a)*(d - c)
print(I)
## [1] 0.7812
#número de pontos gerados
n <- 100000

#limites
a <- 0
b <- 1
c <- 0
d <- 1

soma <- 0
for (i in 1:n) {
  X <- runif(1, a, b)
  Y <- runif(1, c, d)
  Z <- (sqrt(1-X^2) >= Y)
  soma <- soma + Z
}
I <- (b - a)*c + (soma/n)*(b - a)*(d - c)
print(I)
## [1] 0.78737

Podemos observar que quanto maior o número de pontos gerados mais o valor estimado se aproxima do valor real calculado anteriormente, assim como esperado.