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

Leave a Reply

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