O cálculo da integral definida por simulação de Monte Carlo de uma função qualquer \(f(x)\) em um intervalo (0,1) é possível considerando que a função de densidade e probabilidade de uma variável uniforme \(pdf(x)\) é igual a 1. Assim, pode-se multiplicar \(f(x)\) por \(pdf(x)\), sem alterar o valor da integral. Além disso, sendo \(x\) uma variável aleatória uniforme, o valor esperado de \(f(x)\) será dado por \(\int_{0}^{1}f(x)dx\) e que, portanto, poderá ser aproximado por uma amostra de tamanho \(n\) como \(\sum_{i=1}^{n}\frac{f(x_i)}{n}\). O passo a passo para o cômputo da integral de Monte Carlo é dado por:
- Converter o domínio da integral definida \(\int_{lb}^{ub}f(x)dx\) para o intervalo (0,1). Isto pode ser feito transformando a variável \(x\) em uma variável \(y=ax+b\). Assim, resolvendo \(a\) e \(b\) em função de \(lb\) e \(ub\), encontramos que \(x= (ub-lb)y+lb\) e \(dx = (ub-lb)dy\). Portanto, a integral equivalente a ser avaliada é \(\int_{0}^{1}f((ub-lb)y+lb)(ub-lb)dy\).
- Calcular a integral como o valor esperado de \(f(y)\) fazendo \(E(f(y))= \sum_{i=1}^{n}\frac{f(y_i)}{n}\)
- O erro da estimativa é dado como \((V(f(y))/n)^{0,5}\) e usando a propriedade do operador variância \(V(X) = E(X^2) – E(X)^2\), podemos calcular \(V(f(y))= E(f(y)^2)-E(f(y))^2\)
Ainda existem algumas alterações possíveis para aumentar a precisão da estimativa de Monte Carlo. Uma alternativa simples é a técnica da variável antitética em que para cada sorteio de uma variável uniforme \(y\), fixa-se um novo sorteio como \(1-y\). A intuição do método é que a amostra será melhor distribuída em torno da média, facilitando a convergência de \(f(y)\) para a média.
Outra alternativa é a estratificação da amostra, técnica na qual o intervalo (0,1) é dividido em intervalos menores (não necessariamente de igual tamanho) e \(f(y)\) é calculado como a média ponderada das médias de cada um desses intervalos. A intuição do método é que a amostra será melhor distribuída em torno do intervalo (0,1), minimizando a chance de clusters em intervalos específicos.
O código a seguir escrito em Python executa as três rotinas descritas acima (Monte Carlo Simples, Monte Carlo com Variáveis Antitéticas e Monte Carlo com Estratificação) para o cálculo da integral da função \(f(x)=3x^2\) no intervalo (1,2). Neste caso, a integral escolhida possui fórmula fechada fazendo \(F(x) = x^3\) e \(F(2)-F(1)=7\). Em outras integrais, contudo, não é possível o uso de fórmula fechada, sendo, portanto, muito útil a integração por Monte Carlo.
#=====================================================
# Definição de funções
#=====================================================
import random
import statistics as stat
def fy(fx,lb,ub,y):
x=(ub-lb)*y+lb
return fx(x)*(ub-lb)
def MonteCarloIntegral(fx,lb,ub,n,ra=False):
acum=0
runningaverage=[]
for i in range(n):
y=random.uniform(0,1)
acum = acum+fy(fx,lb,ub,y)
runningaverage.append(acum/(i+1))
if(ra==False):
ex=acum/n #E(x)
ans=ex
if (ra==True):
ans =runningaverage
return ans
def MonteCarloIntegralAnthitetic(fx,lb,ub,n,ra=False):
acum=0
runningaverage=[]
for i in range(n):
y=random.uniform(0,1)
acum = acum+fy(fx,lb,ub,y)+fy(fx,lb,ub,1-y)
runningaverage.append(acum/(2*(i+1)))
runningaverage.append(acum/(2*(i+1)))
if(ra==False):
ex=acum/(2*n) #E(x)
ans=ex
if (ra==True):
ans =runningaverage
return ans
def MonteCarloIntegralEstratificado(fx,lb,ub,n,NumeroDeEstratos,ra=False):
acum=0
runningaverage=[]
tamanhoestrato=1/NumeroDeEstratos
k=0
for j in range(NumeroDeEstratos):
for i in range(int(n/NumeroDeEstratos)):
k=k+1
y=random.uniform(0,1)*tamanhoestrato+j*tamanhoestrato
acum = acum+fy(fx,lb,ub,y)
runningaverage.append(acum/(k))
if(ra==False):
ex=acum/n #E(x)
ans=ex
if (ra==True):
ans =runningaverage
return ans
#=====================================================
# Chamada principal
#=====================================================
#Tamanho da amostra
n=1000
#Função a ser integrada
def f(x):
return 3*x**2
lb=1
ub=2
#Valor esperado das integrais e erro padrão das integrais
print('Monte Carlo Simples')
print('..... E(x):')
print(MonteCarloIntegral(f,lb,ub,n,False))
print('..... ep(x):')
temp=list()
for i in range(100):
temp.append(MonteCarloIntegral(f,lb,ub,n,False))
print(stat.stdev(temp))
print('Monte Carlo com Variáveis Antitéticas')
print('..... E(x):')
print(MonteCarloIntegralAnthitetic(f,lb,ub,int(n/2),False))
print('..... ep(x):')
temp=list()
for i in range(100):
temp.append(MonteCarloIntegralAnthitetic(f,lb,ub,int(n/2),False))
print(stat.stdev(temp))
print('Monte Carlo com Estratificação')
print('..... E(x):')
print(MonteCarloIntegralEstratificado(f,lb,ub,n,10,False))
print('..... ep(x):')
temp=list()
for i in range(100):
temp.append(MonteCarloIntegralEstratificado(f,lb,ub,n,10,False))
print(stat.stdev(temp))
#Gráficos de convergência
import matplotlib.pyplot as plt
x = range(1,n+1)
y = MonteCarloIntegral(f,lb,ub,n,True)
plt.axis([1, n, 3, 10])
plt.plot(x, y)
plt.ylabel('Valor esperado da integral')
plt.xlabel('Número de simulações')
plt.title('Monte Carlo Simples')
plt.show()
y = MonteCarloIntegralAnthitetic(f,lb,ub,int(n/2),True)
plt.plot(x, y)
plt.axis([1, n, 3, 10])
plt.ylabel('Valor esperado da integral')
plt.xlabel('Número de simulações')
plt.title('Monte Carlo com Variáveis Antitéticas')
plt.show()
y = MonteCarloIntegralEstratificado(f,lb,ub,n,10,True)
plt.plot(x, y)
plt.axis([1, n, 3, 10])
plt.ylabel('Valor esperado da integral')
plt.xlabel('Número de simulações')
plt.title('Monte Carlo com Estratificação')
plt.show()
Os outputs do código, apresentados a seguir, compreendem o valor esperado da integral, o erro padrão da estimativa e um gráfico de convergência em função do tamanho da amostra.
Monte Carlo Simples
..... E(x):
6.8674639729439635
..... ep(x):
0.08664225839190447

Monte Carlo com Variáveis Antitéticas
..... E(x):
6.999674514060293
..... ep(x):
0.009715483320812318

Monte Carlo com Estratificação
..... E(x):
6.998639691869883
..... ep(x):
0.008343465510077329

Em conclusão, para o exemplo estudado, a técnica de Monte Carlo Simples apresentou uma estimativa de 6,867464 para o valor da integral, com um erro padrão de 0,086642. Usando as técnicas de redução de variância acima descritas o valor da integral foi de 6,999675 (erro padrão de 0,009715) para as variáveis antitéticas e de 6,998639 (erro padrão de 0,008343) para o método da estratificação. Portanto, as técnicas de redução de variância melhoram a precisão da estimativa e reduziram a incerteza associada ao seu valor verdadeiro.
A convergência foi alcançada com poucas simulações com o método das variáveis antitéticas. O método da estratificação, como esperado, precisou percorrer todo o espaço amostral para a obtenção de uma estimativa válida para a integral.