Equilíbrio em duas populações, entendendo o ponto de atração.

Começamos a olhar a alguns posts, como esse aqui sobre interação entre duas espécies. E no último post, esse aqui, terminamos desenhando as isoclinas para as duas espécies e falando que existia um ponto de atração no meio onde elas se cruzavam, mas o que aquilo quer dizer?

Bem primeiro vamos rever qual gráfico estamos falando.

01

Então temos ai, no eixo x o tamanho populacional da espécie 1, no eixo y o tamanho populacional da espécie 2. Nesse caso temos representado as isolinhas das espécies, que é quando zeramos o crescimento de uma espécie em relação a outra, e no final isso so depende das interações, os [latez]\alpha[/latex] das especies.

Os que foram usados para gerar essa figura foram esses:

> a [,1] [,2] [1,] 0.010 0.005 [2,] 0.005 0.010

Representando como uma matriz por fica simples a notação, ja que por exemplo o [latez]\alpha_{2,1}[/latex] ta na linha 2 coluna 1.
Até aqui beleza, mas se você olhou o primeiro post que coloquei, ele tem essa formula:

 \begin{cases}  N_{1,t+1}= N_{1,t}+ N_{1,t}\cdot(r_{1})\cdot (1-{N_{1,t}\cdot \alpha_{11}} - {N_{2,t}\cdot \alpha_{21}}) \\  N_{2,t+1}= N_{2,t}+N_{2,t}\cdot(r_{2})\cdot (1-{N_{2,t}\cdot \alpha_{22}} - {N_{2,t}\cdot \alpha_{12}})  \end{cases}

Vimos como ver dinâmica ao longo do tempo da espécies usando um parzinho de equações diferenciais. Passando para aquele formatinho de dn/dt, podemos usar o pacote deSolve no R para ficar resolvendo ela e avaliar o resultado da seguinte forma.

#Modelo para crescimento continuo
lvcomp2 <- function(t, n, parms) {
  with(as.list(parms), {
    dn1dt <- r1*n[1]*(1-a11*n[1] - a12*n[2])
    dn2dt <- r2*n[2]*(1-a22*n[2] - a21*n[1])
    list(c(dn1dt, dn2dt))
  } )
}
 
#Resolvendo continuamente com deSolve
parms <- c(r1=0.8,r2=0.5,a11=0.010,a21=0.005,a22=0.010,a12=0.005);
initialN<-c(1,1)
out<-ode(y=initialN, times=1:50, func=lvcomp2, parms=parms)

Ai normalmente a gente faz um gráfico para o resultado assim.

02

Bem, fora a matriz de alpha, eu adicionei aqui valores de r para cada espécie, para esse gráfico a espécie um tem r=0.8 e para a especie 2 temos r=0.5. além disso temos que começar de algum lugar, aqui começamos as duas espécies com 1 indivíduo cada.
Nesse gráfico a gente tem a mudança do tamanho das populações, um valor de N1 e N2 para cada valor de tempo, avaliamos aqui 50 unidades de tempo, dias, meses, anos, vai depender do que estamos falando, mas aqui é irrelevante, o que interessa é que passou 50 unidades, e temos os tamanhos populacionais no tempo 1, no tempo 2, 3, 4 … ,50. Agora ao invés de fazer o gráfico assim podemos simplesmente ter no eixo x o tamanho da população 1, N1, no eixo y o tamanho da população 2, N2, e para cada um das 50 unidades de tempo, a gente poe uma bolinha, e para conseguir manter a noção de pra que lado o tempo ta rolando, vamos manter os pontos ligados.
Então desse jeito a figura fica assim:

03

Ae, continuamos tendo os mesmos dados, mas agora estamos vendo como fica o tamanho das populações das 2 espécies ao longo do tempo. Veja que o que ta acontecendo ai é exatamente o mesmo nos dois gráficos, o inicio é com as duas espécies tendo um individuo cada, ai vai mudando e chega uma hora que estabiliza num ponto.

Agora, a gente pode fazer o seguinte. A gente disse que esse ponto de atração não depende do r (taxa de crescimento) e nem da onde começa as populações das espécies, a gente pode tentar o seguinte, vamos sortear dois valores quaisquer de r (taxa de crescimento) e dois tamanhos populacionais iniciais. E vemos o que acontece.

04

Olha ai, começamos em outro ponto, teve uma curvinha diferente, essa curvinha depende das diferenças entre os r, veja denovo a primeira figura e como uma cresce antes e depois diminui, então gera essa curvinha, mas fomos parar no mesmo lugar.

05

Agora como a gente quer verificar se realmente isso sempre ta acontecendo, a gente faz esse mesmo esquema, simula 10 casos, cada caso tem então 2 espécies, cada especie com um r próprio, um tamanho inicial próprio, e isso diferente pra cada caso, e vemos o que acontece, exatamente que esperávamos, todo mundo ta indo pro mesmo lugar, como sendo puxado, por isso chamam esse pontinho onde as isoclinas se cruzam de ponto de atração.

So que lembramos que a espécie 1 tem uma capacidade suporte de 1\over{\alpha_{11}}, então no eixo x, podemos substituir o valor de 100 por 1\over{\alpha_{11}}.

Relembrando os valores de \alpha temos:

> a [,1] [,2] [1,] 0.010 0.005 [2,] 0.005 0.010

E o \alpha_{11} é 0.010, 1/0.010 é 100, logo a gente substitui ali na no eixo do gráfico, e assim temos eixos como o do primeiro gráfico. E agora o mais legal, se fizermos o mesmo procedimento usado no gráfico anterior, mas plotando os pontinhos no primeiro gráfico, agora com as isoclinas, e com mais populações, veremos o seguinte.

06

Tharam, uma convergência sempre para o mesmo ponto, e essa convergência depende simplesmente de como as espécies interagem umas com as outras. Isso é um resultado interessante, ja que começamos a esperar que interações sejam um bocado determinantes para coexistência entre especies, ou ainda que isso seja importante também para pensar em outras características derivadas de coexistência ou competição, como riqueza de especies no ambiente, diversidade biológica, etc.

Agora estamos a um passo para entender essa parte do livro de ecologia do Begon, que sempre foi um mistério 🙂

begonecology

Bem ficamos por aqui, o endereço do repositório no github é esse, com o código para produzir essas figuras, além de aqui em baixo. Até mais.

Referencias
M Henry H Stevens 2011 Primer of Ecology With R. Springer
M. Begon, C. R. Townsend, J. L. Harper 2005 – Ecology: From Individuals to Ecosystems, 4th Edition – Wiley

#########################################################################
#Testando o atrator
#########################################################################
library(deSolve)
 
a <- matrix(c(0.01, 0.005, 0.005, 0.01), ncol = 2, byrow = TRUE)
a
 
### Ambas isoclines no mesmo grafico
plot((a[2,2]-a[1,2])/(a[2,2]*a[1,1] - a[1,2]*a[2,1]),(a[1,1]-a[2,1])/(a[2,2]*a[1,1] - a[1,2]*a[2,1]), cex=2,
     pch=1, col=1, axes=FALSE,ylim=c(0, max(1/a[2,2], 1/a[1,2])), xlim=c(0, max(1/a[1,1], 1/a[2,1])),
     ylab=expression("N"[2]), xlab=expression("N"[1]), pty='s')
curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], add=TRUE, lty=2)
curve(1/a[1,2] - (a[1,1]/a[1,2])*x, 0, 1/a[1,1], add=T)
axis(1, c(0, 1/a[1,1], 1/a[2,1]),c(0, expression(1/alpha[11]),expression(1/alpha[21])))
axis(2, c(0, 1/a[2,2], 1/a[1,2]),c(0, expression(1/alpha[22]),expression(1/alpha[12])))
arrows(c(25,25,25), c(25,25,25), c(25,50,50), c(50,25,50),col=1, length=.12, lwd=c(1,1,2), lty=c(2,1,1))
arrows(c(150,150,150), c(150,150,150), c(150,120, 120), c(120,150,120),col=1, length=.12,  lwd=c(1,1,2), lty=c(2,1,1))
arrows(c(10, 10, 10), c(125,125,125), c(10,30,30), c(105,125,105),length=0.12, col=1,  lwd=c(1,1,2), lty=c(2,1,1))
arrows(c(125,125,125), c(10, 10,10), c(125,105,105), c(30,10,30),length=0.12, col=1,  lwd=c(1,1,2), lty=c(2,1,1))
 
#Modelo para crescimento continuo
lvcomp2 <- function(t, n, parms) {
  with(as.list(parms), {
    dn1dt <- r1*n[1]*(1-a11*n[1] - a12*n[2])
    dn2dt <- r2*n[2]*(1-a22*n[2] - a21*n[1])
    list(c(dn1dt, dn2dt))
  } )
}
 
 
#Resolvendo continuamente com desolve
a
parms <- c(r1=0.8,r2=0.5,a11=0.010,a21=0.005,a22=0.010,a12=0.005);
initialN<-c(1,1)
out<-ode(y=initialN, times=1:50, func=lvcomp2, parms=parms)
 
#
matplot(out[,1],out[,-1],type='l',xlab="Tempo",ylab="Tamanho populacional",
        frame.plot=F,main="Crescimento logístico continuo para duas espécies")
legend("right",c(expression("Espécie 1 "*(alpha[21]==0.005)),
                expression("Espécie 2 "*(alpha[12]==0.005))),
       lty=1:2,bty='n',col=c("black","red"))
 
#Representação das abundancias das duas espécies num plano
plot(1, 1, type = "p", ylim = c(0, 200), xlim = c(0,200),frame=F,
     ylab = expression("N"[2]),xlab=expression("N"[1]),pch=19)
points(out[,2],out[,3],type="b",cex=0.5)
 
#Agora um caso onde as taxas de crescimento são definidas ao acaso,
#assim como a população inicial
plot(1, 1, type = "n", ylim = c(0, 200), xlim = c(0,200),frame=F,
     ylab = expression("N"[2]),xlab=expression("N"[1]),pch=19)
 
parms[1:2]<-runif(2,0,1)
initialN<-runif(2,1,200)
out<-ode(y=initialN, times=1:50, func=lvcomp2, parms=parms)
points(out[,2],out[,3],type="b",cex=0.5,col=1)
 
#O mesmo procedimento anterior varias vezes no mesmo grafico
plot(1, 1, type = "n", ylim = c(0, 200), xlim = c(0,200),frame=F,
     ylab = expression("N"[2]),xlab=expression("N"[1]),pch=19)
 
for(i in 1:10) {
    parms[1:2]<-runif(2,0,1)
    initialN<-runif(2,1,200)
    out<-ode(y=initialN, times=1:50, func=lvcomp2, parms=parms)
    points(out[,2],out[,3],type="b",cex=0.5,col=i)
    }
 
#O mesmo procedimento para o nosso grafico de isoclines
plot((a[2,2]-a[1,2])/(a[2,2]*a[1,1] - a[1,2]*a[2,1]),(a[1,1]-a[2,1])/(a[2,2]*a[1,1] - a[1,2]*a[2,1]), cex=2,
     pch=1, col=1, axes=FALSE,ylim=c(0, max(1/a[2,2], 1/a[1,2])), xlim=c(0, max(1/a[1,1], 1/a[2,1])),
     ylab=expression("N"[2]), xlab=expression("N"[1]), pty='s')
curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], add=TRUE, lty=2)
curve(1/a[1,2] - (a[1,1]/a[1,2])*x, 0, 1/a[1,1], add=T)
axis(1, c(0, 1/a[1,1], 1/a[2,1]),c(0, expression(1/alpha[11]),expression(1/alpha[21])))
axis(2, c(0, 1/a[2,2], 1/a[1,2]),c(0, expression(1/alpha[22]),expression(1/alpha[12])))
arrows(c(25,25,25), c(25,25,25), c(25,50,50), c(50,25,50),col=1, length=.12, lwd=c(1,1,2), lty=c(2,1,1))
arrows(c(150,150,150), c(150,150,150), c(150,120, 120), c(120,150,120),col=1, length=.12,  lwd=c(1,1,2), lty=c(2,1,1))
arrows(c(10, 10, 10), c(125,125,125), c(10,30,30), c(105,125,105),length=0.12, col=1,  lwd=c(1,1,2), lty=c(2,1,1))
arrows(c(125,125,125), c(10, 10,10), c(125,105,105), c(30,10,30),length=0.12, col=1,  lwd=c(1,1,2), lty=c(2,1,1))
 
 
for(i in 1:50) {
    parms[1:2]<-runif(2,0,1)
    initialN<-runif(2,1,200)
    out<-ode(y=initialN, times=1:50, func=lvcomp2, parms=parms)
    points(out[1,2],out[1,3],type="p",pch=19,col=i,cex=0.5)
    points(out[,2],out[,3],type="l",lwd=0.7,col=i)
    }

Um exemplo do algoritimo de Metropolis–Hastings aplicado a distribuição binomial

E la vamos nos denovo olhar o algoritmo Metropolis–Hastings

Bem depois de olhar o exemplo do como estimar o valor de pi usando o método de monte carlo, devemos progredir um pouco no entendimento desse negócio de mcmc.

Ja tivemos posts aqui e aqui, onde olhamos como o algoritmo funcionava para fazer uma regressão linear.

Mas vamos tentar um exemplo mais simples para estimar a probabilidade encontrar ou não aranhas num arbusto.

Certo então existe uma probabilidade \theta (letra grega theta) na qual eu encontro ou não as aranhas nos arbustos.

Então eu fui la e olhei 14 arbustinhos e vi o seguinte:

> meusdados [1] 0 1 1 0 1 0 0 1 0 1 1 1 1 1

Ok 1 significam que nesse arbusto eu encontrei uma aranha, e 0 significa que não tinha aranha.

Eu posso estimar a verosimilhança para várias probabilidades, para vários valores de \theta, ou mais precisamente, o p(D|\theta) onde D são meus dados

likelihood <- function( theta , data ) {
  z <- sum( data == 1 )
  N <- length( data )
  pDataGivenTheta <- theta^z * (1-theta)^(N-z)
  pDataGivenTheta[ theta > 1 | theta < 0 ] <- 0
  return( pDataGivenTheta )
}

Ou seja qual a chance de determinado \theta dado meus dados (p(D|\theta)). Aqui tem uma revisão rápida de probabilidade condicional.

Bom mais fácil que ficar nisso é fazer um gráfico para entender o que essa função ta fazendo.

01

Então temos a como são as chances dado um valor de \theta, ou seja olhando os dados la tivemos 9 arbustos com aranhas, dos 14 que olhamos, sabemos que o maximum likelihood nesse caso é 9/14, ou seja o pico do likelihood é quanto \theta esta num valor igual a 0.64, ou seja um \theta de 0.64 é o mais provável de ter produzido esses resultado, dado os dados que temos, veja que se você pegar um valor de vamos supor 0.2, quase não tem como achar 9 arbustos com aranhas de 14 arbustos verificados, teria que ser uma cagada muito grande.

Mas lembramos que em inferência bayesiana, ponderamos nossos dados, ou seja multiplicamos esse likelihood por um prior.

Então fazemos uma função para gerar esse prior.

prior <- function( theta ) {
  prior <- dbeta(theta,1,1)
  #exemplo de um prior bimodal
  #prior < dbeta( pmin(2*theta,2*(1-theta)),2,2 )
  prior[ theta > 1 | theta < 0 ] <- 0
  return( prior )
}

Aqui, o prior ta definido ali como uma distribuição beta com parâmetros a=1 e b=1. Isso vai dar uma gráfico assim:

02

Veja que todas as possibilidades estão em 1, ou seja eu estou atribuindo uma chance igual para todos os valores possíveis de \theta.

Veja que podemos mudar nossa crença inicial mudando os parâmetros a e b para outros valores, se você não sabe o que acontece com diferentes valores, é mais ou menos isso, veja alguns exemplos:

01

Certo, mas para facilidade do exemplo vamos assumir aquele prior e não ficar discutindo muito isso agora.

Agora se lembrarmos o teorema de bayes usado na inferência bayesiana:

p(\theta|D) = {p(D|\theta) * p(\theta) \over {\int p(D|\theta)p(\theta)}}

Então o que falta é sabe qual o valor desse divisor. Na verdade para esse exemplo da pra resolver analiticamente como feito aqui ou aproximar como aqui, mas vamos usar o algorítimo de Metropolis–Hastings para construir a probabilidade posterior.

A ideia é a mesma de quando estimamos o valor de pi. Relembrando os passos:

1. Defina um domínio de todas as possibilidades 2. Gere aleatoriamente entradas de uma distribuição de probabilidade que cobre o domínio. 3. Faça uma computação determinística das entradas 4. Reuna os resultados

Então o domínio são as probabilidades de achar aranha em um arbusto, que são valores entre 0 e 1.

Eu tenho que gerar valores ao acaso de probabilidade, que cubram todas as possibilidades e fazer uma computação determinística para levar ele em conta ou não, no caso do pi o que a gente gerava um ponto ao acaso dentro de um quadradinho e tínhamos um círculo dentro desse quadrado, ai olhávamos se ele estava ou não dentro do círculo, para no final fazer uma razão.

Aqui o que temos são valores possíveis de \theta e eu quero saber qual gerou os dados, bom isso eu nunca vou saber, mas eu quero ter uma estimativa, uma distribuição que se o valor de \theta seja bem razoável, ele tem que ser alto, senão ele é baixo.
Talvez eu esteja falando bobeira aqui, mas o que acontece é o seguinte, no caso do Metropolis–Hastings, ficamos escolhendo um valor de \theta ao acaso, e colocamos ele ou não na nossa cadeia numa probabilidade que depende de o quão bom o novo \theta é dividido pelo \theta anterior.

É difícil explicar isso, acho que entender também, mas fazemos uma função assim:

targetRelProb <- function( theta , data ) {
  targetRelProb <-  likelihood( theta , data ) * prior( theta )
  return( targetRelProb )
}

Essa função multiplica o likelihood pelo prior, dai geramos uma valor de theta vai ter essa multiplicação do valor novo de theta pela multiplicação com o antigo.

Vamos olhar alguns exemplos, bem eu usei 0.7 como chance para gerar os dados. Então vamos supor que o valor anterior na cadeia é 0.7 e fosse sugerido fosse de 0.4, qual a chance de aceitar esse valor e colocar na cadeia?

> targetRelProb(0.4,meusdados)/targetRelProb(0.7,meusdados) [1] 0.2078775

20%, baixinho certo, ou seja eu quero “visitar” esse valor pouco, quero poucas ocorrências dele na cadeia, assim no histograma da distribuição posterior ele vai ter uma barrinha baixinha, mas impossível ver um valor desse? Não uai, mas é pouco provável.

Agora vamos mais pra perto, se propusermos um valor de 0.5.

> targetRelProb(0.5,meusdados)/targetRelProb(0.7,meusdados) [1] 0.6224313

Olha ai, temos 60% de chance de enfiar esse valor na cadeia, porque ta mais pertinho do valor real. E quanto mais próximo vamos do valor real, agora com 0.55.

> targetRelProb(0.55,meusdados)/targetRelProb(0.7,meusdados) [1] 0.8666388

Temos mais e mais chance de enfiar esses valores pra dentro da cadeia.
Agora até aqui eu sempre tava dividindo com um theta de 0.7, mas e se sei la estamos no começo da cadeia, e estamos em valores bem distantes desses. Por exemplo estamos em 0.2 e temos uma proposta de ir para 0.1

> targetRelProb(0.1,meusdados)/targetRelProb(0.2,meusdados) [1] 0.003519595

Meu, 0.2 ta longe pra caramba das probabilidade real, mas 0.1 é muito mais longe, então a chance de aceitar esse valor é muito muito baixa. Agora num caso ao contrário, estamos em 0.9, que é longe e recebemos uma proposta de 0.8, que é bem mais próximo, qual a chance?

> targetRelProb(0.8,meusdados)/targetRelProb(0.9,meusdados) [1] 11.08606

Puxa, um valor muito maior que 1, veja que isso pode acontecer, na verdade aqui o que fazemos e com certeza aceitar o valor proposto, como se estamos em 0.5 e recebemos a proposta de 0.6

> targetRelProb(0.6,meusdados)/targetRelProb(0.5,meusdados) [1] 1.690757

Mas para fugir desses valores e ficar com valores sempre entre 0 e 1 na cadeia temos que fazer umas correções nos valores, mas basicamente é isso o algoritimo, gera um theta novo e aceita ele ou não, um bocado de vezes, mas um bocado são muitas e muitas vezes mesmo.

Então o algoritimo vai ser isso aqui:

for ( t in 1:(Ncadeia-1) ) {
    valoratual = cadeia[t]
    valorproposto = rnorm( 1 , mean = 0 , sd = 0.1 )
    chanceaceitar = min( 1,
        targetRelProb( valoratual + valorproposto, meusdados ) /
        targetRelProb( valoratual , meusdados )
        )
    if ( runif(1) < chanceaceitar ) {
        cadeia[ t+1 ] <- valoratual + valorproposto
        } else {
            cadeia[ t+1 ] <- valoratual
	}
}

A gente salva numa variável o valor atual, na prática a gente gera uma quantidade, que é o valor proposto, que é uma quantidade para mudar o valor de theta atual, dai calculamos como comentamos acima a chance de aceitar esse novo valor, e se ele for aceito colocamos ele na cadeia, senão repetimos o último valor. E milagrosamente vamos visitar cada valor possível de theta um número de vezes proporcional ao valor de theta que gerou os dados. Na verdade não é milagre, essa estabilidade é prevista para cadeias de markov dado a chance de transição entre estados, mas pra mim é mágico como funciona 🙂

Então ao usarmos nosso algoritimozinho, o que vemos como resultado?

03

Uma boa estimativa de valor real de theta, o intervalo de 95% da distribuição cobre o valor real de theta.

Mas temos uma distribuição continua de thetas, não quero ficar olhando as coisas na forma de histograma, podemos olhar a mesma distribuição na forma de densidade contínua.

04

Lembrando que aqui a gente ja arrancou o início da cadeia, que como começamos num lugar qualquer, traz estimativas ruim. Vamos dar uma olhadinha na cadeia, para o início e um pedaço la na frente.

05

Como começamos num valor arbitrário qualquer, começamos num valor de theta que não tem nada a ver com o valor real e vamos caminhando ao valor próximo do real, ai chegando la ficamos rodeando em volta dele. Veja que o tamanho dos passos aqui depende do valor proposto, se propomos valores pouco diferentes do anterior, vamos andando devagarinho na cadeia, agora poderíamos propor passos mais largos, mais ai também, assim como poderíamos chegar num valor próximo do real mais rápido, poderíamos derrepente sair dele pra muito longe. Esse tipo de ajustamento, como o tamanho do passo, como propor uma esse passo que é a parte difícil, aqui tudo funciona mas em algum momento alguém deve ter ficado testando possibilidades e ficou frustado em como as coisas não funcionavam. Outra coisa é que pode ser proposto valores maiores que 1 ou negativos, menor que zero, veja que ficamos mudando esses valores para os extremos, 0 e 1, mas eu não sei qual o impacto de fazer isso se o valor de theta fosse 1 por exemplo, mas todas essas coisas influenciam no funcionamento e eficácia do uso de MCMC. Nisso pacotes como OpenBugs ou Stan salvam a vida, ao possibilitar a gente abstrair dessa coisas e focar em entender melhor nossas observações da natureza. Mas ainda sim é importante entender pelo menos um pouco do de tudo que estamos fazendo.

Bem acho que ja esta bom por aqui, depois olhamos mais sobre MCMC. Esse script como outros do blog estão sempre disponíveis no repositório do github aqui.

Referência:
John K. Kruschke 2011 Doing Bayesian Data Analysis: A Tutorial with R and BUGS. Academic Press / Elsevier

#Adaptado de uma função do livro Bayesian Data Analysis: A Tutorial with R and BUGS
#Semente para replicar os resultado do post.
set.seed(51)
 
#Gerando uma amostra de dados
probabilidade<-0.70
meusdados<-rbinom(n=14, size=1, prob=probabilidade)
 
#Definindo a função de likelihood binomial.
likelihood <- function( theta , data ) {
  z <- sum( data == 1 )
  N <- length( data )
  pDataGivenTheta <- theta^z * (1-theta)^(N-z)
  pDataGivenTheta[ theta > 1 | theta < 0 ] <- 0
  return( pDataGivenTheta )
}
 
plot(seq(0.01,0.99,0.01),likelihood(seq(0.01,0.99,0.01),meusdados),frame=F,
     xlab="Probabilidade",ylab="Likelihood",type="l")
 
#Função para definir o prior p(D),
prior <- function( theta ) {
  prior <- dbeta(theta,1,1)
  #exemplo de um prior bimodal
  #prior < dbeta( pmin(2*theta,2*(1-theta)),2,2 )
  prior[ theta > 1 | theta < 0 ] <- 0
  return( prior )
}
 
plot(seq(0.01,0.99,0.01),prior(seq(0.01,0.99,0.01)),frame=F,
     xlab="Probabilidade",ylab="Likelihood",type="l")
 
#Função auxiliar para ajudar a calcular a probabilidade relativa
targetRelProb <- function( theta , data ) {
  targetRelProb <-  likelihood( theta , data ) * prior( theta )
  return( targetRelProb )
}
 
# Tamanho da cadeia
Ncadeia = 100000
 
#Inicializando um vetor para salvar os resultados do MCMC
cadeia = rep( 0 , Ncadeia )
 
# Iniciando o MCMC com um valor arbritario
cadeia[1] = runif(1,min=0,max=1)
 
# Especificando um valor de burnning, ou seja o quanto do inicio da cadeia vamos jogar fora
# Aqui vamo jogar fora 10%
burnIn = ceiling( .1 * Ncadeia )
 
# Vamos contar quantas vezes aceitamos ou não os valores propostos
nAceito = 0
nRejeitado = 0
 
# Agora iniciando o "Random Walk"
for ( t in 1:(Ncadeia-1) ) {
    #Valor Atual
    valoratual = cadeia[t]
    #Propomos uma quantidade para alterar esse valor
    valorproposto = rnorm( 1 , mean = 0 , sd = 0.1 )
    #Computamos a chance de aceitar a mudança
    #Essa chance é igual a Likelihood*prior com o novo valor de theta
    #dividido pele Likeliood*prior antigo
    chanceaceitar = min( 1,
        targetRelProb( valoratual + valorproposto, meusdados ) /
        targetRelProb( valoratual , meusdados )
        )
    #Agora aqui a gente ve na sorte se aceitamos ou não o novo valor.
    #Veja que temos uma chance de fazer a mudança ou manter o ultimo valor
    #Se a gente sortei um número qualquer de 0 a 1 e compara se ele é
    #menor que a chance de aceitar, dessa forma mudamos numa probabilidade
    #da chance de mudar
    if ( runif(1) < chanceaceitar ) {
        #Se aceitarmos mudamos ele na cadeia
        cadeia[ t+1 ] <- valoratual + valorproposto
        #So vamos começar a contar depois do burnIn que escolhemos
        #Afinal, com o começo arbitrario, esperamos aceitar bastante
        #a não ser que escolhemos um valor proximo do final no começo
        if ( t > burnIn ) { nAceito <- nAceito + 1 }
 
        } else {
            # Quando não aceitamos o valor proposto mantemos o anterior
            cadeia[ t+1 ] <- valoratual
            # E incrementamos aqui apos o BurnIn, como acima
            if ( t > burnIn ) { nRejeitado <- nRejeitado + 1 }
	}
}
 
#Agora vamos olhar como ficou a cadeia final, ou seja, tirando os valores iniciais
cadeiafinal <- cadeia[ (burnIn+1) : length(cadeia) ]
 
#Alguns Gráficos
 
#-----------------------------------------------------------------------
HDIofMCMC = function( sampleVec , credMass=0.95 ) {
    sortedPts = sort( sampleVec )
    ciIdxInc = floor( credMass * length( sortedPts ) )
    nCIs = length( sortedPts ) - ciIdxInc
    ciWidth = rep( 0 , nCIs )
    for ( i in 1:nCIs ) {
        ciWidth[ i ] = sortedPts[ i + ciIdxInc ] - sortedPts[ i ]
    }
    HDImin = sortedPts[ which.min( ciWidth ) ]
    HDImax = sortedPts[ which.min( ciWidth ) + ciIdxInc ]
    HDIlim = c( HDImin , HDImax )
    return( HDIlim )
}
 
 
#------------------------------------------------------------------------
intconf<-HDIofMCMC(cadeiafinal)
 
hist(cadeiafinal,ylab="Frequência na cadeia final",xlab="Probabilidades",xlim=c(0,1))
lines(c(intconf[1],intconf[2]),c(500,500),lwd=4,lty=3,col="red")
abline(v=mean(cadeiafinal),lwd=4,col="red")
abline(v=probabilidade,lwd=2,col="blue")
 
legend("topleft",lty=c(3,1,1),lwd=c(4,4,2),bty="n",cex=0.8,col=c("red","red","blue"),
       legend=c("Intervalo de confiança","Média das estimativas","Probabilidade usada nas estimativas"))
 
plot(density(cadeiafinal),ylab="Frequência",xlab="Probabilidades",xlim=c(0,1),frame=F)
lines(c(intconf[1],intconf[2]),c(0.25,0.25),lwd=3,lty=2,col="red")
abline(v=mean(cadeiafinal),lwd=1,col="red")
abline(v=probabilidade,lwd=1,col="blue")
 
legend("topleft",lty=c(3,1,1),lwd=c(2,1,1),bty="n",cex=0.8,col=c("red","red","blue"),
       legend=c("Intervalo de confiança","Média das estimativas","Probabilidade usada nas estimativas"))
 
layout(matrix(c(1,2,3,4),ncol=2))
hist(cadeia[1:100],ylab="Frequência na cadeia",xlab="Probabilidades",main="100 primeiros valores",xlim=c(0,1))
plot(0,0,xlim=c(0,100),ylim=c(0,1),type="n",frame=F,main="Cadeia",
     xlab="Posição na cadeia",ylab="Probabilidade")
points(cadeia[1:100],type="l")
 
hist(cadeia[10000:10500],ylab="Frequência na cadeia",xlab="Probabilidades",main="Do 10000 ao 10500",xlim=c(0,1))
plot(0,0,xlim=c(0,500),ylim=c(0,1),type="n",frame=F,main="Cadeia",
     xlab="Posição na cadeia",ylab="Probabilidade",xaxt="n")
axis(1,at=seq(0,500,by=50),labels=seq(10000,10500,by=50))
points(cadeia[10000:10500],type="l")
 
# Evidencia para o modelo, p(D).
 
media<-mean(cadeiafinal)
desvio<-sd(cadeiafinal)
 
a =   media   * ( (media*(1-media)/desvio^2) - 1 )
b = (1-media) * ( (media*(1-media)/desvio^2) - 1 )
 
wtdEvid = dbeta( cadeiafinal , a , b ) /
    ( likelihood( cadeiafinal , meusdados ) * prior( cadeiafinal ) )
pData = 1 / mean( wtdEvid )
 
format(pData,scientific = F)

Problema – Rosalind – Completing a Tree – TREE

Ja faz um tempo que não olhava mais o Rosalind, aqui vai mais um problema.

Nesse problema a gente tem que completar um grafo do tipo arvore.

Bem, a primeira pergunta é o que diabos é um grafo, um grafo é um conjunto de pontos (vértices ou node em inglês) e suas ligações (as arestas ou edges em inglês).

Existe uma disciplina chamada teoria dos grafos que estuda propriedades e características de grafos. Dentro dessa teoria eles são divididos em vários tipos, que apresentam características em comum, um desses tipos é o tipo arvore ou Tree citado no nome do problema.

O grafo de árvore é um grafo conexo (existe caminho entre quaisquer dois de seus vértices) e acíclico (não possui ciclos). Basicamente esses caras são as arvores filogeneticas.
Ele é conexo porque aceitamos, na hora de reconstruir uma arvore qualquer, cada par de especies tem pelo menos um ancestral comum, e as especies, que aqui são as folhas, não se ligam umas nas outras, folhas no nosso caso são as espécies. Para ter uma ideia de uma olhada aqui.

Mas após ver tudo isso, é so olhar as propriedades dos grafos tipo arvore para matar esse problema, todo grafo do tipo arvore terá n-1 arestas, sendo n o número de folhas ou espécies.

Agora o problema fornece as seguintes informações

10 1 2 2 8 4 10 5 9 6 10 7 9

10 é o número vértices, e essas duplas de números são as arestas.

Usando o R e o pacote igraph, a gente pode tentar construir esse grafo de curiosidade da seguinte forma

library(igraph)
plot(graph.formula( 1--2, 2--8, 4--10, 5--9,6--10,7--9,3 ))

Eu coloquei o 3 ali sozinho no final, porque ele não aparece na lista de arestas, então ele não apareceria.

01

Mas é isso ai. Um monte de partezinha perdidas. Para ele ser uma árvore, ele precisa n-1 arestas, olha ai temos 7 linhas, então ta faltando 3 linhas, então se adicionarmos essas linhas respeitando as outras propriedades do grafo para ser uma uma árvore, temos 7 arestas já, 10-1 da 9, temos que ter nove arestas, então ta faltando 3, adicionamos elas e ficamos com o seguinte.

plot(graph.formula( 1--2, 2--8, 4--10, 5--9,6--10,7--9,2--3,10--3,9--3 ))

02

Agora a visualização pode parecer estranha, isso porque para arvores filogenéticas, estamos mais acostumados com algo mais ou menos assim:

arvore-eps

Uma arvorezinha com raiz. você pode arrumar a disposição dos vértices usando o comando tkplot assim:

tkplot(graph.formula( 1--2, 2--8, 4--10, 5--9,6--10,7--9,2--3,10--3,9--3 ))

Então a conta a fazer é bem simples se temos que ter n-1 arestas sempre, então contamos quantas arestas ja temos e diminuímos o número de vértices do número de arestas ja existentes -1. A não ser que as arestas estejam zuadas, teremos uma árvore.

Então em python fazemos a função:

def minarv(n, edges):
    return n - len(edges) - 1

E para os dados do exemplo é só usar a função:

n=10
edges=[(1, 2), (2, 8), (4, 10), (5, 9), (6, 10), (7, 9)]
 
print minarv(n,edges)
>>>3

A diferença para o teste é que temos muito mais vértices, e muito mais arestas, então não da para escrever a mão, logo temos que ler o arquivo txt de dados, basicamente a primeira linha é o número de vértices o todas as outras linhas são arestas.

file = open( "/home/augusto/Downloads/rosalind_tree.txt" )
n = file.readline()
edges = []
 
for linha in file:
    edges.append(linha.split())
 
file.close()
 
print int(n)
print len(edges)
 
print minarv(int(n),edges)

Aqui é só tomar cuidado que quando a gente le o arquivo em python, tudo é lido como string, então para o número de vértices temos que usar int para transformar ele em um número inteiro, para poder fazer continhas.
Como as arestas apenas temos que contar por enquanto, então não precisamos se preocupar por enquanto.
E voala, concluímos mais um problema do Rosalind.

Agora entrando nesses problemas de árvores filogenéticas deve ficar bem legal, essa é uma área legal.

Boa semana 🙂

Equilíbrio em duas populações, desenhando a Isóclina

Isóclina, ou isocline é uma reta ou curva que tem a mesma inclinação quando comparada a uma outra. Em mapas topográficos por exemplo, é aquela linha que junta os pontos que estão em uma mesma altura. No nosso caso aqui é a linha onde onde o crescimento de uma espécies N_i fica zerado, vemos o que acontece com as outras especies do sistema. Da para procurar essa idéia como Zero net growth isoclines (ZNGI) ou Nullcline no wikipedia.

Certo, se esta complicado, vai ficar mais simples quando tentarmos fazer isso no R.
Primeiro lembramos o do ultimo post como podemos representar um sistema com duas espécies usando o crescimento logístico.

 \begin{cases}  N_{1,t+1}= N_{1,t}+ N_{1,t}\cdot(r_{1})\cdot (1-{N_{1,t}\cdot \alpha_{11}} - {N_{2,t}\cdot \alpha_{21}}) \\  N_{2,t+1}= N_{2,t}+N_{2,t}\cdot(r_{2})\cdot (1-{N_{2,t}\cdot \alpha_{22}} - {N_{2,t}\cdot \alpha_{12}})  \end{cases}

A gente pode representar quantas especies quiser assim, mas infelizmente o número de interações, o número de \alpha_{ij} aumenta numa taxa de2^n, onde n é o número de espécies. Ou seja aqui com duas espécies temos 4 \alpha_{ij}, mas com 3 especies temos 8 \alpha_{ij}, com 4 temos 16, e rapidinho fica impossível pensar com esse modelinho. Mas salvo esse problema, vamos avaliar algumas idéias num modelo com duas espécies que é fácil de entender e pode ajudar a conseguir alguns insights quanto a dinâmica da competição entre especies.

Primeiro vamos olhar como é a taxa de mudança no crescimento, e não a mudança nos tamanhos da população, basicamente é a mesma coisa, mas as taxas de crescimento ficam assim:

 \begin{cases}  {dN_{1}\over{dt}} = r_{1}\cdot N_{1}\cdot (1-{N_{1}\cdot \alpha_{11}} - {N_{2}\cdot \alpha_{12}}) \\  {dN_{2}\over{dt}} = r_{2}\cdot N_{2}\cdot (1-{N_{2}\cdot \alpha_{22}} - {N_{1}\cdot \alpha_{21}})  \end{cases}

Certo, é assim que as populações mudam, então para a população 1, a mudança depende do dela mesma, ou seja depende de quanto é a taxa de crescimento dela r_{1}, mas essa taxa de crescimento diminui de acordo com a interação dela com ela mesma, o \alpha_{11} e como ela interage com as outras especies, nesse caso so existem duas, então so tem uma outra interação que é a \alpha_{21}.
A espécie 2 é a mesma coisa. Mas note que todos esses efeitos multiplicam sempre os tamanhos das populações, então o tamanho “atual” influencia no quão forte é esse efeito, e nesse caso continuamente, não em passos discretos.

Agora a gente pode perguntar como fica o equilíbrio dessas populações, a gente ja falou aqui de equilíbrio e estabilidade de uma população. Mas agora cada população não depende somente dela mesmo, depende do estado da outra população. Então depende de como as coisas estão.

Mas podemos traçar a isóclina da população de uma especie em relação a outra para melhor entender o que deve acontecer.

Primeiro vamos pensar o seguinte, se temos duas espécies, então quais são as possibilidades, qualquer combinação de densidade dessas duas especies certo, então o State space é o espaço de possibilidades possiveis, aqui temos o seguinte.

01

Note que eu deixei um pontinho, o que é esse pontinho? É o caso da espécie 1 e a 2 estarem com tamanho populacional 0. Então o estado dessas duas populações estarem em zero seria ali. Mas poderíamos ter esse pontinho em muitos outros lugares, dependendo de como estão num determinado momento do tempo, ja que as especies mudam a sua abundancia de acordo com o tempo.

Certo, agora a isocline é o linha onde o crescimento de uma população é zero. Ou seja todos os pontinhos onde a população tem crescimento zero. Como a gente tem duas populações aqui, uma da especie 1 e outra da especie 2, teremos duas isocline nesse gráfico acima, que são é o state space, todas as possibilidades de abundancias.

Agora como sabemos onde ficam esses pontos, ora, vimos como é o crescimento para as duas populações.

\begin{cases}  {dN_{1}\over{dt}} = r_{1}\cdot N_{1}\cdot (1-{N_{1}\cdot \alpha_{11}} - {N_{2}\cdot \alpha_{12}}) \\  {dN_{2}\over{dt}} = r_{2}\cdot N_{2}\cdot (1-{N_{2}\cdot \alpha_{22}} - {N_{1}\cdot \alpha_{21}})  \end{cases}

Ele é assim certo?
Agora vamos trabalhar primeiro com a especie 2, a população da especie 2.

 {dN_{2}\over{dt}} = r_{2}\cdot N_{2}\cdot (1-{N_{2}\cdot \alpha_{22}} - {N_{1}\cdot \alpha_{21}})

Agora falar que o crescimento é zero é meramente falar que {dN_{2}\over{dt}} = 0. Certo, então como a gente quer saber quando esse zero acontece, a gente substitui ele ali na formula e temos:

 0 = r_{2}\cdot N_{2,t}\cdot (1-{N_{2}\cdot \alpha_{22}} - {N_{1}\cdot \alpha_{21}})

Agora é so resolver essa equação para N2, e saberemos onde N2 estará de acordo com os valores de N1. A primeira é notar que temos uma coisa multiplicando todo um parenteses, que tem que passar para o outro lado dividindo.

 {0\over{r_{2}\cdot N_{2}}} =  1-{N_{2}\cdot \alpha_{22}} - {N_{1}\cdot \alpha_{21}}

Certo, agora zero dividindo qualquer coisa da zero, então podemos voltar o zero ali.

 0 =  1-{N_{2}\cdot \alpha_{22}} - {N_{1}\cdot \alpha_{21}}

Agora a gente traz o N2 pro outro lado, que vem positivo.

 {N_{2}\cdot \alpha_{22}} =  1 - {N_{1}\cdot \alpha_{21}}

Agora pra deixar somente o N2 isolado, a gente manda aquele \alpha_{22} dividindo.

 N_{2} =  {{1 - N_{1}\cdot \alpha_{21}} \over{ \alpha_{22} }}

E separamos os termos, pra facilitar a visualização do que ta acontecendo

 N_{2} =  {1\over{\alpha_{22}}}-{{\alpha_{21}\over{\alpha_{22}}}\cdot N_{1}}

E onde fomos para, na boa e velha equação da reta, a+bx, então a isocline nesse caso é uma reta.
Veja que {1\over{\alpha_{22}}} é um termo sozinho, um intercepto, enquanto -{{\alpha_{21}\over{\alpha_{22}}}\cdot N_{1}} temos um termo multiplicando o tamanho da população da espécie 1, a inclinação.

Agora a gente passa isso pro R.

Primeiro, para estudar isso precisamos de um exemplo, então vamos estabelecer os parametros da população,precisamos chutar valores para alpha para tentar fazer um exemplo, então criamos uma matriz de alphas.

> a [,1] [,2] [1,] 0.010 0.005 [2,] 0.005 0.010

Então temos o elemento a[1,1] que é o \alpha_{11}, e assim para todos os elementos.

Agora podemos descrever aquela equação como uma expressão no R, da seguinte forma.

N2iso <- expression(1/a[2, 2] - (a[2, 1]/a[2, 2]) * N1)

O que é isso?
Bem criamos uma objeto da classe expressão, que podemos resolver para valores de N1

expression(1/a[2, 2] - (a[2, 1]/a[2, 2]) * N1) > class(N2iso) [1] "expression"

Veja que qualquer expressão pode ser avaliada no R, por exemplo dois ao quadrado:

> eval(2^2) [1] 4 > espressao<-expression(2^2) > eval(espressao) [1] 4

Então quando guardamos a expressão, podemos avaliar ela, agora vamos supor que ao inves de somente números temos uma variável, x por exemplo.

> espressao<-expression(x^2) > eval(espressao) Erro em eval(expr, envir, enclos) : objeto 'x' não encontrado

Se a gente tenta resolver com eval, ela não funciona, porque x não tem valor, não existe nenhum valor atribuído a x.

Veja o que temos na memoria até agora:

> ls() [1] "a" "espressao" "N2iso"

Então se atribuirmos um valor a x, ele passa a existir

> x<-2 > ls() [1] "a" "espressao" "N2iso" "x"

E agora conseguimos avaliar o que queremos.

> eval(espressao) [1] 4

Então vamos gerar vários valores de N1, para tentarmos observa como vai ficar essa expressão naquele gráfico de state-space que fizemos ali em cima.

N1 <- 0:200
plot(N1, eval(N2iso), type = "l", ylim = c(0, 200), xlim = c(0,200), ylab = expression("N"[2]))

Então criamos valores para N1, e avaliamos no eixo Y, onde a população da especie 2 deveria estar.

02

E ai temos, com crescimento zero da espécie 2, onde deveríamos ver a espécie 1.

A primeira observação aqui, é ver que quando a espécie 1 esta em zero, a espécie 2 esta na sua capacidade suporte. \alpha_{22} é 0.010. Lembre-se que a capacidade suporte é igual a k = {1\over{\alpha_{22}}} e 1 dividido por 0.01 da 100. Veja no grafico que quando a especie 1 esta ausente, esta com sua densidade em zero, temos a capacidade suporte da especie 2, e conforme vamos aumentando a quantidade da especie 1, temos menos indivíduos da especie 2, se a gente pensar no mundo real, em duas especies competindo por espaço por exemplo, conforme uma consome espaço, sobra menos espaço para a outra, assim para uma aumentar a outra tem que diminuir. Então esta fazendo sentido essa isocline.

Bem agora que ja estamos captando o que esta acontecendo, vamos refazer essa curva, agora usando a função curve do R. O que essa função vai fazer é o mesmo procedimento que acabamos de fazer, so que vai evitar todo esse tramite de escrever a expressão, apenas fornecemos uma função e falamos da onde a onde queremos avaliar ela.

curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], axes=FALSE, lty=2,
      ylim=c(0, max(1/a[2,2], 1/a[1,2])), xlim=c(0, max(1/a[1,1], 1/a[2,1])),
      ylab=expression("N"[2]), xlab=expression("N"[1]), pty="s" )

Note que mandado como primeiro argumento uma função com x como variável, mas podemos colocar qualquer letra ou variável ali, so mudar no argumento xname = “x”, mas vamos deixar x, mas poderíamos colocar N1 ali sem problemas, so mudar o xname. Então colocamos isso e até onde queremos avaliar, o intervalo, que é de zero a 1/a[2,1], o limite da competição, e deixamos a figura sem eixos para mudar a notação.

03

Agora veja que mudamos de números, que são valores que mudariam de caso a caso, para a notação relativa as interações.
Como {1\over{\alpha_{22}}} = 100 então deixar {1\over{\alpha_{22}}} é a mesma coisa, é a capacidade suporte. Agora o outro lado da linha ta batendo onde é o valor máximo de N1 dado que N2 é zero.

Para visualizar isso melhor, vamos fazer também a isocline da especie 2, isso é a isocline quando o crescimento da especie 1 é zerado. E vejamos o que acontece, o processo é exatamente a mesma coisa, alias a expressão vai ser a mesma, so mudando quais alphas usamos.

Plotando as duas isocline, temos o seguinte.

04

Olha ai, elas se cruzam, e note que a capacidade suporte tanto da espécie 1 como da especie 2 esta antes do limite de competição, então para esses valores de alpha…

> a [,1] [,2] [1,] 0.010 0.005 [2,] 0.005 0.010

Temos uma tendencia a coexistência, como sabemos para onde é a tendencia de cada especie, dentro do state-space, podemos somar os atratores e eles vão a direção da bolinha, do ponto onde as linhas se cruzam, ou seja esperamos as duas especies tenham uma abundancia, as duas existam.

Por enquanto vamos ficar por aqui, é meio complicado a primera vez que se ve esse negocio de isocline, mas elas são uteis para ver o deveria acontecer de acordo com as interações que observamos, apesar de um modelo simplista e que não leva em conta o espaço, os insights que vamos conseguindo aqui, são muito uteis para entender modelos mais complexos que representam populações. E como vimos aqui, são de modelos simples assim que vem as ideias de como conciliar de forma eficiente a exploração de especies sem destruir elas.
Bem o proximo post sobre esse assunto, a gente testa essa ideia de que sempre vamos cair naquela bolinha ali, onde as linhas se cruzam e depois quais são as possibilidades de equilíbrio para duas populações. Esse post ja esta muito longo, até a próxima.

O script vai estar no repositório recologia do github além de aqui embaixo.

Referencias
M Henry H Stevens 2011 Primer of Ecology With R. Springer
M. Begon, C. R. Townsend, J. L. Harper 2005 – Ecology: From Individuals to Ecosystems, 4th Edition – Wiley

#State-Space
plot(0, 0, type = "p", ylim = c(0, 200), xlim = c(0,200), ylab = expression("N"[2]),xlab=expression("N"[1]))
 
#alphas
a <- matrix(c(0.01, 0.005, 0.005, 0.01), ncol = 2, byrow = TRUE)
a
 
N2iso <- expression(1/a[2, 2] - (a[2, 1]/a[2, 2]) * N1)
 
N2iso
class(N2iso)
 
eval(2^2)
 
espressao<-expression(2^2)
 
eval(espressao)
 
espressao<-expression(x^2)
 
eval(espressao)
 
ls()
 
x<-2
 
ls()
 
eval(espressao)
 
#Isocline da Especie 1
N1 <- 0:200
plot(N1, eval(N2iso), type = "l", ylim = c(0, 200), xlim = c(0,200), ylab = expression("N"[2]))
arrows(x0 = 90, y0 = 150, x1 = 90, y1 = 80, length = 0.1)
arrows(x0 = 75, y0 = 0, x1 = 75, y1 = 50, length = 0.1)
 
 
#Isocline da Especie 1 usando a função curve
curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], axes=FALSE, lty=2,
      ylim=c(0, max(1/a[2,2], 1/a[1,2])), xlim=c(0, max(1/a[1,1], 1/a[2,1])),
      ylab=expression("N"[2]), xlab=expression("N"[1]), pty="s" )
 
axis(1, c(0, 1/a[1,1], 1/a[2,1]),c(0, expression(1/alpha[11]),expression(1/alpha[21])))
axis(2, c(0, 1/a[2,2], 1/a[1,2]),c(0, expression(1/alpha[22]),expression(1/alpha[12])))
arrows(c(75, 175), c(150,150), lty=2,c(75,175), c(80,45), length=0.12)
arrows(c(25, 125), c(5, 5), lty=2,c(25,125), c(55, 20), length=0.12)
 
### Isocline da especie 2
curve(1/a[1,2] - (a[1,1]/a[1,2])*x, 0, 1/a[1,1],axes=FALSE,ylim=c(0, max(1/a[2,2], 1/a[1,2])),
      xlim=c(0, max(1/a[1,1], 1/a[2,1])),ylab=expression("N"[2]), xlab=expression("N"[1]),pty='s' )
axis(1, c(0, 1/a[1,1], 1/a[2,1]),c(0, expression(1/alpha[11]),expression(1/alpha[21])))
axis(2, c(0, 1/a[2,2], 1/a[1,2]),c(0, expression(1/alpha[22]),expression(1/alpha[12])))
arrows(c(150,150), c(75, 175),c(80,45), c(75,175), length=0.12)
arrows( c(5, 5), c(25, 125),c(55, 20), c(25,125), length=0.12)
 
 
### Ambas isoclines no mesmo grafico
 
plot((a[2,2]-a[1,2])/(a[2,2]*a[1,1] - a[1,2]*a[2,1]),(a[1,1]-a[2,1])/(a[2,2]*a[1,1] - a[1,2]*a[2,1]), cex=2,
     pch=1, col=1, axes=FALSE,ylim=c(0, max(1/a[2,2], 1/a[1,2])), xlim=c(0, max(1/a[1,1], 1/a[2,1])),
     ylab=expression("N"[2]), xlab=expression("N"[1]), pty='s')
curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], add=TRUE, lty=2)
curve(1/a[1,2] - (a[1,1]/a[1,2])*x, 0, 1/a[1,1], add=T)
axis(1, c(0, 1/a[1,1], 1/a[2,1]),c(0, expression(1/alpha[11]),expression(1/alpha[21])))
axis(2, c(0, 1/a[2,2], 1/a[1,2]),c(0, expression(1/alpha[22]),expression(1/alpha[12])))
arrows(c(25,25,25), c(25,25,25), c(25,50,50), c(50,25,50),col=1, length=.12, lwd=c(1,1,2), lty=c(2,1,1))
arrows(c(150,150,150), c(150,150,150), c(150,120, 120), c(120,150,120),col=1, length=.12,  lwd=c(1,1,2), lty=c(2,1,1))
arrows(c(10, 10, 10), c(125,125,125), c(10,30,30), c(105,125,105),length=0.12, col=1,  lwd=c(1,1,2), lty=c(2,1,1))
arrows(c(125,125,125), c(10, 10,10), c(125,105,105), c(30,10,30),length=0.12, col=1,  lwd=c(1,1,2), lty=c(2,1,1))

Geradores de números aleatórios.

Estatística Bayesiana, do jeito que temos tentado usar depende estritamente da geração de números aleatórios. Afinal algorítimos como o de Metropolis–Hastings ou Gibbs sampling depende da produção de números aleatórios dentro de certas distribuições específicas.
Por curiosidade a gente pode dar uma olhada como são produzidos esses números aleatórios.

Acho que o jeito mais fácil de produzir número aleatórios é construir alguma maquinha, do tipo que joga moedas para cima, maquininhas de bingo, sei la essas coisas, se elas não forem viciadas, elas podem ajudar a produzir números aleatórios.

Em livros de estatísticas, as vezes no final, entre as tabelas de valores críticos de alpha para estatística T, F ou chi-square, sempre tinha umas tabelinhas de números aleatórios, onde alguém usou algo assim para produzir um set de números aleatórios.

p06

Existem ainda páginas, como essa aqui, que afirma que pode te fornecer muitos números realmente aleatórios provenientes de ruídos (o erro em modelos probabilísticos) de eventos atmosféricos.

Mas usando MCMC, a gente precisa de muitos números aleatórios, e se for enrolado pra conseguir, com certeza a análise não vai pra frente.

Mas nos computadores, nas verdade, a gente usa números pseudos-aleatórios. Que tem esse nome mais ou menos porque eles apresentam todas as características de uma variável aleatória, mas ele é dependente de um início, um número usado de semente para iniciar o processo.
Por isso no R, quando a gente faz uma simulação mas coloca set.seed(um número), se a gente repetir esse número, geramos os mesmo números através das simulações. O que é muito interessante para reproduzir os resultados dos outros, como normalmente quando eu lembro eu sempre deixo nos scripts o seed que eu usei.

Mas so de curiosidade, vamos dar uma olhada em algumas formas de gerar números aleatórios.

Uma das primeiras formas é conhecida como Método do meio do quadrado(Middle-square method).

Basicamente a gente pega um número grande, olha o “miolo dele”, ai eleva ao quadrado esse miolo que vai dar outro número grande, dai olha o “miolo dele” e assim sucessivamente. Então a gente fornece um número grande para iniciar o processo e vai.

storageArray<-vector()
x<-987678
 
for(i in 1:10){
    newseed = (x / 1000) %% 1000000
    storageArray[i] = newseed
    x = newseed * newseed
    }
> storageArray [1] 987.678000 975.507832 951.615530 905.572116 820.060858 672.499811 452.255995 204.535485 [9] 41.834765 1.750148

Mas aparentemente esse método é uma desgraça, ja que existe uma tendencia dele ir parar no zero e la ficar estacionado, ja que zero ao quadrado é zero e o miolo de zero é sempre zero, então se o número começa a diminuir, aparentemente ja era. Então você tem que ter estratégias para evitar isso, o que a gente não se preocupou aqui, tanto que a sequência é ruim e com uma clara tendência.

Mas uma observação legal é que esse método, apesar de inventado a muito tempo atras, foi descrito por John von Neumann formalmente, o cara dessa frase aqui, em 1949. Isso da o que 64 anos? Então antes disso era sem chance de se ver o uso intensivo de MCMC por não se ter muitas opções para gerar uma quantidade grande de número aleatórios de maneira rápida.

Uma outra estratégia é conhecida como Linear congruential generator. Assim como o primeiro caso, aqui a gente começa de algum ponto, o primeiro número, e através de uma relação de recorrência, onde o próximo número depende do anterior, a gente vai gerando mais números.

Basicamente:

 x_{n+1} = a \cdot x_n + c \mod{m}

No R a gente pode fazer um pequeno teste.

m<-100
a<-2
c<-3
 
x<-vector()
x[1]<-1
 
for(i in 1:9) {
    x[i+1]<-(a*x[i]+c) %% m
}
> x [1] 1 19 25 27 61 39 65 7 21 59

10 números, legal. Mas existe um problema, se tentarmos gerar mais, vamos dizer 100 números desse jeito.

Apesar de estarmos gerando números que tem uma distribuição uniforme.

01

Veremos uma auto-correlação forte.

02

Note as barrinhas estourando forte para fora do limite do que seria esperado para ausência de auto-correlação.

O que esta acontecendo é o seguinte, se olharmos a sequencia que criamos.

> x [1] 1 19 25 27 61 39 65 7 21 59 5 87 81 79 45 67 41 99 85 47 1 19 25 27 61 39 65 7 21 59 5 [32] 87 81 79 45 67 41 99 85 47 1 19 25 27 61 39 65 7 21 59 5 87 81 79 45 67 41 99 85 47 1 19 [63] 25 27 61 39 65 7 21 59 5 87 81 79 45 67 41 99 85 47 1 19 25 27 61 39 65 7 21 59 5 87 81 [94] 79 45 67 41 99 85 47

Note que começamos no 1, ai vai para 19 o segundo número, 25, até 47, depois geramos 1, a partir do momento que geramos o 1, o próximo vai ser 19, porque é o numero que sai ao jogarmos na formula la. Existem para esse exemplo sempre 19 valores antes de começar a repetir, então apesar de ter uma distribuição uniforme, para esses parâmetros e uma semente a gente não gera mais que esse intervalo, denominado de período do gerador de números aleatórios.

É claro que se trocarmos os valores dos parâmetros da relação. Principalmente porque que o período máximo seria o valor do modulo, podemos conseguir um período maiores.

Por exemplo, se usarmos:

m<-1000 a<-17 c<-11 Continuamos com uma distribuição uniforme. 04

Mas conseguimos um período maior.

05

Mas a partir daqui, existem limites físicos, que é a quantidade de memoria que você tem, ou mais especificamente o quanto de memoria você vai gastar guardando um número, um número com mais casas decimais precisa de mais e mais bits, e bits e bits.

Daqui a gente chega no método que é comumente usado no R, que é também uma relação de recursão, mas a partir da ideia de um shift register que muda de acordo com o ultimo estado, o método de Linear feedback shift register.

O método que é empregado no R é chamado de Mersenne twister e é um tipo de Linear feedback shift register.

Se a gente tentar gerar números de uma distribuição uniforme como estávamos fazendo aqui, vemos que conseguimos o resultado desejado.

06

E melhor ainda, sem auto-correlação.

07

Isso porque o período, que estávamos na casa dos 20, 30 nos exemplos acima, para o Mersenne Twister é de 2^{19937} - 1.
Esse período é muito, muito grande, então da para gerar muitos números sem nem ter chance de bater no período dele.

E uma vez que se tenha uma forma de gerar número aleatórios com uma distribuição, como a uniforme, através de uma transformação ou outros métodos, que depois dar uma olhada depois, podemos gerar qualquer distribuição que se possa escrever na forma de PDF.

Se você olhar o comando ?RNG(), verá que “The default is “Mersenne-Twister”.”, esse é o gerador usado em tudo que fazemos, mas tem um monte de outros geradores de números aleatórios, mas isso é bem complicado, é melhor não mexer por enquanto, pelo menos eu não mexo.

Bem é isso, matando a curiosidade apenas sobre da onde vem os números pseudo-aleatórios que a gente usa no computador.

#Middle-square method
storageArray<-vector()
x<-987678
 
for(i in 1:10){
    newseed = (x / 1000) %% 1000000
    storageArray[i] = newseed
    x = newseed * newseed
    }
 
storageArray
 
#Linear congruential generator
m<-100
a<-2
c<-3
 
x<-vector()
x[1]<-1
 
for(i in 1:9) {
    x[i+1]<-(a*x[i]+c) %% m
}
 
x
 
for(i in 1:99) {
    x[i+1]<-(a*x[i]+c) %% m
}
 
plot(density(x),frame=F)
acf(ts(x))
x
 
#Modulo Maior
m<-1000
a<-17
c<-11
 
x<-vector()
x[1]<-2
 
for(i in 1:999) {
    x[i+1]<-(a*x[i]+c) %% m
}
 
plot(density(x),frame=F)
acf(ts(x))
x
 
#Mersenne Twister
set.seed(1)
x<-runif(100000)
 
plot(density(x))
acf(x)
 
?RNG()