Integração de Monte Carlo: teoria e aplicação em R e C++.
A integração de Monte Carlo é um método numérico simples e de fácil entendimento. É uma técnica que se baseia na interpretação de uma integral como sendo a área embaixo de uma curva. Consiste, basicamente, em inscrever a função que se deseja integrar em uma outra forma geométrica conveniente — digamos, um retângulo — depois disso, estima-se a proporção da área abaixo da função em relação à área total do retângulo.
Vejamos o raciocínio seguinte. Pense em uma função qualquer, digamos ƒ(x)=x², com o domínio estabelecido entre x=0 e x=2, sua integral é dada como segue.
Agora, imagine um retângulo no qual poderíamos inscrever — por a função dentro — essa função. Um retângulo com base=2 e altura=4. Dentro do domínio restrito de nossa função, ela atinge um valor máximo em x=2, com ƒ(2) =4, portanto 4 seria uma boa escolha para a altura do retângulo, embora pudesse ser qualquer altura, desde que acima do valor máximo da função. A área total do quadrado, portanto, é de 2x4=8. Ao inscrevermos a função, a proporção da área abaixo dela em relação à área total do retângulo é de (8∕3) ∕ 8=1∕3. Agora pensemos, dado que sabemos a área do quadrado e sabemos que a área embaixo da função representa 1/3 da área do retângulo, basta fazermos (1∕3)×8 e temos a área embaixo da curva.
Ótimo, fizemos o caminho mais fácil. Identificamos o valor da área abaixo da função e calculamos a proporção dessa área em relação ao retângulo. Com base nessa proporção e a área do retângulo pudemos obter a área embaixo da curva, ou seja, valor da integral. Isso mostra o fundamento do método e como ele funciona. É claro que já sabíamos de antemão o valor da integral, mas isso foi um recurso meramente ilustrativo. O que precisamos é justamente de um meio de descobrir o valor da integral de forma numérica, e é aí que o método tem seu valor.
O importante agora é conhecer, dado que temos a função e o retângulo, qual a proporção da área abaixo da função em relação à área do retângulo, com isso o problema vira uma conta de padaria.
A estimativa da proporção buscada é feita de forma bem simples. Observe a figura acima, cada ponto dentro do retângulo é uma coordenada aleatória (x, y) com x ∈ [a, b] e y ∈ [0, H]. Note que podemos preencher todo o retângulo com pontos aleatórios, eventualmente toda a área ficará coberta pelos pontos gerados, depois disso basta “contarmos” quantos pontos ficaram abaixo da curva e dividir esse número pelo total de pontos gerados.
Uma vez que sabemos o número de pontos gerados e a quantidade destes os quais estão embaixo da curva e dado que a área do retângulo é calculada por: (b-a)×H, o valor da integral é aproximadamente:
Dito isso, precisamos gerar as coordenadas e encontrar um meio de decidir quando cada coordenada está abaixo ou acima da curva da função ƒ(x). Primeiro de tudo, precisamos definir a altura do retângulo, esta será equivalente ao valor máximo da função ƒ(x). Para acharmos o valor máximo, derivamos ƒ(x) e igualamos a zero, resolvendo para x. Depois fazemos a segunda derivada e verificamos se é um ponto local máximo dentro do intervalo, enfim, aplicamos os métodos do cálculo. Em resumo, as condições para o ponto de máximo devem obedecer:
Agora, se a função é, respectivamente, monótona crescente ou monótona decrescente, os valores máximo e mínimo são, respectivamente, ƒ(b) e ƒ(a), com a<b; e isso facilita o trabalho.
Feito isso, agora vem a parte interessante, preencher o retângulo com coordenadas aleatórias e verificar se dada coordenada está na curva ou embaixo dela. Para tanto, imagine ƒ(x)=x² e x ∈ [1, 4], com x=2, portanto ƒ(2) = 2² = 4, logo a coordenada (2,4) encontra-se exatamente sobre a curva. Com isso, é fácil perceber que, qualquer coordenada (x, y) que obedeça a restrição abaixo estará embaixo de ƒ(x).
Eis o método. Agora, basta que seja feita a simulação de tantos pontos (x, y) quanto se desejar, lembrando que os pontos devem estar dentro do retângulo escolhido, e verificar quais deles estão abaixo da curva ƒ(x), e então dividir o total de pontos abaixo de ƒ(x) pelo número total de pontos gerados, multiplicando a proporção pela área do retângulo. Para que cada coordenada gerada esteja dentro do retângulo, geramos N y’s uniformes com y ∈ [0, H] e um N x’s uniformes com x∈ [a, b]. Se para dado x, y gerados em uma iteração, y ≤ f(x), então a coordenada (x, y) está embaixo de f(x).
Vamos a um exemplo usando a linguagem R.
Seja a função ƒ(x)=x² para qual desejamos encontrar a integral, com os parâmetros a=0, b=2, H=4.
Vamos investigar agora a convergência da estimativa para o valor teórico da integral, conforme aumentamos o tamanho de pontos gerados.
Conforme aumentamos o número de pontos aleatórios, o valor da integral converge para o valor teórico, com a variância diminuindo. Note que, quanto maior for a área do retângulo escolhido em relação a função, maior a quantidade de pontos necessários para obter uma boa estimativa, portanto, menos eficiente será o programa.
Por fim, é válido notar que, poder-se-ia inscrever ƒ(x) em qualquer outra função com integral conhecida, além de um retângulo, desde que esta tenha uma área maior que a área de ƒ(x). Vejamos, considerando que a convergência do valor estimado se dá conforme aumentamos a quantidade de pontos N, se quiséssemos estimar a integral de uma função ƒ(x) a partir de uma função cs(x) (c vezes s(x)). A proporção estimada seria dada por, considerando n como sendo a quantidade de pontos abaixo de ƒ(x):
A constante c garante que a curva de s(x) esteja sempre acima de ƒ(x). Portanto, o valor estimado da integral seria dado por:
Observe que fizemos uma mudança na notação, mas a ideia é a mesma já apresentada.
A questão agora é que, se por um lado, ao adotarmos um retângulo para estimar a proporção a simulação é bem fácil, quando envelopamos a função por outra função mais complexa, precisamos simular y a partir dessa função, o que pode ou não ser mais complicado. De qualquer forma, não importa quão complicada seja a função ƒ(x), sempre se pode inscrever tal função em um retângulo. Esse método de resolução de integral deriva diretamente de um modelo de simulação chamada aceitação-rejeição, além disso, sobretudo quando se trata de funções de várias variáveis, o método é bastante eficiente, muito além dos métodos tradicionais de integração numérica.
Para aqueles que desejarem simulações mais rápidas, eis o código em C++.
Recomendação de leitura:
Cohen, Harold. Numerical approximation methods. New York: Springer, 2011. Brandt, Siegmund, and S. Brandt. Data analysis. New York etc: Springer, 1999.
