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

Como utilizar Importance Sampling para reduzir a variância em simulações Monte Carlo?

+1 voto
15 visitas
perguntada Jul 3 em Estatística por Pedro Gomes (16 pontos)  
Compartilhe

1 Resposta

0 votos
respondida Jul 3 por Pedro Gomes (16 pontos)  
editado Jul 4 por Pedro Gomes

Essa resposta está baseada em "Monte Carlo Methods - Malvin H. Kalos and Paula A. Whitlock" seção 4.1.

Suponha que temos uma integral G, com uma função f(x) que não é necessariamente a melhor função de densidade de probabilidade para usar em uma simulação Monte Carlo:

\[ G = \int g(X)f(X)dX \]

Pode-se introduzir uma função densidade de probabilidade diferente em G, da seguinte forma:

\[ G = \int [g(X)f(X)/f\'(X)]f\'(X)dX \]

onde f'(x) é maior ou igual a zero e \( \int f'(X)dX = 1 \)
Se o G for fixo, gostaríamos de encontrar o f'(X) que minimiza a variância de G:

\[ var(G) = \int [g^2(X)f^2(X)/f\'^2(X)]f\'(X)dX - G^2 \]

Resolvendo um lagrangiano com a restrição: \( \int f'(X)dX = 1 \) e com g(X) maior ou igual a zero, chegamos a seguinte solução:

\[ f\'(X) = (g(X)f(X))/G \]

Uma simulação Monte Carlo para avaliar a integral seria uma amostra de series de X de f'(X) e construir a soma: \( G'_{N} = \frac{1}{N}\sum_{i=1}^{N} \frac{g(X_{i})f(X_{i})}{f'(X_{i})} \) . Substituindo o resultado obtido na minimização do Lagrangeano: \( G'_{N} = G \) .
Se soubermos a resposta correta, G, previamente, então a simulação Monte Carlo vai entregar uma variância igual a zero.

Exemplo 1:
Considere a equação: \( G = \int_{0}^{1} cos(\frac{\pi x}{2})dx \), sendo \( g(x) = cos(\frac{\pi x}{2}) \) e f(x) = 1.
Fazendo uma simulação Monte Carlo com x uniformemente distribuído entre 0 e 1, a variância será de aproximadamente, 0.0944.
A imagem será apresentada aqui.
Se expandirmos a g(X) numa série de potência, podemos encontrar a função \( f'(x) = \frac{3}{2} (1 - x^2) \) , fazendo a simulação no novo estimador g'(x) = g(x)/f(x), é possível reduzir a variância para 0.0018, com uma só amostra.

Código exemplo 1:

def est_monte_carlo(n):
    x = np.random.uniform(0,1,n)
    a = []
    for i in range(n):    
        h_x = np.cos(np.pi*x[i]/2)
        a = np.append(a,h_x)
    GN = np.mean(a)
    var = np.var(a)

    return var, GN

def imp_samp(n):
    x = np.random.uniform(0,1,n)
    a = []
    for i in range(n):
        f_x = (3/2)*(1 - x[i]**2)
        h_x = np.cos(np.pi*x[i]/2)/f_x
        a = np.append(a,h_x)
    GN = np.mean(a)
    var = np.var(a)

    return var, GN

if __name__ == '__main__':
    numeroAmostras=1000
    n=1000
    vetor_var = np.empty([numeroAmostras])
    vetor_GN = np.empty([numeroAmostras])    

    for i in range(0,numeroAmostras):
        [var, GN] = imp_samp(n)
        vetor_var[i] = var
        vetor_GN[i] = GN

    plt.figure(0)
    plt.hist(vetor_var)
    print('Variância exemplo 1:', var)
    var, gn = imp_samp(1000)
    print('Variância Import Sampling:',var)

Exemplo 2:
Considere: \( \int_{0}^{1} \sqrt[]{(1-x^2)} = \pi/4 \), a variância associada ao estimador dessa integral num processo Monte Carlo é aproximadamente 0.0465.
A imagem será apresentada aqui.
Expandindo numa série de potência encontramos uma função de importância: \( f'(x) = \frac{1-Bx^2}{1-B/3} \)
Escolhendo um B = 0.5 conseguimos reduzir a variância para 0.016.
A imagem será apresentada aqui.

O código utilizado:

def est_monte_carlo2(n):
    x = np.random.uniform(0,1,n)
    a = []
    for i in range(n):    
        h_x = (1 - x[i]**2)**(1/2)
        a = np.append(a,h_x)
    GN = np.mean(a)
    var = np.var(a)

    return var, GN

def imp_samp2(n, B):
    x = np.random.uniform(0,1,n)
    a = []
    for i in range(n):    
        f_x = (1 - B*(x[i]**2))/(1 - B/3)
        h_x = ((1 - x[i]**2)**(1/2))/f_x
        a = np.append(a,h_x)
    GN = np.mean(a)
    var = np.var(a)

    return var, GN

if __name__ == '__main__':
    numeroAmostras=1000
    n=1000
    vetor_var = np.empty([numeroAmostras])
    vetor_GN = np.empty([numeroAmostras])    
    vetor_var1 = np.empty([numeroAmostras])
    vetor_GN1 = np.empty([numeroAmostras])  

    for i in range(0,numeroAmostras):
        [var, GN] = est_monte_carlo2(n)
        [var1, GN1] = imp_samp2(n, 0.5)
        vetor_var[i] = var
        vetor_GN[i] = GN
        vetor_var1[i] = var1
        vetor_GN1[i] = GN1
    plt.figure(0)
    plt.hist(vetor_var)
    print('Variância exemplo 2:', var)
    plt.figure(1)
    plt.hist(vetor_var1)
    print('Variância importance sampling ex2:',var1)

Outra utilidade do Importance Sampling é eliminar os casos com variância infinita. Suponha queremos avaliar a integral \( G = \int_{0}^{1} \frac{dx}{\sqrt{x(1-x)}} \)
Um Monte Carlo direto entrega uma variância infinita, para eliminar isso deve-se introduzir uma função f'(x) com singularidades em 0 e 1 (o método para encontrar f'(x) foi descrito na seção 3.4.4 do livro):
\[ f\'(x) = \frac{1}{4\sqrt{x}} + \frac{1}{4\sqrt{1-x}} \]
O estimador se torna:
\[ g\'(x) = \frac{4}{\sqrt{x} + \sqrt{1-x}} \]
E a variância:
\[ var(g) = \int f\'(x)(g\' - G)^2dx < (4 - \frac{4}{\sqrt{2}})^2 = 1.37 \]

comentou 6 dias atrás por Carlos A T Haraguchi (16 pontos)  
Pedro, o seu post está bem didático. As explicações e os exemplos deixam bem claro o ganho obtido com o uso de uma "importance function" no cálculo de uma integral com simulação Monte Carlo. Há somente um inconveniente (não sei se somente comigo) nas fórmulas: estão aparecendo algumas "sujeiras", como "\" antes do apóstrofo. Isso me confundiu um pouco no início.

Uma questão, porém, não consegui entender. Há alguma técnica para inferir a "importance function"? Entendi que ela deve ser uma função densidade de probabilidade que se pareça, no intervalo considerado, com a função que está sendo integrada. Mas, como a expansão em uma série de potências ajuda a achar a função? Infelizmente, também não encontrei explicação na referência (Kalos & Whitlock, 2008).

Por fim, apenas para acrescentar outro modo de ver o ganho do uso da "importance function" é que, se tomarmos a variância como uma restrição (utilizando métodos adaptativos Monte Carlo, por exemplo), a convergência será mais rápida. E isso é importante para reduzir o esforço computacional.
...