Ao contrário de contagens, as medidas contínuas podem assumir qualquer valor dentro de um intervalo numérico. São exemplos extensões, massa e tempo decorrido. Dizemos então que uma variável é contínua se pode ter qualquer valor no conjunto de números reais, ou um subconjunto deles. As variáveis contínuas necessitam de modelos de distribuição de probabilidade distintos das variáveis discretas.

Nesse tópico apresentamos esses modelos e exploramos suas propriedades básicas, com simulações dos cenários probabilísticos em computador. Também vamos verificar que a soma de muitas fontes de variação no limite cria a distribuição contínua mais conhecida, a Gaussiana.

• Função de densidade probabilística
• Função de probabilidade acumulada
• Esperança e variância de variáveis contínuas
• Teorema Central do Limite
• Distribuições de probabilidade contínuas mais usadas em ecologia

É importante ter clara a diferença entre a função de densidade de probabilidade de uma variável contínua e a probabilidade dessa variável assumir um valor dentro de um intervalo.

A primeira não é uma probabilidade, e pode ter valores maiores do que um. No caso de uma variável normal, isto acontece quando o desvio-padrão é baixo.

A função “dnorm“ no R retorna o valor da função de densidade da normal para um dado valor da variável aleatória. Para uma normal com parâmetros $\mu=2$ e $\sigma=0,1$ os valores da função de densidade para $\mu \pm 1,5\sigma$ são:

> dnorm(2+c(-0.15,0.15),mean=2,sd=0.1)
[1] 1.295176 1.295176

Aplicando a função ”curve“ à função “dnorm” podemos traçar o gráfico da função dentro de um intervalo dado pelos argumentos “from” e “to”. Vamos fazer isto para a variável normal definida acima, no intervalo $\mu \pm 3\sigma$

curve(dnorm(x,mean=2,sd=0.1),from=1.7, to=2.3, 
      xlab="x",ylab="densidade de probabilidade")

Vamos acrescentar neste gráfico os valores que calculamos acima:

points(2+c(-0.15,0.15),dnorm(2+c(-0.15,0.15),mean=2,sd=0.1),
       col="red", pch=19, cex=1.5)
segments(x0=c(2+c(-0.15,0.15,0.15)),
         y0=c(0,0,dnorm(2.15,mean=2,sd=0.1)),
         x1=c(2+c(-0.15,0.15),1.5),
         y1=dnorm(2+c(-0.15,0.15,0.15),mean=2,sd=0.1),
         col="red", lty=2)

Para uma variável aleatória a probabilidade só pode ser definida para um intervalo, e corresponde à área sob a curva que está delimitada por este intervalo. Vamos criar uma função no R para preencher esta área no gráfico da curva normal ¹⁾:

area.norm <- function(from,to,mean,sd,cor="red",shade=1e4,...){
  len <- (to-from)*shade
  x <- seq(from,to,length=len)
  segments(x0=x,x1=x,y0=rep(0,len),y1=dnorm(x,mean=mean,sd=sd),col=cor)
}

E agora aplicá-la para colorir de vermelho a área delimitada por $\mu \pm 1,5\sigma$ no gráfico que criamos:

     area.norm(1.85,2.15,mean=2,sd=0.1)

A área sob uma função é sua integral. Podemos calculá-la com a função de integração do R ²⁾:

     integrate(dnorm,lower=1.85, upper=2.15,mean=2,sd=0.1)

A integral de qualquer função de densidade probabilística no intervalo que contém todos os valores da variável aleatória é um. Verifique isto para a normal, integrando-a de menos infinito ($-\infty$) a mais infinito ($\infty$):

integrate(dnorm,lower=-Inf, upper=Inf,mean=2,sd=0.1)
integrate(dnorm,lower=Inf, upper=Inf) ## normal padronizada

A integral do limite inferior de uma função de densidade até um certo valor é a função de densidade acumulada. Para a normal, a função no R é “pnorm”. A área entre dois valores pode ser obtida subtraindo-se suas probabilidades acumuladas:

<box centered 70% red|A probabilidade de um intervalo é a diferença de suas probabilidades acumuladas> </box>

O que você obtém com o comando:

    pnorm(2.15,mean=2,sd=0.1)-pnorm(1.85,mean=2,sd=0.1)

A funções no R para calcular probabilidades acumuladas, como a “pnorm”, têm o argumento “lower.tail”, que permite calcular a probabilidade a partir de um ponto, isto é, da cauda posterior ao valor da variável (“lower.tail=FALSE”). Com isto, outra maneira de obter a probabilidade para o mesmo intervalo é:

    1-(pnorm(2.15,mean=2,sd=0.1,lower.tail=F)+pnorm(1.85,mean=2,sd=0.1))

Por fim, como o intervalo é simétrico em torno da média, e a curva normal também o é, uma terceira maneira de fazer o cálculo é:

    1-2*pnorm(1.85,mean=2,sd=0.1,lower.tail=TRUE)

Se quiser reproduzir a figura deste tutorial, o código em R é

png("dif_normal.png",height=200,width=600)
par(mfrow=c(1,3))
curve(dnorm(x,mean=2,sd=0.1),from=1.7, to=2.3, 
      xlab="",ylab="",axes="F")
area.norm(1.7,2.15,mean=2,sd=0.1,cor="black")
curve(dnorm(x,mean=2,sd=0.1),from=1.7, to=2.3, 
      xlab="",ylab="",axes="F")
area.norm(1.7,1.85,mean=2,sd=0.1,cor="black")
curve(dnorm(x,mean=2,sd=0.1),from=1.7, to=2.3, 
      xlab="",ylab="",axes="F")
area.norm(1.85,2.15,mean=2,sd=0.1,cor="black")
mtext(c("-","="),at=c(0.6,1.6),padj=1.1,cex=8)
dev.off()

A esperança de uma variável aleatória contínua é

$$ E[X]\ = \ \int x \cdot f(x) \ dx $$

Onde $f(x)$ é a função de densidade probabilística.

Vamos verificar isto numericamente para a variável exponencial, cuja média é o inverso de seu único parâmetro, $\lambda$. Começamos definindo uma função no R que é o produto da densidade pelo valor da variável aleatória:

esper.exp <- function(x,taxa){
  dexp(x,rate=taxa)*x
}

Agora usamos a função “integrate” para fazer a integração numérica desta função no suporte da exponencial (0 a infinito), para uma variável exponencial com $\lambda=0.1$:

 integrate(esper.exp,0,Inf,taxa=1/10) 

Compare o valor obtido com o valor teórico da esperança, $ \lambda^{-1} $ .

Para comparar a distribuição de dados com alguma variável aleatória, podemos usar um histograma re-escalonado para densidade.

Para exemplificar, usaremos os perímetros à altura do peito de fustes de caixeta amostradas em um inventário florestal (caixeta.csv).

Primeiro criamos um objeto com os dados de interesse, que são os fustes da caixeta em Chauás:

caixeta <- read.csv2("caixeta.csv")
chauas <- caixeta[caixeta$local=="chauas"
                  &caixeta$especie=="Tabebuia cassinoides",]

E então criamos um fator para as classes de CAP, com intervalos de 100 mm:

chauas$cap.class <- cut(chauas$cap,breaks=seq(100,800,by=100))

Um gráfico da tabulação deste fator é um histograma de frequências:

barplot(table(chauas$cap.class),space=0)

Dividindo-se os valores pelo total de fustes temos um histograma de frequências relativas:

barplot(table(chauas$cap.class)/length(chauas$cap.class))

Finalmente, dividindo novamente estes valores pelas amplitude do intervalo de classe usado (100 mm) temos um histograma de densidade:

barplot(table(chauas$cap.class)/(length(chauas$cap.class)*100),
        space=0)

Neste histograma, a soma das áreas das barras é um. A função ”hist“ tem um argumento “prob” que converte a escala para densidade:

hist(chauas$cap,breaks=8,prob=T,
     xlab="CAP (cm)",ylab="Densidade")

Nesta escala podemo sobrepor funções de densidade, como a Gaussiana:

 curve(dnorm(x,mean=mean(chauas$cap),
      sd=sd(chauas$cap)),add=T,col="blue") 

As companhias aéreas norte-americanas garantem que 98% de seus passageiros não terão problemas para se acomodar nos assentos da classe econômica. Para isso, usam um modelo que descreve a distribuição das larguras dos quadris dos homens norte-americanos como uma variável normal, com média de 14,4 polegadas, e desvio-padrão de 1,0 polegada.

Faça um gráfico desta distribuição de probabilidades:

curve(dnorm(x,mean=14.4,sd=1),from=8, to=20, 
      xlab="Largura (in)",ylab="densidade de probabilidade")

Para definir a largura do assento, é preciso saber o limite superior do intervalo que contêm 98% da população. Esse é o quantil de 98% da distribuição, no caso uma normal com média de 14,4 e desvio-padrão de 1,0 polegadas.

Calcule o valor deste quantil com o comando:

(larg <- qnorm(mean=14.4,sd=1,p=0.98))

Este é o valor que delimita 98% da área sob a curva normal. Usando a função criada no tutorial anterior podemos visualizar no gráfico a área que vai até o quantil calculado acima:

area.norm(8,larg,mean=14.4,sd=1)

Imagine que o desvio-padrão para o exemplo dos quadris norte-americanos seja 50% maior. O que deve ocorrer com o quantil de 98%? Verifique sua expectativa com o comando:

qnorm(mean=14.4,sd=1.5,p=0.98)

Repita o gráfico do tutorial anterior, e adicione ao gráfico a distribuição de probabilidades com o novo valor do parâmetro desvio-padrão:

curve(dnorm(x,mean=14.4,sd=1),from=8, to=20, 
      xlab="Largura (in)",ylab="densidade de probabilidade")
curve(dnorm(x,mean=14.4,sd=1.5), add=T, col="blue")

Os parâmetros da distribuição normal correspondem à média e ao desvio-padrão. Explore o comportamento da curva com a mudança de seus parâmetros:

cores <- rainbow(n=4)
curve(dnorm(x),from=-8,to=8, col=cores[1])
curve(dnorm(x,sd=2), col=cores[2], add=T)
curve(dnorm(x,mean=2), col=cores[3], add=T)
curve(dnorm(x,mean=2,sd=2), col=cores[4], add=T)
legend(-7.5,0.3,legend=c("média=0,sd=1","média=0,sd=2",
                         "média=2,sd=1","média=2,sd=2"),lty=1,col=cores)

O Teorema Central do Limite (TCL) prova que as médias de amostras independentes tomadas de uma mesma distribuição tendem para uma distribuição normal. A média desta normal é a mesma da distribuição de onde vieram as amostras, e a variância é igual à variância da distribuição original, dividida pelo tamanho das amostras.

Verifique isto com o seguinte tutorial:

Sorteie dez mil números de uma distribuição muito diferente de uma normal, como uma Poisson com média baixa:

vals <- rpois(10000,lambda=0.5)

Agora imagine que estes valores foram obtidos de mil amostras independentes, cada uma de tamanho dez. Vamos criar um fator que rotule as amostras de um a mil:

cod.amostras <- factor(rep(1:1000,each=10))

Agora podemos calcular a média de cada uma dessas amostras:

medias.vals <- tapply(vals,cod.amostras,mean)

E fazer um histograma dessas mil médias:

hist(medias.vals,prob=T)

E sobrepor a função de densidade da normal, com os parâmetros previstos pelo TCL:

curve(dnorm(x,mean=0.5,sd=sqrt(0.5/10)),add=T,col="blue")

Uma outra ótima maneira de avaliar a normalidade é com o gráfico de quantis normais ³⁾:

qqnorm(medias.vals)
qqline(medias.vals)

Repita os procedimentos, com mil amostras de tamanho 50. Quais suas conclusões?

Uma maneira simples de formular o Teorema Central do Limite é que a soma ou a média de muitas variáveis aleatórias tende a uma variável normal. Esta tendência é válida mesmo para variáveis com distribuições muito diferentes da normal e que não tenham os mesmos parâmetros.

Neste tutorial investigaremos se o TCL pode nos ajudar a propor um modelo para a riqueza de espécies em um conjunto de áreas. Vamos simular um cenário em que a ocorrência de cada espécie é uma variável Bernoulli (1 = presente, 0 = ausente), cujo único parâmetro (p a probabilidade de ocorrência) varia entre espécies.

Vamos assumir que as probabilidades de ocorrência também são variáveis aleatórias. Uma escolha natural é a variável Beta, que está restrita a valores entre zero e um, e cuja distribuição é muito flexível.

<box centered 90% red|Função de densidade da distribuição beta usada para representar as probabilidades de ocorrência das espécies. Os parâmetros usados são a=1 e b=5.> </box>

Para simular uma metacomunidade com 300 espécies, sorteamos 300 valores da distribuição beta representada acima:

p.oc <- rbeta(300, shape1=1,shape2=5)

Em seguida criamos uma função que faz N experimentos de Bernoulli, tomando um valor 1 com probabilidade $p$ e valor $0$ com probabilidade 1-p

ocorre <- function(p,N) { sample( c(0,1), size=N, replace=T, prob=c(1-p,p)) }

Aplicando esta função às probabilidades de ocorrência sorteadas no passo anterior, obtemos uma matriz com 100 linhas (locais) e 300 colunas (espécies):

matriz <- sapply(p.oc, ocorre, N=100)

A soma de cada linhas é o número de espécies em cada um dos 100 locais

riquezas <- apply(matriz, 1, sum)

Cuja distribuição pode ser comparada com função de densidade de uma normal de mesma média e desvio-padrão:

hist(riquezas,prob=T,xlab="N de espécies",
     ylab="Densidade de Probabilidade")
curve(dnorm(x,mean=mean(riquezas),sd=sd(riquezas)),
      add=T,col="blue")

Avalie também o ajuste a uma normal com o gráfico de quantis teóricos:

qqnorm(riquezas)
qqline(riquezas)

Se quiser reproduzir o gráfico da variável beta deste tutorial o comando é:

curve(dbeta(x,shape1=1,shape2=5),0,1,
      xlab="x",ylab="Densidade de Probabilidade")

A variável exponencial é a análoga contínua da variável geométrica. Ela pode ser usada para descrever o tempo de espera até a primeira ocorrência de um evento, dado que ele tem uma taxa de ocorrência $\lambda$ constante.

Vamos simular esta situação, imaginando que você acompanha um animal por até duas horas⁴⁾, sempre no mesmo horário. Sua variável de interesse é o tempo até que ele apresente um certo comportamento pela primeira vez. A taxa média com que este comportamento ocorre neste horário é de 0,1/minuto.

A distribuição desta variável aleatória pode ser simulada repetindo o estudo muitas vezes. A cada uma delas haverá um certo número de comportamentos, com uma média de $0.1 \times 120 = 12$ vezes. Supondo que estes eventos estão distribuídos aleatoriamente entre os dias, o número de eventos a cada dia de observação será uma variável Poisson. Para simular 10.000 repetições, sorteamos esta quantidade de valores de uma Poisson com taxa de ocorrência por dia de doze no período de observação :

n <- rpois(10000,12)

O total de eventos será então a soma destes valores sorteados. Em seguida sorteamos este número de valores de uma distribuição uniforme, com o mínimo de 0 e máximo de 120. Estes serão os horários de ocorrência de cada evento, em cada repetição do estudo ⁵⁾.

vals <- runif(sum(n),0,120)

E criamos um fator para identificar cada repetição do estudo, ou seja, cada dia de observação:

cod.amostras <- factor(rep(1:10000,n))

Até aqui temos um vetor com tempo transcorrido até cada um dos eventos que ocorreram em cada dia de observação. Com isto, podemos usar a função ”tapply“ para identificar o menor valor de cada repetição, que é o tempo da primeira ocorrância do evento em cada dia:

primeiro <- tapply(vals,cod.amostras,min)

E agora construímos um histograma de densidade desta variável, e sobrepomos a função de densidade da variável exponencial:

hist(primeiro,prob=T,  xlab="Tempo de Espera (min)", 
     ylab="Densidade Probabilística")
curve(dexp(x,rate=12/120),add=T,col="blue")

Avalie visualmente a concordância com o modelo teórico.

A variável Gama foi criada como uma extensão da exponencial, para descrever o tempo de espera até que um certo número de eventos ocorram, dada uma taxa constante de ocorrência. Devido à sua flexibilidade, posteriormente foi adotada como modelo heurístico para descrever variáveis com distribuições de probabilidade assimétricas.

Usando a simulação do tutorial anterior, podemos obter os tempos de espera até que o animal apresente o comportamento de interesse pela segunda vez. Para isto, criamos uma função que retorna o n-ésimo menor valor de um vetor de números

enesimo <- function(x,n)sort(x)[n]

E aplicamos esta função a cada repetição do estudo, guardada no vetor cod.amostras, do tutorial anterior:

segundo <- tapply(vals,cod.amostras,enesimo,n=2)

Para então fazer um histograma de densidade com estes dados

hist(segundo,prob=T, xlab="Tempo de Espera (min)", 
     ylab="Densidade Probabilística")

Podemos então sobrepor a este histograma a curva da função de densidade probabilística da variável Gama, com os parâmetros teóricos de taxa de ocorrência (“rate”) e número de eventos a esperar (“shape”):

curve(dgamma(x,rate=12/120,shape=2),add=T,col="blue")

A distribuição Weibull também é usada para descrever tempo contínuo de espera até que um evento aconteça. O seu parâmetro de forma determina como a taxa instantânea de mortalidade muda com o tempo decorrido. Se isso acontece, dizemos que a distribuição descreve um processo com memória. Para entender esta propriedade, vamos primeiro verificar que ela não existe na distribuição exponencial.

Suponha um processo de mortalidade em tempo contínuo, a uma taxa constante⁶⁾ de $\lambda = 0.5 \; \text{ano}^{-1}$. Neste caso o tempo de vida segue uma distribuição exponencial.

Neste exemplo, qual será a probabilidade de viver mais do que um ano? Para responder isso, lembre-se que a probabilidade acumulada da exponencial até um valor $x$ é a integral da função de densidade probabilística de zero até $x$:

$$P(t \leq x) \,=\, \int_0^x \lambda e^{- t \lambda} \; dx\, = \, 1-e^{-\lambda x}$$

e portanto a função de probabilidade acumulada a partir deste valor é seu complemento para a unidade:

$$P(t \geq x) \,=\, 1-P(t \leq x) \,=\, e^{-\lambda x}$$

Portanto, em nosso exemplo $P(t \geq 1) = e^{-lambda}$. Podemos obter este valor no R com a função de densidade acumulada pexp, usando o parâmetro lower.tail=TRUE para indicar que queremos a probabilidade acumulada a partir do quantil de $x=1$:

(p1 <- pexp(1, rate=0.5, lower.tail=FALSE))

Qual a probablilidade que esta distribuição atribui a sobreviver um ano e meio ou mais?

(p1.5 <- pexp(1.5, rate=0.5, lower.tail=FALSE))

Com estes dois valores podemos calcular a probabilidade de sobrevivência por um ano e meio ou mais, dado que sobreviveu-se até um ano. Esta é uma probabilidade condicional:

$$P(t\geq1.5 \,| \,t\geq1)$$

Que obtemos divindo a condicionada pela condicionante:

p1.5/p1

Compare este valor com a probabilidade de viver meio ano ou mais:

pexp(.5, rate=0.5, lower.tail=FALSE)

Os dois valores são idênticos! Isso quer mostra que a exponencial não tem memória: a probabilidade de viver por um tempo adicional $\Delta$ não é afetada pelo tempo já vivido. Em linguagem matemática:

$$P(t\geq x + \Delta\,| \,t\geq x) \ = \ P(t\geq \Delta)$$

Como explicamos acima, a probabilidade de se viver mais do que um certo valor $x$ de tempo é

$$P(t \geq x) \,=\, 1-P(t \leq x) \,=\, e^{-\lambda x}$$

e a probabilidade condicional de viver mais que um tempo adicional $\Delta$ dado que viveu mais que o tempo $x$ é a razão entre probabilidade de viver mais que $x+\Delta$ dividido pela probabilidade de viver mais que $x$:

$$ P(t\geq x + \Delta\,| \,t\geq x) \, = \, \frac{P(t \geq x + \Delta)}{P(t \geq x)} $$

Usando nesta razão a expressão da distribuição de probabilidade acumulada da exponencial temos:

$$ \frac{P(t \geq x + \Delta)}{P(t \geq x)} \, = \, \frac{e^{-\lambda (x + \Delta) }}{e^{-\lambda x}} \, = \, e^{-\lambda (x+\Delta)} \times e^{\lambda x } \, = \, e^{-\lambda x - \lambda \Delta + \lambda x} \, = \, e^{- \lambda \Delta} $$

Ou seja, sob a distribuição exponencial, a probabilidade de viver mais que $x + \Delta$ é igual a probabilidade de viver mais que $\Delta$, que é $e^{- \lambda \Delta}$. Neste modelo, a chance de viver um tempo a mais não é afetada pelo quanto já se viveu.

O parâmetro de forma da distribuição Weibull define o tipo de memória que ela pode descrever. Se o parâmetro é igual a um, a Weibull reduz-se a uma exponencial e portanto não tem memória. Verifique com⁷⁾:

## Probabilidade de viver mais de 1 ano
(pw1 <- pweibull(1, shape=1, scale=2, lower=FALSE))
## P de viver mais de 1 ano e meio
(pw1.5<- pweibull(1+ 0.5, shape=1, scale=2, lower=FALSE))
## Probabilidade de viver mais meio ano, dado que viveu um
pw1.5/pw1
## Probabilidade de viver meio ano no inicio da vida
pweibull(0.5, shape=1, scale=2, lower=FALSE)

Quando o parâmetro de forma é maior que um cria-se um efeito de memória que pode ser interpretada como envelhecimento:

(pw1 <- pweibull(1, shape=1.5, scale=2, lower=FALSE))
(pw1.5<- pweibull(1+ 0.5, shape=1.5, scale=2, lower=FALSE))
pw1.5/pw1
pweibull(0.5, shape=1.5, scale=2, lower=FALSE)

Quando o parâmetro de forma é menor do que um temos uma memória que pode descrever mortalidade mais intensa de jovens:

(pw1 <- pweibull(1, shape=0.2, scale=2, lower=FALSE))
(pw1.5<- pweibull(1+ 0.5, shape=0.2, scale=2, lower=FALSE))
pw1.5/pw1
pweibull(0.5, shape=0.2, scale=2, lower=FALSE)

A função de risco, ou taxa de falha, é outra maneira de descrever o tipo de memória de distribuições de tempo. A função de risco é a razão entre a densidade de probabilidade em um tempo e a probabilidade acumulada a partir daquele tempo:

$$h(t) \ = \ \frac{f(t)}{1-F(t)}$$

Esta função expressa como o risco de morte ou falha muda com o tempo já decorrido, ou idade. Vamos fazer um gráfico das funções de risco das três distribuições Weibull usadas acima:

## Funcao de risco
hweib <- function(x, shape, scale){
    dweibull(x, shape, scale)/
        pweibull(x,shape, scale, lower=FALSE)
}
## Grafico
curve(hweib(x, shape=1.5, scale=2),0,3, xlab="Tempo decorrido",
      ylab="Função de risco", col="red")
curve(hweib(x, shape=0.2, scale=2), add=T, col="blue")
curve(hweib(x, shape=1, scale=2), add=T)
legend(x=2.5, y=0.3, c("a > 1", "a < 1", "a = 1"),
       lty=1, col=c("red", "blue", "black"), bty="n")

Assim como a variável normal descreve a soma ou média de muitas variáveis aleatórias, a variável log-normal descreve a multiplicação de variáveis aleatórias.

O exemplo mais conhecido de processo multiplicativo em ecologia é o modelo de crescimento populacional geométrico estocástico.

Neste modelo, o tamanho da população no tempo t+1 é o produto do tamanho da população no tempo anterior, t, multiplicado pela taxa de crescimento populacional R:

$$ N_{t+1} \ = \ R \cdot N_t $$

Se R varia ao longo do tempo, podemos descrevê-la como uma variável aleatória. Portanto, o tamanho da população será o resultado da multiplicação de muitas realizações de uma variável aleatória. Neste cenário, as abundâncias de um conjunto de espécies com dinâmicas populacionais independentes seguirão uma distribuição log-normal.

Para simular esta situação, criamos uma função que reitera a equação de crescimento de uma população time vezes, e repete este procedimento para N populações de tamanho inicial N0, sendo o valor de R uma variável uniforme, com valores entre Rmin e Rmax:

pops <- function(N,N0,Rmin,Rmax,time){
  results <- matrix(nrow=N,ncol=time)
  results[,1] <- N0
  for(i in 1:N){
    for(z in 2:time){
      results[i,z] <- results[i,z-1]*runif(1,Rmin,Rmax)
    }
  }
  results
}

O resultado desta função é uma matriz com N linhas (espécies) e time colunas (os intervalos de tempo). Com ela, simulamos a dinâmica populacional de 200 espécies por 20 intervalos de tempo, com valor médio do R=1, e atribuímos o resultado a um objeto:

abunds <- pops(200,500,Rmin=0.8,Rmax=1.2,time=20)

Com isto podemos fazer o histograma de densidade das abundâncias das 200 espécies no último intervalo de tempo:

h1<-(hist(abunds[,ncol(abunds)], prob=T))

E sobrepor a curva da função de densidade da log-normal, com os parâmetros estimados da amostra de 200 populações:

curve(dlnorm(x,meanlog=mean(log(abunds[,ncol(abunds)])),
                 sdlog=sd(log(abunds[,ncol(abunds)]))),
                        add=T, col="blue")

Se a log-normal é um bom modelo para descrever as abundâncias, seus logaritmos devem seguir uma distribuição normal. Podemos verificar isto visualmente com um gráfico de quantil:

qqnorm(log(abunds[,ncol(abunds)]))
qqline(log(abunds[,ncol(abunds)]))

Um gráfico que pode ajudar a entender a simulação:

nf <- layout(matrix(c(1,2),1,2,byrow=TRUE), c(3,1), TRUE)
par(mar=c(5,4,4,0))
matplot(t(abunds), col=1, lwd=0.5, lty=1, type="l",
        xlab= "Tempo",ylab="N de indivíduos", bty="l", ylim=range(h1$breaks))
par(mar=c(5,0,4,2))
barplot(h1$counts, horiz=T, axes=F, space=0)

Faça os exercícios 202.1 a 202.5 no sistema notaR.

1. A inferência estatística clássica baseia-se fortemente na distribuição Gaussiana. Quais seriam as razões? Quais vantagens e desvantagens você vê nessa estratégia?
2. A Gaussiana é a única distribuição que estudamos nesta unidade que tem a esperança indendente da variância. Você esperaria esta propriedade em medidas tomados de sistemas biológicos? Por que?
3. Um pesquisador descobriu que a distribuição de biomassas das populações das espécies de uma comunidade é muito diferente da família Gaussiana. Ele tentou algumas transformações matemáticas até descobrir que os logarítmos das biomassas podem ser descritas pela família Gaussiana. Qual deve ser a relação entre a esperança dos logarítmos das biomassas e da esperança das biomassas?

• Probability and stochastic distributions for ecological modeling. Capítulo 4, item 4.5.2 de: Bolker, B.M. 2008 Ecological Models and Data in R Princeton: Princeton University Press.

• Randon variables and probability distributions. Cap.2 de: Gotelli, N.J. & Ellison, A.M. 2004. A primer of ecological statistics. Sinauer.

• Probability theory. Cap.3 de: Otto, S. P. & Day, T. (2007). A biologist's guide to mathematical modelling in ecology and evolution. Princenton, Princenton University Press.

• Houle, D., Pélabon, C., Wagner, G. P., & Hansen, T. F. (2011). Measurement and meaning in biology. The Quarterly Review of Biology, 86(1), 3-34. Excelente revisão sobre consequências das transformações matemáticas de dados em biologia. pdf na página do autor.

• Ruel, J. J., & Ayres, M. P. (1999). Jensen’s inequality predicts effects of environmental variation. Trends in Ecology & Evolution, 14(9), 361-366. Introdução à desigualdade de Jensen e suas consequências biológicas. pdf na página do autor.

• Verbete da wikipedia sobre a Lei de Taylor, veja também as referências.

• Newman, M. E. (2005). Power laws, Pareto distributions and Zipf's law. Contemporary physics, 46(5), 323-351. Ótima introdução a uma família de distribuições sem esperança. preprint no repositório arxiv.

• Portal sobre distribuições de probabilidades na Wikipedia: http://en.wikipedia.org/wiki/Probability_distribution

• O Probability Distributions Viewer é um aplicativo em shiny que cria uma página html com gráficos interativos de distribuições de probabilidade. Você pode executá-lo de seu site interativo, ou localmente em seu computador. Instruções para executar localmente aqui

• Virtual Laboratories in Probability and Statistics: portal completo e aprofundado sobre probabilidade, teoria estatística e processos estocásticos. Excelentes textos de estudo, aplicativos, dados, bibliografia e recursos para programação.

• Distribuições interativas on-line do Statistics Online Computational Resource da UCLA: http://www.socr.ucla.edu/htmls/SOCR_Distributions.html

• Capítulo sobre variáveis aleatórias do e-book de Probabilidade e Estatística da UCLA: http://wiki.stat.ucla.edu/socr/index.php/EBook#Chapter_IV:_Probability_Distributions

• Excelente vídeo-aula sobre o conceito de momentos de distribuições probabilísticas, do projeto STEM do MIT: https://ocw.mit.edu/resources/res-tll-004-stem-concept-videos-fall-2013/videos/probability-and-statistics/moments-of-distributions/.

• Um excelente resumo das propriedades e relações entre distribuições de probabilidade:
  • Leemis, L. M., and J. T. McQueston. 2008. Univariate Distribution Relationships. The American Statistician 62:45-53.pdf na página do autor
  • Mapa interativo de relações entre as distribuições de probabilidadesfeito a partir do artigo de Leemis & McQueston (2008).

• Probability cheatsheet: teoria básica da probabilidade em 10 páginas.

• Várias atividades em JAVA com a normal no livro on-line Seeing Statistics : http://www.seeingstatistics.com/seeing1999/normal/normal.html

• Life is log-normal : argumentação provocadora a favor do uso da log-normal ao invés da normal, com link para o artigo e uma demonstração animada.

• ProbFact, um canal no twiter que posta fatos sobre distribuições de probabilidades.

• Sobre o conceito matemático de infinito:
  • dois ótimos vídeos do fantástico canal de popularização da matemática Numberphille:
    • O infinito é maior do que você acredita
    • Paradoxos da infinitude
  • Aline e o Infinito, do projeto Ela está em tudo (Sim, Matemática é coisa de menina) , do IME -USP e UFABC.