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

Generalized additive model, quando não sabemos que modelo usar.

Vamos olhar algo que um amigo meu me perguntou sempre sobre, a algum tempo atrás, o GAM, ou Generalized additive model.

Primeiro, vamos ver qual o interesse em outro método estatístico. Sabemos que existem muitos modelos matemáticos, que explicam partes determinísticas de processos biológicos, como crescimento populacional com o modelo do Lotka-Volterra ou o comportamento de populações no espaço com o modelo de metapopulações do Levins.

Normalmente a gente não coleta dados e ajusta um modelo de crescimento populacional neles. A gente usa coisas como análise de variância e regressão linear entrou outros modelos, porque são modelos mais simples, mais robustos que acabamos por detectar padrões nos dados, e depois vamos mais afundo nesses padrões.

Vamos pensar num exemplo aqui. Vamos gerar dados de um preditor, e duas respostas diferentes.

1
2
3
4
5
###Simulando Dados
set.seed(123)
preditor <- runif(30,0,4*pi)
resposta_pot<-5+2*(preditor)^2+rnorm(30,0,10)
resposta_seno<-5+2*sin(preditor)+rnorm(30)

Seja o preditor a quantidade de açúcar no néctar e a resposta a territorialidade na fonte desse néctar. Suponha que temos algo assim:

Bem a gente pode ajustar uma reta ali, usando uma regressão linear e vai dar um resultado legal.

E teremos uma relação significativa.

Call: lm(formula = resposta_pot ~ preditor) Residuals: Min 1Q Median 3Q Max -43.001 -24.431 -2.866 23.330 53.420 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -52.67 11.82 -4.456 0.000123 *** preditor 26.04 1.47 17.718 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 28.98 on 28 degrees of freedom Multiple R-squared: 0.9181, Adjusted R-squared: 0.9152 F-statistic: 313.9 on 1 and 28 DF, p-value: < 2.2e-16

Com um valor de R quadrado bastante alto, e fim de papo? Pode ser que sim, mas podemos fazer melhor. Veja que nesse caso, temos um modelo exponencial.

Apesar da regressão resolve ali localmente, o primeiro problema é o fato de que não podemos extrapolar baseado nesse modelo de regressão, já que ele vai falhar miseravelmente além do intervalo de dados que temos, é so pensar que a curva vai subir muito mais “rápido” que a reta da regressão para imaginar isso, outro problema é que nosso residuo não tem uma distribuição normal, temos um padrão forte nele, o padrão que vem do determinismo dos dados que não está sendo coberto.

Na figura Residual vs Fitted a gente vê do que estou falando, aqui era para nao ter padrão nenhum, no entanto temos uma curva bem fácil de identificar até no olhômetro.

Mas esse é quando a gente ta com sorte, agora imagine se os dados tem um comportamento todo estranho, um padrão, mas que não é algo com cara de reta.

Ai a regressão não vai dar em nada, no caso anterior a gente poderia resolver com uma transformação, mas nesse caso não.

E aqui a robustez não funciona

Call: lm(formula = resposta_seno ~ preditor) Residuals: Min 1Q Median 3Q Max -2.8625 -1.3596 -0.2053 1.5895 3.9992 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.46404 0.77687 7.033 1.19e-07 *** preditor -0.06799 0.09658 -0.704 0.487 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.904 on 28 degrees of freedom Multiple R-squared: 0.01739, Adjusted R-squared: -0.0177 F-statistic: 0.4956 on 1 and 28 DF, p-value: 0.4873

Mas o padrão está la.

Denovo a gente vê que não deveria acreditar na regressão, a primeira figura do Residual vs Fitted, a gente vê que estamos quebrando premissas, não há uma distribuição normal nos resíduos, tem um padrão neles, mas não conseguimos pegar no modelo.

Então é pra esse tipo de coisa que existe o GAM, a gente sabe que existe um padrão determinístico nos dados, mas não sabemos qual é ele, então ou a gente sai chutando vários modelos, ou a gente uma o Generalized additive model. De certa forma é bem simples usar, ele tem uma implementação no R no pacote mgcv, que é bem simples de usar.

a sintaxe é praticamente a mesma de um glm ou lm, so adicionando aquela função s(), então o ajuste fica assim.

> modelo_gam<- gam(resposta_pot~s(preditor,k=5,fx=T),family = gaussian) > summary(modelo_gam) Family: gaussian Link function: identity Formula: resposta_pot ~ s(preditor, k = 5, fx = T) Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 134.639 1.618 83.22 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Approximate significance of smooth terms: edf Ref.df F p-value s(preditor) 4 4 908.1 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 R-sq.(adj) = 0.992 Deviance explained = 99.3% GCV = 94.22 Scale est. = 78.516 n = 30

E temos um sumario bem parecido com o da regressão linear, o modelo significativo e tal, e podemos plotar o modelo.

E ele consegue capitar bem o modelo original não? Apesar de o intercepto não ter ficado muito correto, bem, existe um padrão nos dados, nos simulamos ele assim, e ele respondeu que existe esse padrão, e deu a forma mais ou menos correta.

Agora vamos ver para o segundo caso, onde a regressão linear não era de muita serventia.

> modelo_gam<- gam(resposta_seno~s(preditor,k=5,fx=T),family = gaussian) > summary(modelo_gam) Family: gaussian Link function: identity Formula: resposta_seno ~ s(preditor, k = 5, fx = T) Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.9750 0.2667 18.66 3.46e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Approximate significance of smooth terms: edf Ref.df F p-value s(preditor) 4 4 5.861 0.00181 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 R-sq.(adj) = 0.401 Deviance explained = 48.4% GCV = 2.5598 Scale est. = 2.1332 n = 30

E de novo conseguimos detectar que algo está acontecendo, esses dados não apenas acaso puro, e vamos ver como é o modelo que ele está propondo gráficamente.

Agora sim não? Mas que magia é essa? Bem o GAM funciona a partir de regressões locais como se medíssemos como os dados mudam ao longo da nossa amostra. Algo como o modelos de Broken-Stick em regressão segmentada ou em piecewise.

Então a é algo como a gente pega uma janelinha nos dados e faz uma regressão.

Depois outra janelinha e faz outra regressão.

E a gente vai ficando com um monte dessas regressões.

Depois é só fazer uma função, tipo piecewise bem grandona e juntar todas elas.

No nosso modelo clássico de regressão, a gente tem

Y_i = \alpha + \beta \cdot X_i + \varepsilon_i

onde o \varepsilon_i=N(0,\sigma^2),

já que falamos que estamos um erro com distrubição normal dos resíduos.

A mudança que o GAM faz é

Y_i = \alpha + f (X_i) + \varepsilon_i

onde o \varepsilon_i=N(0,\sigma^2)

e o  f (X_i) é uma função que descreve os dados, e não mais um coeficiente como era na regressão.

So que você pode pensar, ei, então esse negocio vai ajustar qualquer coisa, qualquer conjunto de dados ele vai achar uma uma função que explica!

Não é tão simples assim, podemos fazer um teste, gerar dados sem padrão nenhum e ver o que ele encontra.

E para um conjunto de dados sem padrão.

Family: gaussian Link function: identity Formula: resposta_aleatoria ~ s(preditor, k = 5, fx = T) Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.2277 0.1509 -1.509 0.144 Approximate significance of smooth terms: edf Ref.df F p-value s(preditor) 4 4 1.592 0.208 R-sq.(adj) = 0.0754 Deviance explained = 20.3% GCV = 0.81982 Scale est. = 0.68318 n = 30

Ele nos diz, não tem padrão nenhum aqui uai. O modelo simplesmente cobre tudo, não existe previsibilidade nenhuma.

Agora nem tudo são flores, e temos que tomar algumas decisões, sempre existe um preço a se pagar por algo tão legal assim, que é como definir a função que descreve os dados, veja que quando estavamo vendo como é feito ela, que é definida naquela função s(), separamos os dados em quadradinhos, fazendo pequenas regressões dentro deles certo, mas qual o tamanho ideal de quadradinho a usar? Isso é o que esta sendo definido naquele argumento k dentro da função s, e esta ligado ao nossos graus de liberdade, veja que quanto menor o quadradinho, menos dados por quadradinho então é mais facil definir errado a inclinação da reta, se tiver somente um ponto num quadradinho, nem tem como definir a inclinação da reta, agora se o quadradinho for muito grande, bem a regressão linear é como se fosse so um quadrado, então nao tem porque desse método, além de que podemos peder mudanças nos dados, então isso pode ser um problema, veja que no caso dos dados oscilando na função seno, se a gente usar um k muito baixo

modelo_gam<- gam(resposta_seno~s(preditor,k=2,fx=T),family = gaussian)

Avaliando os resíduos ainda vemos que algo não está correto.

Agora se aumentarmos o k
modelo_gam<- gam(resposta_seno~s(preditor,k=5,fx=T),family = gaussian)

E podemos ir aumentando até ficar aceitavel a distribuição dos residuos, a partir do momento que estiver aceitavel, aumentar mais é disperdicio pois estamos diminuindo os graus de confiança sem obter nenhum ganho em ter um modelo com melhor ajuste.

E é claro, alguém automatizou isso, então se você coloca -1 no k, a gente usa um algoritimo que define esse valor automaticamente.

Bem é isso ai, uma pequena introdução ao GAM, veja que aquela similaridade com a ideia de regressão não é so para ser bonitinho, isso permite trabalharmos com mais de um termo, interação entre termos, variáveis aleatória, bem tudo que da pra fazer com regressão da para fazer com GAM, sendo que regressão linear pode ser considerado um caso especial de GAM, que é quando a gente cai na função da reta. Eu acho que ando muito desanimado em casa, mas acho que preciso voltar a fazer posts com análise de dados que é muito legal pra ver se animo. 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:

Zuur, AF 2012 Beginner’s Guide to Generalized Additive Models with R

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
147
148
149
150
151
152
153
154
###Simulando Dados
set.seed(123)
preditor <- runif(30,0,4*pi)
resposta_pot<-5+2*(preditor)^2+rnorm(30,0,10)
resposta_seno<-5+2*sin(preditor)+rnorm(30)
 
##
jpeg("01.jpg")
plot(resposta_pot~preditor,frame=F,pch=19,xlab="Preditor", ylab="Resposta")
dev.off()
 
##
jpeg("02.jpg")
plot(resposta_pot~preditor,frame=F,pch=19,xlab="Preditor", ylab="Resposta")
modelo_linear <- lm(resposta_pot~preditor)
abline(modelo_linear, col="red")
dev.off()
 
##
summary(modelo_linear)
 
##
jpeg("03.jpg")
plot(resposta_pot~preditor,frame=F,pch=19,xlab="Preditor", ylab="Resposta")
curve(5+2*x^2,add=T,lwd=2,col="blue")
legend("topleft",col=c("blue","red"),lwd=c(2,1),legend=c("Realidade","Modelo"),bty="n")
 
modelo_linear <- lm(resposta_pot~preditor)
abline(modelo_linear, col="red")
novos_dados <- seq(0,4*pi,0.25)
pred_interval <- predict(modelo_linear, newdata=data.frame(preditor=novos_dados), interval="prediction",level = 0.95)
lines(novos_dados, pred_interval[,2], col="orange", lty=2)
lines(novos_dados, pred_interval[,3], col="orange", lty=2)
dev.off()
 
##
jpeg("04.jpg")
par(mfrow=c(2,2))
plot(modelo_linear)
dev.off()
 
##
jpeg("05.jpg")
plot(resposta_seno~preditor,frame=F,pch=19,ylim=c(0,12),xlab="Preditor",ylab="Resposta")
curve(5+2*sin(x),add=T,lwd=2,col="blue")
legend("topleft",col=c("blue","red"),lwd=c(2,1),legend=c("Realidade","Modelo"),bty="n")
 
modelo_linear <- lm(resposta_seno~preditor)
abline(modelo_linear, col="red")
novos_dados <- seq(0,4*pi,0.25)
pred_interval <- predict(modelo_linear, newdata=data.frame(preditor=novos_dados), interval="prediction",level = 0.95)
lines(novos_dados, pred_interval[,2], col="orange", lty=2)
lines(novos_dados, pred_interval[,3], col="orange", lty=2)
dev.off()
 
##
summary(modelo_linear)
 
##
jpeg("06.jpg")
par(mfrow=c(2,2))
plot(modelo_linear)
dev.off()
 
##Pacote utilizado  para o GAM
##install.packages("mgcv")
library(mgcv)
 
## Ajuste de modelo
modelo_gam<- gam(resposta_pot~s(preditor,k=5,fx=T),family = gaussian)
summary(modelo_gam)
 
## Visualização do modelo
jpeg("07.jpg")
plot(modelo_gam)
dev.off()
 
##
modelo_gam<- gam(resposta_seno~s(preditor,k=5,fx=T),family = gaussian)
summary(modelo_gam)
 
##
jpeg("08.jpg")
plot(modelo_gam)
dev.off()
 
##
jpeg("09.jpg")
plot(resposta_seno~preditor,frame=F,pch=19,ylim=c(0,12))
 
points(x=c(4,4,6,6,4),y=c(2,4,4,2,2),type="l",lty=2)
selecao <- preditor>=4 & preditor<=6
modelo_linear <- lm(resposta_seno[selecao]~preditor[selecao])
curve(coef(modelo_linear)[1]+x*coef(modelo_linear)[2],4,6,add=T,lwd=2,col="red")
dev.off()
 
##
jpeg("10.jpg")
plot(resposta_seno~preditor,frame=F,pch=19,ylim=c(0,12))
 
points(x=c(4,4,6,6,4)+3,y=c(2,4,4,2,2)*2+1,type="l",lty=2)
selecao <- preditor>=7 & preditor<=9
modelo_linear <- lm(resposta_seno[selecao]~preditor[selecao])
curve(coef(modelo_linear)[1]+x*coef(modelo_linear)[2],7,9,add=T,lwd=2,col="red")
dev.off()
 
##
jpeg("11.jpg")
plot(resposta_seno~preditor,frame=F,pch=19,ylim=c(0,12))
 
points(x=c(4,4,6,6,4),y=c(2,4,4,2,2),type="l",lty=2)
selecao <- preditor>=4 & preditor<=6
modelo_linear <- lm(resposta_seno[selecao]~preditor[selecao])
curve(coef(modelo_linear)[1]+x*coef(modelo_linear)[2],4,6,add=T,lwd=2,col="red")
 
points(x=c(4,4,6,6,4)+3,y=c(2,4,4,2,2)*2+1,type="l",lty=2)
selecao <- preditor>=7 & preditor<=9
modelo_linear <- lm(resposta_seno[selecao]~preditor[selecao])
curve(coef(modelo_linear)[1]+x*coef(modelo_linear)[2],7,9,add=T,lwd=2,col="red")
dev.off()
 
##Testando gam contra dados aleatórios
resposta_aleatoria <- rnorm(30)
 
##
jpeg("12.jpg")
plot(resposta_aleatoria~preditor,pch=19,frame=F)
dev.off()
 
modelo_gam<- gam(resposta_aleatoria~s(preditor,k=5,fx=T),family = gaussian)
summary(modelo_gam)
 
##
jpeg("13.jpg")
plot(modelo_gam)
dev.off()
 
##Definindo um valor de K
modelo_gam<- gam(resposta_seno~s(preditor,k=2,fx=T),family = gaussian)
jpeg("14.jpg")
hist(resid(modelo_gam))
dev.off()
 
##
modelo_gam<- gam(resposta_seno~s(preditor,k=5,fx=T),family = gaussian)
jpeg("15.jpg")
hist(resid(modelo_gam))
dev.off()
 
##Auto ajuste de k
modelo_gam<- gam(resposta_seno~s(preditor,k=-1,fx=T),family = gaussian)
jpeg("16.jpg")
hist(resid(modelo_gam))
dev.off()

Criando um gráfico com bandeiras dos países.

Esses tempos, um colega me pediu para fazer alguns gráficos tipo os que tem no google analytics, mais especificamente nesse estilo aqui.

É de uma versão mais antiga dessa figura, mas a ideia é essa, basicamente um gráfico de barras na vertical. A única parte difícil, foi as bandeiras, primeiro foi achar essas bandeiras, depois inserir na figura.

Para achar as figuras, a primeira coisa foi procurar algum pacote que tivesse isso pronto, achei que encontraria, mas não achei nenhum, então ao procurar pelas bandeiras, uma fonte seria baixar diretamente do wikipedia, que tem bandeiras para todos os paises, mas eu achei esse conjunto aqui, que inclusive está em um repositorio do github que podemos clonar, o que facilita bastante o trabalho, ali temos as bandeiras em formato png em vários tamanhos, eu preferi usar o de maior resolução.

então clonado o repositório com as bandeiras, eu fiz a figura com o seguinte código:

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
library(png)
paises<-list.files("64")
 
##Figura country.png
 
country<-c("Brazil","United States","Portugal","Mozambique","Mexico","Germany","Peru","Australia","Cape Verde","Angola")
country_br<-c("Brasil","Estados Unidos","Portugal","Mozambique","México","Alemanha","Perú","Austrália","Cabo Verde","Angola")
sessions<-c(2113,259,85,63,54,49,38,37,36,36)
sessions_percentage<-c(60.39,7.4,2.43,1.8,1.54,1.4,1.09,1.06,1.03,1.03)
 
ampliacao <- 2
jpeg("country.jpg",width = ampliacao*480,height = ampliacao*480,pointsize = ampliacao*12)
escala<-barplot(rev(sessions_percentage),horiz=T,xlim=c(-150,100),xaxt="n",col="#76A7FA")
text(rep(0,10),escala,rev(sessions),pos=2)
text(rev(sessions_percentage)+14,escala,paste(rev(format(round(sessions_percentage, 2), nsmall = 2)),"%",sep=""))
text(rep(-130,10),escala,10:1,pos=2)
text(rep(-120,10),escala,rev(country_br),pos=4)
mtext(c("País","Sessões","% Sessões"),3,0,at=c(-125,-15,35))
 
list.files()
 
for(i in 1:10){
    country[i]
    posicao<-agrep(country[i],paises)
    paises[posicao]
    bandeira<-readPNG(paste("64/",paises[posicao],sep=""))
    rasterImage(bandeira,-133,rev(escala)[i]-0.5,-117,rev(escala)[i]+0.5)
}
dev.off()

Que gera a seguinte figura.

Bem legal eu achei, agora vou fazer alguns comentários do meu código.

Primeiro:

1
2
country<-c("Brazil","United States","Portugal","Mozambique","Mexico","Germany","Peru","Australia","Cape Verde","Angola")
country_br<-c("Brasil","Estados Unidos","Portugal","Mozambique","México","Alemanha","Perú","Austrália","Cabo Verde","Angola")

Bem como são poucos dados, eu digitei o nomes do países em inglês e em português, isso porque o nome em inglês eu uso para achar o nome do arquivo da bandeira, que é o nome do pais, e em português para usar na figura em si.

Depois disso para facilitar a mudança do tamanho, eu costumo guardar uma variável para o tamanho, que multiplica a largura, a altura e o tamanho do ponto, assim só altero ela para mudar o tamanho da figura, a sua resolução.

1
2
3
4
5
6
ampliacao <- 2
jpeg("country.jpg",width = ampliacao*480,height = ampliacao*480,pointsize = ampliacao*12)
 
##Figura
 
dev.off()

A figura em si é um barplot na vertical, a única coisa é que mudo bastante as margens, o tamanho do eixos de forma a caber as outras informações e bandeiras, que são colocadas com a função text para dentro da figura e mtext para a parte de fora. Depois disso, talvez essa parte ainda tenha algo de interesse.

1
2
3
4
5
6
7
for(i in 1:10){
    country[i]
    posicao<-agrep(country[i],paises)
    paises[posicao]
    bandeira<-readPNG(paste("64/",paises[posicao],sep=""))
    rasterImage(bandeira,-133,rev(escala)[i]-0.5,-117,rev(escala)[i]+0.5)
}

Então eu pego o nome do pais em inglês, e uso agrep na lista de nomes de arquivos, agrep, porque é um grep, ou seja acha a palavra na lista, mas permite alguns erros de caracteres, que da para regular com argumentos, mas então com o agrep a gente vai achar o nome do arquivo, com ele a gente le a figura png com o readPNG, que le o raster da figura, então com rasterImage a gente adiciona ela a figura, aqui a gente tem que colocar os quatro pontos da imagem para adicionar ela, veja se da até para colocar a imagem espelhada, feito isso ta pronto. Para economizar, eu fiz um loop para todas as bandeiras.

Bem é isso ai, 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. Um post bem simples, mas é o tipo de código que é legal guardar, pois pode ser útil copiar e colar um dia.

Poisson Point Process

A gente ja falou aqui sobre a distribuição dos pontos nos espaço, mas existe mais coisa a se pensar.

Eu comecei a ler o novo livro do Marc Kery, Applied hierarchical modeling in Ecology, eu ja gostava muito das coisas que ele escrevia desde os dois primeiros livros dele, Introduction to WinBUGS for Ecologists e Bayesian population analysis using WinBUGS, bem se o livro do Ben Bolker Ecological Models and Data in R abria os olhos para um monte de coisa, o Marc Kery detalhava bem modelos bayesianos desde de um começo simples, e me parece que com esse novo livro, vamos um pouquinho mais para a frente, ao mesmo estilo.

Bem uma coisa que comum no livro dele, que fala principalmente de modelos populacionais relacionados a abundância e ocupação, é o processo de pontos de poisson, ou Poisson Point Process. A gente ja falou da distribuição de poisson, que aliás é comum a muitas disciplinas, não so a ecologia. Mas vamos pensar um pouco mais sobre ela.

O Poisson point process (me perdoem usar o termo em inglês, mas as vezes é mais fácil ja ver o termo em inglês, tanto porque é mais fácil para procurar no google, achar na literatura, como muitas vezes eu traduzo as coisas erradas com meu inglês nórdico e depois vira uma confusão) tem a propriedade que cada ponto é estocasticamente independente de todos os outros pontos no processo. Apesar disso, ele é amplamente usado para modelar fenomemos representados por pontos, como fazemos direto ao usar o glm com distribuição de poisson, mas isso implica que ele nao descreve adequadamente fenomenos que tem interação suficientemente forte entre os pontos, o que é bem comum em biologia, e que explica porque a gente ve coisas como modelos quasi-poisson no glm e porque nos preocupamos com overdispersion.

O processo tem seu nome do matemático Siméon Denis Poisson e vem do fato que se uma coleção de pontos aleatórios no espaço formam um processo de poisson, então o número de pontos em uma região finita é uma variável aleatória com uma distribuição de poisson. O processo depende de um único parâmetro, que, dependendo do contexto, pode ser constante ou uma função localmente integrável (na maioria dos nossos casos a+bx ou uma função com várias médias, como numa anova). Esse parâmetro representa o valor esperado da densidade de pontos em um processo de Poisson, também conhecido como intensidade e normalmente representado pela letra grega lambda (\lambda).

A função massa de probabilidade é dada por:

p(n)=\frac{\lambda^n e^{-\lambda}}{n!}

Bem, a integral de p(n) da 1, como toda função massa de probabilidade (isso parece obvio hoje, mas eu demorei so me liguei disso depois de uma disciplina no coursera).

Para um \lambda grande, teremos o equivalente a uma distribuição normal.

Lembrando, a distribuição normal é:

p(x)=\frac{1}{\sqrt{2\sigma^2\pi}}e^{\frac{(x-\mu)^2}{2\sigma^2}}

Tudo isso basta dar ?rnorm e ?rpois, que você vai ver que ambas as distribuições são implementadas assim no R. Mas agora olha so que mágico.

Se você olhar uma parada chamada aproximação de stirling, que é uma formula fechada que mais ou menos da o resultado de um número fatorial para valores grandes, nos temos que

x! \approx \sqrt{2\pi x}e^{-x}x^x \ quando \ x \rightarrow \infty

Então, poisson é uma distribuição discreta e a normal continua, mas seja x=n=\lambda(1+\delta) onde \lambda \geq 10 (essa relação so vale para números grandes) e \delta \leq 1 (poisson é discreto, normal é contínuo). Dai substituindo temos…

p(x)=\frac{\lambda^{\lambda(1+\delta)} e^{-\lambda}}{\sqrt{2\pi \lambda(1+\delta)}e^{-\lambda(1+\delta)}(\lambda(1+\delta))^{\lambda(1+\delta)}}

Primeiro trazemos a potência do e do divisor para cima

p(x)=\frac{ e^{\lambda \delta}}{\sqrt{2\pi \lambda(1+\delta)}(\lambda(1+\delta))^{\lambda(1+\delta) - \lambda(1+\delta)}}

Depois descemos a potência do lambda la de cima

p(x)=\frac{ e^{\lambda \delta}}{\sqrt{2\pi \lambda(1+\delta)}(\lambda(1+\delta))}

simplificamos

p(x)=\frac{ e^{\lambda \delta}}{\sqrt{2\pi \lambda}(\lambda(1+\delta)) (1+\delta)^{1/2}}

ai tiramos o delta de dentro da raiz

p(x)=\frac{ e^{\lambda \delta} \lambda(1+\delta)^{-1/2}}{\sqrt{2\pi \lambda}(1+\delta)}

subimos ele

p(x)=\frac{ e^{\lambda \delta} \lambda^{-1/2}}{\sqrt{2\pi \lambda}}

cortamos o um mais delta

p(x)=\frac{ \lambda e^{-(\lambda \delta)/2}}{\sqrt{2\pi \lambda}}

e arrumamos as potências

E agora tem uma parte que eu não entendo perfeitamente, mas lembre-se que a média e a variância da distribuição de poisson são iguais, então \lambda=\mu=\sigma^2.
Assim, a parte debaixo já está pronta, e a parte de cima da equação, já é praticamente a pdf da distribuição normal, só não sei exatamente o que fazer agora, mas como vemos a relação dela com a distribuição normal nesse diagrama do blog do John Cook, é algo por ae, se alguém quiser comentar o que fazer daqui para frente, eu ficaria feliz 🙂

distribution_chart

Tudo bem que eu não sou muito bom em matemática, mas se a gente fizer um gráfico das pdfs, assim:

1
2
3
plot(1:30,dpois(1:30,10),type="l",lty=3,lwd=3,col="red",frame=F,xlab="x",ylab="Probabilidade")
points(1:30,dnorm(1:30,10,sqrt(10)),type="l",lty=3,lwd=3,col="blue")
legend("topright",legend=c("Poisson","Normal"),lty=3,lwd=3,col=c("red","blue"),bty="n")

A gente vai ver o seguinte:

01

Bem é isso ai, sem script hoje, mas olhe o repositório recologia de qualquer forma, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail, que para esse post em especial, eu tenho certeza que devo ter falado algo errado!!! Então correções são bem vindas.

Referência:

John D. Cook – Diagram of distribution relationships
Marc Kéry & Andy Royle 2015 Applied hierarchical modeling in Ecology Elsevier 183pp

Mapas aleatórios e o número de componentes

Uma forma de entender o ambiente melhor, é tentar ver se apesar de construirmos mapas aleatórios, ainda sim conseguimos enxergar algum padrão nessa aleatoriedade, como o número e tamanho de manchas formados

Um jeito de produzir um mapa aleatório, é simplesmente percorrer uma matriz, e sortear preencher cada célula com 0 ou 1 com uma certa probabilidade.

Podemos fazer uma função simples de ler da seguinte forma:

1
2
3
4
5
6
7
8
9
10
11
gera_mapa<-function(n,p) {
    matriz<-matrix(0,n,n)
    for(i in 1:nrow(matriz)){
        for(j in 1:ncol(matriz)){
            if(runif(1)<=p){
                matriz[i,j]<-1
            }
        }
    }
    return(matriz)
}

Então, pensando em construir mapas quadrados (menos número de linhas e colunas), nos geramos uma matriz de zeros de um tamanho n, então para cada coluna de cada linha (i,j) nos sorteamos um número entre 0 e 1, com distribuição uniforme, e se for menor que o nosso threshold, nos trocamos o zero por 1, senão passamos pro próximo.

Veja que podemos fazer uma imagem com essa matriz, onde quadradinhos pretos são 1 e brancos são 0.

01

Nesse caso, usamos um threshold de 0.5, mas se pegarmos um valor menor, tipo 0.2, temos menos quadradinhos pretos.

02

Agora, agora veja que isso é uma imagem dessa matriz:

> mapa [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 0 0 0 0 0 0 0 0 0 1 [2,] 0 0 1 0 0 1 0 0 0 0 [3,] 0 0 0 0 0 0 1 1 1 0 [4,] 0 0 0 0 0 0 0 0 0 0 [5,] 0 0 1 1 0 0 0 1 0 0 [6,] 0 0 0 0 1 0 0 0 0 0 [7,] 0 0 0 0 0 0 0 0 0 0 [8,] 0 1 0 0 0 0 0 0 0 0 [9,] 0 0 0 0 0 0 0 0 0 0 [10,] 0 0 0 0 0 0 0 0 1 0

Veja que da para ver a semelhança dos pixels pretos com os 1, em relação a localização. Até ai tudo bem. Agora o que seria uma mancha num mapa? Uma mancha seria a união de vários 1 adjacentes.

Mais ou menos assim, pegando um pixel qualquer, se houver algum pixel nas 8 posições em volta desse pixel escolhido, esses 2 pixels são u conjunto, são uma mancha.

Entao temos um pixel e 8 pixels em volta dele

1 2 3 4 x 5 6 7 8

Se em alguma posição de 1 a 8 tem outro 1, unimos os dois, como? Ue, só pensar em um pixel como um vértice e a união então é uma aresta. E assim fazemos nosso grafo.

03

Veja que nessa figura, eu dei nome para os vértices da seguinte forma:

(linha,coluna)

Ou seja, o número de vértices é igual ao número de 1 da matriz, agora, tendo o grafo é facil, basta a gente ver os componentes do grafo. Que no pacote igraph ja tem uma função implementada para fazer isso, mas poderíamos usar o algorítimo union-find como descrito no livro do Robert Sedgewick

Usando a função do igraph, teremos:

> components(grafo) $membership (8,2) (2,3) (5,3) (5,4) (6,5) (2,6) (3,7) (3,8) (5,8) (3,9) (10,9) (1,10) 1 2 3 3 3 4 4 4 5 4 6 7 $csize [1] 1 1 3 4 1 1 1 $no [1] 7

E veja que aqui conseguimos ver tanto o número de componentes formados, como o tamanho de cada componente em vértices, que é a quantidade de pixels do elemento, ou seja, o tamanho da mancha.

Certo agora estamos olhando isso para um mapinha que geramos ao acaso, e se olharmos para muitos em diferentes probabilidades de geração. Bem podemos fazer isso, com 10 repetições para cada probabilidade de geração e em média veremos isso:

04

Primeiro, podemos ver que conforme aumentamos o threshold de geração do mapa, temos menos componente, a não ser quando ele é muito baixo. Mas é so pensar, se for 0, não teremos nenhum componente, e se for 1 teremos exatamente 1 componente, cobrindo todo o mapa. Agora o meio termo é onde as coisas acontecem, rapidamente a gente começa a ter muitos componentes, até um pico la pelo 0.2 e 0.3, variável dependendo de quando fazemos a simulação, mas sempre para por aqui. Depois disso, a quantidade de componentes so diminui.

Agora, avaliando em relação ao mesmo valor de threshold, não é difícil imaginar que, quanto maior o threshold, mais arestas aparecem, mais os pixels acabam unidos, então temos manchas, em média, cada vez maiores. Mas isso demora um tempo para acontecer, tanto em probabilidades de 0.1 como 0.2, temos manchas pequenas, isso porque, a relação não é exatamente linear, so a partir de um momento, que começamos a ter um aumento de probabilidade resultando num aumento do tamanho das manchas.

Bem é isso ai, apesar dos mapas serem gerados aleatoriamente, ainda sim temos algumas certezas, o que permite alguma previsibilidade, mesmo 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:
Robert Sedgewick Kevin Wayne 2011 Algorithms 4th ed Addison-Wesley Professional 992 pp

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
#####################################
## Conectividade
#####################################
library(igraph)
 
##Função para gerar mapas quadrados aleatorios, n é o tamanho e p é a probabilidade de threshold
gera_mapa<-function(n,p) {
    matriz<-matrix(0,n,n)
    for(i in 1:nrow(matriz)){
        for(j in 1:ncol(matriz)){
            if(runif(1)<=p){
                matriz[i,j]<-1
            }
        }
    }
    return(matriz)
}
 
##Função para transformar um mapa em um grafo, com pixels com vértices e pontos adjacentes ligados por arestas
matriz_para_grafo<-function(matriz){
    aresta<-1
    grafo<-data.frame(origem=vector(),destino=vector())
    for(i in 1:nrow(matriz)){
        for(j in 1:ncol(matriz)){        
            if(matriz[i,j]!=0){            
                for(m in -1:1){
                    for(n in -1:1){
                        if(i+m>0 & i+m<=nrow(matriz) &
                           j+n>0 & j+n<=ncol(matriz)){
                            if(matriz[i+m,j+n]!=0){
                                grafo[aresta,]<-c(paste("(",i,",",j,")",sep=""),paste("(",i+m,",",j+n,")",sep=""))
                                aresta<-aresta+1
                            }
                        }
                    }
                }
            }
        }
    }
    indice<-!grafo[,1]==grafo[,2]
    indices<-which(mapa!=0,arr.ind = T)    
    grafo<-simplify(graph.data.frame(grafo[indice,],directed=F,vertices=paste("(",indices[,"row"],",",indices[,"col"],")",sep="")))  
    return(grafo)
}
 
##Função para fazer uma imagem da matriz como vemos ela na tela
plota_mapa<-function(mapa){
    image(t(mapa)[,ncol(mapa):1],col=0:1)
    }
 
 
#####################################
## Mapas aleatorios
#####################################
set.seed(123)
 
##Gera um mapa com threshold de 0.5
mapa<-gera_mapa(10,0.5)
mapa
#jpeg("01.jpg")
plota_mapa(mapa)
#dev.off()
 
##Threshold de 0.2
mapa<-gera_mapa(10,0.2)
#jpeg("02.jpg")
plota_mapa(mapa)
#dev.off()
 
##Transformando o mapa em grafo
grafo<-matriz_para_grafo(mapa)
 
#jpeg("03.jpg")
plot(grafo,vertex.size=20)
#dev.off()
 
##Usando a função do igraph
components(grafo)
 
##Realizando uma simulção com 10 repetições para probabilidades de 0.1 a 0.9
repeticoes<-10
probs<-seq(0.1,0.9,0.05)
tamanho<-15
 
##Matrizes para receber os resultados da simulação
n_compomentes<-matrix(0,ncol=repeticoes,nrow=length(probs),dimnames = list(paste("p=",probs,sep=""),paste("Amostra",1:10)))
maior_componente<-matrix(0,ncol=repeticoes,nrow=length(probs),dimnames = list(paste("p=",probs,sep=""),paste("Amostra",1:10)))
 
##Realizando a simulação
for(i in 1:length(probs)){
    for(j in 1:repeticoes){
        mapa<-gera_mapa(tamanho,probs[i])
        grafo<-matriz_para_grafo(mapa)
        n_compomentes[i,j]<-components(grafo)$no
        maior_componente[i,j]<-max(components(grafo)$csize)
        print(paste("Concluido i=",i," j=",j,sep=""))
    }
}
 
##Avaliando o resultado
#jpeg("04.jpg")
par(mfrow=c(2,1))
plot(probs,rowMeans(n_compomentes),type="b",xlim=c(0,1),ylab="Número médio de componentes",xlab="Threshold geração mapa",frame=F)
plot(probs,rowMeans(maior_componente),type="b",xlim=c(0,1),ylab="Tamanho do maior componente",xlab="Threshold geração mapa",frame=F)
#dev.off()