Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Como usar método de Monte Carlo e valor esperado para reduzir a variância no cálculo de integrais?

0 votos
5 visitas
perguntada 3 dias atrás em Programação Computacional por Athos Carvalho (11 pontos)  
Compartilhe

1 Resposta

0 votos
respondida 3 dias atrás por Athos Carvalho (11 pontos)  

Caso queiramos calcular a integral
\[ G= \int g(x, y) f(x, y) dx dy \]
Onde x e y podem ser vetores multidimensionais. A distribuição marginal de x é definida por:
\[m(x)=\int f(x,y)dy\]
E podemos definir uma quantidade \(h(x)\) como:
\[h(x)=\frac{1}{m(x)} \int g(x,y) f(x,y) dy \]
Assumiremos que as integrais \(m(x)\) e \(h(x)\) podem ser calculadas por métodos que não Monte Carlo. A integral inicial pode ser reescrita como:
\[ G= \int m(x)h(x) dx \]
Onde assumimos também que a ordem de integração é irrelevante. A diferença de variância entre as duas formas para \(G\) são:
\[var(g)-var(h)=E(g^2)-E((E(g|x))^2)=\]
\[E(E(g^2|x))-E((E(g|x))^2)=\]
\[E(E(g^2|x)-(E(g|x))^2)=\]
\[E(var(g|x)) \geq 0 \]
De modo que \(var(g(x,y))-var(h(x)) \geq 0 \).

Testaremos agora com um exemplo. Queremos calcular a integral \(\int_{0}^{1} \int_{0}^{1} g(x, y)\ dx\ dy =\frac{\pi}{4}\), onde:
\[g(x, y)= \begin{cases} 1, x^2 + y^2 \leq 1\\ 0, x^2 + y^2 \gt 1 \end{cases} \]
Se escolhermos \(x\) e \(y\) de maneira uniforme, a variância de um estimador Monte Carlo é:
\[var(g)=\frac{\pi}{4}-(\frac{\pi}{4})^2=0,168\]

Utilizando o método apresentado, utilizaremos a distribuição marginal de \(x\)
\[m(x)= \int_{0}^{1} dy =1\]
E definiremos \(h(x)\) como:
\[h(x)=\frac{1}{m(x)} \int g(x, y) f(x, y) dy = \int_{0}^{\sqrt{1-x^2}} dy = \sqrt{1-x^2}\]

A integral reescrita será:
\[G= \int_{0}^{1}\sqrt{1-x^2} dx \]

Um estimador Monte Carlo com amostragem uniforme de \(x\) no intervalo \((0,1)\) teria a seguinte variância:

\[var(h)= \int_{0}^{1} (1-x^2) dx - (\frac{\pi}{4})^2 =0,05 \leq 0,168 \]

A melhora na variância seria perdida caso o segundo método leve mais tempo computacional (no caso, o cálculo da raiz quadrada).

Podemos observar os resultados já calculados analiticamente de maneira empírica, utilizando o código abaixo:

import numpy as np
def funcx(x):
    return np.sqrt(1 - x**2)

def funcxy(x, y):
    if x**2 + y**2 <= 1:
        r = 1
    else:
        r=0
    return r

def IntegralDMC(n):
    x = np.random.uniform(0,1,n)
    a = []
    for i in range(n):
        for j in range(n):
            h_x = funcxy(x[i],x[j])
            a = np.append(a,h_x)
    GN = np.mean(a)
    var = np.var(a)

    return var, GN

def IntegralMC(n):
    x = np.random.uniform(0,1,n)
    a = []
    for i in range(n):    
        h_x = funcx(x[i])
        a = np.append(a,h_x)
    GN = np.mean(a)
    var = np.var(a)

    return var, GN

Ao rodar o código, obtive os seguintes resultados:

IntegralMC(100000)
(0.05009388415259798, 0.7842652790528831)

e

IntegralDMC(1000)
(0.17217518115900005, 0.778971)

Reduzi \(N\) para o cálculo da dupla integral, pois o tempo de cálculo aumenta em relação a \(N^2\), de modo que valores maiores levariam muito tempo. Ainda assim, é possível observar que os valores calculados se aproximam daqueles encontrados analiticamente.

...