Random Web Surfer, como o google classifica suas páginas

Eu fiz algumas aulas no coursera, e em computação, as que mais me chamaram atenção foram as do Robert Sedgewick, muito famoso acho que pela tese dele sobre a análise de algorítimos do Quicksort, que ele mostra que é em média mais rápido que muitos outros, e acho que o mais comum em uso hoje. Bem também depois da aula eu comecei a olhar os livros deles, que são bem legais, entre eles o Introduction to Programming in Java, que é legal principalmente pelos exercícios que ele propõem.

Um que eu achei muito legal foi o Random Web Surfer, que eu vou colocar aqui, adaptando o código para python.

Nos primórdios da internet, a gente tinha os listão de página, para conseguir acessar, veja o post do pi-hole para mais alguma informação, mas, então, como a gente tinha uma lista, chegou a hora de ordenar ela, ordem alfabética não era muito útil, talvez classificar por tipo de páginas, mas ainda sim fica a pergunta, qual seria a melhor, qual consultar primeiro? Os primeiros web search engines usaram o critérios de quanto mais links apontavam para uma página, melhor ela era, porque todo mundo apontava para ela, logo a gente sempre ia parar nele, logo ela era muito boa. Porém todo mundo começou a burlar isso, se alguém queria ficar pra cima nesse critério, era só pagar para outras páginas colocarem um link para sua página, e rapidamente você estava no topo, ou ainda você poderia criar mil páginas mequetrefes, e colocar o link que quiser nelas, e assim controlar quem está no topo da lista de uma busca da internet, e logo esse artificio começou dificultar a vida dos buscadores de página. Mas o google usou uma forma diferente para classificar páginas.

Bem como funciona a internet, temos uma conjunto de páginas, que tem coneções umas para as outras.

O primeiro passo é como representar essa estrutura, que para nossa sorte pode ser representado facilmente usando um grafo, onde as páginas são vértices e os link entre elas são arestas, um grafo direcional.

Então a entrada pode ser um arquivo onde na primeira linha a gente tem o número de páginas que vamos trabalhar, e depois todos as arestas existentes, qual página tem um link que leva para qual página.

5 0 1 1 2 1 2 1 3 1 3 1 4 2 3 3 0 4 0 4 2

Bem então temos uma representação da internet, o que fazer com ela? De forma a não cair na armadilha anterior, que é contar o quantos links chegam a você?
Uma forma é simular como uma pessoa navega, uma pessoa entra numa página, ai pode clicar num link para outra página, ou digitar o endereço de uma nova página. Mas esses eventos tendem a acontecer com uma distribuição de pareto, a maior parte do tempo, você segue links para as próximas páginas que vai visitar, mas de vez em quanto você simplismente entra numa página por vontade própria, isso numa proporção mais ou menos de 90% do tempo seguindo links e 10% do tempo entrando em novas páginas. Seguindo essa lógica, podemos fazer uma matriz de transição entre páginas, para realizar uma simulação de como o navegador se comporta, e contabilizar quais páginas ele visita mais frequentemente.

Diferente do que uso aqui, numa programação mais em um script continuo realizando as coisas passo a passo, vamos tentar usar um pouco de orientação a objetos, que nada mais é um paradigma de programação para deixar as coisas mais organizadas.

Bem vamos começar criando uma classe em python, e vamos ter um script que vai usar essa classe.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
 
#!/usr/bin/env python
 
import sys
import random
 
class Randomwebsurfer(object):
    #Aqui vai métodos(funções da minha classe)
 
#####################################
## Random Surfer
#####################################
 
if __name__ == '__main__':
    #Aqui eu faço o script, como normalmente fazemos.

Alguns detalhes é que essa linha aqui:

1
#!/usr/bin/env python

Avisa o sistema operacional, quem é o interpretador desse programa, ou seja, se a gente já avisou quem é o interpretador, ao inves de digitar python programa.py, a gente so precisa digitar ./programa.py no terminal.

Bem toda classe tem seus atributos, ou seja a informação dela, nesse caso tem que ser o grafo que representa a relação entre as páginas, que podemos representar como uma matriz, que em python vai ser uma lista de listas, e para não ficar olhando toda hora a matriz para saber quantas páginas temos, vamos anotar o numero de páginas também, como python não precisa realmente declarar os atributos da classe, eu coloquei apenas para não me perder do que tem na classe.

1
2
3
4
5
6
7
8
9
10
11
12
13
class Randomwebsurfer(object):
    n=None
    matriz=None
 
Além disso, toda classe tem que ter um construtor, uma função que vai rodar toda vez que um objeto dessa classe for instanciado, que em classes de python a gente define o __init__
 
class Randomwebsurfer(object):
    n=None
    matriz=None
 
    def __init__(self,n):
        self.n=n
        self.matriz=self.inicia_matriz(n)

Agora toda vez que eu instancio um objeto dessa classe, por exemplo

1
objeto=Randomsurfer(5)

Eu tenho que rodar esse método (para classes a gente chama funções de métodos, mas é a mesma coisa).

Mas o que eu quero fazer, eu quero que alguém me diga quantas páginas existem, para fazer uma matriz quadrada de zeros, para anotar nela os links.

1
2
3
4
5
6
7
8
9
10
11
class Randomwebsurfer(object):
    n=None
    matriz=None
 
    def __init__(self,n):
        self.n=n
        self.matriz=self.inicia_matriz(n)
 
    def inicia_matriz(self,n):
        matriz = [[0 for i in xrange(n)] for j in xrange(n)]
        return matriz

Para tanto eu fiz esse método chamado inicia_matriz, que inicia a matriz de zeros, então dentro da classe eu chamo um método que so existe nela mesma, por isso o “self.”, assim o python sabe onde procurar, nesse caso nela mesma, na classe.

Para ver o que fizemos, podemos fazer um método para imprimir essa matriz, para ver como ela é.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
class Randomwebsurfer(object):
    n=None
    matriz=None
 
    def __init__(self,n):
        self.n=n
        self.matriz=self.inicia_matriz(n)
 
    def inicia_matriz(self,n):
        matriz = [[0 for i in xrange(n)] for j in xrange(n)]
        return matriz
 
    def imprime_matriz(self):
        for linha in self.matriz:
            print ''
            for item in linha:
                print item,
        print ''

Agora temos que ter um método para adiconar links também nessa matriz, porque ela não pode ser zerada. Então fazemos mais um método.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
class Randomwebsurfer(object):
    n=None
    matriz=None
 
    def __init__(self,n):
        self.n=n
        self.matriz=self.inicia_matriz(n)
 
    def inicia_matriz(self,n):
        matriz = [[0 for i in xrange(n)] for j in xrange(n)]
        return matriz
 
    def adiciona_link(self,linha,coluna):
        self.matriz[linha][coluna]+=1
 
    def imprime_matriz(self):
        for linha in self.matriz:
            print ''
            for item in linha:
                print item,
        print ''

Assim adicionamos o método adiciona_link, que recebe o link, da onde ele vem, que é a linha e para onde ele vai, que é uma coluna, ai adicionamos a esse item 1, que é um link, uma aresta, uma conexão entre páginas.

Agora nos vamos precisar calcular a matriz de transição, para conseguir realizar a simulação, a matriz de transição vai ter que incorporar aquelas porcentagens de 10% de ir parar por acaso em qualquer página e 90% de seguir um link na página.

As porcentagens ao acaso é fácil, basta dividir esse 10% pelo número de páginas, que fica uma pouquinho de porcentagem para cada um, se temos 5 páginas, 2% para cara uma. Já os 90% temos que dividir de forma ponderada, pelo número de vezes que outra página é linkada a ela. Para isso podemos calcular o grau de saída da página, dividimos 90% por esse número, que é dividir 90% por quantos links tem nessa página, e ai multiplicamos pela linha da matriz, o valor que esta na linha, se é zero, não há chance, se tem 2 links a chance é maior que uma que tem 1 link.

Então no código fazemos assim, um método para calcular o grau:

1
2
3
4
5
6
    def calcula_grau(self):
        self.grau=[0]*self.n
 
        for i in range(self.n):
            for j in range(self.n):
                self.grau[i]+=self.matriz[i][j]

E um método para iniciar a matriz de transição, onde calculamos o grau, fazemos uma matriz de pulo ao acaso, que é dividir o 10% entre todas as páginas, uma de link que faz a conta que eu disse, divide o 90% pelo grau da página, todos os links que ela tem, e dai multiplica pela quantidade de link para cada local, feito essas duas, a gente simplesmente soma essa duas matrizes que temos a matriz de transição.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
    def inicia_matriz_transisao(self):
        self.calcula_grau()        
        pulo=[[0 for i in xrange(n)] for j in xrange(n)]
        link=[[0 for i in xrange(n)] for j in xrange(n)]
        self.transisao=[[0 for i in xrange(n)] for j in xrange(n)]
 
        for i in range(self.n):
            for j in range(self.n):
                pulo[i][j]=0.1/float(self.n)
 
        for i in range(self.n):
            for j in range(self.n):
                link[i][j]=self.matriz[i][j] * (0.9/float(self.grau[i]))
 
        for i in range(self.n):
            for j in range(self.n):
                self.transisao[i][j]=pulo[i][j]+link[i][j]

Agora veja que ela é a essência da simulação, então ela tem que ser um atributo também, além é claro, do resultado da simulação, então já vamos deixar um frequencia, para guardar o resultado da simulação.

1
2
3
4
5
6
class Randomwebsurfer(object):
    n=None
    matriz=None
    grau=None
    transisao=None
    frequencia=None

Agora vamos finalmente fazer a simulação, que vai ser assim:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
 def simulacao(self,passos):
        self.inicia_matriz_transisao()
 
        if self.frequencia==None:
            self.frequencia=[0]*self.n
 
        pagina=0
        i=0
        while i < passos:
            sorteio=random.random()
            soma=0.0
 
            for j in xrange(self.n):
                soma+=self.transisao[pagina][j]
                if sorteio < soma:
                    pagina=j
                    break
 
            self.frequencia[pagina]+=1
            i+=1

Primeiro, a gente não vai navegar para sempre, então número de passos é o número de páginas que vamos simular, precisamos desse valor para realizar a simulação.
Nesse momento o modelo da internet, as páginas já tem que estar inserido, acho que poderíamos colocar um erro aqui, caso esse não seja o caso, mas depois eu faço as exceções.

A gente usa nosso proprio método de iniciar a matriz de transisão, depois, se não temos um local para guardar resultado, nos o iniciamos, eu fiz dessa forma testando se o resultado ja existia porque imaginei que poderiamos continuar a simulação depois, mas ok.

E temos que ter um ponto de partida, que nesse caso é a página zero, que sempre será a inicial.

Certo então começamos a simulação no while, nele a gente sorteia um numero de zero a um, depois ve em que página caiu. Que é essa parte em espeficio

1
2
3
4
5
            for j in xrange(self.n):
                soma+=self.transisao[pagina][j]
                if sorteio < soma:
                    pagina=j
                    break

suponha que temos três páginas com as porcentagens [20% 30% 50%], eu sorteio um número ao acaso entre zero e um, se caiu 0.97, qual página ele é? A porcentagem pode ser dividida assim:

20% 30% 50% [----------0.2------------0.5------------------------------------1] Página 1 Página 2 Página 3

Então temos intervalos entre zero e um, cada um de tamanho correspondente a porcentagem dada para cada página, e dependendo onde o ponto cai, é qual página a pessoa foi parar.

Por isso o for vai somando as porcentagens e comparando com o ponto, até encontrar a página.

Feito isso a gente faz um update de qual página estamos e somamos 1 a frequência com que a página é visitada, ou seja, é essa simulação se comporta como uma cadeia de markov, é um random walk sem memoria das páginas de internet.

Agora que temos nossa classe arrumadinha, podemos utilizá-la da seguinte forma

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
if __name__ == '__main__':
 
    print "Lendo entrada..."
    arquivo = open(sys.argv[1], 'r')
    n=int(arquivo.readline())
    print "Numero de paginas: " +  str(n)
 
    minharede=Randomwebsurfer(n)
 
    for linha in arquivo:
        vetor=linha.strip('\n')
        vetor=vetor.split(' ')        
        while '' in vetor: vetor.remove('')
 
        for i in range(0,len(vetor),2):
            minharede.adiciona_link(int(vetor[i]),int(vetor[i+1]))
 
 
    print "Matriz de entrada:",
    minharede.imprime_matriz()
 
    print ''
    print 'Rodando simulacao...'
    minharede.simulacao(100)
 
    print ''
    print 'Resultado:'
    minharede.imprime_resultado_simulacao()

A gente faz uma método main, que é o que o python realmente vai rodar, é o nosso script. Bem a entrada vem de um arquivo, naquel formato que falamos, usandos o sys.argv para ler o nome do arquivo que a pessoa fornece ao utilizar o programa.

1
arquivo = open(sys.argv[1], 'r')

Na primeira linha temos o número de páginas

1
2
n=int(arquivo.readline())
print "Numero de paginas: " +  str(n)

A gente le e imprime ele, para ver se ta fazendo correto, lembre-se que lemos tudo como string, então temos que converter para números.
Iniciamos um objeto da nossa classe recem criada

1
minharede=Randomwebsurfer(n)

E vamos depois lendo os pares de inteiro que são os links das páginas e adicionando ao modelo com o método adiciona_link

1
2
3
4
5
6
7
for linha in arquivo:
        vetor=linha.strip('\n')
        vetor=vetor.split(' ')        
        while '' in vetor: vetor.remove('')
 
        for i in range(0,len(vetor),2):
            minharede.adiciona_link(int(vetor[i]),int(vetor[i+1]))

Com isso pronto, a gente pode imprimir a matriz que representa o modelo, pra ver se está correta

1
2
print "Matriz de entrada:",
minharede.imprime_matriz()

Realizamos a simulação, com 100 passos para experimentar

1
2
print 'Rodando simulacao...'
minharede.simulacao(100)

E imprimimos o resultado para olhar

1
2
print 'Resultado:'
minharede.imprime_resultado_simulacao()

Então utilizando nosso programinha, vemos o seguinte:

augusto@augusto-xps:~/git/recologia$ ./randomwebsurfer.py tiny.txt Lendo entrada... Numero de paginas: 5 Matriz de entrada: 0 1 0 0 0 0 0 2 2 1 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 Rodando simulacao... Resultado: 25 30 14 24 7

E pronto, temos uma pequena implementação do Random Surfer.

Acho que depois eu faço um post discutindo esse resultado, porque esse post ta ficando muito grande já, mas é legal programar um pouquinho, eu ando muito preguiçoso e a gente vai esquecendo como fazer as coisas se não treina, mas depois eu continuo esse post. Bem é isso ai, o script vai estar la no repositório recologia, vejam o exercicio orignal aqui , provavelmente com uma explicação melhor que a minha, mas eu queria apensar programar alguma coisa aqui para passar algum tempo, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:

Robert Sedgewick e Kevin Wayne 2007 An Introduction to Programming in Java: An Interdisciplinary Approach – Pearson Addison Wesley 736pp

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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
#!/usr/bin/env python
 
import sys
import random
 
class Randomwebsurfer(object):
    n=None
    matriz=None
    grau=None
    transisao=None
    frequencia=None
 
    def __init__(self,n):
        self.n=n
        self.matriz=self.inicia_matriz(n)
 
 
    def simulacao(self,passos):
        '''
        Esse método realiza a simulacao de navegação por um usuário.
 
        '''
        self.inicia_matriz_transisao()
 
        if self.frequencia==None:
            self.frequencia=[0]*self.n
 
        pagina=0
        i=0
        while i < passos:
            sorteio=random.random()
            soma=0.0
 
            for j in xrange(self.n):
                soma+=self.transisao[pagina][j]
                if sorteio < soma:
                    pagina=j
                    break
 
            self.frequencia[pagina]+=1
            i+=1
 
 
    def inicia_matriz(self,n):
        '''
        Esse método inicia a matriz que representa a rede, mas sem nenhuma conecao entre paginas,
        recebe um argumento que é o numero de paginas do sistema
 
        '''
        matriz = [[0 for i in xrange(n)] for j in xrange(n)]
        return matriz
 
    def adiciona_link(self,linha,coluna):
        '''
        Esse método adiciona uma conecao entre duas paginas, sendo o primeiro argumento a origem e o segundo o destino
 
        '''
        self.matriz[linha][coluna]+=1
 
    def calcula_grau(self):
        '''
        Esse método calcula o grau de cada pagina.
 
        '''
        self.grau=[0]*self.n
 
        for i in range(self.n):
            for j in range(self.n):
                self.grau[i]+=self.matriz[i][j]
 
    def inicia_matriz_transisao(self):
        '''
        Esse método calcula a matriz de transicao, a matriz e links precisam estar todos adicionados
 
        '''
        self.calcula_grau()        
        pulo=[[0 for i in xrange(n)] for j in xrange(n)]
        link=[[0 for i in xrange(n)] for j in xrange(n)]
        self.transisao=[[0 for i in xrange(n)] for j in xrange(n)]
 
        for i in range(self.n):
            for j in range(self.n):
                pulo[i][j]=0.1/float(self.n)
 
        for i in range(self.n):
            for j in range(self.n):
                link[i][j]=self.matriz[i][j] * (0.9/float(self.grau[i]))
 
        for i in range(self.n):
            for j in range(self.n):
                self.transisao[i][j]=pulo[i][j]+link[i][j]       
 
    def imprime_matriz(self):
        '''
        Esse método imprime a matriz de dados que presenta o sistema
 
        '''
        for linha in self.matriz:
            print ''
            for item in linha:
                print item,
                print ''
 
    def imprime_resultado_simulacao(self):
        '''
        Esse método imprime o resultado da simulação
 
        '''
        for i in self.frequencia:
            print i, 
 
 
 
 
#####################################
## Random Surfer
#####################################
 
if __name__ == '__main__':
 
    print "Lendo entrada..."
    arquivo = open(sys.argv[1], 'r')
    n=int(arquivo.readline())
    print "Numero de paginas: " +  str(n)
 
    minharede=Randomwebsurfer(n)
 
    for linha in arquivo:
        vetor=linha.strip('\n')
        vetor=vetor.split(' ')        
        while '' in vetor: vetor.remove('')
 
        for i in range(0,len(vetor),2):
            minharede.adiciona_link(int(vetor[i]),int(vetor[i+1]))
 
 
    print "Matriz de entrada:",
    minharede.imprime_matriz()
 
    print ''
    print 'Rodando simulacao...'
    minharede.simulacao(100)
 
    print ''
    print 'Resultado:'
    minharede.imprime_resultado_simulacao()

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

Snake and Ladders, usando Markov Chain para entender um jogo.

Snake and Ladders, ou Cobras e Escadas em português, é um jogo de tabuleiro originário da índia. Acho que é parecido com a maioria dos jogos de tabuleiros de hoje em dia que so se anda para frente, só mudando a temática, enquanto no Snake and Ladders, você usa uma escada para avançar casas e uma cobra te engole e, bem, você sai pela outra extremidade casa atras, no jogo da vida você pega carona para frente, ou volta no tabuleiro como punição. No wikipedia tem um tabuleiro dele nas antigas.

02

Mas qual o interesse nisso? Olha como o mundo científico é um lugar pequeno. Lembra do senhor John Conway, do jogo da vida que discutimos aqui? Esse mesmo cara escreveu um livro chamado Winning Ways for your Mathematical Plays, e nesse livro ele usa um tipo de cadeia de Markov chamado Absorbing Markov chain para modelar esse joguinho, e aqui vamos olhar esse exemplo.

Bem o jogo é simples, temos esse tabuleiro, mas na nossa simulação vamos usar um tabuleiro um pouco modificado, um com 100 casinhas. Bem para facilitar uma representação no R, vamos usar setinhas ao invés de Escadas e cobras, ou seja se você cair uma no início da seta(lado contrário a ponta) da seta, você vai para onde ela ta apontando, se você cair na ponta da seta, continue sua vida.

Então aqui está nossa representação.

01

Então todo mundo começa com o peão fora do tabuleiro, então joga um dado de 6 faces e vê onde vai parar. Acontece que podemos ver isso como um processo markoviano, em que cada posição no tabuleiro é um estado do sistema, e temos, a cada jogada, as respectivas chances de mudar de estado. So que a gente não vai para qualquer estado, a gente tem uma chance em seis de mover de estado do quadradinho atual para o próximo, uma chance em seis de andar 2 casas, 1 chance em seis de andar 3 casas, até seis.

Isso pode ser representado por uma matriz de estados, com 100 linhas e 100 colunas, uma para cada posição do jogo. 101 na verdade, poque vamos supor que o estado 0 é quando o jogo não começou, o peão não esta ainda no tabuleiro, mas ele tem que começar de algum lugar, que é a posição zero. Assim, a primeira linha vai ser a seguinte.

> M[1,] 0 1 2 3 4 5 6 7 8 9 0.0000000 0.1666667 0.1666667 0.1666667 0.0000000 0.1666667 0.1666667 0.0000000 0.0000000 0.0000000 . . . 90 91 92 93 94 95 96 97 98 99 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 100 0.0000000

Vou suprimir um pouco dos resultado para ficar mais fácil de entender.
Mas veja que essa é a primeira linha da matriz, e ela ta dizendo que uma vez que estamos na posição zero, temos uma chance em seis, 1 divido por 6 que da 0.1666667 arredondando, podemos ir para na casa 1, 2 ,3, 5 ,6. Opa, porque o 4 está com zero de chance? Porque, se você olhar no tabuleiro ali em cima, temos que quando a gente cai no 4, a gente pega a escadinha e vai parar na casa 14.

Podemos colocar esse números na nossa representação gráfica do jogo.

02

Certo, mas a primeira jogada é fácil, mas daqui pra frente que vem a parte bolada.
Se você está na casa 3, e jogar o dado de 6 faces denovo, existem 6 possibilidades ainda.

02

Mas na verdade se a gente jogar o dado duas vezes, podemos saber qual a chance de parar em qualquer lugar ainda, simplesmente multiplicando as chances de cada valor de cada um dos dois dados jogados.

01

Ja vimos isso aqui, em probabilidade condicional. So que a matriz de estados aqui vai respeitar essa propriedade também.

Então podemos fazer uma função para multiplicar linhas dessa matriz.

powermat=function(P,h) {
    Ph=P
    if(h>1) {
        for(k in 2:h) {
            Ph=Ph%*%P
        }
    }
    return(Ph)
}

Essa função multiplica linhas da nossa matriz de transição, e multiplicando as duas primeiras linhas vemos o seguinte.

03

Veja que estamos representando a intensidade das probabilidades com cores, quanto mais vermelho, maior a chance, quanto mais clarinho, menor a chance. Aqui a chance maior é de estarmos na posição sete, já que podemos chegar no sete tirando 1 e 6, 2 e 5, 3 e 4, 4 e 3, opaaaaaa, lembre-se que quando a gente cai no 4, a gente pega a escada para o 14, então com 4 e 3, vamos para o catorze, e depois andamos 3, só que no 17 a gente é engolido pela cobra e volta pro 7, eu diria que isso é uma ironia, já que era para estarmos no 7 mesmo. Mas veja que como temos a matriz, estamos levando em conta as escadas e cobras, o que ignoraríamos, ou teríamos bem mais dificuldade se formos no braços, simulando cada acontecido.

É fácil se perder, se a gente andar 10 vezes, podemos andar de 6 em 6 casas, com muita sorte, você esperaria que ninguém tenha ganho o jogo não? Afinal andando direto, isso da 60 casas, se você começa no zero, so tirando 6 estaria na casa 60, mas lembre-se que quem tem sorte, não vai tirar só 6, vai é pegar as escadas certas, e andar muito mais rápido que isso, veja a distribuição depois de 10 jogadas como fica.

04

As vezes da a impressão que tem algo de especial nessa matriz, mas confira, a gente tem que cada casa só tem chance de ir para 6 outras casas, só voltando em casos de quando encontramos a boca de uma cobra.

> sum(M[2,]>0) [1] 6 > sum(M[10,]>0) [1] 6 > sum(M[50,]>0) [1] 6

Agora vamos supor que depois dessas 10 jogadas, um cara que não sabe perder da um chute na mesa e derruba todos os peões, agora a gente não lembra onde estava, mas podemos pegar a figura acima e podemos então tentar olhar onde você estava, primeiro que sabemos exatamente onde você não pode estar, só ver os lugares com probabilidade zero.

Mais do que isso, se você afirma que estava em algum lugar entre 59, 60 ou 61, podemos ver a probabilidade de estar em cada um desses lugares.

> h=10 > posição<-(initial%*%powermat(M,h))[59:61]/sum((initial%*%powermat(M,h))[59:61]) > posição [1] 0.1597003 0.5168209 0.3234788 > h=10 > posição<-(initial%*%powermat(M,h))[59:61]/sum((initial%*%powermat(M,h))[59:61]) > sum(posição) [1] 1 > posição [1] 0.1597003 0.5168209 0.3234788

Veja que a gente multiplica onde podemos estar depois de 10 jogadas, olhas a probabilidade de cada uma dessas 3 posições, divido pela soma dessas 3 posições, e temos nossa distribuição para a sua afirmação de estar nas casas 59, 60 ou 61.

Veja que podemos seguir a mesma linha, e perguntar quando ainda estaremos jogando, se a gente olhar a distribuição para cada jogada, é só ver que a chance da gente estar jogando ainda, é a chance de não estarmos na posição 100. Veja que uma vez que caímos no estado “posição 100”, ele é estacionário, e todo o jogo tende pra ele, se a gente chegar la, a gente não sai mais, isso é o jogo acabou, então podemos olhar a chance de não estarmos no 100 jogada por jogada e construir uma curva de probabilidade para isso.

05

Após 50 jogadas, é praticamente certeza que o jogo acabou. Veja a distribuição de chances depois de 50 jogadas.

06

Tem que ser muito pé frio para não ter terminado o jogo depois de 50 jogadas, mas veja que isso é uma situação pouco provável, mas não impossível.

Mas em média, vão levar 32 jogadas para finalizar o jogo.

> sum(1-game) [1] 32.16499

Veja que 1 é 100% de chance do jogo ter terminado, quando estamos diminuindo 100% de chance, o que está sobrando é a chance de ainda estarmos jogando. Quando o jogo acabou com certeza, temos 1, ou 100% de chance do jogo ter acabado, então 100% – 100% da zero, a chance da zerada, então estamos somando ai em cima todas as chances de estarmos jogando, para a distribuição ali em cima. Assim talvez seja mais difícil pensar.

Para ficar mais simples, podemos perguntar, quando na distribuição, temos menos de 50% de chance de estar jogando, 50% de chance do jogo ter acabado.

> max(which(1-game>.5)) [1] 29

São vinte e nove jogadas, e já temos 50% de chance do jogo ter acabado.

07

Agora tudo isso que estamos falando, é tendo em mente apenas um jogador. E com dois jogadores, mudaria muito?

Bem, é intuitivo pensar que, quanto mais gente jogando, mais rápido o jogo deve acabar. Se temos dois jogadores, é mais difícil exatamente termos dois jogadores azarados.

08

É só fazer a potência da distribuição anterior, elevar ela pelo número de jogadores, no caso dois, que temos a nova distribuição de quando o jogo deve acabar, podemos calcular também qualquer uma das estatísticas anteriores.

Mais do que isso, podemos usar isso para qualquer quantidade de jogadores.

09

Bem, para quem não pegou ainda a ideia das cadeias de markov, esse foi mais um exemplinho bem legal.
Esse não é o primeiro exemplo de Markov Chain que a gente olha, ja vimos aqui como prever a chuva usando Markov Chain, ou ainda como podemos prever a chuva, sem ver a chuva, apenas pelo comportamento de alguém, nesse exemplo está aqui. Isso pode parecer todo, mas o crescimento populacional pode ser pensado dessa forma, o que a gente ve pode ser apenas uma parte do processo que realmente esta acontecendo, como vimos no crescimento de passarinhos aqui.

Outro processo que temos muito interesse e que é assim é a evolução, a gente viu como simular cadeias de DNA com cadeias de markok aqui. E a evolução é mais ou menos isso, é uma passagem quase sempre, so de ida, não tem volta, a gente dificilmente volta para um estado anterior.

Bem, acho que é isso ai por hoje. Ficamos por aqui. Os scripts sempre no repositório recologia e é isso, até o próximo post. A, e se você curtir esse post, de uma olhada num blog chamado Freakonometrics, foi la que eu vi pela primeira vez essa ideia do Snake and Ladders

#Codigo original
#http://freakonometrics.blog.free.fr/index.php?post/2011/12/20/Basic-on-Markov-Chain-(for-parents)
 
#Matrix de transição
n=100
M=matrix(0,n+1,n+1+6)
rownames(M)=0:n
colnames(M)=0:(n+6)
 
for(i in 1:6) {
    diag(M[,(i+1):(i+1+n)])=1/6
}
 
M[,n+1]=apply(M[,(n+1):(n+1+6)],1,sum)
M=M[,1:(n+1)]
starting=c(4,9,17,20,28,40,51,54,62,64,63,71,93,95,92)
ending  =c(14,31,7,38,84,59,67,34,19,60,81,91,73,75,78)
 
for(i in 1:length(starting)) {
    v=M[,starting[i]+1]
    ind=which(v>0)
    M[ind,starting[i]+1]=0
    M[ind,ending[i]+1]=M[ind,ending[i]+1]+v[ind]
}
 
#Olhando o que criamos
nrow(M)
M[1,]
M[2,]
 
sum(M[2,]>0)
sum(M[10,]>0)
sum(M[50,]>0)
 
#Função para multiplicar a matriz de transição
powermat=function(P,h) {
    Ph=P
    if(h>1) {
        for(k in 2:h) {
            Ph=Ph%*%P
        }
    }
    return(Ph)
}
 
#Função para plotar o calculo das posições
initial=c(1,rep(0,n))
COLOR=rev(heat.colors(101))
u=1:sqrt(n)
boxes=data.frame(index=1:n,ord=rep(u,each=sqrt(n)),abs=rep(c(u,rev(u)),sqrt(n)/2))
 
#
position=function(h){
    D=initial%*%powermat(M,h)
    plot(0:10,0:10,col="white",axes=FALSE,xlab="",ylab="",main=paste("Posição após",h,"jogadas"))
 
    for(i in 1:n) {
        polygon(boxes$abs[i]-c(0,0,1,1),boxes$ord[i]-c(0,1,1,0),col=COLOR[min(1+trunc(500*D[i+1]),101)],border=NA)
    }
 
    segments(c(0,10),rep(0,2),c(0,10),rep(10,2))
    segments(rep(0,2),c(0,10),rep(10,2),c(0,10))
    segments(0:10,rep(0,11),0:10,rep(10,11))
    segments(rep(0,11),0:10,rep(10,11),0:10)
 
    for(i in 1:length(starting)) {
        arrows(x0=boxes$abs[starting[i]]-0.5,
               y0=boxes$ord[starting[i]]-0.5,
               x1=boxes$abs[ending[i]]-0.5,
               y1 =boxes$ord[ending[i]]-0.5,
               lty=3,length=0.10,col="darkgray",lwd=2)
    }
    text(boxes$abs-.5,boxes$ord-.5,boxes$index,cex=.7)
}
 
#figuras
 
position(1)
position(2)
position(10)
 
#Se quiser ver uma animação olhe esse codigo, no Sys.sleep da para configurar
#uma pausa entre os plots, deixei em 1/4 de segundo por jogada.
#for (i in 1:100) {
#                  position(i)
#                  Sys.sleep(0.25)
#                  }
 
 
#Probabilidade parcial
h=10
posição<-(initial%*%powermat(M,h))[59:61]/sum((initial%*%powermat(M,h))[59:61])
sum(posição)
posição
 
#Distribuição de jogadas
distrib<-initial%*%M
game<-rep(NA,1000)
 
for(h in 1:length(game)){
    game[h]<-distrib[n+1]
    distrib<-distrib%*%M
}
 
#Figuras de distribuição
plot(1-game[1:200],type="l",lwd=2,col="red",ylab="Probabilidade de ainda estar jogando",xlab="Jogadas")
position(50)