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 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:
- 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.
- 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.