Método de Monte Carlo e transformações

Transformação linear

Sorteando um número aleatório ${\textstyle x\in \left[0,1\right)}$ então fazemos uma transformação para obter um número ${\textstyle y\in \left[a,b\right)}$. Isto é, obtemos a seguinte transformação ${\textstyle f:\left[0,1\right)\rightarrow \left[a,b\right)}$ da seguinte forma:

$$y\left(x\right)=a+\left(b-a\right)x$$

Se todos os números entre ${\textstyle 0}$ e ${\textstyle 1}$ tinham igual probabilidade de serem sorteados, após a transformação todos os números entre ${\textstyle a}$ e ${\textstyle b}$ também possuem igual probabilidade, pois ${\textstyle y}$ varia com ${\textstyle x}$ de forma linear, isto é, a distribuição uniforme de números é mantida.Por sua vez, a distribuição uniforme significa que a probabiliade de obter um número entre ${\textstyle x}$ e ${\textstyle x+dx}$ é dada pela função densidade de probabilidade da seguinte forma:

$$p\left(x\right)dx={\begin{cases}dx,&0<x<1\\0,&{\text{outra forma}}\end{cases}}$$

Sendo ${\textstyle p\left(x\right)=p}$ constante, temos que:

$$\int _{-\infty }^{+\infty }dx=p\int _{0}^{1}dx=p=1$$

Pela normalização. Mas se ampliarmos o intervalo dos números possíveis para entre ${\textstyle a}$ e ${\textstyle b}$:

$$\int _{-\infty }^{+\infty }p'dx=\int _{a}^{b}p'dx=p'\int _{a}^{b}dx=p'\left(b-a\right)=1$$

Então agora ${\textstyle p\left(y\right)=p'=1/\left(b-a\right)}$. Isto é distribuição de probabilidade continua constante, mas com uma menor probabilidade de sortear um número ${\textstyle y}$ qualquer, quando em comparação de sortear um ${\textstyle x}$ qualquer.

Transformação não-linear

O mesmo não ocorre com uma transformação não linear. Por exemplo se ${\textstyle y\left(x\right)=4x^{2}}$, derivando temos que:

$${\frac {dy}{dx}}=8x\quad \longrightarrow \quad dy=8xdx$$

Diferente do caso anterior que tínhamos apenas a transformação linear ${\textstyle dy=\left(b-a\right)dx}$. Podemos ver ainda usando a própria definição de derivada:

$${\begin{aligned}y\left(x+dx\right)-y\left(x\right)&=4\left(x+dx\right)^{2}-4x^{2}\\&=4x^{2}+4xdx+dx^{2}-4x^{2}\\&=4\cdot 2xdx+dx^{2}\end{aligned}}$$

E sendo ${\textstyle dx}$ um diferencial então ${\textstyle dx^{2}\approx 0}$, logo ${\textstyle dy=8xdx}$. Agora a distribuição de probabilidade é alterada com a transformação. Para uma transformação ${\textstyle y\left(x\right)}$ qualquer, como a probabilidade se conserva, ainda temos:

$$\left|p\left(y\right)dy\right|=\left|p\left(x\right)dx\right|$$

Considrando que ${\textstyle y=y\left(x\right)}$ tem inversa, então ${\textstyle x=x\left(y\right)}$ então:

$${\begin{aligned}\left|p\left(y\right)dy\right|=&\left|p\left(x\right)dx\right|\\p\left(y\right)\left|dy\right|=&p\left(x\right)\left|dx\right|\\p\left(y\right)=&p\left(x\right)\left|{\frac {dx}{dy}}\right|\end{aligned}}$$

Sendo ${\textstyle y\left(x\right)=-\ln \left(x\right)}$ então reescrevendo ${\textstyle x=e^{-y}}$, logo ${\textstyle {\frac {dx}{dy}}=-e^{y}}$ Então temos:

$$p\left(y\right)=p\left(x\right)e^{-y}$$

E sendo nossa ditribuição em ${\textstyle x}$ uniforme com ${\textstyle p\left(x\right)=p=1}$ como vimos anteriormente, ficamos com:

$$p\left(y\right)=e^{-y}$$

Para a transformação ${\textstyle y=-\ln \left(x\right)}$. O mais comum é que saibamos a distribuição de probabilidade ${\textstyle p\left(y\right)}$ que queremos, e uma vez que:

$$p\left(y\right)=\left|{\frac {dx}{dy}}\right|$$

E integramos então para encontrar ${\textstyle x\left(y\right)}$.

Método da rejeição

Nem sempre o ${\textstyle p\left(y\right)}$ desejado é fácil de definir matematicamente. O método da rejeição é um método rústico para obtermos ${\textstyle p\left(y\right)}$.

• Desenhar a função ${\textstyle p\left(y\right)}$ desejada dentro dos limites ${\textstyle a\leq y<b}$ e ${\textstyle 0\leq p\left(y\right)\leq p_{max}}$.
• Geramos um ponto aleatório ${\textstyle \left(y,p\left(y\right)\right)}$, se estiver abaixo da curva desejada é aceito.

Se gerarmos pontos aleatórios em grande quantidade, a razão de pontos aceitos para cada ${\textstyle y}$ em relação à todos o pontos aleatórios gerados neste ${\textstyle y}$, nos dá uma estimativa de ${\textstyle p\left(y\right)}$ neste ponto, em relação do ${\textstyle p_{max}}$. Isto é, se para um ${\textstyle y}$ qualquer, metade dos pontos gerado aleatoriamente foram válidos, então ${\textstyle p\left(y\right)\approx {\frac {p_{max}}{2}}}$.

Cálulo de integrais definidas

Em uma ideia bastante análoga à anterior, aqui utilizamos a ideia que a integral definida é cálculo da área sob a curva. Definindo novamente limites ${\textstyle x_{min}\leq x\leq x_{max}}$ e ${\textstyle y_{min}\leq y\leq y_{max}}$ no qual a região que queremos calcular está contida, geramos uma série de pontos aleatórios, se gerado pontos suficientes, a razão de pontos que foram gerados dentro da área que queremos calcular em relação ao total de pontos gerados, será a mesma razão da área que queremos calcular em relação à área total definida pelos limites.

#Cálculo da área de um círculo
import matplotlib.pyplot as plt   #Plotar gráfico
import numpy as np                #Funções matemáticas
import random                     #Números aleatórios

r=1                               #Raio do círculo
Np=1000000                        #Número de pontos
hitx=[];hity=[];missx=[];missy=[] #Gurdar resultados
for n in range(Np):
  x=random.random()*2-1           #Pontos aleatórios gerados entre -1 e +1 
  y=random.random()*2-1           #Pontos aleatórios gerados entre -1 e +1
  if(np.sqrt(x**2+y**2)<=r):      #Se o ponto caiu dentro do círculo
    hitx.append(x);hity.append(y)
  else:
    missx.append(x);missy.append(y)

#Plotar o círculo
x=np.arange(-1,1,1E-5)
plt.plot(x,-np.sqrt(r-x**2),'-k')
plt.plot(x,np.sqrt(r-x**2),'-k')
#Plotar pontos
plt.plot(hitx,hity,'or')
plt.plot(missx,missy,'ob')

#Cálcular a área
a_ret=2*2                         #Área do retângulo: geramos pontos entr -1 e +1 na duas dimensões, então é um quadrado de lado 2
a_cir=a_ret*len(hitx)/Np          #A_ret *(hit/N)
print('A área do círculo estimada é: {:.2f}'.format(a_cir))
print('A área do círculo exata    é: {:.2f}'.format(np.pi*r**2))