Estimando o valor de pi usando o método de Monte Carlo

Ae, denovo com algum tempo livre para mais alguns posts 🙂

Eu já fiz alguns posts sobre Markov Chain, para pelo menos termos uma ideia do que se tratava, mas tudo que a gente faz em estatística Bayesiana é sempre usando o Monte Carlo Markov Chain (MCMC), a gente ja deve ter alguma noção do que se trata o Markov Chain, mas e o tal do Monte Carlo?

Bem eu ja comentei em algum post que não me lembro, que Monte Carlo é um bairro famoso de Mônaco por seus cassinos, e cassino é o lugar da probabilidade. Basicamente, o método de Monte Carlo é um algorítimo para se obter soluções numéricas de optimização, integral e amostrar distribuições probabilísticas.

O método de Monte Carlo mais ou menos é o seguinte:

1. Defina um dominio de todas as possibilidades 2. Gere aleatoriamente entradas de uma distribuição de probabilidade que cobre o dominio. 3. Faça uma computação deterministica das entradas 4. Reuna os resultados

Certo, num entendi porra nenhuma. Eu também não, mas é tentando com algo fácil que sabemos o resultado que a gente aprende.

Então vamos calcular a area de um circulo:

Então desenhamos um circulo no R

01

Certo, podemos ver que o raio desse circulo é 1. E sabemos que a área do círculo é

 Area = \pi \cdot r^{2}

Então a área desse exemplo tem que ser:

> #Area do circulo > pi*1^2 [1] 3.141593

Certo, mas como aplicar aquela ideia do Monte Carlo?
Bem o primeiro passo é definir um domínio de todas as possibilidades, vamos supor que a gente não soube-se o valor de \pi(pi), a gente tava ferrado porque não tinha como calcular a área, mas uma figura que é mais simples de deduzir qual a área é o quadrado. Então a gente desenha um quadrado em volta desse círculo.

02

Sabemos aqui que o lado desse quadrado é 2 e a área do quadrado é:

 Area = lado^{2}

Então a área do quadrado é:

> #Area do Quadrado > 2^2 [1] 4

Pode parecer idiota, mas se o nosso objetivo aqui é saber a área do circulo dentro do quadrado, podemos dizer que ele contem toda a área do circulo e mais um pouquinho.

O primeiro passo ta pronto. Definimos algo que contém todas as possibilidades.

Agora vamos para o passo dois, gere entradas de uma distribuição que cobre todo o domínio.

Vamos la, pense que se quisermos colocar um pontinho nesse quadrado num lugar qualquer, se sortearmos um ponto no eixo x ao acaso de uma distribuição uniforme dos valores -1 até o valor 1, e o mesmo para o eixo y, um ponto ao acaso de uma distribuição uniforme dos valores -1 até o valor 1, teremos um ponto x e y , que com certeza vai cair dentro desse quadrado.

Podemos simular isso no r. Então sorteamos um ponto que são duas coordenadas, x e y, de uma distribuição uniforme que vai de -1 a 1.

Então basicamente estamos fazendo é isso.

03

Bang, jogamos um pontinho ali dentro, e ele caiu dentro do circulo legal. Se funcionou, vamos gerar uns 500 pontos desses dai e ver o que vai acontecer

04

Olha que legal, os pontinhos vermelhos cairam dentro do círculo e os azuis cairam fora do circulo. Tem muito mais pontos vermelhos que azuis certo, isso porque a área do circulo é grande, alias, ela é uma proporção da área do quadrado. Proporção é a palavra, se o quadrado é o total, então a parte que queremos dividido pelo total é uma proporção que nos permite saber da area total quanto é área do circulo, ou seja, essa proporção vezes a área do quadrado é a área do circulo. Mas como a gente vê essa proporção a partir desse monte de ponto sorteados ao acaso, vamos fazer uma simulação maior, cinco mil pontos, então o número de pontos que ficaram dentro do circulo dividido pelo total de pontos que simulamos no quadrado é a proporção da área entre os dois.

Então fazemos a simulação da seguinte forma:

#Fazendo essa simulação 5 mil vezes e olhando a proporção de pontos dentro do circulo
n<-5000
x.pos <- runif(n, min=-1, max=1)
y.pos <- runif(n, min=-1, max=1)
local.pos <- ifelse(x.pos^2 + y.pos^2 <= 1, TRUE, FALSE)
dentro <- length(which(local.pos == TRUE))

Escolhemos um tamanho n para a simulação, cinco mil, dai sorteamos 5 mil pontos para x e para y, veja que é exatamente o que tamo fazendo no gráfico ali visualmente, só que 5 mil pontos ia deixar poluído o gráfico.

Dai o passo três é fazermos uma computação determinística, olhamos se caiu dentro ou fora do circulo. Se esta difícil pegar o que esta quarta linha esta fazendo, so lembrar pra que serve o famoso  c^2 = a^2 + b^2 , conhecido como teorema de pitágoras, mas eu também apanhei para entender essa parte, mas aqui tem um tutorialzinho legal para entender como fazer um gráfico de um circulo.
Por ultimo, agregamos os resultados, vemos a razão de pontos dentro do circulo do total de pontos que simulamos.

Temos nesse caso que:

> #Propoção > dentro/n [1] 0.786

Mas a gente sabe qual deve ser essa proporção, com infinitos pontos, a area do circulo dividido pela area do quadrado que é o total.
E para o estamos, essa proporção da.

> #Exato > (pi*1^2) / (2^2) [1] 0.7853982

Puxa vida, como diabos esse treco funciona heheheh, é impressionante. Se antes tava complicado de entender o algorítimo de Metropolis-hasting, olha ai como é bem parecido, ficamos gerando alguma coisa, temos um função pra aceitar ou não essas coisas que geramos e depois vemos uma distribuição.

Agora se multiplicarmos essa proporção pela área do quadrado, que no exemplo é 4, redescobriremos a área do circulo, que é o mesmo valor de \pi nesse caso, estimamos o valor de \pi, que legal eeeee.

> 4*(dentro/n) [1] 3.144

Bem mas podemos ainda ver o seguinte, com 1 pontinho apenas, a estimativa seria péssima, porque seria ou 1/1 ou 0/1, que daria 1 ou 0, que multiplicando por 4, daria ou 4 ou 0 a nossa estimativa de pi, mas com 2 pontos, ja temos 5 possibilidades 4/4, 3/4, 2/4, 1/4 e 0/4, note que quanto mais pontos sorteamos, mais possibilidades de valores, mas o que vai acontecendo de acordo com que sorteamos mais pontos, de acordo com que iteramos mais nessa ideia, nesse algorítimo?

pi

Olha ai, no começo o valor oscila muito, e vai de ruim a pior, mas conforme usamos mais pontos, mais iterações, a estimativa fica melhor e melhor, e começa a variar cada vez menos, ficamos mais e mais próximos do valor de pi, mas vamos chegar no valor de pi?
Ai eu ja não entendo, pi não é um número real, então num tem como saber exatamente, poxa mas a estimativa é muito boa. O valor de pi mesmo que temos no R tem um número finito de casas decimais, enquanto pi realmente tem um número infinito de casas decimais. Mas isso é uma discussão complexa e recorrente em fóruns, mas tem mais haver com computação que matemática nesse caso, mas tem um comentário interessante no faq do R aqui.

Mas beleza, agora a gente ja tem uma ideia melhor do que Monte Carlo, do que é Markov Chain, então da para tentar progredir pro MCMC talvez entendo melhor o que ta acontecendo.

Bem segue o script, ele esta no repositório do github também e até o próximo post 🙂

set.seed(171)
t <- seq(0.2*pi,length=1000)
coords <- t(rbind(0+sin(t)*1,0+cos(t)*1))
 
#desenhando um circulo
#jpeg("01.jpg")
plot(0,0,type="n",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5),frame=F,xlab="Eixo x",ylab="Eixo y")
points(coords,cex=0.1)
abline(h=c(0),v=c(0),lty=3,lwd=0.5)
lines(c(0.01,0.99),c(0.01,0.01),lwd=2,lty=2,col="red")
text(0.5,0.1,"Raio=1")
#dev.off()
 
#Desenhando um quadrado sobre o circulo
#jpeg("02.jpg")
plot(0,0,type="n",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5),frame=F,xlab="Eixo x",ylab="Eixo y")
points(coords,cex=0.1)
abline(h=c(0),v=c(0),lty=3,lwd=0.5)
lines(c(-1,-1,-1,1,1,1,1,-1),c(-1,1,1,1,1,-1,-1,-1),lwd=2,col="darkgray")
lines(c(-1.01,1.01),c(-1.01,-1.01),lwd=2,lty=2,col="red")
text(0,-1.08,"Lado=2")
#dev.off()
 
#Area do circulo
pi*1^2
 
#Area do Quadrado
2^2
 
#proporção
(pi*1^2) / (2^2)
 
#Jogando um ponto ao acaso dentro do quadrado
#jpeg("03.jpg")
plot(0,0,type="n",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5),frame=F,xlab="Eixo x",ylab="Eixo y")
points(coords,cex=0.1)
abline(h=c(0),v=c(0),lty=3,lwd=0.5)
lines(c(-1,-1,-1,1,1,1,1,-1),c(-1,1,1,1,1,-1,-1,-1),lwd=2,col="darkgray")
 
x.pos <- runif(1, min=-1, max=1)
y.pos <- runif(1, min=-1, max=1)
local.pos <- ifelse(x.pos^2 + y.pos^2 <= 1, TRUE, FALSE)
points(x.pos,y.pos,pch=19,cex=0.5,col=ifelse(local.pos,"red","blue"))
#dev.off()
 
#Jogando 500 pontos ao acaso dentro do quadrado
#jpeg("04.jpg")
plot(0,0,type="n",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5),frame=F,xlab="Eixo x",ylab="Eixo y")
points(coords,cex=0.1)
abline(h=c(0),v=c(0),lty=3,lwd=0.5)
lines(c(-1,-1,-1,1,1,1,1,-1),c(-1,1,1,1,1,-1,-1,-1),lwd=2,col="darkgray")
 
for(i in 1:500) {
x.pos <- runif(1, min=-1, max=1)
y.pos <- runif(1, min=-1, max=1)
local.pos <- ifelse(x.pos^2 + y.pos^2 <= 1, TRUE, FALSE)
points(x.pos,y.pos,pch=19,cex=0.5,col=ifelse(local.pos,"red","blue"))
}
#dev.off()
 
#Fazendo essa simulação 5 mil vezes e olhando a proporção de pontos dentro do circulo
n<-5000
x.pos <- runif(n, min=-1, max=1)
y.pos <- runif(n, min=-1, max=1)
local.pos <- ifelse(x.pos^2 + y.pos^2 <= 1, TRUE, FALSE)
dentro <- length(which(local.pos == TRUE))
 
#Propoção
dentro/n
#Exato
(pi*1^2) / (2^2)
 
#Multiplicando pera area do quadrado para termos o valor de pi
4*(dentro/n)
#Valor de pi
pi
 
#Guardando as estimativas num vetor
estimativa<-vector()
 
for(i in 1:5000) {
    estimativa[i] <- length(which(local.pos[1:i] == TRUE)) / (i)
}
 
#Grafico animado
 
passos<-round(seq(1,5000,length=200))
estimativa<-estimativa*4
 
for(i in passos) {
jpeg(sprintf("pi%05d.jpg",i), width = 300, height = 300)
plot(c(1:5000),estimativa[1:5000]*4,type="n",ylim=c(min(estimativa),max(estimativa)),xlab="Iterações",
     ylab="Valor de Pi",main=paste("Iteração",i),frame=F)
points(c(1:i),estimativa[1:i],type="l",col="blue",lwd=2)
abline(h=pi,lty=3,col="red")
dev.off()
}
system("convert pi*.jpg -delay 10 pi.gif")
system("rm pi*.jpg")

6 thoughts on “Estimando o valor de pi usando o método de Monte Carlo

  1. Muito bom todos os seus posts. Parabéns. Seus posts tem me ajudado muito. Gostaria de saber como fazer para programar a posteriori do local de cada aranha do arbusto do exemplo.
    suponha que eu coloquei coordenadas e fiz tipo um plano cartesiano das árvores e quero saber como posso repartir o espaço(plano) que contém aranhas em quantidade alta, média ,baixa e sem aranhas, por exemplo. Como achar a probabilidade de dividir esta área e a probabilidade de cada parte da área? como achar estas probabilidades identificadas a localização? Estou precisando conseguir isso e não sei como fazer.

Leave a Reply

Your email address will not be published. Required fields are marked *