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'))

Leave a Reply

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