Question

Como calcular a integral definida de uma função utilizando a técnica de simulação de Monte Carlo? Como é possível aumentar a precisão da estimativa utilizando o método das variáveis antitéticas e o método da estratificação.

Answer 1

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:

1. 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\).
2. Calcular a integral como o valor esperado de \(f(y)\) fazendo \(E(f(y))= \sum_{i=1}^{n}\frac{f(y_i)}{n}\)
3. 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.

Comments

Também é possível combinar o método de variáveis antiéticas com o método de estratificação. Selecionando Xi de forma:

\[ \alpha R_{i} = \int_{-\inf}^{X_{i}}f(x)dx \]

e

\[ \alpha + (1-\alpha) R_{i} = \int_{-\inf}^{X'_{i}}f(x)dx \]

um estimador não-viesado seria:

\[ 1/N\sum \alpha g(X_{i}) + (1 - \alpha)g(X'_{i}) \]