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
99 visitas
perguntada Nov 22, 2020 em Programação Computacional por Athos Carvalho (31 pontos)  
Compartilhe

1 Resposta

0 votos
respondida Nov 22, 2020 por Athos Carvalho (31 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.

comentou Dez 9, 2020 por JOAO PAULO (26 pontos)  
Athos, sua resposta foi completa e concisa. Você utilizou o exemplo da secção do livro para mostrar como reduzir a variância no cálculo de integrais utilizando Monte Carlo. Eu teria duas observações somente:
(1)    a primeira delas se refere a forma que você apresentou a aproximação do cálculo da variância em relação a seu valor empírico encontrado. Na resposta você mostra apenas um ponto específico (n=100 ou n=1000). Para que tenhamos uma ideia de como ocorre essa aproximação ficaria mais interessante mostrar a convergência da variância (e da média também) quando aumentamos N. No código abaixo mostro como criar um gráfico que ajuda a visualizar essa convergência. Também, de acordo com esse gráfico, é possível notar como a convergência via integral simples é muito menos volátil do que via integral dupla. Eu fiz uma pequena alteração no seu código para mostrar esse gráfico de convergência de valores para um N a escolha.

(2)    A segunda se refere ao tempo despendido entre as duas integrais. Eu calculei o tempo entre elas e a integral dupla demorou muito mais (o que obviamente, já era esperado) em relação a integral única. Na resposta você mencionou essa demora e afirmou que seria N2. Nos meus cálculos para n=100, enquanto a integral simples demorou 0,8 segundos aproximadamente, a integral dupla demorou mais de 4 segundos. Pode ter havido alguma alteração no meu computador no momento, mas acho que essa demora pode ser maior do N2.
Segue o código do item (1) que mencionei. Abs

import matplotlib.pyplot as plt
import numpy as np
import time



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

if __name__ == '__main__':
#    IntegralDMC(100)
#    IntegralMC(100)
    start_time = time.time()
    n = 100
    maximo = np.arange(n)
    funcao = [IntegralMC(i) for i in maximo]
    plt.plot(maximo, funcao)
    plt.legend(['Variância','Média'])
    plt.xlabel('número de interações')
    plt.ylabel('valor')
    plt.show()
    print("--- %s seconds ---" % (time.time() - start_time))
...