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

Teoria neutra na evolução e na ecologia

A teoria neutra, na evolução, começou em 1960 com Motto Kimura, basicamente diz que a frequência de alguns alelos podem variar por acaso devido a esses alelos terem um sucesso reprodutivo, um fitness neutro, ambos tem o mesmo valor de fitness.

Bastante simples, mas com implicações fortes. Lembrando que alelos são as formas de expressão de um gene.
Se a gente pensar na natureza e na espécie modelo mais comum em genética, a mosquinha Drosophila melanogaster que tiver, uma mutação, por exemplo, asas curtas.

08
08

Uma mosca, na natureza, sem a capacidade de voar está ferrada, ela vai morrer de fome. Ela vai ter um fitness muito menor que a mosca com asas normais e a frequência desse gene vai tender a sumir, isso é conhecido como seleção natural.

Agora talvez, outra mutações não tenham impacto direto no fitness. Estou supondo aqui, a titulo de exemplo, ja que eu não entendo de Drosophilas, mas a cor dos olhos, se não houver seleção natural qualquer, as fêmeas não verem moscas de olhos vermelhos como mais bonitas que as de olhos laranjas, o que aconteceria?

08
08

Bem essa mutação poderia ser considerada neutra, ou seja ela num melhora nem piora o fitness, mas isso não implica em nenhum momento que se você acompanhar essa população, a frequência desse alelo sera constante. Isso porque temos uma população, na natureza sempre finita, a a frequência é sempre o número de indivíduos que exibe uma característica, o alelo em questão dividido pelo tamanho total da população.

Nesse caso que temos 2 alelos em um gene (situação fictícia, eu estou assumindo que essa mutação e neutra e não existe mais variantes, sei que tem, mas é só para o exemplo ser fácil de entender), os alelos são, olhos vermelhos e olhos laranjas.

Mas esse exemplo é só para visualizar a situação, a ideia é que se esse gene é neutro, existe um processo onde se você começa cada alelo presente em 50% da população, 50 indivíduos, vamos supor, a qualquer momento que alguém usar um mata moscas aleatoriamente e mate, por cagada, 9 moscas de olhos vermelhos e 1 mosca da olhos laranjas, sobrara 41 moscas de olhos vermelhos a 49 de olhos laranjas, e se cada uma tem sei la, 1 filho, na próxima geração, veremos 45,5% de moscas de olhos vermelhos e 54.4% de moscas de olhos laranjas.

O problema é que quando eu bato o mata moscas, eu tento mata moscas, eu não fico olhando qual mosca eu vou matar, e pior que isso, nenhum tipo de moscas é mais rápido, ou mais ágil para fugir nesse exemplo, o que gera uma situação onde essa mudança na frequência dos alelos é completamente imprevisível, pra que direção ela vai.

Essa é um processo conhecido como deriva genética, ou genetic Drift, um efeito da teoria neutra.

Existem várias formas de simular a deriva genética, lembrando que a herança pode ser variável mesmo na gente. Lembre-se que mitocôndrias a gente so ganha da nossa mãe enquanto para genes nucleares alelo vem do nosso pai e outro alelo vem da nossa mãe, e tem crossing-over. Se começarmos a pensar muito nisso a gente não vai pra frente. Mas vamos usar aqui só pra visualizar o modelo do Wright–Fisher, que considera gerações discretas, indivíduos diploides, e cada nova geração tem um tamanho fixo, basicamente cada nova geração a gente pega um número aleatória de indivíduos da frequência da geração anterior.

01

Veja que cada linha é uma simulação, no eixo x temos as gerações, e no eixo x estamos representando a frequência de um alelo, lembrando que as duas frequências vão dar 1, então se um alelo tem 100% o outro tem 0, e vice versa.

Veja que em alguns casos, laranja e verde por exemplo, esse alelo domina, e o outro some, em outras simulações, ele some e só fica o outro alelo na população.

Apesar de imprevisível, a gente pode dizer que a deriva genética acaba sempre mudando as frequências de alelos, ainda mais, quase sempre sumiu com algum dos alelos, ela tirou alguém da jogada, e assim diminuiu a diversidade genética, já que ter uma população com 2 alelos é mais diversa que uma população com 1 alelo, seja ele qual for.

Certo, mas ai vendo tudo isso, um cara chamado Stephen P. Hubbell pensou se isso não poderia acontecer em outra escala, na escala de espécies. Será que existe uma teoria neutra para riqueza de espécies? Então ele pensou na teoria neutra unificada da biodiversidade e biogeografia (the unified neutral theory of biodiversity and biogeography), mas vamos dizer teoria neutra que é menor, ninguém tem saco de escrever esse nome por muito tempo.
A sua importância é que ela possibilita mais que muitos modelos de comunidade, ela faz predições claras em muitos níveis de organização, tanto de processos evolucionários como ecológicos.

Comunidades e genes são computacionalmente e matematicamente relacionadas, geralmente usando modelos iguais aos usados na genética, logo as lições que aprendemos sobre a deriva genética se aplicam a modelos de comunidade. A dinâmica ecológica numa escala local é conhecida de como deriva ecológica, e populações mudam ao acaso, devido a processos estocásticos como o nosso mata-moscas.

Hubbell propôs a teoria neutra como um modelo nulo para a dinâmica de indivíduos, para explicar a abundância de árvores tropicais ao na costa rica, barro colorado. No pacote vegan, se você abrir o pacote e digitar data(BCI), você pode ver parte dos dados que o Hubbell trabalhava, as abundâncias na floresta da reserva do Barro Colorado na costa rica, alias, como curiosidade, se você digitar data(), sem nenhum argumento, o R lista todos os data-set de exemplos disponíveis, de todos os pacotes abertos.

Certo, mas agora o controverso é que o Hubbell propôs, diferente do que comumente associamos a modelos nulos (que são como as coisas acontecem ao acaso normalmente, para teste de hipoteses), que é realmente assim que as coisas acontecem na dinâmica de comunidades.

Mas isso é controverso? Porque?

Porque estamos constantemente avaliando fatores que influenciam na riqueza e diversidade de espécies, e vem um cara e fala que tudo isso num é assim, que quem define quais espécies estão num local é o acaso. Realmente uma afirmação que vai contra muitos trabalhos que estamos acostumados a ver.

Mas se a gente não for muito extremista (e começar a achar que tudo tem que ser ou 100% estocásticos ou 100% determinísticos, sem acordos), a palavra chave para a teoria neutra do Hubbell é metacomunidade, uma coleção de comunidades similares conectadas por dispersão.

A metacomunidade é povoada inteiramente por indivíduos que são funcionalmente idênticos, segundo a teoria neutra. A teoria neutra fica então uma teoria de dinâmica de indivíduos, modelando nascimentos, morte, migração e mutação desses indivíduos, sendo a riqueza de espécies uma medida derivada desses parametros. Assumimos que dentro de uma guilda, como a de árvores do final de sucessão em floresta tropical, as espécies são essencialmente neutras, todo mundo tem um fitness igual, equivalente. Isso significa que probabilidades de nascimento, morte, mutação e migração são próximas ou equivalentes para todos os indivíduos. Assim, mudanças nos tamanhos das populações ocorrem devido a “random walks”, eventos estocásticos. Como a gente viu nas frequências gênicas. A gente pode simular praticamente igual fizemos ao gene, mas vamos usar aqui o pacote untb, que já tem essa simulação pronta.

02

Mas calma, isso não implica necessariamente na ausência de competição ou outros efeitos indiretos mediados pela dependência de densidade. A competição é vista como difusa e igual entre indivíduos.
Efeitos dependentes da densidade aparecem através de formas especificas no número total de indivíduos da comunidade, ou como características dos indivíduos relacionado as probabilidades de natalidade, mortalidade e especiação.

Um paradoxo gerado na teoria neutra de Hubbell é que na ausência de migração ou mutação, a diversidade declina para zero, riqueza para um, uma monodominância. Uma caminha aleatória (random walk), devido a equivalência de fitness, vai eventualmente resultar na perda de todas as espécies, exceto uma. Entretanto a perda de diversidade em uma comunidade local é esperada de forma muito, muito devagar, e é contra-balanceada por imigração e especiação. Assim, espécies não aumentam deterministicamente quando raras, isso faz o conceito de coexistência diferente do conceito que a gente vê nos modelos de Lotka-Volterra. Coexistência no caso da teoria neutra é um efeito estocástico balanceado pelo aparecimento de novas especies localmente. Lembrando que estamos pensando dentro de uma guilda, um nicho, não estamos comparando gramas com árvores aqui, ou onças com aranhas.

Mas se todas as comunidades fazem um random walk para a monodominâncias, como a diversidade é mantida em qualquer região em particular? Dois fatores mantém espécies em qualquer comunidade local. Primeiro, imigração para a comunidade local de indivíduos da metacomunidade. Mesmo com cada comunidade sofrendo um random walk para a monodominância, cada comunidade pode vir a ser dominada por uma espécies diferente, formando o pool de espécies dessa metacomunidade, já que todas as comunidades apresentam fitness diferentes. Assim comunidades separadas são preditas de se tornarem dominadas por espécies diferentes, e essas diferenças entre espécies locais ajudam a manter a diversidade na paisagem ocupada pela metacomunidade.

Um segundo fator, numa escala de tempo muito maior, é a especiação que dentro da metacomunidade mantém a biodiversidade. Mutação e consequentemente, especiação prove uma forma de introdução de diversidade na metacomunidade. Random walks em direção a extinção em comunidades que são tao lentas que são contrabalanceadas por especiação.

E agora temos uma pequena ideia do que se trada a teoria neutra da ecologia, e como ela veio da teoria neutra da genética. As simulações foi mais pra fazer alguma figurinha que realmente explicar algo, mas post sem figurinha é chato, e ainda as figuras foram feitas no ggplot2, que faz figuras bem bonitas, mas que até hoje eu não aprendi a usar direito. o script está no repositório do github, além de aqui em baixo, e precisamos definitivamente de muitos posts ainda pra entender a teoria do Hubbell, até onde ela vai, e como estimar parâmetros em metacomunidades para confrontar dados com teoria nesse caso aqui, mas paciência é uma virtude, sempre.

Referencia:
M Henry H Stevens 2011 A primer of ecology with R. Springer

library(untb)
library(ggplot2)
library(reshape)
library(RColorBrewer)
set.seed(123)
 
###################################################
### Deriva Genética
###################################################
#Original - http://statisticalrecipes.blogspot.com.br/2012/02/simulating-genetic-drift.html
 
#Parametros da Simulação de Deriva Genetica
N = 20           # Número de individuos diploides
N.chrom = 2*N    # Número de cromossomos
N.gen = 100      # Número de gerações
N.sim = 10       # Número de simulações
p = 0.2
q = 1-p
 
# Simulation
saida <- matrix(0,nrow=N.gen,ncol=N.sim,dimnames = list(paste("Geraçao",1:N.gen),paste("S",1:N.sim)))
saida[1,] = rep(N.chrom*p,N.sim) # Inicializando o numero de alelos na primeira geração
for(j in 1:N.sim){
  for(i in 2:N.gen){
    saida[i,j]<-rbinom(1,N.chrom,prob=saida[i-1,j]/N.chrom)
    }
  }
#Transformando em frequencia
saida<-data.frame(saida/N.chrom)
 
# Modificando os dados para a estrutura que da para plotar
sim_data <- melt(saida)
ggplot(sim_data, aes(x = rep(c(1:100), N.sim), y = value, colour = variable)) +
    geom_line() +
    scale_colour_brewer(palette="Set3") +
    ggtitle("Simulação de deriva genética") +
    xlab("Geração") + ylab("Frequencia alélica") +
    ylim(0,1) +
    labs(colour = "Simulações")
#ggsave("01.jpg")
 
###################################################
### Deriva Ecologica
###################################################
saida <- untb(rep(1:10, 90), prob=0, gens=5000, D=9, keep=TRUE)
 
sim_data <- melt(data.frame(species.table(saida)))
 
ggplot(sim_data, aes(x = rep(c(1:5000), 10), y = value, colour = variable)) +
    geom_line() +
    scale_colour_brewer(palette="Set3") +
    ggtitle("Simulação de deriva populacional") +
    xlab("Geração") + ylab("Abundância da espécie") +
    labs(colour = "Simulações")
#ggsave("02.jpg")

Olhando o MCMC mais de perto…

Algumas vezes pode ser frustrante tentar usar estatística bayesiana.
Mas vamos continuar usando o Metropolis-Hasting para um modelo de regressão linear. E para tanto simulamos mais alguns dados

01

Primeiro vamos ver o seguinte, eu alterei o prior para sempre usar uma distribuição normal agora

prior <- function(param){
alfa = param[1]
beta = param[2]
sd = param[3]
alfaprior = dnorm(alfa, sd=2, log = T)
betaprior = dnorm(beta, sd=1, log = T)
sdprior = dnorm(sd, sd=5, log = T)
return(alfaprior+betaprior+sdprior)
}

E alterei um pouquinho a função que propõe um novo passo.

proposalfunction <- function(param){
return(rnorm(3,mean = param, sd= c(1,0.1,1)))
}

Aqui reside a seguinte decisão, se você quer passos grandes ou passos pequenos. Como assim? A cada novo passo, uma condição é testada, se vamos ficar no mesmo lugar ou vamos para esse novo lugar, se esse novo lugar tem um likelehood ruim, ficamos onde estamos, lembre-se que tem um if que faz essa restrição, sempre é sorteado uma chance de aceitar ou não o novo passo, mas quanto melhor posicionado estamos perto daquilo que deve ser a verdade, menor a chance de aceitar um passo para muito longo.
O ponto também são que usamos alfa e beta, os parâmetros do modelo, o anterior da cadeia, então ali nos desvios, quando estamos falando sd=número, é o quão longe podemos sair da onde estamos.

metropolis-hastings

Mas pagaremos um preço se dermos passos curtos, que é de ter que caminhar mais, meio óbvio não.

Mas vamos olhar num plano, como mudam os valores de intercepto e inclinação ao longo da cadeia. Então o modelo sempre tem um ponto que ele esta , e vamos unindo os pontos com linhas, então a linha vermelha é a união de como esta indo a cadeia. e vamos apagar os pontos mais antigos para facilitar a visualização, clareando eles com o tempo. Como é uma simulação, vamos deixar um pontinho azul que representa a “verdade”, que infelizmente nunca saberemos numa análise com dados reais, mas como aqui o intuito é entender melhor como a análise ta se comportando, temos esse privilegio de saber a verdade.

E vemos o seguinte

rw

Se alguém lembrou da figurinha do random walk não é por acaso, só que aqui a probabilidade da uma direção, e o random walk é em torno dos parâmetros reais.

Olhe como começamos de um ponto bem longe e a cadeia caminha para perto do ponto azul, que são os parâmetros reais usados na simulação, e la começa a ficar rodando. Usamos uma cadeia de 50 mil passos aqui, e fui colocando exponencialmente os passos no gráfico, cada vez mais e mais passos de uma vez.
Outra coisa é que não importa da onde começamos, a cadeia sempre vai convergir para a “verdade”. Por isso vai ser comum ao invés de usar um ponto inicial fico, começar de um ponto aleatório, já que se o ponto inicial influenciar no resultado final, algo não deve estar certo.

Mas vamos olhar como esta a distribuição desses 50 mil pontos ai nesse plano dos dois parâmetros (lembrando que a cadeia tem 3 partes, existe o desvio que não esta ai representado).

03

Fazendo curvas de nível para representar as densidades, igual curvas de altitude em mapinhas, vemos que o pico esta bem próximo da verdade.

Outra coisa que a gente pode ver é como parece que um parâmetro compensa o outro.
Veja a cadeias em relação aos parâmetros reais.

02

Veja como os valores de inclinação ficam um pouco abaixo da linha vermelha (o parâmetro real) e o intercepto faz o contrario, fica acima, mas quando um muda outro muda junto, vemos isso na figura anterior também, já que a figura parece uma bolinha achatada, mas que tem uma inclinação.

Por isso é importante olhar como ficam as cadeias no final de uma análise e não somente computar alguma estatística do tipo aceitação (acceptance), que apesar de ficar alta (alto ou baixo sempre são conceitos abstratos sem algo para comparar), quase 55%, ainda podemos ver esse pequeno desviozinho.

Realmente todo esse processo é bem difícil de fazer, e ficar mudando os algoritmos para qualquer alteraçãozinha, mas agora nos sabemos (apesar que bugs usa um esquema ligeiramente diferente, o gibbs sampler) o que o BUGS vai fazer com os dados, ele pula toda essa parte sofrida de fazer o algoritmos, e só vamos falar o likelehood e o prior e dar como queremos as cadeias (aqui só estávamos fazendo uma, mas o ideal é fazer mais de uma) e ele faz essa parte chata, e nos preocupamos em interpretar os resultados.

Bem por último, para aqueles que tem dificuldade de entender as curvas de nível, podemos montar aquele gráfico em perspectiva, onde a montanha é o empilhamento dos pontos e quanto mais alto a montanha, maior era a densidade de pontos ali.
Gráficos assim são ruins de visualizar, principalmente representar numeração em 3d é, mas fica tão legal ele girando que não quis jogar ele fora.

per

Agora podemos voltar a usar o Openbugs ou Winbugs com um pouco mais de conforto quanto ao que diabos esta acontecendo. Ou mesmo o stan 🙂

Aqui também estão o blog original onde vi a ideia de como fazer no R o Metropolis-hasting, o Theoretical Ecology, e se alguém quer ver um gibbs sampler, olhem o Darren Wilkinson’s research blog. É isso 🙂

###################################
# Metropolis-Hastings MCMC Denovo #
###################################
#dados de exemplo
set.seed(5)
 
alfa.verdadeiro <- 4
beta.verdadeiro <- 8
desv.verdadeiro <- 6
tamanho.amostra <- 30
 
# Criando valores do preditor
x <-runif(30,0,10)
# Criando valores dependentes a partir de "beta*x + alfa + N(0,sd)"
y <-alfa.verdadeiro + beta.verdadeiro*x + rnorm(
n=tamanho.amostra,mean=0,sd=desv.verdadeiro)
 
#Grafico
plot(x,y, main="Dados Simulados",xlab="Preditor",ylab="Resposta",
frame.plot=F,pch=19)
abline(lm(y~x))
 
#Likelihood
likelihood <- function(param){
alfa = param[1]
beta = param[2]
sd = param[3]
pred = alfa + beta*x
singlelikelihoods = dnorm(y, mean = pred, sd = sd, log = T)
sumll = sum(singlelikelihoods)
return(sumll)
}
 
#Prior distribution (Agora com distribuição normal para todos os priors)
prior <- function(param){
alfa = param[1]
beta = param[2]
sd = param[3]
alfaprior = dnorm(alfa, sd=2, log = T)
betaprior = dnorm(beta, sd=1, log = T)
sdprior = dnorm(sd, sd=5, log = T)
return(alfaprior+betaprior+sdprior)
}
 
posterior <- function(param){
return (likelihood(param) + prior(param))
}
 
######## Metropolis algorithm ################
#Note que aqui voce vai configurar o tamanho dos passos.
proposalfunction <- function(param){
return(rnorm(3,mean = param, sd= c(1,0.1,1)))
}
 
run_metropolis_MCMC <- function(startvalue, iterations){
chain = array(dim = c(iterations+1,3))
chain[1,] = startvalue
for (i in 1:iterations){
proposal = proposalfunction(chain[i,])
probab = exp(posterior(proposal) – posterior(chain[i,]))
if (runif(1) < probab){
chain[i+1,] = proposal
}else{
chain[i+1,] = chain[i,]
}
}
return(chain)
}
 
#vamos começar de algum lugar aleatorio agora
#mas cuidado, um start value ruim, pode não ser bom
startvalue <- runif(3,0,10)
chain <- run_metropolis_MCMC(startvalue, 50000)
burnIn <- 25000
 
#Agora a cadeia é maior, mas temos um acceptance menor
#Mas lembre-se que o acceptance vai depender do prior
#e do proposal
acceptance <- 1-mean(duplicated(chain[-(1:burnIn),]))
acceptance
 
#olhando estabilidade das cadeias
par(mfrow = c(3,1))
plot(chain[-(1:burnIn),1], type = "l",
xlab="" , main = "Cadeia dos valores de intercepto", )
abline(h=alfa.verdadeiro,lwd=2,col="red")
plot(chain[-(1:burnIn),2], type = "l",
xlab="" , main = "Cadeia dos valores da inclinação", )
abline(h=beta.verdadeiro,lwd=2,col="red")
plot(chain[-(1:burnIn),3], type = "l",
xlab="" , main = "Cadeia dos valores do desvio", )
abline(h=desv.verdadeiro,lwd=1,col="red")
 
#Grafico animado do MCMC para os parametros
#intercepto e inclinação
passos<-c(round(seq(1,223,length=200)^2),50000)
 
for(i in 1:200) {
jpeg(sprintf("RW%03d.jpg",i), width = 300, height = 300)
plot(0,0,type="n",xlab="Intercepto",ylab="Inclinação",
xlim=c(-2,12),ylim=c(0,10),main=paste("Passo",passos[i+1]))
cores<-rev(heat.colors(i, alpha =0.5))
for(j in 1:i) {
lines(chain[passos[j]:passos[j+1],1:2],type="l",col=cores[j])
}
points(alfa.verdadeiro,beta.verdadeiro,pch=19,col="blue")
dev.off()
}
 
system("convert RW*.jpg -delay 5 rw.gif")
system("rm RW*.jpg")
 
#preparando a densidade dos pontos
library(MASS)
z<-kde2d(chain[,1],chain[,2])
 
#contorno da densidade de pontos (passos)
plot(chain[,1],chain[,2],col=heat.colors(nrow(chain),alpha=0.4))
contour(z,lwd=3,add=T)
points(alfa.verdadeiro,beta.verdadeiro,pch=19,col="blue")
 
#perspectiva animada da densidade de pontos
theta<-seq(0,360,2)
for(i in 1:length(theta)) {
jpeg(sprintf("per%03d.jpg",i), width = 350, height = 350)
persp(z,phi = 15,theta=theta[i],xlab="Intercepto",
ylab="Inclinação",zlab="Densidade")
dev.off()
}
system("convert per*.jpg -delay 30 per.gif")
system("rm per*.jpg")
 
#lembre-se que trabalhar em outro diretorio apenas para os graficos
#pode evitar transtornos.

Random Walk ou a caminhada aleatória

rw

Random walk, caminhada ao acaso ou talvez caminhada aleatória é um termo recorrente em biologia, mas que eu achei bem difícil entender.

Ele aparece e reaparece dentro da biologia sempre não importa para onde vamos.

A primeira vez que ouvi menção sobre ele foi ligado a teoria unificada da biodiversidade (UNTB) que diz que dentro dos níveis tróficos, as espécies são equivalentes, e que mudanças na abundância, composição e riqueza de espécies variam ao acaso, como um random walk. Se a gente vê uma espécie ali é porque ela teve sorte e não porque é melhor adaptada. Isso pode parecer bater de frente com a evolução, mas talvez tenha uma ponta de verdade já que como a evolução esta sempre em ação, realmente quem a gente esta vendo agora é porque sobreviveu a seleção natural, ou seja pelo menos adaptado ao momento de agora a espécie esta até certo ponto, então assumir alguma equivalência nessa adaptação entre espécies me parece razoável. Mas caso a idéia seja olhar melhor para a UNTB, o R tem o pacote UNTB para melhor entender como ela funciona e suas previsões.

Outro lugar onde a palavra Random Walk aparece é na genética, na teoria neutra, que alias diferente da UNTB, não tem muita discussão sobre sua ocorrência. Então temos alguns pedaços do DNA, sequências que não influenciam na nossa sobrevivência, logo não estão sobre influencia da seleção natural, logo elas ficam mudando e acumulando mudança, caminhando aleatoriamente sobre o plano das possibilidades. Proposta pelo senhor Kimura, tem desdobramentos muito importantes, como a teoria da coalescência e os relógios moleculares.

Mas então, Random Walks são somente uma coleção de passos ao acaso, sejam esses passos as mudanças das bases em uma sequência de DNA, mudanças da abundância das espécies em uma floresta ou pecinhas andando em um joguinho de ludo.

Podemos simular isso fácil no R simplesmente sorteando valores para coordenadas. Por exemplo, pra desenhar um ponto precisamos de cordenadas no eixo x e y, para o eixo X, a cada geração ou passo a gente sorteia entres valores 1,0 ou -1, a gente pode ir para frente, para traz ou ficar parado, fazemos o mesmo para o eixo Y e temos uma caminhada em duas dimensões. Mas ao invés de eixos e esses valores poderíamos colocar pares de bases, abundância das espécies ou ainda a diferenciação entre sequências e ver o que acontece.

Vamos simular aqui 6 populações, 6 caminhadas diferentes e ver o que acontece, vamos começar todas do centro do gráfico e para todas e andar mil passos.

01

Após caminhados 1000 passos, cada caminhada acabou em um lugar. Eu coloquei cada população com uma cor, é legal como a vermelha e a preta acabaram perto, mas cada uma das outras acabaram em um extremo do gráfico.

Eu vi a ideia desse script pela primeira vez no livro do Michael Crawley e sempre achei legal poder pensar como com o passar do tempo é assim que podem estar acontecendo as mudanças no nosso DNA e nas espécies das florestas que vemos aqui e ali.

Nesse script também tem o exemplo do inicio do post, o gráfico animado e como ele foi feito.
Eu não sei exatamente o ponto que gostaria de chegar nesse post, mas é legal ver a animação do random walk.

#Random Walk
 
#criando uma lista para depositar o resultado
random.walks<-list()
 
#simulação, o j são quantas populações eu quero, aqui são 6
#dai eu ponho no item j da lista uma matrix para receber o resultado
#ai no segundo loop interno, eu simulo o random walk, aqui para mil passos
#lembrando que tem 1001 linhas na matrix, pois para todo mundo a linha 1 é 50,50
#o meio do grafico
for (j in 1:6) {
random.walks[[j]]<-matrix(NA,ncol=2,nrow=1001,dimnames=list(1:1001,c("x","y")))
random.walks[[j]][1,]<-c(50,50)
for (i in 1:1000) {
xi<-sample(c(1,0,-1),1)
yi<-sample(c(1,0,-1),1)
random.walks[[j]][i+1,"x"]<-random.walks[[j]][i,"x"]+xi
random.walks[[j]][i+1,"y"]<-random.walks[[j]][i,"y"]+yi
}
}
 
#depois faço um grafico basico
plot(0:100,0:100,type="n",xlab="",ylab="")
 
#coloco as 6 linhas
for(i in 1:6){
lines(random.walks[[i]],col=i)
}
#uma bolinha no ultimo ponto para identificarmos onde acabou a caminhada
#apos os mil passos
for(i in 1:6){
points(random.walks[[i]][1001,"x"],random.walks[[i]][1001,"y"],col=i,pch=19,cex=1.5)
}
 
#Fazendo 200 figura jpg, organizadas por numero
for(i in 1:200){
jpeg(sprintf("RW%03d.jpg",i), width = 300, height = 300)
plot(0:100,0:100,type="n",xlab="",ylab="",main=paste("Passo",i*5))
lines(random.walks[[1]][1:(i*5),])
dev.off()
}
 
#depois é so converter para um gif animado usando o comando
#convert do imagemagik aqui no linux
system("convert RW*.jpg -delay 5 rw.gif")
#e apagar as 200 imagens
system("rm RW*.jpg")
#por segurança, é bom fazer isso em um diretorio separado