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.
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.
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.
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 \]