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.

Biogeografia de ilhas

MacArthur e Wilson propuseram a teoria de biogeografia de ilhas onde o número de espécies de uma ilha oceânica seria uma função da taxa de imigração e extinção de espécies. O número de espécies em qualquer tempo é um equilíbrio dinâmico, resultando de ambos os lentos processos de extinção e continua chegada de novas espécies. Assim, é esperado que a composição de espécies pode mudar ao longo do tempo, gerando um “turnover” de espécies.

Vamos pensar da taxa de imigração, y, como uma função do número de espécies, x, que ja está em uma ilha. Essa relação tem uma inclinação negativa, porque conforme o número de espécies aumenta, a chance de um novo individuo chegar ser de uma nova espécie fica mais baixa, lembre-se de como é a distribuição das abundancias das espécies aqui. Assim, a imigração vai ter seu maior valor quando a ilha está vazia, quando o x é igual a zero, e vai cair para zero quando todas as espécies possíveis tiverem chegado la. Além disso, a inclinação deve desacelerar porque algumas espécies são muito mais prováveis que outras de imigrar. Isso significa que a taxa de imigração decai abruptamente assim que as espécies mais prováveis de imigrar chegam, e depois só aquelas que tem pouca chance de dispersar, ou baixas abundancias estão faltando chegar. Denovo pensando nas distribuições de abundancias de espécies, como Preston observou, algumas espécies são simplesmente mais abundantes que outras, produzindo mais indivíduos para dispersar, sendo “melhor dispersores”
Ou seja, esperamos algo desse jeito.

01

Agora vamos imaginar a taxa de extinção, y, como uma função do número de espécies, x. Essa relação vai ter uma relação positiva com o número de espécies, de tal forma que a probabilidade de extinção aumenta com o número de espécies, fácil pensar nisso, quando mais espécies, mais a gente vai adicionando a chance de extinção delas a taxa geral de extinção.

Essa relação também é predita de ter uma inclinação cada vez maior, mais ou menos pelos menos motivos do declínio cada vez maior da imigração. Algumas espécies são menos comuns que outras, e assim mais fácil de se tornarem extintas devido a estocasticidade demográfica e ambiental, tipo um drift genético, e segundo algumas espécies vão ter um menor fitness por muitas razões. Com o acumulo de espécies, maior sera o número de especies propicias a extinção (em relação a baixa abundancia e fitness) presentes, e que vão se extinguir. Lembrando que estamos falando de extinção localmente. Porque lembro que uma vez que apresentei um seminário sobre meta-populações,a palavra extinção soava agressiva para as pessoas, mas veja que extinção aqui não significa o fim da especie, só o fim local.
Então localmente, a taxa de extinção vai ser algo assim:

02

Assim, a taxa de mudança do número de especies em uma ilha vai ser \Delta R, que vai ser a diferença entre imigração I e extinção E, ou

\Delta R = I-E

Quando \Delta R=0, nos temos um equilíbrio. Se soubermos a forma quantitativa da imigração e extinção, nos podemos resolver para o equilíbrio. Esse equilíbrio vai ser um ponto no eixo x, onde as duas taxas se cruzam.

Na teoria do MacArthur e Wilson de biogeografia de ilhas, essas taxas podem ser determinadas pelo tamanho das ilhas, onde:

  • Ilhas maiores tem menor taxa de extinção por causa do maior tamanho médio das populações
  • Ilhas maiores tem maior taxa de colonização porque são alvos maiores para a dispersão

A distância entre as ilhas e as fontes de propágulos também influência essas taxas da seguinte forma:

  • Ilhas próximas do continente tem maiores taxas de colonização de novas espécies porque mais propágulos podem chegar, eles não morrem na jornada por exemplo
  • O efeito da área seria mais importante em ilhar longe do continente do que ilhas perto do continente, para reverter o continuo processo de extinção

Derivar uma curva de imigração ou extinção é um bocado complicado. Mas para ilustrar a teoria, a gente pode assumir que a taxa de imigração I, pode ser representada como uma função exponencial negativa e^{I_0-iR}, onde I_0 é a taxa de imigração em uma ilha vazia e i é o efeito negativo por espécies, que multiplica R que é quantas espécies tem la.

Então a gente pode criar um exemplo, com o número de espécies na ilha:

1
exemplo<-data.frame(especies=seq(0.1,50,0.2))

E ver a taxa de acordo com o número de espécies:

1
2
3
4
##Imigração
IO<-log(1)
b<-0.1
exemplo$imigracao<-exp(IO-b*exemplo$especies)

Agora para a extinção, E, está precisa ser zero se não houver espécies presentes. Imaginando que a extinção é uma função da densidade e que a densidade média declina conforme o número de espécies aumenta, então \bar{N}=\frac{1}{R}. então a gente vai ter algo mais ou menos assim

e^{dR}-1

Nos subtraimos esse 1 para ter certeiza que E=0 quando R=0.
O número de espécies, R, vai resultar de \Delta R = 0 = I-E, o ponto onde as linhas se cruzam.

I=e^{I_0-iR}
E=e^{dR}-1
\Delta R = 0 = I-E

03

O exato ponto onde essas linhas se cruzam podem ser encontrada de muitas formas, a gente tem visto isso em evolução aqui, mas vamos criar uma função no R para minimizar, nos vamos minimizar (I-E)^2, porque elevando ao quadrado a diferença nos temos uma quantidade com a conveniente propriedade que o mínimo ira ser aproximado tanto pelo lado positivo como o negativo.

1
2
3
deltaR<-function(R){
    (exp(IO-b*R) - (exp(d*R)-1))^2
}

Nos usamos essa função de um parâmetro no optimize, e especificamos que queremos saber o ótimo em um intervalo de 1 a 50.

> optimize(f=deltaR,interval = c(1,50)) $minimum [1] 16.91321 $objective [1] 2.813237e-13

A saída nos diz que o mínimo é algo em torno de 16.9, e podemos colocar isso na figura

04

Agora suponha uma ilha um pouco mais longe, do mesmo tamanho, mas mais longe do continente, veja que a taxa de imigração diminuirá, porque pensa, as espécies vão morrer antes de alcançar a ilha e isso vai resultar um outro ponto mínimo, outro ponto de equilíbrio.

05

A beleza dess teoria é que ela foca nossa atenção em processos a nível de paisagem, geralmente fora dos limites espaciais e temporais das nossas capacidades de amostragem. Ele demostra que qualquer fator que ajude a determinar taxas de imigração ou extinção, incluindo a area de uma ilha ou a sua proximidade a uma fonte de propágulos, vai alterar o equilíbrio do número de espécies. Mas temos que lembrar que a identidade das espécies pode mudar ao longo do tempo, isso é, um turnover de espécies, mas um turnover devagar, porque as espécies que tem uma grande chance de imigrar e menor chance de extinguir vão se manter no ambiente muito provavelmente.

Bem é isso ai, incrível como as teorias mais legais são extremamente simples, mas com um efeito profundo na forma como pensamos sobre as coisas, 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:
M.H.H. Stevens 2011, A Primer of Ecology 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
library(ggplot2)
 
exemplo<-data.frame(especies=seq(0.1,50,0.2))
##Imigração
IO<-log(1)
b<-0.1
exemplo$imigracao<-exp(IO-b*exemplo$especies)
 
ggplot(data=exemplo,aes(x=especies,y=imigracao))+geom_line()+
    xlab("Número de espécies") +
    ylab("Taxa de imigração")
ggsave("01.jpg")
##curve(exp(IO-b*x),0,50,xlab="n. de espécies",ylab="Taxa Imigração")
 
##Extinção
d<-0.01
exemplo$extincao<-exp(d*exemplo$especies)-1
ggplot(data=exemplo,aes(x=especies,y=extincao))+geom_line()+
    xlab("Número de espécies") +
    ylab("Taxa de extinção")
ggsave("02.jpg")
##curve(exp(d*x)-1,0,50,xlab="n. de espécies",ylab="Taxa Extinção")
 
 
##Os dois juntos
exemplo2<-stack(exemplo[,c("imigracao","extincao")])
colnames(exemplo2)<-c("valores","tipo")
exemplo2$especies<-rep(exemplo$especies,2)
 
ggplot(data=exemplo2,aes(x=especies,y=valores,group=tipo,colour=tipo))+geom_line()+
    scale_colour_discrete(name="Taxas",labels=c("Extinção","Imigração"))+
    xlab("Número de espécies") +
    ylab("Taxa")
ggsave("03.jpg")
##
deltaR<-function(R){
    (exp(IO-b*R) - (exp(d*R)-1))^2
}
estavel<-optimize(f=deltaR,interval = c(1,50))
segmento<-data.frame(x = c(estavel$minimum,estavel$minimum),y=c(0,max(apply(exemplo[,2:3],1,min))))
 
ggplot(data=exemplo2,aes(x=especies,y=valores,group=tipo,colour=tipo))+geom_line()+
    scale_colour_discrete(name="Taxas",labels=c("Extinção","Imigração"))+
    geom_line(aes(x = x, y = y,group="Estável", colour = "Estável"),data = segmento,colour="black")+
    geom_point(aes(x = x, y = y,group="Estável", colour = "Estável"), data = segmento[2,],colour="black")+
    xlab("Número de espécies") +
    ylab("Taxa")
ggsave("04.jpg")
 
exemplo<-data.frame(especies=seq(0.1,50,0.2))
exemplo$isp1<-exp(IO-0.1*exemplo$especies)
exemplo$isp2<-exp(IO-0.14*exemplo$especies)
exemplo$extincao<-exp(d*exemplo$especies)-1
exemplo2<-stack(exemplo[,c("isp1","isp2","extincao")])
colnames(exemplo2)<-c("valores","tipo")
exemplo2$especies<-rep(exemplo$especies,3)
 
##
deltaR<-function(R,b){
    (exp(IO-b*R) - (exp(d*R)-1))^2
}
b<-0.1
estavel<-optimize(f=deltaR,interval = c(1,50),b=b)
segmento1<-data.frame(x = c(estavel$minimum,estavel$minimum),y=c(0,max(apply(exemplo[,c(2,4)],1,min))))
b<-0.14
estavel<-optimize(f=deltaR,interval = c(1,50),b=b)
segmento2<-data.frame(x = c(estavel$minimum,estavel$minimum),y=c(0,max(apply(exemplo[,c(3,4)],1,min))))
 
ggplot(data=exemplo2,aes(x=especies,y=valores,group=tipo,colour=tipo))+geom_line()+
    scale_colour_discrete(name="Taxas",labels=c("Extinção","Imigração Sp1","Imigração Sp2"))+
    geom_line(aes(x = x, y = y,group="Estável", colour = "Estável"),data = segmento1,colour="black")+
    geom_point(aes(x = x, y = y,group="Estável", colour = "Estável"), data = segmento1[2,],colour="black")+
    geom_line(aes(x = x, y = y,group="Estável", colour = "Estável"),data = segmento2,colour="black")+
    geom_point(aes(x = x, y = y,group="Estável", colour = "Estável"), data = segmento2[2,],colour="black")+
    xlab("Número de espécies") +
    ylab("Taxa")
ggsave("05.jpg")

Rosalind – Introduction to Random Strings

Nesse problema, a gente tem que simular strings baseado na quantidade de GC dela. Uma das relevâncias do conteúdo de GC, se eu não estiver falando muita bobeira, é que G e C se ligam fazendo 3 pontes de hidrogênio, o que confere estabilidade a Molécula de DNA, mas tem muitos motivos diversos para medir a quantidade de GC.

Ok, mas sem devagar do problema, nossa temos uma sequência de DNA e um vetor de valores entre 0 e 1, que são proporções do conteúdo de GC de uma molécula de DNA. E temos que calcular qual a probabilidade de, ao acaso, a sequência acima ter sido gerada a partir desse conteúdo de GC.
Isso é muito parecido com sortear caras e coroas a partir de um likelihood (aqui, um post de 2012, esse é velho hehe).

Nossa entrada é:

ACGATACAA 0.129 0.287 0.423 0.476 0.641 0.742 0.783

Como são quatro bases, temos 4 possibilidades para cada posição, então temos 4\cdot4\cdot4\cdot4\cdot4\cdot4\cdot4\cdot4\cdot4=4^9=262144 possibilidades de moléculas de DNA. Agora suponto que a gente atribua uma probabilidade p de 0.129 para sair C ou G e 1-p, 0.871 para A,T, qual a chance de sair uma molécula com a mesma quantidade de GC que estamos vendo.

Primeiro vamos colocar esses dados no python

1
2
3
dna='ACGATACAA'
prob_gc='0.129 0.287 0.423 0.476 0.641 0.742 0.783'       
prob_gc = map(float, prob_gc.split())

A função map é para transformar string da leitura em um tipo que nos interessa, nesse caso float. Dai contamos a quantidade de GC

1
2
3
4
5
6
7
conteudo_gc = 0
conteudo_at = 0
for codon in dna:
	if codon in ['C', 'G']:
		conteudo_gc += 1
        else:
                conteudo_at += 1

Agora, pensando em uma quantidade apenas, para cada base, ela pode ser GC ou não, então é a chance da base ser GC ou não, coincidir com o que esperamos e a próxima também, e a próxima também. E em probabilidade, um ‘e’ significa que multiplicamos probabilidades.

1
2
3
4
5
6
7
gc=0.129
total=1.0
for base in dna:
        if base in ["G","C"]:
                total*=  gc*0.5
        else:
                total*=  (1-gc)*0.5

Mas veja que como falamos de uma base no conjunto ATCG, estamos falando da metade da probabilidade, dae multiplicamos por meio, porque a outra metade não ia resolver, certo? Mas então somamos tudo isso e é só tirar o log10

1
print log10(total)

Mas ai entramos em um problema de computação, que é o underflow, veja que a gente ta multiplicando e multiplicando um número, e ele ta ficando pequeno e pequeno até não ser possível mais representá-lo, e é ai que o log começa a valer a pena.

Veja que logs tem a propriedade de

log(x\cdot y) = log(x)+log(y)

E somando, a gente não vai encolhendo os números, então resolve o problema de precisão. Facilita bem pelo menos. Então podemos modificar o código acima para

1
2
3
4
5
6
7
8
total=0.0
for base in dna:
        if base in ["G","C"]:
                total+=  log10(gc*0.5)
        else:
                total+=  log10((1-gc)*0.5)
 
print total

Agora é só fazer isso, conteúdo de GC por conteúdo de GC

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
from math import log10
 
with open('/home/augusto/Downloads/rosalind_prob.txt') as input_data:
    dna, prob_gc = input_data.readlines()
 
prob_gc = map(float, prob_gc.split())
 
#Conta o número de GC e AT
conteudo_gc = 0
conteudo_at = 0
for codon in dna:
	if codon in ['C', 'G']:
		conteudo_gc += 1
        else:
                conteudo_at += 1
 
prob_lista  = []
for prob in prob_gc:        
	prob_final = conteudo_gc*log10(0.5*prob) + (len(dna)-conteudo_gc) * log10(0.5*(1-prob))
	prob_lista.append(str(prob_final))
 
print ' '.join(prob_lista)

Outra coisa, é que como estamos falando de quantidade de GC, não interessa a ordem das bases, e por último, temos a intuição, que o maior likelihhod vai ser no conteúdo de GC mais próximo do valor encontrado na sequência mesmo. E o join no string final é só para não imprimir a resposta desorganizada, como quando a gente imprime uma lista no python.

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:
Rosalind

Á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