Question

Como calcular a integral definida de uma função utilizando o método Quase-Monte Carlo?

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 possibilidade interessante é a de substituir o gerador de números aleatórios por uma sequência predeterminada, um método conhecido como Quase-Monte Carlo.

A sequência de Van der Corput pode ser utilizada para a amostragem dos valores \(x_i)\) para o cálculo da integral. A sequência de \(n\) números de Van der Corput na base \(k\) é gerada convertendo a sequência dos naturais \(1...n\) pelo seguinte roteiro:

1. Para cada um dos números na sequência dos naturais, reescrever o número correspondente por meio da mudança de base de notação.
2. Para cada um dos \(m\) algarismos do número calculado no passo anterior, calcular o número de Van der Corput correspondente como \(Algarismo 1\times Base^{-1}, ..., Algorismo m\times Base^{-m}\).

Por exemplo, o número 121 pode ser escrito na base 5 como 441, pois \(4\times5^2+4\times5^1+1\times5^0\) é igual a 121. O número de Van der Corput será dado então por \(1\times5^{-1}+ 4\times5^{-2}+ 4\times5^{-3}= 0,392\).

Os gráficos abaixo ilustram como o quadrado unitário é preenchido com a técnica de Monte Carlo Simples e com a sequência de Van der Corput nas bases 2 e 3. Como é possível observar, o preenchimento é muito mais uniforme com a técnica dos números quase-aleatórios.

O código abaixo implementa a rotina de geração da sequência de Van der Corput e a utiliza 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\).

#=====================================================
# 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 Base10ToBaseN(DecimalNumber,BaseN):
   n=DecimalNumber
   ans=''
   while (n>0):
       resto = n % BaseN
       quociente = n //BaseN
       n=quociente
       ans=(str(resto))+ans
   return ans   
def ReverteOrdemDaString(String):
   #ReverteOrdemDaString('abacaxi')
   #Retorna: 'ixacaba'
   ans=''
   for i in String:
       ans = i+ans
   return ans
def ComputaVanDerCorputDaStringReversa(StringReversa,BaseN):
   #ComputaVanDerCorputDaStringReversa('01',2)
   #Retorna: 0.25, ou 0*2^-1+1*2^-2
   ans=0
   k=0
   for i in StringReversa:
     k=k+1
     ans=ans+float(i)*BaseN**(-k)
   return ans
def SequenciaDeVanDerCorput(n,BaseN):
    ans=list()
    for i in range(n):
        StringBaseN=Base10ToBaseN(i+1,BaseN)
        StringReversa=ReverteOrdemDaString(StringBaseN)
        ans.append(ComputaVanDerCorputDaStringReversa(StringReversa,BaseN))
    return ans  

def MonteCarloIntegralVanDerCorput(fx,lb,ub,n,BaseN,ra=False):
    acum=0
    runningaverage=[]
    sequencia=SequenciaDeVanDerCorput(n,BaseN)
    for i in range(n):
      y=sequencia[i]
      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 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

#=====================================================
# Chamada principal
#=====================================================

#Tamanho da amostra
n=1000
# Base para Van der Corput
baseN=4
#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 Sequência de Van der Corput')
print('.....  E(x):')
print(MonteCarloIntegralVanDerCorput(f,lb,ub,n,baseN,False))

#Gráficos comparação: aleatório x quase-aleatório
import matplotlib.pyplot as plt
x1=list()
y1=list()

for i in range(1000):
    x1.append(random.uniform(0,1))
    y1.append(random.uniform(0,1))

x2=SequenciaDeVanDerCorput(1000,2)
y2=SequenciaDeVanDerCorput(1000,3)


plt.scatter(x1, y1)
plt.title('Monte Carlo Simples')
plt.show()

plt.scatter(x2, y2)
plt.title('Sequência de Van der Corput - Base 2 e 3')
plt.show()

#Gráficos de convergência

x = range(1,n+1)

y = MonteCarloIntegral(f,lb,ub,n,True)
plt.axis([1, n, 5, 8])
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 = MonteCarloIntegralVanDerCorput(f,lb,ub,n,baseN,True)
plt.plot(x, y)
plt.axis([1, n, 5, 8])
plt.ylabel('Valor esperado da integral')
plt.xlabel('Número de simulações')
plt.title('Monte Carlo com Sequência de Van der Corput')
plt.show()

Os resultados da rotina são apresentados a seguir:

Monte Carlo Simples
..... E(x): 7.046164638542711
..... ep(x): 0.08635636479076761

Monte Carlo com Sequência de Van der Corput
..... E(x): 6.986629786491394

Em conclusão, o método de Quase-Monte Carlo foi eficiente para melhoria da precisão das estimativas e aumento da velocidade de convergência para o valor verdadeiro.