A evolução da cooperação e o torneio de Axelrod

A primeira vez que li o livro O gene egoísta, do Richard Dawkins, apesar de muito legal o livro, foi meio triste. Porque com todas aquelas teorias de evolução, do gene egoísta e tal, da a impressão que que evoluímos para todos terem um comportamento egoísta, a luta pela sobrevivência. A gente acaba vendo poucos exemplos onde ser “legal” é eficiente na natureza. Segundo a teoria neodarwinista.

Antes de ler os próximos livros eu fiquei preocupado com isso. No fim do livro o Dawkins, ele ainda fala que uma das características que nos distingue das outras espécies é podermos escapar do “Gene egoísta” porque pensamos, mas ainda sim, isso não era muito reconfortante, isso não explicava muito bem porque vemos tanta gente legal no mundo, essas pessoas não podem ser o “fitness” ruim.

Mas então, lendo os outros livros, como O Relojoeiro Cego, eu li a primeira vez sobre o torneio de Axelrod, e pela primeira vez uma explicação que fazia sentido, em como ser legal pode ser evolutivamente estável.

Muitas perguntas em biologia, como em outras ciências, tem o problema de serem muito difíceis de se fazer hipóteses testáveis. So pensar em evolução, como diabos se testa algo em evolução, separar grupos tratamento e controle e esperar organismos evoluírem nunca é uma proposta viável, dado o tempo que leva para algo acontecer, nossa vida é curta para essa escala de tempo. Talvez de para acompanharmos alguns organismos, como bactérias, para algumas perguntas específicas, mas como testar algo como cooperação, evolução da cooperação.

Um cara chamado John Maynard Smith teve uma sacada. Ele resolveu usar na biologia a teoria dos jogos para testar coisas. Quase todo mundo ja ouviu falar de teoria dos jogos, é o que aquele cara do filme uma mente brilhante estuda, o John Nash.

Mas então, dento desse cenário, o Robert Axelrod, que é um cientista politico e o William D. Hamilton publicaram um artigo na revista science que deu o melhor exemplo, de como a cooperação aparece. Mas não a cooperação como em insetos eu-sociais, tipo formigas, abelhas, etc. Cooperação entre indivíduos independente do seu parentesco.

O que o Axelrod fez, e eu não sei exatamente aonde o W.D. Hamilton entrou na jogado, acho que deve ter entrado para escrever o artigo, analisar os dados, já que o Axelrod fez o Axelrod Tournament. Em teoria dos jogos, as pessoas quantificam jogos, cada jogada tem um custo ou retorno aos seus participantes. Talvez o jogo mais famoso e primeiro exemplo da teoria dos jogos é o dilema do prisioneiro.

O Dilema do prisioneiro é relativamente simples, existem dois caras presos, o cara A e o cara B. Eles são colocados separados em celas diferentes. Dai vem um policial para cada um, e pergunta se ele quer confessar o crime ou não. Então três coisas podem acontecer (na verdade quatro, mas ja vemos isso).

Uma primeira possibilidade é todo mundo é dedo duro, o A entrega o B, e o B entrega o A. Agora o policial tem a confissão dos dois e manda os dois para a cadeia para cumprir dois anos cada.

A segunda possibilidade é quando os dois ficam em silencio, se ambos ficam em silencio, o policial não tem nada contra eles, a não ser as provas que já tinha antes e só consegue colocar cada um por um ano na prisão.

A terceira possibilidade é quanto temos um dedo duro e um cara legal, ai o policial vai ter as provas que tem mais o testemunho do parceiro de crimes, colocando quem foi dedurado por três anos na prisão e liberando o dedo duro, porque ele ajudou a criar provas contra o outro prisioneiro. Mas aqui tanto o A pode ser o dedo duro e o B ficar em silencio como o B ser o dedo duro e o A ficar em silencio, em ambos os casos um se ferra e o outro se da bem.

Agora veja, podemos resumir essas possibilidades em uma matriz onde as linhas são as ações de um prisioneiro e colunas as ações do outro prisioneiro.

 \begin{array}{ccc}  & {Prisioneiro\ B\ fica\ em\ silencio} & {Prisioneiro\ B\ entrega\ o\ parceiro\ A} \\ {Prisioneiro\ A\ fica\ em\ silencio} & {1/1} & {3/0} \\ {Prisioneiro\ A\ entrega\ o\ parceiro\ B} & {0/3} & {2/2} \\ \end{array}

Esses valores podem variar de acordo como a gente quiser, mas a ideia do jogo é essa, temos participantes, aqui o prisioneiros, e cada decisão deles leva a custo, ou ganho, só que esse custo/ganho não depende da ação de um dos prisioneiros apenas, o custo/ganho sempre vai ser dado pela combinação da ação dos dois.

Agora vamos pensar um pouco, se a gente pega somente o prisioneiro A. Vamos supor que o Prisioneiro B ja tenha tomado sua decisão de ficar em silencio, aqui a decisão do prisioneiro A é entre também ficar em silencio e ficar um ano preso ou entregar o B e não ficar nenhum ano preso. Agora vamos supor que diferente disso, o Prisioneiro B tenha escolhido entregar A, agora A enfrenta a decisão de passar três anos presos ao ficar em silencio, ou 2 anos se denunciar B. Veja que para cada um dos participantes, se a gente pensar independente do outro prisioneiro, compensa sempre denunciar o outro prisioneiro, ser o dedo duro, sempre você é o dedo duro, você passa menos tempo na cadeia, é a escolha entre 1 ano ou nenhum ano ou a decisão entre 3 anos ou 2 anos.
Para esses valores, compensa sempre ser dedo duro, nunca cooperar, quem coopera com os outros sempre pode se ferrar bonito, passar 3 anos na cadeia se ficar em silencio, ou se der sorte 1 se o outro prisioneiro também ficar em silencio, mas se cada prisioneiro está em uma cela separada, incomunicável, quem garante que o outro prisioneiro não vai te sacanear? Sendo que ele so tem a “ganhar”.

Seguindo essa analise, nunca esperaríamos que o altruísmo aparece nas espécies. Para testar essa hipótese, o Axelrod teve uma idéia, ele bolou um torneio de computador, onde as pessoas submetiam programas de computador para jogar o dilema do prisioneiro, mas de forma iterada. Você não joga somente uma vez, você joga varias “partidas”, e ao invés de anos de cadeia a gente contabiliza pontos seguindo a seguinte matriz.

 \begin{array}{ccc}                      & {B\ Coopera}   & {B\ nao\ coopera} \\   {A\ Coopera}       & {1/1}          & {2/0} \\   {A\ nao\ Coopera}  & {0/2}          & {-1/-1} \\   \end{array}

Agora a gente tira ponto quanto os dois não cooperam, da um ponto pra cada quando os dois cooperam, e se somente um cooperar e o outro não cooperar, que não cooperou ganha 2 pontos e o outro não ganha nada.

E o jogo vai ser jogado por varias partidas, só que cada jogador sabe o historio do outro jogador com quem já interagiu, então a primeira vez que dois jogadores, vamos chamar de agentes, a primeira vez que eles interagem é novidade, a segunda cada um já lembra o resultado da primeira interação.

Mas além disso, teremos vários agentes, com diferentes estratégias. A gente vai fazer essa simulação do torneio do Axelrod.

Então a gente vai fazer varias funções, que serão os agentes, a usaremos TRUE e FALSE para cooperar ou não cooperar respectivamente.

Mas são essas funções? Vamos pegar um exemplo, a estratégia do ditador agressivo:

# Ditador Agressivo
ag.ditador <- function(h) {
    if (length(h)==0) {
        return(FALSE)
    } else {
        return(mean(h,na.rm=T)<=.5)
    }
}

O ditador agressivo começa o jogo não cooperando.

> ag.ditador(c()) [1] FALSE

Se ele nunca interagiu com você, ele solta logo um “não coopero” (FALSE). Mas ao longo do jogo, se em média você começar a nunca mais cooperar com ele.

> ag.ditador(c(T,T,F,F,F)) [1] TRUE

Veja que o outro jogador cooperou as duas primeiras rodadas e parou de cooperar, ai o ditador agressivo volta atras, e começa a cooperar com você pra ganhar sua amizade denovo.

> ag.ditador(c(T,T,F,F,F,F,T,T,T)) [1] FALSE

Mas ai a hora que você perdoou o ditador agressivo, e começou a cooperar com ele, ele vai e te ferra denovo e volta a não cooperar, e assim vai.

Essa é uma estratégia, podemos fazer funções para as estratégias mais simples possíveis, como sempre cooperar, ou sempre não cooperar, como estratégias elaboradas, para tentar prever o padrão do outro agente e tomar proveito disso, porque sempre que você não coopera e o outro coopera, você ganha 2 pontos.

E na nossa simulação, como eu fiz poucas estratégias, vamos colocar mais de um agente com a mesma estratégia, vamos usar 20 agentes de 11 estratégias diversificadas, e vamos jogar por 500 partidas.
Basicamente cada partida a gente faz 10 pares de agentes, sorteando os pares ao acaso, e fornece para função o histórico de interação entre esses agentes, sendo somente o historio entre esses agentes, e não todas as interações do agente.

Depois disso transformamos as interações em pontos, e podemos ver a pontuação ao longo do tempo.

01

Veja que na minha simulação, eu tive três agentes usando a estrategia tit4tat, essa estrategia foi a ganhadora do torneio, nas duas edições que aconteceram. Quando o Axelrod fez o torneio, ele pediu para um monte de cientista da área de teoria dos jogos, optimização, para enviarem programas com estratégias pro campeonato, as que estão nessa simulação, não são nem de longe tão complexas como a que as pessoas tentaram. Tinha de tudo, no entanto, a estrategia extremamente simples chamada tit4tat venceu. Além disso, quando acabou a primeira edição, o Axelrod publicou o resultado do torneio e deixou publico todos os programas enviados. Mas na segunda edição, mesmo todo mundo podendo estudar os programas dos outros, estudar o resultado do jogo, podendo bolar agora maneiras de como explorar as estratégias das outras pessoas, o tit4tat venceu denovo na segunda edição.

Mas como funciona o tit4tat?
O Tit for Tat foi feito por um russo chamado Anatol Rapoport. Basicamente ele começa cooperando, e sempre coopera enquanto o outro agente também cooperar, se em alguma rodada o outro agente não cooperar, então ele começa a não cooperar também, mas ele só vai não cooperar até que o outro agente volte a cooperar, a partir dai ele volta a cooperar também. Simples assim.

Vamos lembrar o esquema de pontuação.

 \begin{array}{ccc}                      & {B\ Coopera}   & {B\ nao\ coopera} \\   {A\ Coopera}       & {1/1}          & {2/0} \\   {A\ nao\ Coopera}  & {0/2}          & {-1/-1} \\   \end{array}

Agora pense nessa estratégia contra as mais simples, que é de nunca cooperar e sempre cooperar.
Quando o tit4tat encara um agente que nunca coopera, como na primeira rodada ele vai cooperar e o outro agente não vai cooperar, o outro agente vai ganhar 2 pontos e o tit4tat não ganha nada, dai pra frente, os dois nunca mais vão cooperar entre si, sempre ganhando -1 ponto cada um e a pontuação vai somente despencar, mas ainda sim o outro agente sempre estará 2 pontos a frente do tit4tat, ou seja, sozinho o tit4tat nunca ganha desse agente.
Por outro lado se o tit4tat for jogar com alguém que sempre coopera, os dois vão sempre cooperar e cada um vai ganhar sempre um ponto, por quantas rodadas forem o jogo, isso faz com que o tit4tat consiga empatar com esse agente, mas nunca ganhar.

Mais ou menos da para ver que tit4tat nunca ganha, no melhor caso ele empata, no entanto os três agentes deles estão na frente na nossa simulação.

Acontece que temos que pensar no resultados em grupos, o tit4tat não esta ganhando na verdade, ele começa a se dar bem, quando interage com outros agentes que também cooperam, ganhar um ponto é sempre melhor que perder 1 ponto. Veja que algumas estratégias, pouco depois do inicio começam a despencar na pontuação, por exemplo o laranja, que é o sempre trapaceia cai fácil, porque como os programas tem a memoria de como foram as interações anteriores, alguns agentes são fácil de decidir que nunca vão cooperar, sempre vão tentar passar a perna, então quando todo mundo começa a não cooperar com ele, ele sempre perde pontos.
Por outro lado algumas estratégias são mais confusas. Pegue o ditador agressivo por exemplo, como ele muda de estratégia de acordo com a sua, sempre tentando te convencer a voltar a ser cooperativo e então trapaceando denovo, ele acaba se mantendo bem, mas veja de temos 3 agentes dele também, mas somente 2 continuam bem no final, um começa a descambar.

No artigo do Axelrod e do Hamilton, apesar do tit4tat ter sido a estratégia vitoriosa, eles tenta ver o que fez algumas estrategias serem melhores que muitas outras. E o que observaram é que cooperar pode não ser bom, quando se está no meio de um monte de trapaceiros, mas quando você encontra outros agentes que também cooperam, esses dois começam a melhorar a pontuação juntos.
É claro que se você ficar cooperando com todo mundo, você vai se ferrar, como a gente o agente que sempre coopera, ele não vai bem, mas se você cooperar com quem coopera, e trapacear quem trapaceia, logos os trapaceiros começam a não conseguir mais somar pontos, e somente os agentes cooperando entre si conseguem pontos.
So que as vezes alguns agente deixam de cooperar, como o agente ditador agressivo, mas quando ele volta a cooperar, se você perdoa-lo, vocês podem voltar a a fazer pontos juntos.

O que o Axelrod e Hamilton viram foi isso, características como sempre cooperar com quem coopera com você, trapacear os trapaceiros e perdoar quem quer mudar de estrategia gerava a melhor pontuação. Logo ser legal com pessoas legais pode ser uma estratégia evolutivamente estável.

Essa conclusão é muito legal né. Claro que comportamento em espécies é muito mais complexo que essas regrinhas simples, mas por exemplo em primatas, que vivem em grupos e a gente ve um tirando os piolhos dos outros é um exemplo. Talvez seria melhor deixar alguém tirar seu parasitas mas não perder tempo fazendo isso pra ninguém. Mas logo ninguém mais vai querer interagir com você no grupo e você vai encher de piolhos, logo todo mundo tem que ser legal um com os outros, se for querer viver sem piolhos.

Outro exemplo é você pensar em dois insetos, onde cada um poe um ovo, se os dois fizerem um ninho, eles podem revesar, onde um cuida dos ovos e o outro come, e depois eles invertem os papeis, mas se na sua vez de revesar, você abandonar o ovo do seu parceiro, ele não deixara descendentes mas você, nessa geração você se deu bem, mas logo na próxima não haverá ninguém que cooperara. Assim apesar de em uma geração você ter uma grande vantagem, na próxima você se ferra. Como no jogo, no inicio todos os agentes estão embolados, mas no final muitos deles despencam. A estratégia que sempre trapaceia é funciona assim. Veja que a linha laranja, até a rodada 50 por ai, estava entre as melhores pontuações, somente depois que ele começa a cair, ai sem ninguém para trapacear, so perde pontos.

Isso tudo é muito legal, vermos como uma simulação pode trazer tantas conclusões assim. O W.D. Hamilton tem modelos muito legais, além de ser uma personalidade muito interessante. No obituário que o Dawkins escreveu quando o Hamilton morreu, ele comenta que o Hamilton, professor titular de Cambridge ia de bicicleta trabalhar todo dia e não ligava muito para muitas coisas do nosso mundo comum, como dinheiro ou luxo. Muito legal.

Existe muita coisa para escrever sobre esse tópico, mas uma hora o post tem que acabar, e melhor acabar por aqui. O script para a simulação está embaixo, é legal ver ela varias vezes, variando outras caracteristicas como números de agentes, ou mesmo escrever outras estrátegias, e ver como as vezes o tit4tat chega a perder, mas quase sempre ele está bem colocado, e aumentando os pontos. Tanto o número de agentes como a constituição deles influencia muito o resultado, tente colocar muitos ditadores ou ainda muitos agentes imprevisiveis, que da resultado bem legais para se pensar. Mas uma evidencia para a cooperação, para um porque as pessoas serem boas é muito legal para dormir tranquilo, pensando sempre que o mundo não deve ter somente pessoas ruins, ja vivemos com tantos problemas na vida, ler esse artigo traz algum conforto na vida, pelo menos pra mim. Os scripts sempre no repositório recologia e eu ja tinha pego um script parecido de onde copiei essa simulação a muito tempo atrás, mas nem lembro mais o link, mas foi outra pessoa que bolou esse esquema para a simulação, so não lembro o nome para dar o devido crédito, peço desculpa, mas tenho certeza que todo mundo fica feliz em ver seus script e funções rodando por ai. E é isso, desculpem os erros de portugues, ja que escrevi muito rapido esse texto.
Uma ultima coisa, eu ouço um podcast sobre ciência chamado fronteiras da ciência, o episodio 19 da quinta temporada é sobre a teoria dos jogos, são fisicos falando disso, que devem entender teoria dos jogos melhor que eu, então vale muito a pena ouvir, caso tenha mais interesse no assunto, e tem esse programa que eu ja vi aqui, que também reproduz o torneio e é bem legal, para quem quer só brincar e não mexer com script, e até o próximo post.

# Simulação do torneio de Axelroad
set.seed(1415)
 
# Retorno para:
# Cooperação mutua
coopera.coopera <- 1
# Você trapaceia mas seu oponente coopera
trapaceia.coopera <- 2
# Você coopera mais seu oponente trapaceia
coopera.trapaceia <- -2
# Todo mundo trapaceia
trapaceia.trapaceia <- -1
 
# Função que retorna o valor de resultado da interação para o participante
retorno <- function(voce,oponente) {
    saida <- 0*voce
    saida[(voce==1)&(oponente==1)] = coopera.coopera
    saida[(voce==0)&(oponente==1)] = trapaceia.coopera
    saida[(voce==1)&(oponente==0)] = coopera.trapaceia
    saida[(voce==0)&(oponente==0)] = trapaceia.trapaceia
    return(saida)
}
 
# Exemplo da função para tres jogos
retorno(c(0,0,0),c(1,1,0))
 
 
#Os tipos de agentes.
 
#Os agentes simples:
sempre.coopera   <- function(h) { return(TRUE) }
sempre.trapaceia <- function(h) { return(FALSE) }
 
 
# Imprevisivel
imprevisivel <- function(h) { return(1==rbinom(1,1,.5)) }
 
#tit4tat
 
tit4tat <- function(h) {
    if (length(h)==0) {
        return(TRUE)
    } else {
        return(rev(h)[1]==1)
    }
}
 
 
# tit4tat inicialmente agressivo
ag.tit4tat <- function(h) {
  if (length(h)==0) {
      return(FALSE)
  } else {
      return(rev(h)[1]==1)
  }
}
 
# tit4tat oportunista
op.tit4tat <- function(h) {
  if (length(h)==0) {
      return(TRUE)
  } else {
      return(((n<n.iteracoes-n.agentes)&(rev(h)[1]==1)))
  }
}
 
 
# fairbot
fairbot <- function(h) {
  if (length(h)==0) {
      return(TRUE)
  } else {
      return(mean(h,na.rm=T)>=.5)
  }
}
 
# fairbot inicialmente agressivo
ag.fairbot <- function(h) {
    if (length(h)==0) {
        return(FALSE)
    } else {
        return(mean(h,na.rm=T)>=.5)
    }
}
 
# Ditador covarde
ditador <- function(h) {
    if (length(h)==0) {
        return(TRUE)
    } else {
        return(mean(h,na.rm=T)<=.5)
    }
}
 
# Ditador Agressiveo
ag.ditador <- function(h) {
    if (length(h)==0) {
        return(FALSE)
    } else {
        return(mean(h,na.rm=T)<=.5)
    }
}
 
# Frenemy
frenemy <- function(h) {
    if (length(h)==0) {
        return(TRUE)
    } else {
        return(mean(rev(h)[1:2],na.rm=T)==1)
    }
}
 
ag.ditador(c())
 
ag.ditador(c(T,T,T,T,T,T))
ag.ditador(c(T,T,F,F,F))
ag.ditador(c(T,T,F,F,F,F,T,T,T))
 
 
#Parametros da simulação
n.iteracoes <- 500
n.agentes   <- 20
 
lista.agentes <- c("sempre.coopera", "sempre.trapaceia", "tit4tat" , "ag.tit4tat",
               "fairbot"    , "ag.fairbot"  , "ditador", "ag.ditador",
               "op.tit4tat" , "imprevisivel"  , "frenemy")
 
# Selecionando agentes para a simulação
agentes <- sort(sample(lista.agentes, n.agentes, replace=T))
 
# Matriz da historia de cada jogador
agentes.numero <- matrix(NA, ncol=n.agentes, nrow=n.iteracoes)
 
# Matriz de ações dos agentes
agentes.acoes.recebidas <- matrix(NA, ncol=n.agentes, nrow=n.iteracoes)
agentes.acoes.aplicadas <- matrix(NA, ncol=n.agentes, nrow=n.iteracoes)
 
# Pontuação historicas
agentes.pontos <- matrix(NA, ncol=n.agentes, nrow=n.iteracoes)
 
for (n in 1:n.iteracoes) {
  while (sum(is.na(agentes.numero[n,]))>0) {
  sempar <- 1:n.agentes
    while (length(sempar)>1) {
      i <- sempar[1]
      sempar <- sempar[sempar!=i]
      i.adversario <- sample(rep(sempar,2), 1)
      sempar <- sempar[sempar!=i.adversario]
      agentes.numero[n,i] <- i.adversario
      agentes.numero[n,i.adversario] <- i
    }
  }
  for (i in 1:n.agentes) {
    a.encarado <- agentes.numero[n,i]
    a.numeros <- agentes.numero[,i]
    a.acao <- agentes.acoes.recebidas[,i]
    h <- a.acao[(a.numeros==a.encarado)&(n>1:n.iteracoes)]
    response <- get(agentes[i])(h)
    agentes.acoes.aplicadas[n,i] <- response
    agentes.acoes.recebidas[n,a.encarado] <- response
  }
}
 
# Calculando os scores por rodada
placar.interacao <- retorno(agentes.acoes.aplicadas, agentes.acoes.recebidas)
placar.acumulado <- apply(placar.interacao, 2, cumsum)
 
 
#Figura 1
cores <- rainbow(length(lista.agentes))
plot(x=c(1,n.iteracoes+25),y=c(min(placar.acumulado),max(placar.acumulado)),type="n",frame=F,
     main="Performance ao longo do tempo", xlab="Encontro", ylab="Pontuação")
for (i in 1:n.agentes) {
    lines(1:n.iteracoes, placar.acumulado[,i], col=cores[agentes[i]==lista.agentes], lwd=2)
}
legend("topleft",lty=1,legend=lista.agentes,col=cores,bty="n",lwd=2)

Regressão linear ajustada com lm, bbmle, OpenBugs, Jags e Stan no R

No último post a gente falou do pacote sads do R. Muito legal, e ele faz vasto uso do pacote bbmle.

O pacote bbmle é bem legal, ele permite você escrever uma função de verosimilhança qualquer e então ele ajusta ela pra você. A gente tem um monte de post, por exemplo aqui, aqui e aqui onde o que a gente faz no fim das contas é escrever a função de verosimilhança e ajustar. Mas a gente sempre tinha ajustado usando estatística bayesiana, usando o OpenBugs pra ser mais preciso, e para conferir usando estatística frequentista usando o lme4, ou algum outro pacote que fazia o modelo pra gente.

Mas o bblme vai funcionar mais ou menos na mesma estrategia da linguagem bugs, escreve seu modelo e o bbmle ajusta pra você.

Vamos pegar um modelo simples aqui, a boa e velha regressão linear, e criar alguns dados.

01

Beleza, o mais simples de fazer é usar a função lm que ajusta modelos lineares no R a ajustar o modelo e olhar o resultado e pronto, simples assim.

> fit.lm <- lm(resposta~preditor) > summary(fit.lm) Call: lm(formula = resposta ~ preditor) Residuals: Min 1Q Median 3Q Max -1.0381 -0.4135 0.1270 0.3737 0.7836 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.41354 0.32834 7.351 5.28e-08 *** preditor 0.47595 0.01535 31.004 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.491 on 28 degrees of freedom Multiple R-squared: 0.9717, Adjusted R-squared: 0.9707 F-statistic: 961.3 on 1 and 28 DF, p-value: < 2.2e-16

Certo, agora vamos lembrar como é a formula da reta da regressão.
O que a gente fez é:

resposta = \alpha + \beta \cdot preditor + erro

E esse erro é a parte onde entra a distribuição normal, gaussiana ou sei la, mas o erro então é o seguinte.

erro = Normal(0,\sigma^2)

A distribuição com média zero e desvio sigma.

Daria na mesma o que a gente escreveu acima estar tudo junto:

resposta = \alpha + \beta \cdot preditor + Normal(0,\sigma^2)

Ou ainda:

resposta = Normal(\alpha + \beta \cdot preditor,\sigma^2)

Ou as vezes ao invés daquele escrito normal acho que a gente vê escrito p ou P, mas as pessoas usam p pra tanta coisa em estatística que não é a toa que tudo fica confuso, e não sou só eu que fica confuso, veja so aqui.

Mas a gente consegue perceber, pelo gráfico, que a resposta e o preditor são nossos dados, que temos valores, e o que sobra são os chamados parâmetros nessa formula toda. Na regressão, a gente está lidando com três deles, dois que determinam a reta, o alpha (\alpha) que é onde ela cruza o eixo de resposta e o beta (\beta) que é a inclinação. Além desses dois caras temos o desvio padrão, que é o sigma (\sigma), esse cara diz o quanto afastado daquela reta as nossas amostras podem estar, não exatamente isso mas mais ou menos isso, a dispersão dos pontos em relação a reta.

Ok, por mais confuso que isso possa parecer, se esse é o modelo, e resposta e preditor nos temos, sempre esses três caras tem que sair do modelo, que no caso a gente otimiza para conseguir o melhor ajuste.

Então o bbmle, assim como o OpenBugs e companhia vai precisar de um modelo construído, que diga pra ele como esta a qualidade do ajuste, verosimilhança do modelo, e ele te diz quais os melhores valores para esses três parâmetros.

Na saida da regressão ali em cima, o alpha é o intercept (2.41354), o beta é o que ele chamou de preditor(0.47595) e por último, o sigma é o Residual standard error (0.491). Os parâmetros alpha e beta não se tem como ter certeza absoluta, so quem tem certeza é quem gerou os dados, como na natureza é ela que gera os dados, então a gente nunca tem certeza, mas a gente acha que esses são os melhores valores (lembrando que cada valor tem seu erro aceitável também, ali do lado).

Mas então, o lm faz todo esse processo pra gente, e não tem porque não usar ele para modelos assim, mas só para treinar, ou curiosidade, vamos ver como ficaria esse modelo no bbmle.

No bbmle a gente tem que fazer nossa própria função de verosimilhança.

funcao.likelihood <- function(alpha=1,beta=1,desvio=1) {
    -sum(dnorm(resposta,mean=alpha+beta*preditor,sd=desvio,log=TRUE))
}

Antes de pensar que essa função faz muita coisa, vamos só usar ela, pra ver como é mais simples que parece, ela simplesmente calcula um número, que é a verosimilhança do modelo, na verdade o negative log-likelihood, veja que tem um menos antes do sum ali.

Então você fala um valor pra cada parâmetro e a função fala se o quão bom é o ajuste, em termos de verossimilhança.

> funcao.likelihood(alpha=1,beta=1,desvio=1) [1] 1488.255

Se a gente muda esses valores, o ajuste muda

> funcao.likelihood(alpha=1.5,beta=1,desvio=1) [1] 1632.549

Quanto menor melhor, quanto mais a gente aproxima todos os valores dos valores que usamos para gerar os dados, menor ele vai ficar.

> funcao.likelihood(alpha=1.5,beta=0.9,desvio=1) [1] 1038.321

E se a gente usar parâmetros muito nada a ver, esse número fica gigante.

> funcao.likelihood(alpha=8,beta=0.9,desvio=1) [1] 3195.411

E isso é tudo que o bbmle precisa. Agora para ajustar os parâmetros é so usar essa função de verosimilhança na função mle do pacote bbmle.

> fit.mle <- mle2(funcao.likelihood) Houve 11 avisos (use warnings() para vê-los) > summary(fit.mle) Maximum likelihood estimation Call: mle2(minuslogl = funcao.likelihood) Coefficients: Estimate Std. Error z value Pr(z) alpha 2.413540 0.317207 7.6087 2.768e-14 *** beta 0.475945 0.014830 32.0925 < 2.2e-16 *** desvio 0.474314 0.061234 7.7459 9.487e-15 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 -2 log L: 40.38304

E pimba, estimamos tudo do modelo. O desvio, ou sigma aqui fica menor que na função lm, porque o método usado é diferente, mas basicamente temos o mesmo resultado. Recebemos um aviso também, que acho que porque pode ser tentando valores negativos de desvio na função, mas não existe desvio negativo, dai a função produz um NaN, tipo assim.

> funcao.likelihood(alpha=8,beta=0.9,desvio=-1) [1] NaN Mensagens de aviso perdidas: In dnorm(resposta, mean = alpha + beta * preditor, sd = desvio, : NaNs produzidos

Mas o mle acha o caminho certo e ajusta. Ele usa o método do Nelder and Mead como default, a gente falou do método do Newton aqui, e vimos o gradient descent aqui, esse é outro método que serve para achar os melhores valores para os parâmetros.

Mas o bbmle é legal que da para usar outros métodos de otimização, só citar qual método você quer usar.

> fit.mle2<-mle2(funcao.likelihood,method="L-BFGS-B",lower=c(alpha=0.0001,beta=0.0001,desvio=0.0001)) > summary(fit.mle2) Maximum likelihood estimation Call: mle2(minuslogl = funcao.likelihood, method = "L-BFGS-B", lower = c(alpha = 1e-04, beta = 1e-04, desvio = 1e-04)) Coefficients: Estimate Std. Error z value Pr(z) alpha 2.413540 0.317207 7.6087 2.768e-14 *** beta 0.475945 0.014830 32.0924 < 2.2e-16 *** desvio 0.474314 0.061234 7.7459 9.488e-15 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 -2 log L: 40.38304

Ou ainda usar outra função para otimizar, que não o optmin que é o otimizador default do R, eu acho.

> fit.mle3<-mle2(funcao.likelihood,optimizer="nlminb",lower=c(alpha=0.0001,beta=0.0001,desvio=0.0001)) > summary(fit.mle3) Maximum likelihood estimation Call: mle2(minuslogl = funcao.likelihood, optimizer = "nlminb", lower = c(alpha = 1e-04, beta = 1e-04, desvio = 1e-04)) Coefficients: Estimate Std. Error z value Pr(z) alpha 2.413539 0.317206 7.6087 2.768e-14 *** beta 0.475945 0.014830 32.0925 < 2.2e-16 *** desvio 0.474313 0.061234 7.7460 9.486e-15 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 -2 log L: 40.38304

Veja que esses outros métodos nos permitiram fazer mais coisas, uma delas foi delimitar limites para os parâmetros. No caso colocamos que tudo tem que ser maior que 0.0001, isso impossibilita um valor negativo pro desvio e veja que depois disso a gente não recebe mais warnings.
Eu acho que a parte boa do bbmle, é que se amanha a gente lê um artigo, com um modelo muito louco que alguém acabou de inventar, se a gente conseguir passar ele para uma função de verossimilhança, a gente consegue ajustar valores pra ele.

Alias o cara que escreveu esse pacote, o Ben Bolker, que é um cara que parece ser, além de inteligente, bem legal, já que ele responde todo mundo nos fóruns de R, e em tudo que é fórum da internet, bem ele escreveu um livro chamado Ecological Models and Data in R, que é muito legal, pois ao invés de ensinar sobre testes e mais testes com seus nomes complicados, ele ensina que existe funções determinísticas, distribuições estocásticas,e você juntas essas coisas e descreve o mundo, ai depois de muitos capítulos ele mostra que você continua fazendo a mesma coisa que fazia nas anovas da vida, mas usando o modelo que você pensa e não pensando o que pode testar com os dados que tem.

Mas já que estamos aqui, veja que na estatística bayesiana, apesar de partir de pressupostos diferentes, paradigmas diferente e tal, no final a gente quer descrever a mesma coisa, e vai ser a mesma coisa na pratica, escrever um modelo, ajustar os dados, essa mesma regressão é feita da seguinte forma no OpenBugs

Escrevemos o modelo

sink("linreg.txt")
cat("
model {
# Prior
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
 sigma ~ dunif(0, 100)
 tau <- 1/ (sigma * sigma)
 
# Verossimilhança
 for (i in 1:n) {
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- alpha + beta*x[i]
 }
}
",fill=TRUE)
sink()

Depois temos que preparar o terreno, preparar os dados para exportar, criar uma função para inicializar a cadeia, definir os parâmetros da cadeia e então otimizar.

> bugs.data <- list(x=preditor,y=resposta,n=length(resposta)) > inits <- function(){ list(alpha=rnorm(1), beta=rnorm(1), sigma = rlnorm(1))} > params <- c("alpha","beta","sigma") > nc = 3 #Numero de cadeias > ni=1200 #Tamanho da cadeira > nb=200 #Numero de simulação que serão descartadas > nt=1 #Thinning rate > modelo1.bugs <-bugs(data = bugs.data, inits = inits, parameters = params,model = "linreg.txt",n.thin = nt,n.chains = nc,n.burnin = nb, n.iter = ni) > print(modelo1.bugs, dig = 3) Inference for Bugs model at "linreg.txt", Current: 3 chains, each with 1200 iterations (first 200 discarded) Cumulative: n.sims = 3000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha 2.460 0.354 1.812 2.223 2.438 2.680 3.230 1.013 170 beta 0.474 0.017 0.436 0.464 0.475 0.485 0.504 1.017 140 sigma 0.514 0.071 0.398 0.462 0.506 0.559 0.674 1.002 1600 deviance 43.640 2.655 40.630 41.700 42.960 44.780 50.750 1.003 1600 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 = Dbar-Dhat) pD = 2.9 and DIC = 46.5 DIC is an estimate of expected predictive error (lower deviance is better).

Aqui o sigma fica um pouquinho maior, mas ainda temos basicamente o mesmo resultado. O bbmle é um pacote do R, aqui é meio diferente, porque o OpenBugs é um programa, de vida própria, que do R a gente um pacote chamado R R2OpenBUGS que tudo que faz é mandar todas essas informações que separamos no R para o OpenBugs, o OpenBugs processa elas e depois devolve bonitinho pro R, tudo organizadinho para então a gente interpretar, e fazer o que mais quiser como gráficos.

O OpenBugs, que na verdade vem da linguagem Bugs, de fazer inferência bayesiana, no caso Bugs significa Bayesian inference Using Gibbs Sampling, que usa um tipo de amostrador para tentar descobrir os valores reais dos parâmetros, e chegar ao melhor chute pelo menos. No final a gente queria saber o valor real do parâmetro, mas nunca saberemos, tem um post sobre Gibbs Sampler aqui, se a curiosidade bater.

Mas o OpenBugs não é o único a usar a linguagem Bugs para inferência, o jeitão de escrever o modelo etc, temos o WinBugs, que não esta mais em desenvolvimento, e o Jags, Just another Gibbs Sampler, que está a todo vapor eu acho. Na verdade ajustar modelos com o Jags, que é um programa externo também, é bem parecido, a gente pode ajustar o mesmo modelo denovo, da seguinte forma.

> jags <- jags.model(file='linreg.txt',data=bugs.data,inits=inits,n.chains = 3,n.adapt=100) Compiling model graph Resolving undeclared variables Allocating nodes Graph Size: 130 Initializing model |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% > update(jags, 1000) |**************************************************| 100% > jags.samples(jags,c('alpha', 'beta','sigma'),1000) |**************************************************| 100% $alpha mcarray: [1] 2.380524 Marginalizing over: iteration(1000),chain(3) $beta mcarray: [1] 0.4774447 Marginalizing over: iteration(1000),chain(3) $sigma mcarray: [1] 0.5097927 Marginalizing over: iteration(1000),chain(3)

E temos nossos três parâmetros denovo, veja que usamos tudo igual ao que usamos no OpenBugs. O modelo, o jeito de mandar os dados, numa lista e tal. A gente ouve muito falar do Jags, eu acho que uso mais o OpenBugs porque os livros que olho sempre o pessoal ta usando ele, mas migrar para os modelos pro Jags é relativamente fácil, como acabamos de ver.

Agora uma outra nova opção que surgiu foi o Stan. Stan é um programa para ajustar modelos bayesianos, e mais novo, acho que tem mais suporte a computação paralela, mas ao contrario do OpenBugs e do Jags que usam o Gibbs Sampler, ele usa um negocio chamado Hamiltonian Monte Carlo, esse eu não tenho nem ideia de como funciona, mas a gente pode ajustar a mesma regressão com ele, o que não é tão difícil, e ver que funciona igual.

O jeito de construir os modelos nele é bem diferente, essa nossa regressão fica algo assim.

linreg_stan <- '
data{
  int n;
  real y[n];
  real x[n];
}
 
parameters{
  real alpha;
  real beta;
  real<lower=0, upper=100> sigma;
}
 
transformed parameters{
    real mu[n];
    for(i in 1:n){
    mu[i] <- alpha + beta*x[i];
  }
}
 
model{
  //Priors
  alpha ~ normal(0, 100);
  beta ~ normal(0, 100);
  sigma ~ uniform(0, 100);
 
  //Likelihood
  for(i in 1:n){
    y[i] ~ normal(mu[i], sigma);
}
}
'

O que no final das contas nem é tão diferente assim do Bugs, é definir a função de verosimilhança, como é o framework bayesiano, a gente tem os prior, só que além disso la no inicio, a gente tem que inicializar todos os valores que vamos usar, como no c++, linguagem que ele é feito. Dai é so preparar os dados para mandar, numa lista igual ao bugs e mandar ajustar.

> stan.data <- list(n = length(resposta), x = preditor, y = resposta) > fit.stan<- stan(model_code=linreg_stan,data=stan.data,iter = 1000,chains = 3,thin = 1) TRANSLATING MODEL 'linreg_stan' FROM Stan CODE TO C++ CODE NOW. COMPILING THE C++ CODE FOR MODEL 'linreg_stan' NOW. SAMPLING FOR MODEL 'linreg_stan' NOW (CHAIN 1). Iteration: 1 / 1000 [ 0%] (Warmup) Iteration: 100 / 1000 [ 10%] (Warmup) Iteration: 200 / 1000 [ 20%] (Warmup) . . . Um monte de mensagens . . . Iteration: 900 / 1000 [ 90%] (Sampling) Iteration: 1000 / 1000 [100%] (Sampling) # Elapsed Time: 0.166985 seconds (Warm-up) # 0.131985 seconds (Sampling) # 0.29897 seconds (Total)

E podemos olhar nossos resultado de forma parecida ao output do OpenBugs

> print(fit.stan,digits=3) Inference for Stan model: linreg_stan. 3 chains, each with iter=1000; warmup=500; thin=1; post-warmup draws per chain=500, total post-warmup draws=1500. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat alpha 2.420 0.020 0.369 1.668 2.184 2.419 2.656 3.127 324 1.005 beta 0.476 0.001 0.017 0.442 0.465 0.475 0.487 0.511 331 1.003 sigma 0.520 0.004 0.077 0.395 0.464 0.513 0.567 0.693 454 1.003 lp__ 4.950 0.078 1.449 1.411 4.339 5.278 5.954 6.517 347 1.002 Samples were drawn using NUTS(diag_e) at Thu Jul 17 23:01:30 2014. For each parameter, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat=1).

E denovo os mesmo valores de alpha, beta e sigma.
Bem a gente começou olhando o bbmle, usando o lm pra ver se estávamos chegando a valores razoáveis para os parâmetros da regressão e acabamos olhando também como ajusta usando o OpenBugs, Jags e Stan.

A primeira pergunta agora seria, qual caminho seguir?

A minha melhor resposta é:
Eu não faço a menor ideia, mas é legal pra caralho chegar quase no mesmo valor de alpha, beta e sigma de todos os jeitos.
Na pior da hipóteses, se o modelo que estivermos tentando ajustar for certo, vai dar certo. Claro que para uma regressão linear a gente vai usar o lm na verdade, acho que não faz sentido adicionar toda a complexidade desses outros métodos para algo simples. Mas para muitos outros modelos de nosso interesse, que não da para ajustar com lm, ou mesmo glm, o jeito é começar a embrenhar nesses métodos.

Claro que se da para encontrar as coisas prontas, como os modelos prontinhos no pacote sads, vamos usar la e não tentar reinventar a roda. Mas as vezes a gente quer mudar algo um pouquinho, e ai saber como tudo funciona é legal. Além disso, como no livro do Bolker, um comentário que acho muito interessante, não interessa o método que você usa, é sempre bom saber um pouquinho de cada método para conseguir ler melhor a literatura, já que sempre cada pessoa vai fazer as coisas de um jeito, e não vamos deixar de ler um artigo porque ele usa estatística bayesiana ou frequentista.

Bem e esse post é legal, porque eu sempre quis ver alguém ajustando algo que eu entenda, como uma regressão linear com todos esses programas e pacotes, como ninguém fez isso pra mim ler, eu mesmo fiz ^^. Agora eu tenho uma referencia quando quiser usar algo diferente do OpenBugs.

E é isso. o script vai estar la no repositório recologia, se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo, ou mande um e-mail e ficamos por aqui.

set.seed(171)
preditor<-runif(30,10,30)
resposta<-rnorm(30,2+0.5*preditor,0.5)
 
plot(resposta~preditor)
 
############################
#lm do stats
############################
fit.lm <- lm(resposta~preditor)
summary(fit.lm)
 
############################
#bbmle
############################
library(bbmle)
 
funcao.likelihood <- function(alpha=1,beta=1,desvio=1) {
    -sum(dnorm(resposta,mean=alpha+beta*preditor,sd=desvio,log=TRUE))
}
 
fit.mle <- mle2(funcao.likelihood)
summary(fit.mle)
warnings()
 
fit.mle2<-mle2(funcao.likelihood,method="L-BFGS-B",lower=c(alpha=0.0001,beta=0.0001,desvio=0.0001))
summary(fit.mle2)
 
fit.mle3<-mle2(funcao.likelihood,optimizer="nlminb",lower=c(alpha=0.0001,beta=0.0001,desvio=0.0001))
summary(fit.mle3)
 
############################
#OpenBugs
############################
library(R2OpenBUGS)
 
sink("linreg.txt")
cat("
model {
# Prior
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
 sigma ~ dunif(0, 100)
 tau <- 1/ (sigma * sigma)
 
# Verossimilhança
 for (i in 1:n) {
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- alpha + beta*x[i]
 }
}
",fill=TRUE)
sink()
 
bugs.data <- list(x=preditor,y=resposta,n=length(resposta))
 
# Função de inicialização
inits <- function(){ list(alpha=rnorm(1), beta=rnorm(1), sigma = rlnorm(1))}
 
# Escolha os parametros que quer estimar
params <- c("alpha","beta","sigma")
 
# Caracteristicas do MCMC
nc = 3 #Numero de cadeias
ni=1200 #Tamanho da cadeira
nb=200 #Numero de simulação que serão descartadas
nt=1 #Thinning rate
 
# Inicie o Amostrador
 
modelo1.bugs <-bugs(data = bugs.data, inits = inits, parameters = params,model = "linreg.txt",n.thin = nt,
                     n.chains = nc,n.burnin = nb, n.iter = ni)
 
print(modelo1.bugs, dig = 3)
 
############################
#Jags
############################
library(rjags)
 
jags <- jags.model(file='linreg.txt',data=bugs.data,inits=inits,n.chains = 3,n.adapt=100)
update(jags, 1000)
jags.samples(jags,c('alpha', 'beta','sigma'),1000)
 
 
############################
#Stan
############################
library(rstan)
 
linreg_stan <- '
data{
  int n;
  real y[n];
  real x[n];
}
 
parameters{
  real alpha;
  real beta;
  real<lower=0, upper=100> sigma;
}
 
transformed parameters{
    real mu[n];
    for(i in 1:n){
    mu[i] <- alpha + beta*x[i];
  }
}
 
model{
  //Priors
  alpha ~ normal(0, 100);
  beta ~ normal(0, 100);
  sigma ~ uniform(0, 100);
 
  //Likelihood
  for(i in 1:n){
    y[i] ~ normal(mu[i], sigma);
}
}
'
 
stan.data <- list(n = length(resposta), x = preditor, y = resposta)
 
fit.stan<- stan(model_code=linreg_stan,data=stan.data,iter = 1000,chains = 3,thin = 1)
print(fit.stan,digits=3)

Distribuição de abundância de espécies

Um dos problemas em análises de dados, ao meu ver, é a constante discrepância entre o que as pessoas realmente pensam e as medidas usadas. Muitas vezes as pessoas tem ideas e descrições elaboradas dos processos que estão acontecendo na natureza, mas depois que os dados são coletados, e vamos confrontar esses dados com os processos descritos, parte da informação é perdida devido a escolha das medidas para representar esses dados.

Uma estatística qualquer é um valor que tenta resumir os dados, podemos pensar no chi-quadrado, valor F, etc. Esses valores resumem uma certa informação, que tem sua utilidade, mas o principal é que sabendo o valor e um grau de liberdade, mais ou menos conseguimos reconstruir a informação original, não perfeitamente, mas ter uma ideia de como eram os dados que geraram essa estatística, esse resumo dos dados.

Em ecologia de comunidades temos também varias formas de resumir os dados, de calcular estatísticas da comunidade. Mas alguns desses valores, dessas estatísticas são bem destrutivas sozinhas. Destrutivas no sentido que somente com o valor calculado, é quase impossível reconstruir os dados originais. Vamos por exemplo pegar a riqueza de espécies. Saber que uma comunidade tem Riqueza 5, 5 especies é uma informação bem ambígua, porque a gente pode ter uma espécie dominante e 4 raras, um deserto com 5 espécies raras, até dizer que uma espécie é rara é complicado, porque ser rara no fim das contas depende da abundância das outras espécies. Índices de diversidade vão nos levar a problemas semelhantes, mas não vamos ficar aqui só reclamando.

Pensando nisso, pessoas bolaram outras formas de descrever comunidades de forma a englobar melhor as informações do ambiente, temos a famosa estatísticas multivariadas, que tem bastante implementações prontinhas no pacote vegan, prontas para usar e, o que vamos ver agora, curvas de distribuição de abundância de espécies.

Como isso funciona? Bem como cada espécie tem uma abundância na comunidade, podemos fazer um histograma das abundâncias das espécies. Vamos fazer um histograma de densidade de árvores nos dados do Hubbell, que comentamos no post anterior, cada espécie tem sua própria densidade. Isso nos mostra de cara uma informação bem comum, a maioria das espécies é rara.

01

Mas a maioria das comunidades são assim, a maioria das espécies são raras. Um cara chamado Frank Preston propôs que nós descrevemos as comunidades utilizando os logaritmos de abundâncias das espécies. Isto frequentemente revela que uma comunidade pode ser descrita aproximadamente com a distribuição normal aplicada aos dados transformados em logaritmo.

02

Veja que com um valor de média e desvio, a gente poderia, depois de descrito, reconstruir melhor os dados originais, melhor que um valor de riqueza ou diversidade por exemplo. Mas também podemos visualizar a comunidade como uma distribuição rank-abundância, atribuindo as espécies mais abundantes como rank = 1, ea menos abundante tem rank = R, em uma amostra de espécies R, então traçamos uma curva log-abundância versus rank.

03

Agora temos vários tipos de curvas que podem descrever uma comunidade seguindo essa linha de pensamento. O legal é que esse tipo de curva capta informações de riqueza e abundâncias, e certos processos, como a teoria neutra do Hubbell pode ser simulada para chegarmos a uma aparência esperada em uma curva dessa, ai começamos a confrontar dados com teoria.

Pena que diferente da diversidade, que tem funções prontas no vegan, distribuições de abundância de espécies não tem funções … , opaaaa, tem um pacote no R chamado sads que faz tudo isso pra gente.

Basta a gente converter as abundâncias das espécies na classe que o pacote usa.

> N.oc<- octav(N,preston=T) > N.oc Object of class "octav" octave upper Freq 1 0 1 9.5 2 1 2 16.0 3 2 4 18.0 4 3 8 19.0 5 4 16 30.0 6 5 32 35.0 7 6 64 31.0 8 7 128 26.5 9 8 256 18.0 10 9 512 13.0 11 10 1024 7.0 12 11 2048 2.0 13 12 4096 0.0

E os plotar os gráficos. Lembrando que o intervalo que cada classe cobre aqui é diferente. A função hist tem um método próprio de classificação de classe, os intervalos aqui são calculados de forma diferente, se eu não estou entendendo errado a função.

04

Veja que esse é equivalente ao nosso segundo gráfico, lembrando que no eixo x aqui estamos olhando valores brutos de abundância, mas a escala é logaritima. O nosso segundo gráfico acabava em 8 e e^8 vale:

> exp(8) [1] 2980.958

Perto da maior abundância ali.

E para traçar curva log-abundância versus rank, a gente ajusta a curva log-normal.

> N.ln <- fitsad(N,"lnorm") Mensagens de aviso perdidas: In (function (x, trunc, start.value, trueLL = TRUE, dec.places = 0, : trueLL used, check if the precision in your data matches the dec.places argument > summary(N.ln) Maximum likelihood estimation Call: mle2(minuslogl = LL, start = list(meanlog = meanlog, sdlog = sdlog), data = list(x = x)) Coefficients: Estimate Std. Error z value Pr(z) meanlog 3.142695 0.119146 26.377 < 2.2e-16 *** sdlog 1.787195 0.084249 21.213 < 2.2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 -2 log L: 2314.499

A gente recebe uma mensagem de aviso, mas sobre precisão. No R a gente pode receber mensagem de erros, que vai fazer a função falhar, e mensagens de aviso, onde a função funciona, mas podem haver problemas, aqui temos um aviso falando sobre precisão. Esse aviso é da propria função do pacote fitlnorm, que é chamada quando a gente usa o argumento “lnorm”. Mas tudo bem, podemos fazer os gráficos de diagnostico do modelo e temos e depois pesquisar melhor o que estamos fazendo de errado até aqui. Mas Basicamente calculamos a média e o desvio padrão da distribuição dos dados, espero que se tivessemos usado mean() o resultado deveria ser o mesmo nesse ajuste, como vimos aqui. Mas a função é importante para possibilitar o ajuste de outros modelos, que não são tão simples de fechar uma formula analiticamente. Mas depois os podemos ver os gráficos num diagnostico do modelo.

05

Todos os gráficos que estávamos fazendo e mais. Então ta tudo prontinho no pacote sads.

Muito legal, agora temos uma teoria, a teoria neutra do Hubbel, temos uma forma legal de representar de comunidade, prontinho na mão e agora o próximo passo é fazer uma ligação entre esses dois tópicos.
O script está no repositório do github, além de aqui em baixo, e até o próximo post.

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

library(sads)
library(vegan)
 
data(BCI)
N <- sort(colSums(BCI), decr = TRUE)
 
#Figura 01
hist(N, main = "Classes de densidade de espécies",ylab="Frequência",xlab="Classes de Abundância")
 
#Figura 02
hist(log(N), xlab = "Log das Classes de Abundância", main = NULL)
m.spp <- mean(log(N))
sd.spp <- sd(log(N))
R <- length(N)
curve(dnorm(x, m.spp, sd.spp) * R, 0, 8, add = T,col="red",lty=2,lwd=2)
 
#Figura 03
plot(log(N), pch=19, ylim = c(0, 8), main = NULL,xlab = "Species Rank", ylab = "Log Density",frame=F)
ranks.lognormal <- R * (1 - pnorm(log(N), m.spp,sd.spp))
lines(ranks.lognormal, log(N),lwd=3,lty=2,col="red")
legend("topright",legend="Curva Log-Normal",bty="n", log(N),lwd=3,lty=2,col="red")
 
 
N.oc<- octav(N,preston=T)
N.oc
 
#Figura 04
plot(N.oc)
 
exp(8)
 
N.ln <- fitsad(N,"lnorm")
summary(N.ln)
 
#Figura 05
par(mfrow=c(2,2))
plot(N.ln)

Teoria neutra na evolução e na ecologia

A teoria neutra, na evolução, começou em 1960 com Motto Kimura, basicamente diz que a frequência de alguns alelos podem variar por acaso devido a esses alelos terem um sucesso reprodutivo, um fitness neutro, ambos tem o mesmo valor de fitness.

Bastante simples, mas com implicações fortes. Lembrando que alelos são as formas de expressão de um gene.
Se a gente pensar na natureza e na espécie modelo mais comum em genética, a mosquinha Drosophila melanogaster que tiver, uma mutação, por exemplo, asas curtas.

08
08

Uma mosca, na natureza, sem a capacidade de voar está ferrada, ela vai morrer de fome. Ela vai ter um fitness muito menor que a mosca com asas normais e a frequência desse gene vai tender a sumir, isso é conhecido como seleção natural.

Agora talvez, outra mutações não tenham impacto direto no fitness. Estou supondo aqui, a titulo de exemplo, ja que eu não entendo de Drosophilas, mas a cor dos olhos, se não houver seleção natural qualquer, as fêmeas não verem moscas de olhos vermelhos como mais bonitas que as de olhos laranjas, o que aconteceria?

08
08

Bem essa mutação poderia ser considerada neutra, ou seja ela num melhora nem piora o fitness, mas isso não implica em nenhum momento que se você acompanhar essa população, a frequência desse alelo sera constante. Isso porque temos uma população, na natureza sempre finita, a a frequência é sempre o número de indivíduos que exibe uma característica, o alelo em questão dividido pelo tamanho total da população.

Nesse caso que temos 2 alelos em um gene (situação fictícia, eu estou assumindo que essa mutação e neutra e não existe mais variantes, sei que tem, mas é só para o exemplo ser fácil de entender), os alelos são, olhos vermelhos e olhos laranjas.

Mas esse exemplo é só para visualizar a situação, a ideia é que se esse gene é neutro, existe um processo onde se você começa cada alelo presente em 50% da população, 50 indivíduos, vamos supor, a qualquer momento que alguém usar um mata moscas aleatoriamente e mate, por cagada, 9 moscas de olhos vermelhos e 1 mosca da olhos laranjas, sobrara 41 moscas de olhos vermelhos a 49 de olhos laranjas, e se cada uma tem sei la, 1 filho, na próxima geração, veremos 45,5% de moscas de olhos vermelhos e 54.4% de moscas de olhos laranjas.

O problema é que quando eu bato o mata moscas, eu tento mata moscas, eu não fico olhando qual mosca eu vou matar, e pior que isso, nenhum tipo de moscas é mais rápido, ou mais ágil para fugir nesse exemplo, o que gera uma situação onde essa mudança na frequência dos alelos é completamente imprevisível, pra que direção ela vai.

Essa é um processo conhecido como deriva genética, ou genetic Drift, um efeito da teoria neutra.

Existem várias formas de simular a deriva genética, lembrando que a herança pode ser variável mesmo na gente. Lembre-se que mitocôndrias a gente so ganha da nossa mãe enquanto para genes nucleares alelo vem do nosso pai e outro alelo vem da nossa mãe, e tem crossing-over. Se começarmos a pensar muito nisso a gente não vai pra frente. Mas vamos usar aqui só pra visualizar o modelo do Wright–Fisher, que considera gerações discretas, indivíduos diploides, e cada nova geração tem um tamanho fixo, basicamente cada nova geração a gente pega um número aleatória de indivíduos da frequência da geração anterior.

01

Veja que cada linha é uma simulação, no eixo x temos as gerações, e no eixo x estamos representando a frequência de um alelo, lembrando que as duas frequências vão dar 1, então se um alelo tem 100% o outro tem 0, e vice versa.

Veja que em alguns casos, laranja e verde por exemplo, esse alelo domina, e o outro some, em outras simulações, ele some e só fica o outro alelo na população.

Apesar de imprevisível, a gente pode dizer que a deriva genética acaba sempre mudando as frequências de alelos, ainda mais, quase sempre sumiu com algum dos alelos, ela tirou alguém da jogada, e assim diminuiu a diversidade genética, já que ter uma população com 2 alelos é mais diversa que uma população com 1 alelo, seja ele qual for.

Certo, mas ai vendo tudo isso, um cara chamado Stephen P. Hubbell pensou se isso não poderia acontecer em outra escala, na escala de espécies. Será que existe uma teoria neutra para riqueza de espécies? Então ele pensou na teoria neutra unificada da biodiversidade e biogeografia (the unified neutral theory of biodiversity and biogeography), mas vamos dizer teoria neutra que é menor, ninguém tem saco de escrever esse nome por muito tempo.
A sua importância é que ela possibilita mais que muitos modelos de comunidade, ela faz predições claras em muitos níveis de organização, tanto de processos evolucionários como ecológicos.

Comunidades e genes são computacionalmente e matematicamente relacionadas, geralmente usando modelos iguais aos usados na genética, logo as lições que aprendemos sobre a deriva genética se aplicam a modelos de comunidade. A dinâmica ecológica numa escala local é conhecida de como deriva ecológica, e populações mudam ao acaso, devido a processos estocásticos como o nosso mata-moscas.

Hubbell propôs a teoria neutra como um modelo nulo para a dinâmica de indivíduos, para explicar a abundância de árvores tropicais ao na costa rica, barro colorado. No pacote vegan, se você abrir o pacote e digitar data(BCI), você pode ver parte dos dados que o Hubbell trabalhava, as abundâncias na floresta da reserva do Barro Colorado na costa rica, alias, como curiosidade, se você digitar data(), sem nenhum argumento, o R lista todos os data-set de exemplos disponíveis, de todos os pacotes abertos.

Certo, mas agora o controverso é que o Hubbell propôs, diferente do que comumente associamos a modelos nulos (que são como as coisas acontecem ao acaso normalmente, para teste de hipoteses), que é realmente assim que as coisas acontecem na dinâmica de comunidades.

Mas isso é controverso? Porque?

Porque estamos constantemente avaliando fatores que influenciam na riqueza e diversidade de espécies, e vem um cara e fala que tudo isso num é assim, que quem define quais espécies estão num local é o acaso. Realmente uma afirmação que vai contra muitos trabalhos que estamos acostumados a ver.

Mas se a gente não for muito extremista (e começar a achar que tudo tem que ser ou 100% estocásticos ou 100% determinísticos, sem acordos), a palavra chave para a teoria neutra do Hubbell é metacomunidade, uma coleção de comunidades similares conectadas por dispersão.

A metacomunidade é povoada inteiramente por indivíduos que são funcionalmente idênticos, segundo a teoria neutra. A teoria neutra fica então uma teoria de dinâmica de indivíduos, modelando nascimentos, morte, migração e mutação desses indivíduos, sendo a riqueza de espécies uma medida derivada desses parametros. Assumimos que dentro de uma guilda, como a de árvores do final de sucessão em floresta tropical, as espécies são essencialmente neutras, todo mundo tem um fitness igual, equivalente. Isso significa que probabilidades de nascimento, morte, mutação e migração são próximas ou equivalentes para todos os indivíduos. Assim, mudanças nos tamanhos das populações ocorrem devido a “random walks”, eventos estocásticos. Como a gente viu nas frequências gênicas. A gente pode simular praticamente igual fizemos ao gene, mas vamos usar aqui o pacote untb, que já tem essa simulação pronta.

02

Mas calma, isso não implica necessariamente na ausência de competição ou outros efeitos indiretos mediados pela dependência de densidade. A competição é vista como difusa e igual entre indivíduos.
Efeitos dependentes da densidade aparecem através de formas especificas no número total de indivíduos da comunidade, ou como características dos indivíduos relacionado as probabilidades de natalidade, mortalidade e especiação.

Um paradoxo gerado na teoria neutra de Hubbell é que na ausência de migração ou mutação, a diversidade declina para zero, riqueza para um, uma monodominância. Uma caminha aleatória (random walk), devido a equivalência de fitness, vai eventualmente resultar na perda de todas as espécies, exceto uma. Entretanto a perda de diversidade em uma comunidade local é esperada de forma muito, muito devagar, e é contra-balanceada por imigração e especiação. Assim, espécies não aumentam deterministicamente quando raras, isso faz o conceito de coexistência diferente do conceito que a gente vê nos modelos de Lotka-Volterra. Coexistência no caso da teoria neutra é um efeito estocástico balanceado pelo aparecimento de novas especies localmente. Lembrando que estamos pensando dentro de uma guilda, um nicho, não estamos comparando gramas com árvores aqui, ou onças com aranhas.

Mas se todas as comunidades fazem um random walk para a monodominâncias, como a diversidade é mantida em qualquer região em particular? Dois fatores mantém espécies em qualquer comunidade local. Primeiro, imigração para a comunidade local de indivíduos da metacomunidade. Mesmo com cada comunidade sofrendo um random walk para a monodominância, cada comunidade pode vir a ser dominada por uma espécies diferente, formando o pool de espécies dessa metacomunidade, já que todas as comunidades apresentam fitness diferentes. Assim comunidades separadas são preditas de se tornarem dominadas por espécies diferentes, e essas diferenças entre espécies locais ajudam a manter a diversidade na paisagem ocupada pela metacomunidade.

Um segundo fator, numa escala de tempo muito maior, é a especiação que dentro da metacomunidade mantém a biodiversidade. Mutação e consequentemente, especiação prove uma forma de introdução de diversidade na metacomunidade. Random walks em direção a extinção em comunidades que são tao lentas que são contrabalanceadas por especiação.

E agora temos uma pequena ideia do que se trada a teoria neutra da ecologia, e como ela veio da teoria neutra da genética. As simulações foi mais pra fazer alguma figurinha que realmente explicar algo, mas post sem figurinha é chato, e ainda as figuras foram feitas no ggplot2, que faz figuras bem bonitas, mas que até hoje eu não aprendi a usar direito. o script está no repositório do github, além de aqui em baixo, e precisamos definitivamente de muitos posts ainda pra entender a teoria do Hubbell, até onde ela vai, e como estimar parâmetros em metacomunidades para confrontar dados com teoria nesse caso aqui, mas paciência é uma virtude, sempre.

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

library(untb)
library(ggplot2)
library(reshape)
library(RColorBrewer)
set.seed(123)
 
###################################################
### Deriva Genética
###################################################
#Original - http://statisticalrecipes.blogspot.com.br/2012/02/simulating-genetic-drift.html
 
#Parametros da Simulação de Deriva Genetica
N = 20           # Número de individuos diploides
N.chrom = 2*N    # Número de cromossomos
N.gen = 100      # Número de gerações
N.sim = 10       # Número de simulações
p = 0.2
q = 1-p
 
# Simulation
saida <- matrix(0,nrow=N.gen,ncol=N.sim,dimnames = list(paste("Geraçao",1:N.gen),paste("S",1:N.sim)))
saida[1,] = rep(N.chrom*p,N.sim) # Inicializando o numero de alelos na primeira geração
for(j in 1:N.sim){
  for(i in 2:N.gen){
    saida[i,j]<-rbinom(1,N.chrom,prob=saida[i-1,j]/N.chrom)
    }
  }
#Transformando em frequencia
saida<-data.frame(saida/N.chrom)
 
# Modificando os dados para a estrutura que da para plotar
sim_data <- melt(saida)
ggplot(sim_data, aes(x = rep(c(1:100), N.sim), y = value, colour = variable)) +
    geom_line() +
    scale_colour_brewer(palette="Set3") +
    ggtitle("Simulação de deriva genética") +
    xlab("Geração") + ylab("Frequencia alélica") +
    ylim(0,1) +
    labs(colour = "Simulações")
#ggsave("01.jpg")
 
###################################################
### Deriva Ecologica
###################################################
saida <- untb(rep(1:10, 90), prob=0, gens=5000, D=9, keep=TRUE)
 
sim_data <- melt(data.frame(species.table(saida)))
 
ggplot(sim_data, aes(x = rep(c(1:5000), 10), y = value, colour = variable)) +
    geom_line() +
    scale_colour_brewer(palette="Set3") +
    ggtitle("Simulação de deriva populacional") +
    xlab("Geração") + ylab("Abundância da espécie") +
    labs(colour = "Simulações")
#ggsave("02.jpg")

Integral pela regra do trapézio

Gostando de cálculo ou não, muitas vezes é necessário calcular integrais definidas de alguma função \int_a^b \! f(x) \, \mathrm{d}x. dada uma função f. Eu acho integral um nome difícil para área embaixo de uma curva, que é o que é a integral.

A partir do teorema fundamental do cálculo, nos sabemos que se podemos achar a antiderivativa ou integral indefinida, como {F^\prime(x)}={\frac{d}{dx}}{F(x)}={f(x)}, então \int_a^b \! f(x) \, \mathrm{d}x. vai ser F(b)-F(a). Lembrando que quando isso começou, as contas eram na mão, num tinha computador, então se isso parece difícil agora, e o método que veremos mais fácil, no passado era o contrário, porque quanto menos conta melhor, eu acho.

Mas o duro é que para muitas funções f esse esquema não é possível, chegar a uma antiderivativa numa formula fechada. Isso é, nos não temos uma formula finita para F. Nesses casos (além dos normais) nos podemos integrar numericamente para uma integral definida aproximada.

Quando estamos analisando dados, usamos isso um monte de vezes sem nos dar conta, pelo menos eu demorei a visualizar isso, mas nos comumente usamos a integral da densidade da distribuição normal, entre outras distribuições.

Para a tabela de valores z, que nada mais é que a distribuição normal com média zero e desvio padrão 1, temos:

\phi(z)=\int_{-\inf}^z \frac{1}{\sqrt{2\pi}}e^\frac{-x^2}{2}dx

Nos sabemos que \phi(0)=\frac{1}{2} e \phi(\inf)=1 mas para todos os outros valores de z, integração númerica é usada, inclusive, esses tempos, o pessoal da lista brasileira de R estava discutindo isso aqui, bem interessante a discussão.

Mas continuando, existem vários meios de se aproximar um valor de integral, e um método para integração numérica é a aproximação trapeizodal, que é obtida aproximando a área de y=f(x) sobre um subintervalo \{x_i,x_{i+1}\} de um trapézio, ou seja, a função f(x) é aproximada por uma linha reta no intervalo \{x_i,x_{i+1}\}, a largura do trapézio é h, o lado direito tem altura f(x_i) e o lado direito tem altura f(x_{i+1}). A área do trapézio é então.

\frac{h}{2}(f(x_i)+f(x_{i+1}))

Quem não se lembra do trápezio, é essa figurinha aqui:

Trapezio

Lembre-se que estamos pensando nela tipo em pé, girando pra esquerda.

Certo, agora a gente soma todos os subintervalos juntos e temos nossa aproximação da integral da função f(x) no intervalo que queríamos \int_a^b \! f(x) \, \mathrm{d}x.

Regra do trapézio.
T=\frac{h}{2}(f(x_0)+2f(x_1)+2f(x_2)+,\dots,+2f(x_{n-1})+f(x_n))

Observe que, para i=1,\dots,i=n-1, f(x_i) contribui para a área do trapézio para a esquerda de x_i e para a direita de x_i e assim aparece multiplicada por 2 na regra do trapézio. Em contraste f(x_0) e f(x_n) contribuem apenas para a área do primeiro e última trapézio, respectivamente.

Então vamos tentar visualizar melhor a situação.
Temos então uma função f qualquer, vamos supor essa aqui.

f(x)= 4 \cdot x^3

No R a gente vai escrever essa função assim então:

ftn <- function(x) return(4 * x^3)

Ai a pra ficar mais fácil de entender, a gente consegue agora colocar ela num gráfico, usando a função curve() do R.

01

Agora o que é uma integral? Integral é calcular uma área. A gente sabe calcular área de alguma figuras, sei la, um quadrado é o lado ao quadrado, um retângulo tem área base vezes altura e podemos ficar olhando o wikipedia por exemplo, várias figuras geométricas e ficar vendo formulas fechadas para as suas áreas.

Mas e se a gente quiser calcular a área sobre essa curva no intervalo de 0 a 1. Vamos colorir a figura com a área desejada pra ver melhor.

02

Não tem uma formula fechada pra isso, só que integrar é fazer isso, é calcular a área embaixo dessa curva, dessa função, que começou com ou o Newton, ou o Leibniz, ja que existe uma controvérsia ai que não é o caso agora.

Mas como a gente disse ali, a gente pode usar a regra do trapézio, que é só desenhar um monte de trapézios embaixo da curva, calcular a área de cada um e somar todas as áreas, claro que precisamos sistematizar como criar esses trapézios.

No caso da curvinha ali em cima:

03

Veja que eu desenhei 5 trapézios ali embaixo da curva. Bem simples, as duas alturas você consegue resolvendo a função f(x) para dois pontos onde começa e termina cada trapézio, a base vai ser a diferenças entre esses dois pontos e pronto (sempre usamos o mesmo valor de h, partimos a área em trapézios de h iguais, uma sequencia de a a b indo em unidades de h) e assim temos todos os ingredientes para a área de cada trapézio a mão.

{A}={\frac{(B+b)}{2} \cdot h}

Então podemos fazer vários trapézios sistematicamente sobre a curva, e a partir da regra do trapézio, conseguir a área, ou integral, que nada mais vai ser que calcular a área de todos esses trapézios e somar.

Então fazemos nossa função para tal:

trapezoid <- function(ftn, a, b, n) {
  h <- (b-a)/n
  x.vec <- seq(a, b, by = h)
  f.vec <- sapply(x.vec, ftn)
  T <- h*(f.vec[1]/2 + sum(f.vec[2:n]) + f.vec[n+1]/2)
  return(T)
}

Então a entrada vai ser uma função, do tipo que você entra um valor e sai o resultado, aqui a gente ta usando a f(x)= 4 \cdot x^3 (diga-se de passagem que da pra calcular a antiderivativa aqui, mas, vamos usar uma função simples ajuda a entender que ta acontecendo).

Ai a gente usa a função e pimba, temos a área aproximada para o intervalo de 0 a 1.

> trapezoid(ftn, 0, 1, n = 5) [1] 1.04

Veja que a função é bem simples, todos os trapézios vão ficar com o mesmo h, que vai depender de quantos trapézios estamos usando.

  h <- (b-a)/n

Dai a gente faz uma sequência de valores começando em a, terminando em b, numa distância simétrica de h, que vai dar n+1 valores, e resolve a função f(x) pra todos eles. Veja que a gente ta resolvendo a função usando o sapply, da família dos apply do R, para discutir, melhor fazer um post a parte.

  x.vec <- seq(a, b, by = h)
  f.vec <- sapply(x.vec, ftn)

E por último é só lembrar a formula la em cima da regra do trapézio, é só implementar ela, soma o primeiro o último e todos outros valores, vezes o h.

  T <- h*(f.vec[1]/2 + sum(f.vec[2:n]) + f.vec[n+1]/2)

Claro que esse não é um método exato, então conforme a gente aumenta o número de trapézios, aumentamos nossa precisão, diminuímos o erro em relação ao número real, mas aumentar os trapézios implica em fazer mais contas, o que demanda mais do computador, claro que nesse exemplo é trivial, mas dependendo da função pode não ser tão trivial assim.

04

Então se a gente quer mais precisão a gente usa mais trapézios.

> trapezoid(ftn, 0, 1, n = 2) [1] 1.25 > trapezoid(ftn, 0, 1, n = 40) [1] 1.000625 > trapezoid(ftn, 0, 1, n = 60) [1] 1.000278 > trapezoid(ftn, 0, 1, n = 1000) [1] 1.000001

Veja que o resultado está tendendo a estacionar em 1.
Claro que chega uma hora que mesmo aumentando o número de trapézios começamos a mudar muito pouco o valor, que é a hora de parar, se tivéssemos a possibilidade de usar infinitos trapézios, teríamos o valor exato (isso fica no mundo da imaginação, e do calculo), mas como mesmo os números guardados no computador tem um limite de precisão, chega uma hora que é inútil aumentar o número de trapézios. Podemos ver no gráfico abaixo o efeito de aumentar o número de trapézios no resultado, como vemos acima, que fica cada vez mais perto de 1.

05

Mas não vamos nos preocupar com esses pontos, o legal é ver que funciona.
No R existe a função chamada integrate(), que calcula integrais, funciona de forma semelhante a nossa função aqui, mas ele usa outros métodos de integrar numericamente, e escolhe o melhor método pra você, mas vamos usar ela pra conferir se não estamos fazendo nada errado, se estamos no caminho certo.

> integrate(ftn,0,1) 1 with absolute error < 1.1e-14

E pimba, veja que no nosso gráfico la em cima, quanto mais trapézios, a gente estava tendendo a estacionar no valor de 1, e a função do R integra exatamente para o valor de 1, olha ai nossa lógica não esta tão errada assim.

Mas essa função, como a maioria dos exemplos pra tudo, é bem nada a ver, não tem sentido nenhum. Pelo menos sentido que a gente vá usar em biologia, mas essa aqui a gente vê direto.

06

Essa curva nada mais é que essa formula:

\phi(z)=\int_{-\inf}^z \frac{1}{\sqrt{2\pi}}e^\frac{-x^2}{2}dx

Que a gente pode implementar no R assim:

ftn <- function(x) return(1/sqrt(2*pi)*exp(-x^2/2))

Se você for olhar a formula da densidade da distribuição normal, lembre-se que ai em cima, o desvio padrão é 1, logo 1 vezes alguma coisa é a própria coisa, então os valores de 1 do desvio padrão estão omitidos, o \sigma da formula no wikipedia aqui.

Outra detalhe, vamos chamar de ftn denovo pela falta de criatividade para nomes.

Mas então podemos perguntar qual o a área entre um intervalo, vamos supor -1.96 e 1.96.

07

E a regra do trapézio continua valendo, é só a gente usar nossa nova função.

08

> trapezoid(ftn, -1.96, 1.96, n = 20) [1] 0.9492712

Olha ai, quase 0.95, que é o famoso intervalo de confiança de 95% da distribuição normal. Por um tempo, quando ia fazer aqueles gráficos de barras, e ia colocar o tamanho das barras, eu nunca entendia porque usava 1.96 e não 2 na hora de multiplicar o desvio padrão. Mas a área de 95% fica sobre -1.96 e 1.96, sendo que assim a área é 95%, 0.95, e o meio é exatamente o zero.
A gente pode conferir com a função integrate do R

> integrate(ftn,-1.96,1.96) 0.9500042 with absolute error < 1e-11

Além disso, essa função normal já esta escrita no R, a função dnorm nada mais que a função acima, mas bem mais caprichada, mas o default dela é media zero e desvio 1, como usamos acima, podemos conferir integrando ela nesse intervalo e vendo o resultado.

> integrate(dnorm,-1.96,1.96) 0.9500042 with absolute error < 1e-11

Alias essa é uma característica legal de toda distribuição de probabilidade, todas, sempre integram pro valor de 1, 100% de chance, embaixo da curva sempre tem todas as possibilidades, e a gente sempre fica olhando a área que nos interessa embaixo dessas curvas.

> integrate(dnorm,-Inf,Inf) 1 with absolute error < 9.4e-05

Bem, acho que é isso ai por hoje. Ficamos por aqui. Os scripts sempre no repositório recologia e aqui nesse link tem um reference card bem legal de estatística, com as principais formulas, um bom exercício é ficar implementando e plotando eles, e é isso, até o próximo post, se você leu até aqui e me viu escrevendo alguma abobrinha, por favor não deixe de me avisar do meu erro por e-mail ou fazendo um comentário.

Referencia:
Owen Jones, Robert Maillardet & Andrew Robinson 2009 – Introduction to Scientific Programming and Simulation Using R. CRC Press 499pp

trapezoid <- function(ftn, a, b, n = 100) {
  h <- (b-a)/n
  x.vec <- seq(a, b, by = h)
  f.vec <- sapply(x.vec, ftn)
  T <- h*(f.vec[1]/2 + sum(f.vec[2:n]) + f.vec[n+1]/2)
  return(T)
}
 
ftn <- function(x) return(4 * x^3)
 
#figura 1
curve(ftn(x),-2,2)
 
#figura 2
curve(ftn(x),-1,1.5,lwd=2)
polygon(c(seq(0,1,length=1000),rev(seq(0,1,length=1000))),c(ftn(seq(0,1,length=1000)),rep(0,1000)),col="red",lwd=0.5)
abline(h=0,lty=2,lwd=2)
 
#figura 3
a=0
b=1
n=5
h <- (b-a)/n
x.vec <- seq(a, b, by = h)
f.vec <- sapply(x.vec, ftn)
curve(ftn(x),-1,1.5,lwd=1,lty=3,add=T)
for(i in 2:(n+1)) {
    polygon(c(x.vec[i],x.vec[i],x.vec[i-1],x.vec[i-1]),c(0,f.vec[i],f.vec[i-1],0),lwd=2)
}
abline(h=0,lty=2,lwd=2)
 
#figura 4
n=20
h <- (b-a)/n
x.vec <- seq(a, b, by = h)
f.vec <- sapply(x.vec, ftn)
curve(ftn(x),-1,1.5,lwd=1,lty=3,add=T)
for(i in 2:(n+1)) {
    polygon(c(x.vec[i],x.vec[i],x.vec[i-1],x.vec[i-1]),c(0,f.vec[i],f.vec[i-1],0),lwd=2)
}
abline(h=0,lty=2,lwd=2)
 
#
trapezoid(ftn, 0, 1, n = 2)
trapezoid(ftn, 0, 1, n = 40)
trapezoid(ftn, 0, 1, n = 60)
trapezoid(ftn, 0, 1, n = 1000)
 
#figura 5
plot(0,0,type="n",xlim=c(0,14),ylim=c(0.9,1.3),frame=F,xaxt="n",xlab="Número de Trapézios",ylab="Valor total aréa")
pontos<-c(2,5,10,seq(100,1000,by=100))
for(i in 1:13) {
    points(i,trapezoid(ftn, 0, 1, n = pontos[i]))
}
axis(1,at=1:13,label=pontos)
 
#
integrate(ftn,0,1)
 
#figura 6
curve(1/(1*sqrt(2*pi))*exp(-x^2/(2*1^2)),-4,4,add=T,lty=2,lwd=2)
 
#figura 7
ftn <- function(x) return(1/sqrt(2*pi)*exp(-x^2/2))
curve(ftn(x),-4,4,lwd=2)
polygon(c(seq(-1.96,1.96,length=1000),rev(seq(-1.96,1.96,length=1000))),c(ftn(seq(-1.96,1.96,length=1000)),rep(0,1000)),col="red",lwd=0.5)
abline(h=0,lty=2,lwd=2)
 
#figura 8
a=-1.96
b=1.96
n=20
h <- (b-a)/n
x.vec <- seq(a, b, by = h)
f.vec <- sapply(x.vec, ftn)
T <- h*(f.vec[1]/2 + sum(f.vec[2:n]) + f.vec[n+1]/2)
curve(ftn(x),-4,4,lwd=2,lty=3,col="red")
for(i in 2:n) {
    polygon(c(x.vec[i],x.vec[i],x.vec[i-1],x.vec[i-1]),c(0,f.vec[i],f.vec[i-1],0))
}
 
trapezoid(ftn, -1.96, 1.96, n = 20)
integrate(ftn,-1.96,1.96)
integrate(dnorm,-1.96,1.96)
integrate(dnorm,-Inf,Inf)