A integração de Monte Carlos é um método estatístico usado para se estimar o valor de integrais definidas. No trabalho em questão, avaliaremos a integração pela média dos valores obtidos numa função unidimensional. Seja \(g(x)\) uma função definida no intervalo real \([a,b]\) a qual queremos integrar. Consideremos \(X\) uma variável aleatória cujo suporte \(\chi\) seja o próprio conjunto \([a,b]\) e com função de densidade \(f(x)\). Teremos então: \[E(g(X)) = \int_{[a,b]} g(x)f(x)dx\]. Tomando-se amostras aleatórias de \(X\) (\(X_{1},...,X_{n}\)independentes e identicamente distribuídas), sabemos pela Lei de Kolmogorov dos grandes números que \(\frac{\sum_{i=1}^{n}g(x_{i})}{n}\) convergirá quase certamente para \(E(g(X))\) à medida que \(n\rightarrow\infty\). Consideremos então \(X\) uma variável aleatória com distribuição uniforme em \([a,b]\), \(X\sim U[a,b]\), teremos então: \[\frac{\sum_{i=1}^{n}g(x_{i})}{n}\underset{q.c}{\rightarrow} E(g(X)) = \int_{a}^{b}\frac{g(x)}{b-a}dx \] Ou seja: \[(b-a)\frac{\sum_{i=1}^{n}g(x_{i})}{n}\underset{q.c}{\rightarrow} \int_{a}^{b}g(x)dx \] Em outras palavras; seleciona-se, segundoa distribuição uniforme, pontos aleatórios do domínio \([a,b]\) a se integrar. Como mostrado anteriormente, a média das imagens dos pontos selecionados convergirá com probabilidade \(1\) para o valor da integral definida. Considerando agora o exemplo da função \(g(x) = \sqrt{1-x^{2}}\) definida em \([-1,1]\). A estimativa da integral em questão será uma estimativa do valor de \(\frac{\pi}{2}\) (área de um semi-círculo de raio \(1\)). Selecionando uma amostra aleatória de N pontos segundo a distribuição uniforme no conjunto \([-1,1]\) temos:

Implementação em R

Para \(10.000\) pontos:

set.seed(123)
N<-10000
x<-runif(N,-1,1)
pi.hat <- mean(sqrt(1-x^2))
print(2*pi.hat)
## [1] 1.577113

Agora para \(N = 100.000\)

set.seed(123)
N<-100000
x<-runif(N,-1,1)
pi.hat <- mean(sqrt(1-x^2))
print(2*pi.hat)
## [1] 1.571857

E para \(N =1.000.000\):

set.seed(123)
N<-1000000
x<-runif(N,-1,1)
pi.hat <- mean(sqrt(1-x^2))
print(2*pi.hat)
## [1] 1.570367

Estimativas estas cada vez mais precisas do valor de \(\frac{\pi}{2}\).