Método de Monte Carlo e transformações

Transformação linear

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

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

Se todos os números entre $0$ e $1$ tinham igual probabilidade de serem sorteados, após a transformação todos os números entre $a$ e $b$ também possuem igual probabilidade, pois $y$ varia com $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 $x$ e $x+dx$ é dada pela função densidade de probabilidade da seguinte forma:

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

Sendo $p\left(x\right)=p$ constante, temos que: $$\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}dx=p\underset{0}{\overset{1}{\int }}dx=p=1$$

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

$$\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}{p}^{\prime }dx=\underset{a}{\overset{b}{\int }}{p}^{\prime }dx=p^{\prime }\underset{a}{\overset{b}{\int }}dx=p^{\prime }(b-a)=1$$

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

Transformação não-linear

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

$$\frac{dy}{dx}=8x⟶dy=8xdx$$

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

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

E sendo $dx$ um diferencial então $d{x}^2\approx 0$, logo $dy=8xdx$. Agora a distribuição de probabilidade é alterada com a transformação. Para uma transformação $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 $y=y\left(x\right)$ tem inversa, então $x=x\left(y\right)$ então:

$$\begin{array}{rl}\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{array}$$

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

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

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

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

Para a transformação $y=-\ln \left(x\right)$. O mais comum é que saibamos a distribuição de probabilidade $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 $x\left(y\right)$.

Método da rejeição

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

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

Se gerarmos pontos aleatórios em grande quantidade, a razão de pontos aceitos para cada $y$ em relação à todos o pontos aleatórios gerados neste $y$, nos dá uma estimativa de $p\left(y\right)$ neste ponto, em relação do $p_{max}$. Isto é, se para um $y$ qualquer, metade dos pontos gerado aleatoriamente foram válidos, então $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 $x_{min}\le x\le x_{max}$ e $y_{min}\le y\le 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=[] #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 curvas
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')

a_ret=2*2 #Área do retângulo
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))