Insertionsort em R e C++ usando o pacote rcpp

Bem, depois de vermos o algorítimo mais simples que tem para ordenar um vetor, o bubble sort. Agora a gente vai olhar outro um pouco mais legal, esse se chama insertion sort.

A ideia é a seguinte:

Então a gente tem um vetor qualquer, com n elementos.
 C = \{X_1 , X_2 , X_3 , \dots , X_n  \}

E a gente quer que esses elementos fiquem assim:

 Regra = \{x_1 \leq x_2 \leq x_3 \leq \dots \leq x_n  \}

Como a gente faz?
Vamos supor um conjunto aqui.

 C = \{6, 5, 3, 1, 8, 7, 2, 4\}

Ele esta desorganizado porque ele não respeita nossa Regra_1. So olhar o primeiro elemento que é 6.

So que a estratégia aqui muda. O que a gente faz?
Primeiro pega um sub-vetor, um pedaço desse la em cima.

 C = \{6\}

Esse vetor está organizado? Sim, não temos ninguém para comparar. Um vetor de 1 elemento sempre está organizado.

Agora pegamos o segundo elemento, o 5.
Ai fazemos a pergunta onde ele deveria está?
Antes ou depois do 6?
Bem ele é menor que o 6 então tem que estar depois, então movemos o 6 um índice pra frente e dai colocamos o 5.

 C = \{5 ,6\}

Agora pegamos outro elemento, o 3, e repetimos o os passos anteriores.

Perguntamos, onde é a posição do 3?
Ele é menor que o 6? Sim então movemos o 6 um índice para a frente.
Ele é menor que o 5? Sim, então movemos o 5 um índice para a frente, como estamos na primeira casa do vetor, paramos. E ficamos com o seguinte.

 C = \{3 ,5 ,6\}

Notem que o que temos é que o inicio do vetor sempre esta organizado dessa forma, e sempre adicionamos um elemento a um vetor organizado. Mas se o elemento que formos adicionar deve ficar na ultima posição, a gente não vai mover nada de lugar, ou seja, vai estar tudo certo e teremos organizado com menos trocas de posição de elementos. Diferente do bubblesort, que tínhamos que comparar o vetor inteiro para garantir que o ultimo elemento é o maior, mas depois temos que refazer varias comparações para garantir que o segundo maior elemento é o penúltimo, ou seja nunca tinha como escapar de comparar tudo com tudo, aqui tem.

Mas talvez vendo uma animaçãozinha fica mais fácil de entender.

Insertion-sort-example-300px

Apesar de parecer complicado, esse é o jeito que a maioria das pessoas organiza cartas na mão, quando a gente quer deixar em alguma ordem. Veja so.

Ok, mas então como fica esse processo no R?

insertionsort<-function(vetor){
    n<-length(vetor)
 
    for(i in 2:n) {
        aux=vetor[i]
        j=i-1
        while(j>=1 && vetor[j]>aux) {
            vetor[j+1]<-vetor[j]
            j=j-1
            }
        vetor[j+1]=aux
        }
    return(vetor)
    }

Bem simples, um loop usando i diz qual o tamanho do vetor que já está organizado, veja que a gente começa do 2, ja que o índice 1 já está organizado, e ai é só salvar a ultima posição, usamos um while para quanto o auxiliar for menor que o valor no vetor ir empurrando ele para a frente e quando isso não é verdade, colocamos o valor que salvamos no auxiliar na posição desejada.

E esse processo funciona que uma beleza.

> vetor<-sample(100) > vetor [1] 78 2 25 43 10 56 8 64 49 41 32 60 59 5 3 38 69 46 13 77 70 57 4 [24] 26 51 45 53 73 85 55 36 66 83 80 61 79 48 50 62 42 35 23 97 19 27 20 [47] 65 94 54 81 1 82 75 93 86 24 71 63 100 34 39 28 6 7 30 33 31 11 72 [70] 18 99 16 17 96 95 52 76 14 89 22 84 98 91 67 37 87 40 88 21 58 92 29 [93] 12 68 15 90 74 47 9 44 > insertionsort(vetor) [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 [24] 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 [47] 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 [70] 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 [93] 93 94 95 96 97 98 99 100

Podemos fazer isso também usando o Rcpp, mas ai é preciso lembrar, que na linguagem C e C++, diferente do R, os índices começam no 0 e não em 1, então o primeiro elemento de um vetor está na posição vetor[1] no R, mas esta na posição vetor[0] no C++.
Com isso em mente, o código é praticamente o mesmo.

library(Rcpp)
cppFunction("
    NumericVector insertionsortC(NumericVector vetor) {
        int n = vetor.size();
 
        double aux;
        int i , j;
 
        for(i=1;i<n;i++) {
            aux=vetor[i];
            j=i-1;
            while(j>=0 && vetor[j]>aux) {
                vetor[j+1]=vetor[j];
                j=j-1;
                }
            vetor[j+1]=aux;
            }
        return vetor;
        }
")

E funciona legal também.

> insertionsortC(vetor) [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 [24] 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 [47] 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 [70] 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 [93] 93 94 95 96 97 98 99 100

Bem a gente vai ver daqui dois algorítimos de ordenação, se eu continuar fazendo isso, que o que a gente manda para o Rcpp, o código que está entre aspas, pode ser um arquivo texto com o código. E poderíamos mandar mais funções juntas, mas no quicksort a gente vai ver como organizar isso, por enquanto deixa desse jeito que é mais simples.

Mas outra forma ainda de implementar isso é com recursão. A gente ja usou recursão para resolver problemas no Rosalind. Mas basicamente é so chamar a mesma função dentro dela mesma. Um exemplo bem simples é o fatorial.

Fatorial é aquele número que multiplicamos por ele menos 1 até chegar a 1.
Assim,  3!=3 \cdot 2 \cdot 1 e 2 fatorial seria  2!=2 \cdot 1

Então outra forma de escrever isso é  3!=3 \cdot 2!. Fazendo isso até chegar no 1 fatorial, logo no R isso fica assim.

fatorial<-function(n) {
    if(n==1) {
        return(1)
        } else {
            return(n*fatorial(n-1))
            }
    }

Que da.

> fatorial(3) [1] 6

Ou seja, temos um caso base, que é o 1 e vamos fazendo alguma coisa até chegar nessa base.
O insertionsort pode seguir esse raciocínio também da seguinte forma:

cppFunction("
    NumericVector insertionsortRC(NumericVector vetor, int n) {
 
        double aux;
        int i;
 
        if(n>1) {
            insertionsortRC(vetor,n-1);
            aux=vetor[n-1];
            i=n-1;
            while(vetor[i-1]>aux && i>=0 ) {
                vetor[i]=vetor[i-1];
                i--;
                }
            vetor[i]=aux;
            }
 
        return vetor;
        }
    ")

Que funciona também

> insertionsortRC(vetor,length(vetor)) [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 [24] 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 [47] 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 [70] 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 [93] 93 94 95 96 97 98 99 100

Agora eu so conseguir fazer o insertionsort recursivo fazendo uma função que recebe dois argumentos. Eu até perguntei no stack-overflow aqui se tinha como fazer de outra forma, mas não sei se é possivel, mas tudo funciona de qualquer forma.
Que descobri que não sabia fazer foi fazer o insertion sort de forma recursiva no R. Tem uma pergunta aberta aqui na lista brasilera de r r-br. Se alguém puder me explicar como faz, eu fico grato.

E bem eu acho que é isso, agora temos dois algorítimos de ordenação. Vamos adicionar mais alguns a essa lista e depois será divertido tentar analisar o quão bom cada um é, e quanto…

Bem o script está aqui embaixo além do repositório recologia no github. Até a próxima.

vetor<-sample(100)
vetor
 
insertionsort<-function(vetor){
    n<-length(vetor)
 
    for(i in 2:n) {
        aux=vetor[i]
        j=i-1
        while(j>=1 && vetor[j]>aux) {
            vetor[j+1]<-vetor[j]
            j=j-1
            }
        vetor[j+1]=aux
        }
    return(vetor)
    }
 
insertionsort(vetor)
 
library(Rcpp)
 
cppFunction("
    NumericVector insertionsortC(NumericVector vetor) {
        int n = vetor.size();
 
        double aux;
        int i , j;
 
        for(i=1;i<n;i++) {
            aux=vetor[i];
            j=i-1;
            while(j>=0 && vetor[j]>aux) {
                vetor[j+1]=vetor[j];
                j=j-1;
                }
            vetor[j+1]=aux;
            }
        return vetor;
        }
")
 
 
 
insertionsortC(vetor)
 
cppFunction("
    NumericVector insertionsortRC(NumericVector vetor, int n) {
 
        double aux;
        int i;
 
        if(n>1) {
            insertionsortRC(vetor,n-1);
            aux=vetor[n-1];
            i=n-1;
            while(vetor[i-1]>aux && i>=0 ) {
                vetor[i]=vetor[i-1];
                i--;
                }
            vetor[i]=aux;
            }
 
        return vetor;
        }
    ")
 
vetor
insertionsortRC(vetor,length(vetor))
 
fatorial<-function(n) {
    if(n==1) {
        return(1)
        } else {
            return(n*fatorial(n-1))
            }
    }
 
fatorial(3)

Estimando colonização e extinção de uma espécie usando OpenBugs e Unmarked

Bem, nesse post aqui, a muito tempo atrás, vimos como estimar parâmetros do modelo de metapopulações. Mais especificamente dentro do nosso interesse, colonização e extinção.
Para relembrar um pouco sobre metapopulações, podemos dar uma olhada nesse post aqui, para ter uma ideia do modelo.

A ideia aqui é so começar a olhar, como estimar, os mesmos parâmetros, como no unmarked, para colonização e extinção. O mesmo que faz a função colext() desse pacote.Mas levando em conta também a detecção.

Pra falar a verdade, vamos primeiro olhar aqui como fica o modelo no na linguagem bugs, e depois em outros posts olhamos melhor pedacinho por pedacinho, mas vamos fazer ele e deixar funcionando para começar a realizar alguns testes.

Bem primeiro vamos olhar como são os dados, então temos a ocupação ao longo do tempo. Vamos simular alguns dados e obtemos o seguinte.

01

Veja que a ocupação que observamos é menor que a ocupação real do ambiente, a realizada realmente. Isso porque a detecção não é 100%, ou seja, se pensarmos em aranhas em arbustos, não é porque eu não vi a aranha que é 100% garantido que o arbusto esta vazio, e pior que isso, pode ser um processo estocástico, quero dizer, não é garantido que sempre eu “não consigo ver igual” a aranha, numa coleta num dia de vento pode ser pior que uma coleta num dia sem vento, no dia de vento a aranha se escondeu e ai eu posso achar que existem menos aranhas no ambiente que realmente tem. E assim vai.
Então ali a linha azul é como se fosse a ocupação (pode chamar de prevalência também) das aranhas se eu visse todas. Mas temos a linha vermelha pontilhada que é a chance de eu ver essas aranhas, ou seja quanto menor ela esta, menos aranha eu vejo, e no final o que eu realmente vejo no campo é a linha preta. Veja bem que em dados reais, tudo que a gente ve é a linha preta, as outras você nunca vê, a gente vê aqui porque simulou os dados, mas elas são ocultas, e é isso que tentamos estimar. Todas essas probabilidades.

E por curiosidade como são esses dados?

> data$y[,,1] [,1] [,2] [,3] [1,] 0 1 0 [2,] 1 0 0 [3,] 0 0 0 [4,] 0 0 0 . . . [27,] 0 1 1 [28,] 0 1 0 [29,] 0 1 0 [30,] 0 0 0

É tipo uma matriz de 3 dimensões, onde as linhas são os meus arbustos, eu tenho 30 nesse caso, eu ia três vezes ao campo todo ano, olha se minha arainha estava la ou não, e por vários anos, anos ficam como a terceira dimensão dessa matriz, ou array seria um termo melhor, da para organizar de outras formas os dados, mas assim fica mais simples de trabalhar com os índices, que ficam [amostras, repetições , anos]

Certo então vamos definir nosso modelo na linguagem bugs.

model {
 
        # Priors
        psi1 ~ dunif(0, 1)
        for (k in 1:(nyear-1)){
            phi[k] ~ dunif(0, 1)
            gamma[k] ~ dunif(0, 1)
            p[k] ~ dunif(0, 1)
        }
        p[nyear] ~ dunif(0, 1)
 
        # Modelo ecologico, definindo estado real
        for (i in 1:nsite){
            z[i,1] ~ dbern(psi1)
            for (k in 2:nyear){
                muZ[i,k]<- z[i,k-1]*phi[k-1] + (1-z[i,k-1])*gamma[k-1]
                z[i,k] ~ dbern(muZ[i,k])
            } #k
        } #i
 
        # Modelo referente as observações
        for (i in 1:nsite){
            for (j in 1:nrep){
                for (k in 1:nyear){
                    muy[i,j,k] <- z[i,k]*p[k]
                    y[i,j,k] ~ dbern(muy[i,j,k])
                } #k
            } #j
        } #i
 
# Parametros derivados
        psi[1] <- psi1
        n.occ[1]<-sum(z[1:nsite,1])
        for (k in 2:nyear){
            psi[k] <- psi[k-1]*phi[k-1] + (1-psi[k-1])*gamma[k-1]
            n.occ[k] <- sum(z[1:nsite,k])
            growthr[k-1] <- psi[k]/psi[k-1]
            turnover[k-1] <- (1 - psi[k-1]) * gamma[k-1]/psi[k]
        }
    }

E a parte mais importante é aqui, onde realmente as coisas estão acontecendo

 
        # Modelo ecologico, definindo estado real
        for (i in 1:nsite){
            z[i,1] ~ dbern(psi1)
            for (k in 2:nyear){
                muZ[i,k]<- z[i,k-1]*phi[k-1] + (1-z[i,k-1])*gamma[k-1]
                z[i,k] ~ dbern(muZ[i,k])
            } #k
        } #i

Veja que aqui a gente ta trabalho com o locais e as repetições dentro de anos.
Usamos distribuição binomial, ou seja nada mais é que um tipo de regressão logístico.
E se a gente lembrar do modelo do Levins:

{dN\over{dt}}={c \dot N \dot (1-N)-e\dot N}

Temos alguns comentários de muito tempo atras aqui.

Veja que isso é o que estamos usando nessa linha aqui:

muZ[i,k]<- z[i,k-1]*phi[k-1] + (1-z[i,k-1])*gamma[k-1]

So trocar o gamma por c, o z por N e o phi por e.

Mas essa quantidade Z pode estar subestimada, pode estar contata a menos, e isso a gente corrige aqui.

        # Modelo referente as observações
        for (i in 1:nsite){
            for (j in 1:nrep){
                for (k in 1:nyear){
                    muy[i,j,k] <- z[i,k]*p[k]
                    y[i,j,k] ~ dbern(muy[i,j,k])
                } #k
            } #j
        } #i

Juntando todos os dados no Bugs, temos o resultado:

> print(out1, dig = 2) Inference for Bugs model at "Dynocc.txt", Current: 3 chains, each with 2500 iterations (first 500 discarded), n.thin = 4 Cumulative: n.sims = 6000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff psi[1] 0.89 0.08 0.69 0.84 0.90 0.95 0.99 1.00 6000 psi[2] 0.83 0.07 0.68 0.78 0.83 0.87 0.94 1.00 2600 ... psi[9] 0.58 0.12 0.35 0.50 0.58 0.67 0.83 1.00 720 psi[10] 0.55 0.11 0.33 0.47 0.54 0.62 0.78 1.00 2400 phi[1] 0.87 0.07 0.71 0.82 0.88 0.93 0.99 1.00 6000 ... phi[9] 0.68 0.15 0.37 0.57 0.68 0.78 0.96 1.00 1200 gamma[1] 0.45 0.27 0.02 0.22 0.43 0.66 0.95 1.00 730 ... gamma[9] 0.35 0.19 0.03 0.21 0.34 0.47 0.73 1.00 1900 p[1] 0.30 0.05 0.20 0.26 0.29 0.33 0.41 1.00 6000 ... p[9] 0.40 0.09 0.24 0.34 0.40 0.47 0.60 1.00 6000 p[10] 0.44 0.09 0.27 0.37 0.43 0.50 0.62 1.00 5700 n.occ[1] 27.45 1.92 23.00 26.00 28.00 29.00 30.00 1.00 4200 ... n.occ[10] 16.58 2.71 13.00 14.00 16.00 18.00 23.00 1.00 1200 growthr[1] 0.94 0.11 0.75 0.87 0.93 0.99 1.17 1.00 2800 ... growthr[9] 0.97 0.29 0.53 0.77 0.93 1.13 1.67 1.00 680 turnover[1] 0.06 0.07 0.00 0.01 0.04 0.09 0.25 1.00 2200 ... turnover[9] 0.28 0.16 0.02 0.15 0.27 0.39 0.61 1.00 1100 deviance 771.25 20.97 732.50 756.60 770.60 785.02 813.20 1.00 1100 For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC info (using the rule, pD = var(deviance)/2) pD = 219.5 and DIC = 990.7 DIC is an estimate of expected predictive error (lower deviance is better).

Cada parâmetro por ano, lembrando que quando vemos por exemplo 9 parâmetros gamma, é porque ele é usado na transição de um tempo para o próximo, por isso temos 10 anos de dados simulados, mas 9 parâmetros estimados, não temos como estimar o gamma para a primeira observação.

Podemos testar a ocupação estimada com as dos dados originais e ver se obtivemos algum sucesso.

02

Bem o vermelhinho é o que estimamos, o pretinho são os dados que tinhamos, e o azul são os valores que fomos tirar la da simulação, o real, utilizado na simulação, com exceção do primeiro ano, o intervalo de confiança cobriu os valores reais.

É, parece que esse tipo de modelagem funciona, mas é preciso uma quantidade grande de dados, e aqui melhorar os dados é determinado por mais replicas, mais locais e mais tempo. Somente um desses três não melhora muito a parametrização.
Mas agora a gente ja sabe mais ou menos o que a função colext esta tentando fazer. Isso ja ajuda a pensar um pouco.

Agora é continuar estudando. Esse script, esta presente aqui, junto com outros varios exemplos legais. Existem também esse site referente ao pacote unmarked e um google group de discussão aqui. Esses são os pontos primeiros para buscar informação e exemplos para entender melhor tudo isso.

#install.packages("unmarked")
library(unmarked)
 
data.fn<-function(R=30,J=3,K=10,psi1=0.4,range.p=c(0.2, 0.4),range.phi=c(0.6,0.8),range.gamma=c(0, 0.1)) {
 
    # Iniciando os vetores
    site <- 1:R # Locais
    year <- 1:K # Anos
    psi <- rep(NA, K) # Probabilidade de ocupação
    muZ <- z <- array(dim = c(R, K)) # Ocorrencia esperada e realizada
    y <- array(NA, dim = c(R, J, K)) # Historico de deteção
 
    # Determinando a ocupação inicial e seus parametros
    psi[1] <- psi1 # Probabilidade de ocupação inicial
    p <- runif(n = K, min = range.p[1], max = range.p[2])
    phi <- runif(n = K-1, min = range.phi[1], max = range.phi[2])
    gamma <- runif(n = K-1, min = range.gamma[1], max = range.gamma[2])
 
    # Estados latentes de ocupação
    # Primeiro ano
    z[,1] <- rbinom(R, 1, psi[1])
    # Estado de ocupação inicial
    # Proximos anos
    for(i in 1:R){
        # Loop sobre locais
        for(k in 2:K){
            # Loop sobre anos
            muZ[k] <- z[i, k-1]*phi[k-1] + (1-z[i, k-1])*gamma[k-1] # Prob for occ.
            z[i,k] <- rbinom(1, 1, muZ[k])
            }
        }
    # Grafico da ocupação realizada
    plot(year,apply(z,2,mean),type="l",xlab="Ano",ylab="Ocupação",col="blue",xlim=c(0,K+1),
         ylim=c(0,1),lwd=2,lty=1,frame.plot=FALSE,las = 1)
    lines(year,p,type="l",col="red",lwd=2,lty=2)
    # Gerando dados de detecção/não detecção
    for(i in 1:R){
        for(k in 1:K){
            prob <- z[i,k] * p[k]
            for(j in 1:J){
                y[i,j,k] <- rbinom(1, 1, prob)
                }
            }
        }
    # Computando ocupação anual
    for (k in 2:K){
        psi[k] <- psi[k-1]*phi[k-1] + (1-psi[k-1])*gamma[k-1]
        }
    # Plotando ocupação aparente
    psi.app <- apply(apply(y,c(1,3), max), 2, mean)
    lines(year, psi.app, type= "l", col = "black", lwd = 2)
    legend("topright",lty=c(1,2,1),lwd=2,col=c("blue","red","black"),
           legend = c("Ocupação Verdadeira","Detecçao"," Ocupação Observada"))
 
 
# Retornando dados
    return(list(R=R,J=J,K=K,psi=psi,psi.app=psi.app,z=z,phi=phi,gamma=gamma,p=p,y=y))
    }
 
set.seed(123)
jpeg("01.jpg")
data<-data.fn(R=30,J=3,K=10,psi1=0.6,range.p=c(0.1,0.9),range.phi=c(0.7,0.9),range.gamma=c(0.1,0.5))
str(data)
dev.off()
 
yy <- matrix(data$y, data$R, data$J* data$K)
str(yy)
head(yy)
 
year <- matrix(c('01','02','03','04','05','06','07','08','09','10'),
nrow(yy), data$K, byrow=TRUE)
head(year)
 
simUMF <- unmarkedMultFrame(y = yy,yearlySiteCovs = list(year = year),numPrimary=data$K)
summary(simUMF)
 
fm1 <- colext(psiformula = ~1,      #Ocupação no primeiro ano
              gammaformula = ~ 1,   #Colonização
              epsilonformula = ~ 1, #Extinção
              pformula = ~ 1,       #Detecção
              data = simUMF)
 
summary(fm1)
system.time(summary(fm2 <- colext(psiformula=~1, gammaformula=~year-1,epsilonformula=~1, pformula =~1, data = simUMF,
                                  control=list(trace=TRUE, REPORT=1, maxit = 500), se = TRUE)))
 
system.time(summary(fm3 <- colext(psiformula=~1, gammaformula=~1,epsilonformula=~year-1, pformula =~1, data = simUMF,
                                  control=list(trace=TRUE, REPORT=1, maxit = 500), se = TRUE)))
 
system.time(summary(fm4 <- colext(psiformula=~1, gammaformula=~1,epsilonformula=~1, pformula =~ year-1, data = simUMF,
                                  control=list(trace=TRUE, REPORT=1, maxit = 500), se = TRUE)))
 
system.time(summary(fm5 <- colext(psiformula=~1, gammaformula=~year-1,epsilonformula=~year-1, pformula =~1,
                                  data = simUMF,control=list(trace=TRUE, REPORT=1, maxit = 500), se = TRUE)))
 
system.time(summary(fm6 <- colext(psiformula=~1, gammaformula=~year-1,epsilonformula=~1, pformula =~year-1,
                                  data = simUMF,control=list(trace=TRUE, REPORT=1, maxit = 500), se = TRUE)))
 
system.time(summary(fm7 <- colext(psiformula=~1, gammaformula=~1,epsilonformula=~year-1, pformula =~year-1,
                                  data = simUMF,control=list(trace=TRUE, REPORT=1, maxit = 500), se = TRUE)))
 
system.time(summary(fm8 <- colext(psiformula=~1, gammaformula=~year-1,epsilonformula=~year-1, pformula =~year-1,
                                  data = simUMF,control=list(trace=TRUE, REPORT=1, maxit = 500), se = TRUE)))
 
 
modelos <- fitList(
'psi(.)gam(.)eps(.)p(.)' = fm1,
'psi(.)gam(Y)eps(.)p(.)' = fm2,
'psi(.)gam(.)eps(Y)p(.)' = fm3,
'psi(.)gam(.)eps(.)p(Y)' = fm4,
'psi(.)gam(Y)eps(Y)p(.)' = fm5,
'psi(.)gam(Y)eps(.)p(Y)' = fm6,
'psi(.)gam(.)eps(Y)p(Y)' = fm7,
'psi(.)gam(Y)eps(Y)p(Y)' = fm8)
modSel(modelos)
 
# Bundle data
bugs.data <- list(y = data$y, nsite = dim(data$y)[1], nrep = dim(data$y)[2],nyear = dim(data$y)[3])
 
data$y[,,1]
 
# Specify model in BUGS language
sink("Dynocc.txt")
cat("
    model {
 
        # Priors
        psi1 ~ dunif(0, 1)
        for (k in 1:(nyear-1)){
            phi[k] ~ dunif(0, 1)
            gamma[k] ~ dunif(0, 1)
            p[k] ~ dunif(0, 1)
        }
        p[nyear] ~ dunif(0, 1)
 
        # Modelo ecologico, definindo estado real
        for (i in 1:nsite){
            z[i,1] ~ dbern(psi1)
            for (k in 2:nyear){
                muZ[i,k]<- z[i,k-1]*phi[k-1] + (1-z[i,k-1])*gamma[k-1]
                z[i,k] ~ dbern(muZ[i,k])
            } #k
        } #i
 
        # Modelo referente as observações
        for (i in 1:nsite){
            for (j in 1:nrep){
                for (k in 1:nyear){
                    muy[i,j,k] <- z[i,k]*p[k]
                    y[i,j,k] ~ dbern(muy[i,j,k])
                } #k
            } #j
        } #i
 
# Parametros derivados
        psi[1] <- psi1
        n.occ[1]<-sum(z[1:nsite,1])
        for (k in 2:nyear){
            psi[k] <- psi[k-1]*phi[k-1] + (1-psi[k-1])*gamma[k-1]
            n.occ[k] <- sum(z[1:nsite,k])
            growthr[k-1] <- psi[k]/psi[k-1]
            turnover[k-1] <- (1 - psi[k-1]) * gamma[k-1]/psi[k]
        }
    }
    ",fill = TRUE)
sink()
 
# Valores iniciais da cadeia
zst <- apply(data$y, c(1, 3), max) # Observações ocorridas para iniciar Z
inits <- function(){ list(z = zst)}
# Parametros que queremos salvar
params <- c("psi", "phi", "gamma", "p", "n.occ", "growthr", "turnover")
# Definições do MCMC
ni <- 2500
nt <- 4
nb <- 500
nc <- 3
# Chamando Openbus (Pode demorar um pouco)
library(R2OpenBUGS)
out1 <- bugs(bugs.data,inits,params,"Dynocc.txt",n.chains=nc,n.thin=nt,n.iter=ni,n.burnin=nb)
 
# Sumario dos resultados
print(out1, dig = 2)
 
str(data)
jpeg("02.jpg")
plot(1:data$K,data$psi,type="l",xlab="Ano",ylab="Probabilidade de ocupação",col="blue",xlim=c(0,data$K+1),ylim=c(0,1),
     lwd=2,lty=1, frame.plot=FALSE,las = 1)
lines(1:data$K, data$psi.app, type = "l", col = "black", lwd = 2)
points(1:data$K, out1$mean$psi, type = "l", col = "red", lwd = 2)
segments(1:data$K, out1$summary[1:data$K,3], 1:data$K,out1$summary[1:data$K,7], col = "red", lwd= 1)
dev.off()

Bubblesort em R e C++ usando o pacote rcpp

Opa, a algum tempo atras eu estava olhando como usar funções escritas em C do R, aqui.

Mas como a gente viu la, é um tramite danado tudo isso, desde a parte de escrever o programa, compilar, etc.

Mas o mais complicado, é conseguir fazer a linguagem C entender como está organizado os dados que vem do R, e vice versa. Mas quem viu aquele post, também viu uma alma caridosa que postou no comentario sobre como fazer isso usando o pacote rcpp.

Bem, dando uma olhada como fazer as coisas, parece que nada é tão impraticável assim. Então como uma primeiro olhada, vamos olha um algorítimo chamado bubblesort feito direto no R e usando o rcpp.

Esse algorítimo entra numa categoria que acho que todo livro de algorítimos básico aborda, e realiza uma tarefa básica essencial, organiza um vetor em alguma ordem, aqui a ordem numérica.

Então a gente tem um vetor qualquer, com n elementos.
 C = \{X_1 , X_2 , X_3 , \dots , X_n  \}

E a gente quer que esses elementos fiquem assim:

 C = \{x_1 \leq x_2 \leq x_3 \leq \dots \leq x_n  \}

Ou seja, uma ordem crescente, eles vão ficar organizados dentro do vetor.

Como a gente consegue isso? A gente vai comparando o primeiro com o segundo, se o segundo for menor a gente troca os dois de lugares, senão a gente mantem eles nos mesmo lugares, depois comparamos o segundo com o terceiro, se for diferente trocamos de lugares esses dois, senão mantemos nos mesmo lugares, fazemos isso até chegar na comparação do ultimo com o penúltimo. Depois dessa comparação e troca se necessário, a gente pode fazer uma afirmação quanto a esse vetor, podemos dizer que o ultimo elemento é o maior elemento do vetor, ja que ele foi empurrado la pro fundo, a partir daqui da para pensar que se continuamos fazendo isso agora até o penúltimo, depois até o antepenúltimo, ao chegar na comparação que vai até o segundo termo, temos certeza que todo o vetor está organizado depois disso.

No fim do post tem uma animação para entender melhor, mas vamos ver como fica o código na linguagem R.

bubblesortR<-function(vetor) {
    aux<-numeric()
    for(i in 1:(length(vetor)-1)) {
        for(j in 1:(length(vetor)-i)) {
            if(vetor[j]>vetor[j+1]) {
                aux<-vetor[j]
                vetor[j]<-vetor[j+1]
                vetor[j+1]<-aux
                    }
                }
            }
        return(vetor)
 
    }

Então a gente recebe o vetor, ai veja que o primeiro loop, que usa i, ele vai até o fim do vetor, depois ele vem diminuindo, isso é o que diminuímos no segundo loop que usa j, e dentro desses dois loops tudo que a gente faz é comparar os pares de valores, e se j for menor que o próximo, a gente troca, senão não fazemos nada.

Mas mais fácil é testar a função, então se temos uma vetor desorganizado, com 100 números assim:

> vetor<-sample(100) > vetor [1] 56 89 96 69 37 65 45 74 9 55 84 7 31 92 14 47 41 21 28 11 42 39 19 [24] 64 76 38 22 24 3 85 63 20 53 88 70 59 32 83 82 18 77 57 91 95 73 81 [47] 27 33 17 10 34 100 71 16 23 29 87 35 13 78 43 46 36 52 58 90 60 40 15 [70] 1 4 68 61 94 50 48 67 86 97 66 72 5 6 62 26 99 79 2 54 8 51 98 [93] 93 75 25 80 49 12 30 44

A função vai organizar eles e deixar assim:

> bubblesortR(vetor) [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 [24] 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 [47] 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 [70] 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 [93] 93 94 95 96 97 98 99 100

Legal né, acontece que para fazer isso então, quanto mais números tivermos, mais operações temos que fazer e mais isso demora.

A gente pode na verdade medir o tempo no R usando a função system.time, fazendo algo como um experimento, sorteamos um vetor de números inteiros do tamanho que queremos e ordenamos ele, marcando quanto tempo demora, e a gente pode fazer isso para diferentes tamanhos e ver o que acontece, nesse caso a gente ve o seguinte.

> system.time(bubblesortR(sample(100))) usuário sistema decorrido 0.048 0.000 0.045 > system.time(bubblesortR(sample(1000))) usuário sistema decorrido 4.684 0.004 4.701 > system.time(bubblesortR(sample(5000))) usuário sistema decorrido 116.784 0.172 117.215

Mais ou menos, conforme aumentamos o número de elementos, aumentamos o tempo de execução, mas não linearmente, porque olha so, ao aumentarmos de 100 para 1000, aumentamos uma ordem de magnitude o tempo ali, mas quanto aumentamos em 5 vezes o número de elementos, aumentamos muito mais que 5 vezes o tempo.
Legal né, mas esse não é o foco ainda. Vamos olhar como fazer a mesma coisa no rcpp, usando código na linguagem c++.

library(Rcpp)
 
cppFunction("
    NumericVector bubblesortC(NumericVector vetor) {
        int n = vetor.size();
 
        double aux;
        int i , j;
 
        for(i=0;i<n-1;i++) {
            for(j=0;j<n-i-1;j++) {
                if(vetor[j]>vetor[j+1]) {
                    aux=vetor[j];
                    vetor[j]=vetor[j+1];
                    vetor[j+1]=aux;
                    }
                }
            }
        return vetor;
        }
")

Bem, o código é bem parecido, so a sintaxe do for que é diferente um pouco e temos que usar uma função para obter o tamanho do vetor, igual quando usamos length no R.
Se usarmos ele no mesmo vetor que criamos antes teremos o mesmo efeito

> bubblesortC(vetor) [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 [24] 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 [47] 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 [70] 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 [93] 93 94 95 96 97 98 99 100

Certo, mas vamos olhar os tempos de execução para matar a curiosidade.

> system.time(bubblesortC(sample(100))) usuário sistema decorrido 0 0 0 > system.time(bubblesortC(sample(1000))) usuário sistema decorrido 0.004 0.000 0.004 > system.time(bubblesortC(sample(5000))) usuário sistema decorrido 0.084 0.000 0.087

São bem menores, isso porque é outra linguagem que é mais rápida que o R, mas como nada na vida é de graça, a gente paga o preço de tudo ser mais complexo de fazer. Complexo no sentido que, pelo menos eu, acho mais fácil escrever as coisas direto no R.

Mas quanto o tempo começar a se mostrar como um empecilho, vale a pensa pelo menos ter idea de soluções possíveis como essa. Existe um tutorial bem bacana aqui sobre como usar o Rcpp, além de um livro do próprio autor que está disponível, e foi publicado esse ano.

Outra coisa é que esse algorítimo é bem lento, mas é interessante pensar como ele funciona, como a descrição sempre é complexa para mim, eu fiz essa animação bem comum em todo post, livro ou sei la quem que está tentando entender o que esta acontecendo numa ordenação.

bubblesort

Bem o script esta aqui embaixo e no repositório recologia, e espero continuar lendo mais um pouco sobre Rcpp e colocando mais coisas aqui, mas aqui ja deu para ter uma ideia, e ver que é bem legal, ja que poupa todo o tramite como visto ao tentarmos usar a função .C
Além de quebrar o recesso sem postar nada, muitas provas ultimamente 🙁

library(Rcpp)
 
cppFunction("
    NumericVector bubblesortC(NumericVector vetor) {
        int n = vetor.size();
 
        double aux;
        int i , j;
 
        for(i=0;i<n-1;i++) {
            for(j=0;j<n-i-1;j++) {
                if(vetor[j]>vetor[j+1]) {
                    aux=vetor[j];
                    vetor[j]=vetor[j+1];
                    vetor[j+1]=aux;
                    }
                }
            }
        return vetor;
        }
")
 
bubblesortR<-function(vetor) {
    aux<-numeric()
    for(i in 1:(length(vetor)-1)) {
        for(j in 1:(length(vetor)-i)) {
            if(vetor[j]>vetor[j+1]) {
                aux<-vetor[j]
                vetor[j]<-vetor[j+1]
                vetor[j+1]<-aux
                    }
                }
            }
        return(vetor)
 
    }
 
vetor<-sample(100)
 
vetor
bubblesortR(vetor)
bubblesortC(vetor)
 
system.time(bubblesortR(sample(100)))
system.time(bubblesortR(sample(1000)))
system.time(bubblesortR(sample(5000)))
 
system.time(bubblesortC(sample(100)))
system.time(bubblesortC(sample(1000)))
system.time(bubblesortC(sample(5000)))
 
 
vetor<-sample(8)
 
#Animação
num<-1
aux<-vector()
resposta<-logical()
 
for(i in 1:(length(vetor)-1)) {
    for(j in 1:(length(vetor)-i)) {
        num<-num+1
        jpeg(sprintf("bubble%04d.jpg",num), width = 350, height = 350)
        cores<-rep("grey",8)
        cores[c(j,j+1)]<-"blue"
        resposta<-vetor[j]>vetor[j+1]
        barplot(vetor,names.arg = vetor,yaxt="n",col=cores,main=paste(vetor[j],"maior que", vetor[j+1],"?"))
        dev.off()
 
         if(vetor[j]>vetor[j+1]) {
             aux<-vetor[j]
             vetor[j]<-vetor[j+1]
             vetor[j+1]<-aux
             }
 
        num<-num+1
        jpeg(sprintf("bubble%04d.jpg",num), width = 350, height = 350)
        if(resposta==TRUE) {
            cores<-rep("grey",8)
            cores[c(j,j+1)]<-"red"
            barplot(vetor,names.arg = vetor,yaxt="n",col=cores,main=paste("Trocamos",vetor[j],"e", vetor[j+1],"de lugar"))
            }else{
                cores<-rep("grey",8)
                cores[c(j,j+1)]<-"blue"
                barplot(vetor,names.arg = vetor,yaxt="n",col=cores,main="Não trocamos nada")
                }
        dev.off()
        }
    }
 
system("convert -delay 90 bu*.jpg  bubblesort.gif")
system("rm bu*.jpg")

Outros casos de estabilidade em populações.

Certo, a gente começou a olhar aqui, como modelar duas espécies usando a ideia do crescimento logístico de Lotka-Volterra, depois aqui como avaliar o resultado da competição entre duas desenhando isolinhas, olhamos aqui um pouco melhor como elas funcionam, mas sempre usamos o mesmo exemplo, os mesmos valores para as interações, que resultavam em coexistência. Mas agora vamos olhar algumas das outras possibilidades.

Como vimos nas análises anteriores, o fator determinante parece ser como são as interações entre as espécies, que no caso da coexistência entre duas espécies tinhamos:

 (\alpha_{12} < \alpha_{22}) \land  (\alpha_{21} < \alpha_{11}) [/latex]  Caso alguém não conheça esse simbolo, "[latex]\land[/latex]", esse é o "e" da matemática, vem da <a href="http://pt.wikipedia.org/wiki/%C3%81lgebra_booleana">álgebra booleana</a>, que merece um post aqui cedo ou tarde, mas estamos querendo dizer aqui que isso [latex] (\alpha_{12} < \alpha_{22})[/latex] e [latex] (\alpha_{21} < \alpha_{11}) [/latex] tem que acontecer ao mesmo tempo, não da para tirar uma conclusão olhando uma das comparações sozinhas, temos que olhar as duas e então chegamos a uma conclusão.  Veja que isso pode parecer complexo, mas traduzindo, temos que para as espécies coexistirem estavelmente, os efeitos delas, nelas mesmas precisam ser maiores que os efeitos na outra espécie.  E o resultado é esse gráfico que já vimos várias vezes.  <a href="http://recologia.com.br/wp-content/uploads/2013/10/01.jpg"><img src="http://recologia.com.br/wp-content/uploads/2013/10/01.jpg" alt="01" width="480" height="480" class="aligncenter size-full wp-image-1694" srcset="http://recologia.com.br/wp-content/uploads/2013/10/01.jpg 480w, http://recologia.com.br/wp-content/uploads/2013/10/01-150x150.jpg 150w, http://recologia.com.br/wp-content/uploads/2013/10/01-300x300.jpg 300w" sizes="(max-width: 480px) 100vw, 480px" /></a>   Mas então como são os outros casos?   Quando a gente ver…  [latex] (\alpha_{12} < \alpha_{22}) \land  (\alpha_{21} > \alpha_{11})

A espécie 1 pode invadir quando rara, mas a espécie 2 não pode.
Isso leva a exclusão competitiva pela espécie 1, sempre. A espécie 1 vence sempre. Isso também é citado como equilíbrio de borda, porque uma das espécies ocorre somente na borda do gráfico, do state-space.
A gente não ta falando que a espécie 1 ganha sempre? Mas lembre-se que se ela for zero, só quando algum indivíduo dela invadir que a população vai começar a crescer, ou seja, só na borda de N1=0 que a espécie 2 está na jogada. O equilíbrio onde N1>0 é chamado de equilíbrio interno.

No gráfico vemos:

02

Agora quando a gente ve…

 (\alpha_{12} > \alpha_{22}) \land  (\alpha_{21} < \alpha_{11}) [/latex]  A espécie 1 não consegue invadir o sistema, mas a espécie 2 consegue. O que gera o inverso do caso anterior, aqui a espécie 2 vence, a espécie 2 exclui a 1 por competição. Mas lembre-se que essas equações e suas relações são derivações de simples equações da reta e tem que ser pensado a limitações delas, por exemplo não usar pontos negativos. Mas de qualquer forma temos podemos extrair hipóteses testáveis daqui.  O gráfico dela fica assim, o inverso do anterior...  <a href="http://recologia.com.br/wp-content/uploads/2013/10/03.jpg"><img src="http://recologia.com.br/wp-content/uploads/2013/10/03.jpg" alt="03" width="480" height="480" class="aligncenter size-full wp-image-1697" srcset="http://recologia.com.br/wp-content/uploads/2013/10/03.jpg 480w, http://recologia.com.br/wp-content/uploads/2013/10/03-150x150.jpg 150w, http://recologia.com.br/wp-content/uploads/2013/10/03-300x300.jpg 300w" sizes="(max-width: 480px) 100vw, 480px" /></a>   Agora partimos pro caso legal.  [latex] (\alpha_{12} > \alpha_{22}) \land  (\alpha_{21} > \alpha_{11})

Porque legal? Porque aqui a gente cria um equilíbrio interno instável, a exclusão de alguma das espécies vai ocorrer, mas qualquer uma pode vencer. Isso normalmente é chamado de controle do fundador, porque quem chega primeiro ta na vantagem, a decisão depende em parte das abundância iniciais que observamos. Aqui o equilíbrio cria uma saddle ou sela no “state-space” ou espaço de estado, espaço de possibilidades. Podemos dizer que de algumas direções, a sela funciona como um atrator, como na coexistência, mas de outro lado, ela age como um repulsor, como existe nos casos de exclusão competitiva. Mas para ver isso a gente tem que dar um passo além e usar matrizes jacobianas.

Mas antes vamos olhar o gráfico como ficou.

04

Veja para onde os vetores estão apontado, em cada parte do gráfico existe uma direção diferente, isso que estamos falando quando falamos de cela, é saddle em inglês, espero não ter traduzido errado.

Então podemos olhar agora nossos casos principais, coexistência, exclusão competitiva e controle do fundador, onde não temos um vencedor definido, depende das condições que começamos a observar o sistema.

05

Agora vale lembrar que usamos valores para as figuras ficarem “bonitinhas”. Mas em qualquer valor veremos essas coisas acontecendo. Saia inventado valores para as interações das espécies e veja o que acontece.

06

Bem pode parecer um bocado de trabalho, para chegar nessas conclusões não? Mas pense que no final vale a pena, agora a gente finalmente entende essa figura do livro de ecologia.

begonecology

Agora temos como fazer hipóteses subsidiadas em alguma teoria, perfeita não é, mas os pressupostos são simples, a salvo por alguns desvios, essas funções ajustam bem, e funcionam no mundo real. Sabemos que não é perfeito, e que em alguns casos elas falham, mas ainda assim ela nos trazem informação, já que sabemos como as coisas deveriam ser e se existe um desvio do que esperamos, podemos procurar o que esta acontecendo com alguma direção, já que se sabemos como o que acontece no mundo real é diferente desses casos, podemos ver também o que devemos mudar nesses modelos para tentar entender melhor o mundo. Ou seja, apesar do sacrifício de ter que entender todas esses gráficos e contas, no final ganhamos um melhor entendimento de como diabos o mundo funciona. E para entender essas coisas não existe muita mágica, se ficar batendo cabeça, cedo ou tarde a gente começa a entender.

Mas é isso, o script está no repositório do github, além de aqui em baixo, e logo vamos começar a tentar entender o que diabos é a tal da matriz jacobiana.

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

###################################################
### Casos de equilibrio para duas populações 
###################################################
 
### Coexistencia estavel
 
 
layout(matrix(c(1:4),ncol=2,byrow=T))
 
a <- matrix(c(0.01,0.005,0.005,0.01),ncol=2,byrow=TRUE)
 
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])))
box()
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))
 
### Espécie 2 pode invadir, a espécie 1 não
 
a <- matrix(c(0.01, 0.02,0.005, 0.01), ncol=2, byrow=TRUE)
 
curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], lwd=1, lty=2,  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]))
 
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]), lwd=1,c(0, expression(1/alpha[22]),expression(1/alpha[12])))
box()
xs <- c(1, 1, .4, .4) * a[1,1]^-1
ys <- c(.2, .4, .4, .7)*a[2,2]^-1
 
arrows(xs[1:3], ys[1:3], xs[2:4], ys[2:4],lty=2:1, lwd=c(1,1), length=0.12)
points(0,1/a[2,2], cex=2)
 
### Espécie 1 pode invadir, espécie 2 não
 
a <- matrix(c(0.01, 0.005,0.02,0.01), ncol=2, byrow=TRUE)
curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], lwd=1, lty=2,  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]))
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])))
 
ys <- c(1, 1, .4, .4) * a[1,1]^-1
xs <- c(.2, .4, .4, .7)*a[2,2]^-1
 
arrows(xs[1:3], ys[1:3], xs[2:4], ys[2:4],lty=1:2, lwd=c(1,1), length=0.12)
points(1/a[1,1],0, cex=2)
 
### Nenhuma especie pode invadir, com equilíbrio instável
 
a <- matrix(c(  0.01, 0.02,0.02, 0.01), ncol=2, byrow=TRUE)
curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], lty=2,  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]))
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])))
 
pxys <- matrix(c(.13,.1,.25,.1,
                 .1,.13,.1,.25,
                 .115,.115,.25,.25,
                 .7,.8,.4,.8,
                 .8,.7,.8,.4,
                 .75,.75,.4,.4,
                 .25,.4,.1, .6,
                 .4,.25,.6, .1),
              nc=4, byrow=T)
xys <- pxys
xys[,1:2] <-  a[1,1]^-1 * pxys[,1:2]
xys[,3:4] <-  pxys[,3:4] * a[2,2]^-1
arrows(xys[,1], xys[,2], xys[,3], xys[,4],lty=c(1,2,3,1,2,3,2,1), lwd=c(1), length=0.12)
n2star <- (a[1,1] - a[2,1])/(a[1,1]*a[2,2] - a[1,2]*a[2,1])
n1star <- (a[2,2] - a[1,2])/(a[1,1]*a[2,2] - a[1,2]*a[2,1])
points(n1star,n2star, cex=2)
 
#Interações ao acaso
par(mfrow=c(3,3))
 
for(i in c(1:9)) {
a <- matrix(runif(4,0,2),ncol=2,byrow=TRUE)
 
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])))
legend("topright",legend=paste(c("1,1=","1,2=","2,1=","2,2="),round(a,digits=2)),bty="n")
}