Método de Monte Carlo e o experimento da Agulha de Buffon

Desde que eu comecei a ler sobre estatística bayesiana (nossa já faz alguns anos), eu comecei a sempre olhar como funcionava os algorítimos dela, que sempre caia no MCMC. Um MC desses é do Monte Carlo, que é um método de amostragem bem legal, já falamos dele aqui quando estimamos o valor de pi, mas agora vamos ver outro experimento com ele bem legal, o experimento das Agulhas de Buffon.

Bem biólogos sempre estão enfiados no meio da análise de dados, desde o Fisher fazendo anova a muitos outros nomes, e não é que no Monte Carlo a gente vê isso também. O Georges-Louis Leclerc, conde de Buffon era um botânico que propôs e estudou esse experimento.

O experimento é o seguinte, ele jogava várias agulhas no chão de madeira.

Então, depois contava quantas dessas agulhas encostavam nas divisões das madeiras.
Como ninguém tem muitas agulhas e as vezes não tem piso de madeira, podemos simular o experimento no R. Fazemos nosso piso.

E agora é so jogar Agulhas nele.

Jogamos uma agulhinha ali, e ela caiu no meio do piso, ou seja, ela não acertou a junção do piso. Mas se jogarmos várias.

Algumas vão acertar, e outras não. E agora é só contar. As agulhas ficaram legais né, são uma reta com um pontinho na ponta hehehe. Mas como eu fiz as agulhas? Foi com o seguinte código.

1
2
b=1
a=0.3

Bem primeiro definimos duas coisas, o tamanho das tábuas e o tamanho das agulhas, as tábuas aqui são de 1 unidade de tamanho e a agulha de 0.3, chamamos de unidade porque se for metro ou centímetro não tem real interesse aqui na simulação. Depois jogamos a agulha.

1
2
3
4
5
6
7
8
9
10
11
    xcenter <- runif(1,0,4*b)
    ycenter <- runif(1,0,4*b)
 
    phi <- runif(1,0,2*pi)
    x0 <- xcenter-(a/2)*cos(phi)
    y0 <- ycenter-(a/2)*sin(phi)
    x1 <- xcenter+(a/2)*cos(phi)
    y1 <- ycenter+(a/2)*sin(phi)
 
    points(x=c(x0,x1),y=c(y0,y1),type="l")
    points(x=x0,y=y0,type="p",pch=19,cex=0.6)

Eu sorteio onde é o meio da agulha, depois calculo onde ela começa no x0 e y0 e onde ela termina, o x1 e y1, que é bem fácil usando seno e cosseno, dai no gráfico eu ploto uma linha e uma bolinha na ponto, para dar um charme. Agora para testar se a agulha encosta borda da do piso, podemos começar a pensar de forma mais simplificada o problema. Veja que temos três variáveis que realmente interessam.

Uma variável é onde estão as ranhuras da madeira, que como a tábuas tem tamanho b, será de b em b na nossa área de amostragem, e onde as agulhas caem, que pode ser representado pelo ponto central onde ela caiu e o ângulo \phi dela em relação ao chão, e o tamanho da agulha que é a.

As agulhas não interagem umas com as outras, nem rolam para dentro das junções do piso, além que de as tábuas são simétricas, então x_{center} \leftrightarrow b-x_{center} e \phi \leftrightarrow -\phi.

Então de forma geral, podemos considerar um plano reduzido de amostragem, simplificado na verdade que é:

 0  \textless \phi   \textless \frac{\pi}{2}

e,

 0  \textless x_{center}   \textless \frac{b}{2}

Dessa forma, podemos testar se a agulha encosta no canto da madeira se x_{tip} \textless 0 onde x_{tip} = x_{center} - \frac{a}{2} cos (\phi).

Confesso que fiquei meio perdido quanto a isso no começo, mas veja que o centro da agulha sempre vai estar entre  0  \textless x_{center}   \textless \frac{b}{2} e temos a fração de b porque passou de \frac{b}{2} temos o outro lado da madeira, lembre-se da simetria x_{center} \leftrightarrow b-x_{center}. O cosseno vem da inclinação da agulha, estamos tipo projetando ela no eixo x, porque o eixo y não interessa aqui. Assim como temos sempre algo entre 0 e \frac{b}{2}, a “sombra” da agulha no eixo x vai ser negativo quando tivermos passado o canto do piso.

No caso das agulhas que estamos fazendo no meu código, podemos sempre avaliar se o tamanho dela ta certinho com:

1
sqrt((x1-x0)^2+(y1-y0)^2)

E para fazermos o teste dessa forma, usamos o modulo em b e em pi para ficarmos nos intervalos acima, e usar essa simetria

1
xcenter%%(b/2)-(a/2)*cos(phi%%(pi/2))

Agora vamos fazer um experimento dessa forma mais simples, duas mil vezes, para ver algo legal, e contar quantas vezes acertamos os cantos do piso

1
2
3
4
5
6
7
8
9
10
11
12
n <- 2000
b <- 1
a <- 0.7
nhits <- 0
for(i in 1:n){
    xcenter <- runif(1,0,b/2)
    phi <- runif(1,0,pi/2)
    xtip <- xcenter-(a/2)*cos(phi)
    if(xtip<0){
        nhits <- nhits+1
        }
}

Que vai resultar em:

> nhits [1] 877

Nesse caso acertamos 877 vezes, que da

> nhits/n [1] 0.4385

Mais ou menos 44 porcento, que é mais ou menos

> ((a/b)*(2/pi)) [1] 0.445633

Hã? O que está acontecendo? Bem, não é de se admirar que o valor encontrado está relacionado a pi, mas veja que a é o tamanho da agulha, e intuitivamente sabemos que quanto mair agulha, mais chance ela tem de encostar no canto do piso certo? Mas quanto maior o piso, menor a chance da agulha encostar no canto, pois tem mais piso para ela cair sem encostar em nada. Veja que podemos olhar, a partir da nossa fórmula para testar se caiu ou no canto, o valor de N_{hit} o que ta acontecendo quando variamos essas quantidades.

Primeiro da agulha em relação a N_{hit} para vários tamanhos de pisos.

E depois o contrário

Veja que dependemos dessas duas variáveis, na verdade, podemos plotar todos os casos juntos

Então ((a/b)*(2/pi)) nada mais é que essa área azul, a integral disso, infelizmente eu não sei cálculo muito bem, ainda mais com duas variáveis, então nem vou tentar, mas a intuição é essa do resultado do experimento, se essa figura é o total das possibilidades, a área azul é a porcentagem desse total que são quando a agulha encosta no canto. Da até um pouco de vergonha que em 1700 o cara já dava conta de calcular isso, sem computador, e hoje eu apanho para entender. Pior ainda o fato do cosseno, que devia ser bem complicado sem computador.

Bem é isso ai, MCMC, como quase tudo em estatística, tem sempre uns biólogos perdidos no meio, acho que principalmente porque precisamos de problemas para resolver, e biólogo sempre observa mais problema na natureza para resolver o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail e bom natal a todos 🙂

Referência:
Werner Krauth 2006 – Statistical Mechanics: Algorithms and Computations Oxford University Press 342pp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
##Ilustrando o problema de buffon
b=1
a=0.3
 
 
plot(0,0,type="n",xlim=c(0,4*b),ylim=c(0,4*b),frame=F,xlab="",ylab="",main="Piso",xaxt="n",yaxt="n")
abline(v=seq(0,3*b,b))
polygon(x=c(0,1,1,0),y=c(-1,-1,5,5),col="brown3")
polygon(x=c(2,3,3,2),y=c(-1,-1,5,5),col="brown3")
 
 
for(i in 1:100){
    xcenter <- runif(1,0,4*b)
    ycenter <- runif(1,0,4*b)
 
    phi <- runif(1,0,2*pi)
    x0 <- xcenter-(a/2)*cos(phi)
    y0 <- ycenter-(a/2)*sin(phi)
    x1 <- xcenter+(a/2)*cos(phi)
    y1 <- ycenter+(a/2)*sin(phi)
 
    points(x=c(x0,x1),y=c(y0,y1),type="l")
    points(x=x0,y=y0,type="p",pch=19,cex=0.6)
}
 
##Tamanho das agulhas
sqrt((x1-x0)^2+(y1-y0)^2)
 
##Testando se ela encostou no canto
xcenter%%(b/2)-(a/2)*cos(phi%%(pi/2))
 
 
##Experimento simplificado
n <- 2000
b <- 1
a <- 0.7
nhits <- 0
for(i in 1:n){
    xcenter <- runif(1,0,b/2)
    phi <- runif(1,0,pi/2)
    xtip <- xcenter-(a/2)*cos(phi)
    if(xtip<0){
        nhits <- nhits+1
        }
}
 
##Porcentagem de hits
nhits/n
 
##Estimativa de quanto deve ser a porcentagem de hits
((a/b)*(2/pi))
 
 
##Avaliando a relação entre tamanho do piso, tamanho da agulha e hits
b_vetor<- seq(0.1,b,length.out = 6)
a_vetor<- seq(0,1,0.001)
 
par(mfrow=c(3,2))
for(i in 1:length(b_vetor)){
    plot(a_vetor, ((a_vetor/b_vetor[i])*(2/pi)),type="l",xlab="Valores de a",ylab="N hits",main=paste("b=",round(b_vetor[i],2),sep=""),ylim=c(0,5))
}
 
b_vetor<- seq(0.1,b,0.001)
a_vetor<- seq(0.1,1,length.out = 6)
par(mfrow=c(3,2))
for(i in 1:length(a_vetor)){
    plot(b_vetor, ((a_vetor[i]/b_vetor)*(2/pi)),type="l",xlab="Valores de b",ylab="N hits",main=paste("a=",round(a_vetor[i],2),sep=""),ylim=c(0,5))
}
 
 
##Plotando as duas variáveis juntas
angulo <- seq(0,pi/2,length.out = 500)
vetor <- seq(0,b/2,length.out = 500)
casos <- expand.grid(angulo=angulo,b=vetor)
casos$cor <- ifelse(casos$b<a/2 & casos$angulo < acos(casos$b/(a/2)),"lightblue","brown1")
 
 
plot(casos$b,casos$angulo,col=casos$cor,pch=19,cex=0.4,xaxt="n",yaxt="n",xlab="",ylab="",frame=F)
axis(1,at=c(0,a/2,b/2),labels=c(0,"a/2","b/2"))
axis(2,at=c(0,pi/4,pi/2),labels=c(0,expression(phi),expression(frac(pi,2))),las=2)
text(a/4,pi/8,expression('N'['hits']*'=1'))
text(a/2,pi/3.5,expression('N'['hits']*'=0'))

Método de Monte Carlo

De forma geral, o método de Monte Carlo é um método estatístico, experimental de calcular integrais usando posições aleatórias, chamadas amostras, cuja distribuição tem que ser escolhidas com muito cuidado. Aleatória, ou “Random” vem de uma palavra francesa antiga randon (to run around, andar em volta) e “sample”, amostras é derivada do latin exemplum, exemplo, um pouco de cultura inútil para alegrar o dia.

Eu as vezes estava dando uma olhada no livro Statistical Mechanics. Algorithms and Computation, mais para olhar o primeiro capitulo, que tem uma introdução muito legal sobre MCMC.

Primeiro ele fala do exemplo que já vimos a alguma tempo atrás aqui.
Mas ele explica isso como um jogo de crianças.

pebbles

Então o joguinho tem um circulo desenhado dentro de um quadrado, e os guris vão jogando pedras, e vão anotando quantas pedras caem dentro do circulo. Ai basicamente ele traz o seguinte algorítimo.

mcalg.

Podemos implementar ele no R, para testar.

1
2
3
4
5
6
7
8
9
10
11
direct_pi<- function(n=4000) {
    n_acertos<-0
    for(i in 1:n){
        x<-runif(n=1,min=-1,max=1)
        y<-runif(n=1,min=-1,max=1)
        if(x^2 + y^2 < 1){
            n_acertos<-n_acertos+1
        }
    }
    return(n_acertos)
}

E fazer alguns experimentos.

> direct_pi() [1] 3152

Se usarmos nossa função, temos o número de acertos, relativos a 4 mil tentativas, 4 mil jogadas de pedras no chão. E se repetirmos vezes o suficiente, que podemos fazer com a função replicate do R, que replica coisas.

> repeticoes<-replicate(100,direct_pi()) > mean(4*(repeticoes/4000)) [1] 3.14076 > pi pi>min(4*(repeticoes/4000)) [1] TRUE

Vemos que fazendo isso 100 vezes, a média é muito próxima do valor de pi real, também que pi está contido dentro do menor e maior valores encontrados.

Massa, mas até aqui estamos falando de Monte Carlo, agora entra o Markov Chain, que é a tal da cadeia. Aqui, ele propõe outro jogo, pensando que os adultos foram jogar, só que eles foram jogar num heliporto. A primeira grande diferença aqui é que agora não da para jogar uma pedra, e cobrir toda a área, porque agora o espaço é muito grande. Mas o que da para fazer é ir andando pela área. Então a gente joga uma pedra, vai onde ela caiu, e joga outra pedra, de olhos fechados e assim vai, devagarinho, cobrindo toda a área do heliporto.

mcmcheli

Agora, o acerto é estar dentro do circulo, então todo mundo vai com um saquinho cheio de pedras, e toda lugar onde joga sua pedra, deixa uma pedrinha onde ela caiu para contar depois.
Mas existe um problema, se estamos jogando a pedra, indo onde ela está e jogando denovo, pode acontecer de eventualmente, chegarmos na borda do heliporto e jogar a pedra pra fora. E agora, a gente vai buscar a pedra? Continua jogando ela la fora até voltar pro heliporto? Bom, uma solução, é contar denovo, deixar 2 pedrinhas no mesmo lugar para contar depois, e jogar denovo do mesmo lugar, a sua pedra de arremessar.

Nisso a gente vai ter o seguinte

mcmcheli2

Agora a gente tem uma pilhas de pedras nas beiradas, e só pedras espalhadas pelo meio. Meio louco, mas esse joguinho funciona (alias, isso é o metropolis-hasting).
E temos um algorítimo para isso também.

markov_pi

E podemos implementar ele no R

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
markov_pi<-function(n=4000,passo=0.1){
    n_acertos<-0
    x<-1
    y<-1
    for(i in 1:n){
        delta_x<-runif(n=1,min=-passo,max=passo)
        delta_y<-runif(n=1,min=-passo,max=passo)
        if(abs(x+delta_x)<1 & abs(y+delta_y)<1){
            x<-x+delta_x
            y<-y+delta_y
        }
        if(x^2 + y^2 < 1){
            n_acertos<-n_acertos+1
        }
    }
    return(n_acertos)
}

Veja que agora a gente começa de um lugar (x=1,y=1), e vai andando no heliporto, andando na área, então, o antes, a gente jogava as pedrinhas, simplesmente sorteava números, agora a gente sorteia uma distância, para adicionar ao passo, que são os deltas.
Outra coisa, veja que se a gente olhar agora, antes de adicionar um acerto, se a pedra está dentro do heliporto, se estiver fora, a gente não soma os deltas, ou seja o local onde a pedra está continua igual, o que é adicionado aos acertos, ou não.

Podemos testar essa nova função, e vemos o seguinte.

> markov_pi() [1] 3397

Temos a contagem de passos.

> repeticoes<-replicate(100,markov_pi()) > mean(4*(repeticoes/4000)) [1] 3.11109

E o na média, se realizamos esse procedimento várias vezes, o valor é próximo de pi, sendo que o valor real está contido entre o maior e menor caso.

> pi pi>min(4*(repeticoes/4000)) [1] TRUE

Mas diferente do Monte Carlo, veja que aqui, nos não pegamos as amostras diretamente, como fizemos antes, aqui uma amostra depende da anterior, depende de onde você está no espaço amostral, e isso pode gerar uma certa dependência a curto prazo (que é que a gente tenta resolver quanto usa o thin no openbugs), mas ainda sim funciona.
Outra coisa é que veja que dependendo do tamanho do passo, a gente pode não explorar totalmente o espaço.

Se aumentarmos, ainda sim conseguimos andar bem em tudo

> repeticoes<-replicate(100,markov_pi(passo=0.5)) > mean(4*(repeticoes/4000)) [1] 3.13863

Mas se os passos forem muito pequenos, simplesmente podemos nunca chegar na borda, simplesmente demorar muito mais para percorrer todo o espaço.

> repeticoes<-replicate(100,markov_pi(passo=0.01)) > mean(4*(repeticoes/4000)) [1] 0.83486

Isso implica que para o MCMC, quando não retiramos nossas amostras diretamente da distribuição de interesse, temos que nos preocupar com esse ajuste do algorítimo, uma regra consensual, é que a rejeição, ou aceitação de novas amostras deve ficar em torno de 50%, o que é possível ser avaliado, é só contar o quanto estamos rejeitando amostras, para dar um parecer.

Bem é isso ai, so queria compartilhar essa explicação sobre MCMC, que achei muito legal, eu vi o cara do livro num curso do Coursera na verdade, aqui e depois eu consegui uma cópia do livro para dar uma olhada, o que é bem legal também. Alias o livro tem uma excelente página wiki, de acompanhamento, com os algorítimos em tudo que é linguagem, menos no R hehe, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:

Werner Krauth 2006 Statistical Mechanics – Algorithms and Computations

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
#Monte Carlo
direct_pi<- function(n=4000) {
    n_acertos<-0
    for(i in 1:n){
        x<-runif(n=1,min=-1,max=1)
        y<-runif(n=1,min=-1,max=1)
        if(x^2 + y^2 < 1){
            n_acertos<-n_acertos+1
        }
    }
    return(n_acertos)
}
 
#Testes
direct_pi()
 
repeticoes<-replicate(100,direct_pi())
mean(4*(repeticoes/4000))
pi<max(4*(repeticoes/4000))
pi>min(4*(repeticoes/4000))
 
#MCMC
markov_pi<-function(n=4000,passo=0.1){
    n_acertos<-0
    x<-1
    y<-1
    for(i in 1:n){
        delta_x<-runif(n=1,min=-passo,max=passo)
        delta_y<-runif(n=1,min=-passo,max=passo)
        if(abs(x+delta_x)<1 & abs(y+delta_y)<1){
            x<-x+delta_x
            y<-y+delta_y
        }
        if(x^2 + y^2 < 1){
            n_acertos<-n_acertos+1
        }
    }
    return(n_acertos)
}
 
 
#Teste
markov_pi()
 
repeticoes<-replicate(100,markov_pi())
mean(4*(repeticoes/4000))
pi<max(4*(repeticoes/4000))
pi>min(4*(repeticoes/4000))
 
 
repeticoes<-replicate(100,markov_pi(passo=0.5))
mean(4*(repeticoes/4000))
 
repeticoes<-replicate(100,markov_pi(passo=0.01))
mean(4*(repeticoes/4000))
pi<max(4*(repeticoes/4000))
pi>min(4*(repeticoes/4000))