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

Duas espécies em metapopulações de Levins

Ja vimos metapopulações para uma espécie aqui, aqui e aqui, mas será que poderíamos adicionar mais uma espécie ao modelo de Levins?

Primeiro, vamos lembrar como é a metapopulação de levins.

\frac{dO}{dt}=cO(1-O)-eO

onde:

  • O – Ocupação
  • c – taxa de colonização
  • e – taxa de extinção

Então basicamente temos dois termos,

cN(1-N)

que faz a população crescer a uma taxa c, mas quanto maior a ocupação, mais propagulos são emitidos, mas veja que quanto maior a ocupação, menor o crescimento, que é a multiplicação por (1-N), conforme o N cresce, o valor tende a zero, porque cada vez teremos 1 menos alguma coisa, que da menos que um, estaremos multiplicando por uma fração menor.

eN

Cada mancha ocupada tem uma chance de extinção e, logo quanto mais manchas ocupadas, mais manchas veremos sendo extintas.

Agora vamos colocar outra espécie, ela vai ocupar as manchas vagas pela primeira espécie, ou seja, vagou, ela tem chance de entrar ocupar, mas a se a espécie um ocupar novamente, ela perde a competição.

Como são duas espécies, precisamos agora acompanhar a mudança na ocupação de cada uma, então precisamos de duas equações

E espécie 1 continua como a equação anterior
\frac{dO_1}{dt}=c_1O_1(1-O_1)-eO_1
mas a segunda fica
\frac{dO_2}{dt}=c_2O_2(1-O_1-O_2)-c_1*O_1*O_2-eO_2

Ou seja, a metapopulação 2 expande sua ocupação, diminuindo a expansão de acordo com o o quanto ja está ocupado pela ocupação da espécie 1 e dela mesmo, e tiramos as manchas que ela perde na competição, além do que ela perde com a extinção.

Agora, quem lembra do pacote desolve, vamos usar ele para ver como isso fica ao longo do tempo, então temos que transformar numa função para usar com a função ode.

1
2
3
4
5
6
7
8
9
10
11
metpop2sp <- function(t, y, p) {
    {
        o1 <-y[1] 
        o2 <-y[2]
    }
    with(as.list(p), {
        do1.dt <- c1*o1*(1-o1)-e*o1
        do2.dt <- c2*o2*(1-o1-o2)-c1*o1*o2-e*o2
        return(list(c(do1.dt, do2.dt)))
    })
}

Agora temos que determinar por quanto tempo queremos observar a ocupação e os paramtros definir os parametros de ocupação e extinção

1
2
3
4
5
6
tempo <- seq(0, 200, by =5)
c1 = 0.4
c2 = 0.5 
o1 = 0.1
o2 = 0.4 
e = 0.25

E vemos o seguinte

01

A espécie um rapidamente domina o ambiente, retirando a espécie 2 da jogada.

Mesmo quando mudamos os parametros um pouco, deixando ambas com a mesma taxa de colonização, maior a taxa de extinção, lembrando que para uma espécie, quando a extinção ficava maior que a colonização, ela tendia a extinção, sendo a diferença o quão rápido isso ia acontecer.

1
2
3
4
5
c1 = 0.1 
c2 = 0.1
o1 = 0.05 
o2 = 0.05 
e = 0.05

02

Agora se mudarmos o jogo para a espécie 1 extinguir, a espécie 2 consegue dominar, salvo que ela também não tenha a colonização menor que a extinção

1
2
3
4
5
c1 = 0.05
c2 = 0.07 
o1 = 0.05 
o2 = 0.05
e = 0.06

03

E existem também os casos onde podemos ter uma coexistencia

1
2
3
4
5
c1 = 0.1
c2 = 0.5
o1 = 0.05 
o2 = 0.05
e = 0.05

04

Agora, as possibilidades de combinações de parâmetros são grandes, tornando mais complexo definir quando uma ou outra espécie ganha, ou quando temos a cooexistencia, mas podemos ter um feeling da coisa já. além disso, nossas equações garantem a vitoria da espécie 1, mas isso não precisa ser assim.

Bem é isso ai, lembrando um pouco de equações diferencias, que não falamos a um bom tempo aqui e vendo a parte essa expansão para a metapopulação, 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:
A. A. de Oliveira e P. I. Prado – Coexistência em Metapopulações – Roteiro no EcoVirtual
Robert May & Angela McLean 2007 – Theoretical Ecology: Principles and Applications. Oxford University 272 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
##
##install.packages("deSolve")
library(deSolve)
rm(list=ls())
 
##Função
metpop2sp <- function(t, y, p) {
    {
        o1 <-y[1] 
        o2 <-y[2]
    }
    with(as.list(p), {
        do1.dt <- c1*o1*(1-o1)-e*o1
        do2.dt <- c2*o2*(1-o1-o2)-c1*o1*o2-e*o2
        return(list(c(do1.dt, do2.dt)))
    })
}
 
##Tempo
tempo <- seq(0, 200, by =5)
 
##Caso 1
c1 = 0.4
c2 = 0.5 
o1 = 0.1
o2 = 0.4 
e = 0.25
 
out.ode<-ode(c(o1,o2), tempo, metpop2sp, list(c1=c1,c2=c2,e=e))
 
##jpeg("01.jpg")
matplot(out.ode[,1], out.ode[, -1], type = "l", lty = 1:2,frame=F,xlab="Tempo",ylab="N",lwd=2,col=c("red","blue"))
legend("right", c("Sp1", "Sp2"), lty = 1:2, col = c("red","blue"), bty = "n",lwd=2)
##dev.off()
 
##Caso 2
c1 = 0.1 
c2 = 0.1
o1 = 0.05 
o2 = 0.05 
e = 0.05
 
out.ode<-ode(c(o1,o2), tempo, metpop2sp, list(c1=c1,c2=c2,e=e))
##jpeg("02.jpg")
matplot(out.ode[,1], out.ode[, -1], type = "l", lty = 1:2,frame=F,xlab="Tempo",ylab="N",lwd=2,col=c("red","blue"))
legend("right", c("Sp1", "Sp2"), lty = 1:2, col = c("red","blue"), bty = "n",lwd=2)
##dev.off()
 
##Caso 3
c1 = 0.05
c2 = 0.07 
o1 = 0.05 
o2 = 0.05
e = 0.06
 
out.ode<-ode(c(o1,o2), tempo, metpop2sp, list(c1=c1,c2=c2,e=e))
##jpeg("03.jpg")
matplot(out.ode[,1], out.ode[, -1], type = "l", lty = 1:2,frame=F,xlab="Tempo",ylab="N",lwd=2,col=c("red","blue"))
legend("right", c("Sp1", "Sp2"), lty = 1:2, col = c("red","blue"), bty = "n",lwd=2)
##dev.off()
 
##Caso 4
c1 = 0.1
c2 = 0.5
o1 = 0.05 
o2 = 0.05
e = 0.05
 
 
out.ode<-ode(c(o1,o2), tempo, metpop2sp, list(c1=c1,c2=c2,e=e))
##jpeg("04.jpg")
matplot(out.ode[,1], out.ode[, -1], type = "l", lty = 1:2,frame=F,xlab="Tempo",ylab="N",lwd=2,col=c("red","blue"))
legend("right", c("Sp1", "Sp2"), lty = 1:2, col = c("red","blue"), bty = "n",lwd=2)
##dev.off()

Existem infinitos números primos!

Essa é a pergunta que foi respondida pelo Euclides, cerca de 350 antes de Cristo, e hoje, até hoje a gente (pelo menos eu) apanha para entender.

Vamos olhar

Teorema: Existem infinitos números primos.
Suponha que exista uma quantidade finita de números primos. Seja a lista desses números p_1,p_2,p_3,\dots,p_n a lista de todos os números primos e um m um número tal que m = p_1 p_2 p_3 \dots p_n + 1. Note que m não é divisivel por p_1 ja que dividir m por p_1 vai dar um quociente de p_2 p_3 \dots p_n e um resto de 1. Similarmente, m não é divisivel por p_2,p_3,\dots,p_n.

Agora considerando o fato que todo número maior que 1 ou é um número primo, ou pode ser escrito como um produto de números primos. É claro que m é maior que 1, já que ele é alguma coisa mais 1, então m é primo ou deve ser um produto de primos.

Suponha que m seja primo, note que m é maior que todos os números na lista p_1,p_2,p_3,\dots,p_n, então nos encontramos um número que é primo que não está na lista, contradizendo nossa suposição original de que tinha a lista de todos os números primos.

Agora suponha que por outro lado m não é primo, então ele tem que ser um produto de primos. Suponha que q é um dos primos desse produto. Então m tem que ser divisível por q. Mas a gente já viu que m não é divisível por nenhum número na lista p_1 , p_2,\dots , p_n, então denovo nos temos uma contradição com a suposição de que temos uma lista de todos os números primos.
Já que a suposição de que existem uma quantidade finita de números primos leva a uma contradição, tanto considerando m como primo ou não primo, que são as duas únicas possibilidades, então deve existir uma quantidade infinita de números primos.

Se a nossa lista de primos for 2 e 3, então 2 vezes 3 da 6, mais 1 da 7 que também é primo. Se a lista fosse 2,3 e 5, multiplicando da 30, mais 1 da 31, que também e primo, ou seja, a gente pode ir gerando primos assim. A gente não gera o próximo primo da lista, mas gera um número primo, ou seja, a sempre tem mais um número primo.

Bem é isso ai, esse negócio de prova é muito divertido, 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:
Daniel J. Velleman 2006 – How to Prove It: A Structured Approach 2ed Cambridge University Press 384pp.

Álgebra de Matrizes, a matriz inversa.

Operações com matrizes é uma atividade básica utilizada em muitas das operações de estatística que vemos. Como ja vimos aqui, por exemplo para determinar o resultado de um modelo linear, os valores preditos, agente pode multiplicar uma matriz pelos valores dos parâmetros, matriz essa que é construído pela função formula, o tilzinho (“~”) que a gente usa nas funções de modelagem.

Certo, mas vamos ver uma propriedade legal da multiplicação de matriz, que é resolver conjuntos de equações.

Por exemplo, vamos pegar o sistema de equações abaixo.

  a+b+c=6\\  3a-2b+c=2\\  2a+b-c=1

Esse sistema pode ser escrito e resolvido usando álgebra de matrizes.

   \left(\begin{array}{ccc} 1&1&1 \\ 3&-2&1 \\ 2&1&-1  \end{array} \right) \cdot   \left(\begin{array}{c} a \\ b \\ c  \end{array} \right) =   \left(\begin{array}{c} 6 \\ 2 \\ 1  \end{array} \right)

Se realizar a multiplicação da matriz pelo vetor, a gente tem o sistema acima.
E ele ele pode ser resolvido da seguinte forma, para encontrar os valores de a, b e c

   \left(\begin{array}{c} a \\ b \\ c  \end{array} \right) =  \left(\begin{array}{ccc} 1&1&1 \\ 3&-2&1 \\ 2&1&-1  \end{array} \right)^{-1} \cdot   \left(\begin{array}{c} 6 \\ 2 \\ 1  \end{array} \right)

Ou seja, é so multiplicar a matriz inversa pelos resultados das equações, que temos os valores de a, b e c. Agora o problema é achar a matriz inversa, essa matriz elevado a -1, normalmente denotada como X^{-1}.
A matriz inversa tem outra propriedade, que quando multiplicada uma matriz pela sua matriz inversa, temos a matriz identidade.

Ou seja:
X \cdot X^{-1} = I

Veja que apesar da complicação dessas operações, tem algoritimos prontos para resolver tudo :).
Para achar a matriz inversa no R, existe a função solve, então pegando o nosso primeiro caso, se transformamos ele em matriz

> x <- matrix(c(1,3,2,1,-2,1,1,1,-1),3,3) > y <- matrix(c(6,2,1),3,1) > x [,1] [,2] [,3] [1,] 1 1 1 [2,] 3 -2 1 [3,] 2 1 -1 > y [,1] [1,] 6 [2,] 2 [3,] 1

Primeiro podemos testar a propriedade de que X \cdot X^{-1} = I:

> solve(x)%*%x [,1] [,2] [,3] [1,] 1.000000e+00 2.775558e-17 -2.775558e-17 [2,] 1.110223e-16 1.000000e+00 -2.775558e-17 [3,] -1.110223e-16 -5.551115e-17 1.000000e+00

Claro que a solução é feito com algum método númerico, como aqui, então a aproximação pode parecer um número esquisito, mas so olhar o número de casas decimais, se a gente arredonda…

> round(solve(x)%*%x) [,1] [,2] [,3] [1,] 1 0 0 [2,] 0 1 0 [3,] 0 0 1

A solução para o sistema de equações é a multiplicação da inversa pelo vetor de respostas

> solve(x)%*%y [,1] [1,] 1 [2,] 2 [3,] 3

Veja que isso é o mesmo que usar o segundo argumento da função, e lembrando que %*% são para operaçoes matriciais no R, usar somente o * para multiplicar vai fazer uma operação element-wise, elemento por elemento

> solve(x,y) [,1] [1,] 1 [2,] 2 [3,] 3

É realmente é milagroso tudo isso. Mas ai bate uma curiosidade, como diabos isso é feito? Bem existem muitos métodos, mas vamos dar uma olhadinha em um chamado Gaussian Elimination.

Basicamente:

  1. Expandir a matriz de coeficientes com mais uma coluna, o vetor de respostas
  2. Transformar a matriz de coeficientes numa matriz triangular
  3. Resolver as linhas de baixo para cima

No R funciona assim

A primeira parte é simples, é so unir a matriz e as respostas, ai ja contamos o número de linhas

1
2
expandido <- cbind(x,y)
p <- nrow(x)

Dai para transforma em matriz triangular, é so realizar essa operação descrita no wolfram:

x_i=\frac{1}{a_{ii}'} (b_i' -\sum\limits_{j=i+1}^k a_{ij}'x_j)

Mas a partir da segunda linha é mais fácil de pegar o que está sendo feito

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
##Transformando para uma matriz triangular
expandido[1,] <- expandido[1,]/expandido[1,1]
 
##A partir da segunda linha
i <- 2
while (i <= p) {
    j <- i
    ## para todas as linhas diminua depois dela da linha de cima, multiplicando pelo elemento da coluna anterior, veja
    ##que toda vez que a gente faz isso, a gente zera uma coluna, da linha 2 para a frente, ou seja, a gente ta transformando
    ##na matriz triangular
    while (j <= p) {
        expandido[j, ] <- expandido[j, ] - expandido[i-1, ] * expandido[j, i-1]
        j <- j+1
    }
    ##depois disso a gente troca as linhas de lugar
    while (expandido[i,i] == 0) {
        expandido <- rbind(expandido[-i,],expandido[i,])
    }
    ##e multiplica por aquel 1/a da formula
    expandido[i,] <- expandido[i,]/expandido[i,i]
    i <- i+1
    ##veja que a cada iteração aqui, uma coluna é zerada, e a gente passa a operar uma submatriz abaixo, ou seja, na
    ##primeira rodada da gente opera da linha 2 pra baixo, na segunda iteração, a da linha 2 pra baixo
}

E veja como ficou

> expandido [,1] [,2] [,3] [,4] [1,] 1 1 1.0 6.0 [2,] 0 1 0.4 3.2 [3,] 0 0 1.0 3.0

Agora é so ir substituindo de baixo para cima. por exemplo, so de olho a gente sabe que a última equação da 3. Aliás, é a divisão ali que faz com que a matriz identidade apareça ali no meio, não exatamente, mas é so substituir o valor de baixo na linha de cima que a gente resolve a linha de cima. Então resolvendo de baixo pra cima…

1
2
3
4
5
6
7
for (i in p:2){
    for (j in i:2-1) {
        print(expandido)
        expandido[j, ] <- expandido[j, ] - expandido[i, ] * expandido[j, i]
 
    }
}

Temos a solução

> expandido [,1] [,2] [,3] [,4] [1,] 1 0 0 1 [2,] 0 1 0 2 [3,] 0 0 1 3

É impressionante como essas coisas funcionam, e como as pessoas pensam nessas soluções. Claro que isso so funciona para matrizes que são invertíveis, por exemplo

1
2
x <- matrix(c(1,1,1,1,1,1,1,1,1),3,3)
y <- matrix(c(1,2,3),3,1)

Essa matriz não tem como ser invertida, porque se você pensa no sistema de equações, nao da para somar a, b e c e dar 1 para uma equação e dois para outra, não existem valores de a, b e c que façam isso, logo, não da para achar uma matriz inversa.

> solve(x,y) Error in solve.default(x, y) : Lapack routine dgesv: system is exactly singular: U[2,2] = 0

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.

Referência:
PH525x series – Biomedical Data Science
Wolfram MathWorld

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
##Exemplo
x <- matrix(c(1,3,2,1,-2,1,1,1,-1),3,3)
y <- matrix(c(6,2,1),3,1)
x
y
 
 
##Matriz identidade
diag(3)
 
##Matriz inversa de x
solve(x)
 
##Matriz multiplicado pela matriz inversa
solve(x)%*%x
round(solve(x)%*%x)
 
##Solução do sistema de equações
solve(x)%*%y
solve(x,y)
 
##Expandir matriz
p <- nrow(x)
expandido <- cbind(x,y)
 
##Transformando para uma matriz triangular
expandido[1,] <- expandido[1,]/expandido[1,1]
 
##A partir da segunda linha
i <- 2
while (i <= p) {
    j <- i
    ## para todas as linhas diminua depois dela da linha de cima, multiplicando pelo elemento da coluna anterior, veja
    ##que toda vez que a gente faz isso, a gente zera uma coluna, da linha 2 para a frente, ou seja, a gente ta transformando
    ##na matriz triangular
    while (j <= p) {
        expandido[j, ] <- expandido[j, ] - expandido[i-1, ] * expandido[j, i-1]
        j <- j+1
    }
    ##depois disso a gente troca as linhas de lugar
    while (expandido[i,i] == 0) {
        expandido <- rbind(expandido[-i,],expandido[i,])
    }
    ##e multiplica por aquel 1/a da formula
    expandido[i,] <- expandido[i,]/expandido[i,i]
    i <- i+1
    ##veja que a cada iteração aqui, uma coluna é zerada, e a gente passa a operar uma submatriz abaixo, ou seja, na
    ##primeira rodada da gente opera da linha 2 pra baixo, na segunda iteração, a da linha 2 pra baixo
}
expandido
 
 
##Resolvendo o sistema
for (i in p:2){
    for (j in i:2-1) {
        print(expandido)
        expandido[j, ] <- expandido[j, ] - expandido[i, ] * expandido[j, i]
 
    }
}
expandido
 
x <- matrix(c(1,1,1,1,1,1,1,1,1),3,3)
y <- matrix(c(1,2,3),3,1)
solve(x,y)

Provas matemáticas

https://xkcd.com/435/

Se existe uma coisa que é confusa na ciência, é a ideia de verdade. Talvez muito porque a própria mídia, noticiários e revista em geral confundem bastante o que é ciência empírica (baseada em evidências). Mas também são tantos termos, lei, teoria, hipótese… e assim vai.
Quem trata realmente dessas coisas é a epistemologia, eu tentei explicar uma visão de conhecimento nesse post aqui, a definição Platão é bem legal, mas apesar que como eu comentei, mesmo ela tem problemas.

Mas na matemática, a parada é “simples” (veja que eu coloquei aspas, porque simples é de acreditar, não de colocar em prática, que nem eu sei direito), apesar de discussão sobre o que é ciência, que os filósofos Karl Popper e antes dele, o Kant tem, bem, olhar como as coisas são feitas é legal, e pode ajudar a quando for escrever um manuscrito (não em matemática, mas sobre discutir qualquer coisa), refletir melhor ao ter uma visão de como o conhecimento é feito.

Vamos ver uma coisa legal, sabemos que um número primo é um valor que é somente divisível por 1 e ele mesmo. É algo testável, logo, podemos fazer uma função que fala se é verdadeiro ou falso, que um número seja primo, por exemplo assim.

1
2
3
4
5
6
7
8
9
primo <- function(numero) {
   if (numero == 2) {
      TRUE
   } else if (any(numero %% 2:(numero-1) == 0)) {
      FALSE
   } else { 
      TRUE
   }
}

O any é so para ver se tem algum TRUE, dessa forma a gente evita usar loop de for, que deixa lento a função. Mas pra nossa necessidades, assim ta bom. Agora vamos usar essa função, para os valores de 2 a 10 para testar.

> exemplo<-data.frame(numero=2:10) > exemplo$primo<-sapply(exemplo$numero,primo) > exemplo numero primo 1 2 TRUE 2 3 TRUE 3 4 FALSE 4 5 TRUE 5 6 FALSE 6 7 TRUE 7 8 FALSE 8 9 FALSE 9 10 FALSE

Certo, dois é o único número par primo, 3 é primo por que so é divisível por 2 e 3, quatro não é, porque 4 é igual a 2 vezes 2, 5 é primo, 6 é igual a 2 vezes 3, ou seja, é divisível por 2, bem e assim vai, nada de mágico até agora. Mas olha que interessante, veja o que acontece se para cada número desse, a gente elevar ele a 2^n -1 e testar também se ele é primo.

> exemplo$pot<-(2^exemplo$numero)-1 > exemplo$primo_pot<-sapply(exemplo$pot,primo) > exemplo numero primo pot primo_pot 1 2 TRUE 3 TRUE 2 3 TRUE 7 TRUE 3 4 FALSE 15 FALSE 4 5 TRUE 31 TRUE 5 6 FALSE 63 FALSE 6 7 TRUE 127 TRUE 7 8 FALSE 255 FALSE 8 9 FALSE 511 FALSE 9 10 FALSE 1023 FALSE

Veja que podemos ver um padrão aqui, sempre que um número, que vamos chamar de n, então, sempre que n é primo, 2^n-1 também é primo.

Será que esse padrão continua para sempre? Da uma vontadinha de olhar esse resultado e dizer que sim não? Isso é o que os matemáticos chamam de conjecturas. Viu alguma coisa legal, ai você acha que ela sempre ocorre, ai escreve isso com conjecturas.

Conjectura 1: Suponha que n é um número inteiro maior que 1 e que n é primo. Então 2^n-1 é primo também.

Conjectura 2: Suponha que n é um número inteiro maior que 1 e que n não é primo. Então 2^n-1 não é primo também.

Então a partir desse padrão, conseguimos gerar 2 conjecturas, mas conjecturas são só chutes, mas podemos ir além. Nesse caso alias, se aumentarmos em mais uma linha nossa tabela.

> exemplo numero primo pot primo_pot 1 2 TRUE 3 TRUE 2 3 TRUE 7 TRUE 3 4 FALSE 15 FALSE 4 5 TRUE 31 TRUE 5 6 FALSE 63 FALSE 6 7 TRUE 127 TRUE 7 8 FALSE 255 FALSE 8 9 FALSE 511 FALSE 9 10 FALSE 1023 FALSE 10 11 TRUE 2047 FALSE

Veja que a conjectura 1 vai pro beleleu, 11 é primo, mas 2^{11}-1=2047=23 \cdot 89, então 2^{11}-1 não é primo, ou seja, ele é um contra exemplo, e só precisamos de um contra exemplo para ver que a nossa conjectura 1 não está correta. Mas nesse caso, é so avançar ainda mais um pouco que veremos muito outros contra exemplos, indo até 30 temos o 23, que é primo mas 2^{23}-1=8388607=47 \cdot 178481 e temos o 29 que é primo mas 2^{29}-1=536870911=2089 \cdot 256999.

Triste, a conjectura 1 já era, mas veja que até o número 30, a conjectura 2 está firme e forte.
Encontrar um contra exemplo da conjectura nos permitiu afirmar que ela está errada, mas até agora não achamos contra exemplo para a conjectura 2, mas isso não garante que ela está certa. Pois a gente sempre pode aumentar a nossa tabela, e em algum momento aparecer um contra exemplo. Ou seja, avaliando exemplo nunca da para ter certeza das coisas, pois sempre podem existir mais casos a se avaliar, a única forma de ter certeza é provando. Então vamos tentar fazer uma prova aqui.

Prova da conjectura 2: Como n não é primo, devem existir números inteiros positivos a e b

tal que
a \textless n
b \textless n
n=a \cdot b

Seja
x=2^b-1
y=1+2^b+2^{2b}+2^{3b}+\dots+2^{(a-1)b}

Então
xy=(2^b-1)\cdot(1+2^b+2^{2b}+2^{3b}+\dots+2^{(a-1)b})
realizando a multiplicação
xy=2^b(1+2^b+2^{2b}+2^{3b}+\dots+2^{(a-1)b})-(1+2^b+2^{2b}+2^{3b}+\dots+2^{(a-1)b})
agora a gente faz so a multiplicação do 2^b pelos termos do parenteses.
xy=(2^b+2^{2b}+2^{3b}+\dots+2^{ab})-(1+2^b+2^{2b}+2^{3b}+\dots+2^{(a-1)b})
so lembrando que tem mais um cara escondidinho ali, para ficar mais fácil de ver como o conteúdo dos parenteses são iguais
xy=(2^b+2^{2b}+2^{3b}+\dots+2^{(a-1)b}+2^{ab})-(1+2^b+2^{2b}+2^{3b}+\dots+2^{(a-1)b})
se são iguais e um parenteses é positivo e o outro negativo, a gente pode cortar alguns termos, como o 2^b
xy=(2^{2b}+2^{3b}+\dots+2^{(a-1)b}+2^{ab})-(1+2^{2b}+2^{3b}+\dots+2^{(a-1)b})
agora, se sairmos cortando tudo, ficamos com
xy=(2^{ab})-(1)
que é
xy=2^{ab}-1

Como b \textless n, nós podemos concluir que x=2^b-1 \textless 2^n -1. Também, como ab = n \textgreater a, isso segue que b \textgreater 1. Desta forma, x=2^b-1 \textgreater 2^1-1=1, então y \textless xy = 2^n -1. Veja que nos mostramos que 2^n-1 pode ser escrito como o produto de dois inteiros positivos x e y, ambos os quais são menores que 2^n-1, assim 2^n-1 não é primo.

Agora que a conjectura 2 foi provada, a gente pode chamar ela de um teorema. Apesar dessa prova ser meio misteriosa, existe uma estrutura na forma como ela é escrita. Mas a gente consegue entender que se n é um inteiro não primo maior que 1, ele pode ser escrito como um produto de dois inteiros positivos a e b, então a prova da um método de como escrever 2^n-1 como um produto de dois inteiros menores positivos x e y. E dada a definição de número primo, se n não é primo, então 2^n-1 também não vai ser primo. Por exemplo, pegamos n=12, então 2^n-1=4096. Como 12=3\cdot4, nos podemos pegar a=3 e b=4 na prova. Então de acordo com as formulas para x e y dados na prova, nos teríamos x=2^b-1=2^4-1=15, e y=1+2^b+2^{2b}+\dots+2^{(a-1)b} = 1+2^4+2^8=273. E assim como é predito na formula da prova, nos temos que xy=15\cdot273=4095=2^{12}-1.

Existem outras formas de fatorar 12 em um produto de dois inteiros menores, e esses podem levar a outras formas de fatorar 4095. Por exemplo, como 12=2\cdot6,nos podemos usar valores de a=2 e b=6. Logo x=2^b-1=2^6-1=63, e y=1+2^b+2^{2b}+\dots+2^{(a-1)b} = 1+2^{(2-1)\cdot6}=1+2^6=65. E assim como é predito na formula da prova, nos temos que xy=63\cdot65=4095=2^{12}-1.
Mas uma forma de fatorar já demonstra que o número não é primo.

Bem é isso ai, algo que eu gostaria de aprender melhor, ai comecei a tentar ler um livro sobre o assunto e escrever alguma coisa, o script, apesar de não ter muita coisa, vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, o que é bem provável, pois não sou muito bom de matemática, deixe um comentário corrigindo ou mande um e-mail.

Referência:
Daniel J. Velleman 2006. How to Prove It – A Structured Approach 2ed

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
primo <- function(numero) {
   if (numero == 2) {
      TRUE
   } else if (any(numero %% 2:(numero-1) == 0)) {
      FALSE
   } else { 
      TRUE
   }
}
 
 
exemplo<-data.frame(numero=2:10)
exemplo$primo<-sapply(exemplo$numero,primo)
exemplo
 
exemplo$pot<-(2^exemplo$numero)-1
exemplo$primo_pot<-sapply(exemplo$pot,primo)
exemplo
 
exemplo<-data.frame(numero=2:11)
exemplo$primo<-sapply(exemplo$numero,primo)
exemplo$pot<-(2^exemplo$numero)-1
exemplo$primo_pot<-sapply(exemplo$pot,primo)
 
exemplo