====== Simulações ======


O codigo do PI parece estar funcionando bem, inclui apenas como atributo do objeto de saída os valores de argumentos de entrada na simulação. Isso facilita a vida quando estamos apresentando os resultados e fazendo gráficos. Estou pensando também em criar uma versão onde o objeto de saída trás apenas as informações que usamos. Rodei quatro simulações e o workspace está demorando para subir. Isso pode complicar nossa vida depois. Uma alternativa viável é usar os dados intermediários dos ciclos ao invés de colocar tudo no arquivo de saída.
  - arquivo com a função {{simula.neutra_19jul2010.r|}}
  - entre na página para ver os resultados da simulacões: [[.:result]]

===== Parâmetros da Simulação =====

==== Definidos ====
  * $S\ = \ $ número de espécies da comunidade  : //inteiro positivo//
  * $J\ = \ $ número de indivíduos da comunidade (tamanho da comunidade): //inteiro positivo//
  * $x_i\ = \ $ número de propágulos que o indivíduo $i$ produz por intervalo
  * $c_v\ = \ $ proporção de variação na produção de propágulos em relação a sua esperança: // real positivo //
  * $X\ = \ $ total de propágulos que um indivíduo produz (esforço reprodutivo total): //inteiro positivo//
  

==== Deduzidos ====
   * Dado que o esforço reprodutivo total é o mesmo para todos os indivíduos ($X$), o tempo médio de vida de cada indivíduo é $$E[t_i]\ = \ X/x_i $$ 
   * A cada intervalo cada indivíduo tem uma probabilidade fixa $p_i$ de morrer. Portanto, a probabilidade de sobreviver $t$ intervalos é dada por uma [[http://en.wikipedia.org/wiki/Geometric_distribution|distribuição geométrica]]((Para tutorial em R sobre esta distribuição veja [[http://cmq.esalq.usp.br/BIE5781|Wiki da BIE5781]])) com suporte $0,\infty$. A esperança desta distribuição é $E[t_i]=p_i^{-1}$, portanto:
                  * $$p_i \ = \ x_i/X$$
  * Dado $s$ igual para todas as espécies, e uma distribuição gaussiana de caracteres contínuos, o número de sementes por intervalo $x_i$ da prole de um indivíduo pode ser definido como uma variável normal, com parâmetros 
                  * $$\mu \ = \ x_i$$
                  * $$\sigma \ = \ c_v * mu\ $$ 

===== Pseudocódigo =====

  - Defina $J$, $S$, $X$, $c_v$ 
  - Estabeleça as condições iniciais de acordo com o modelo de Hubbell:
          * Comece com o mesmo número de indivíduos por espécie, que será $J/S$.  
((vamos usar número de espécies e indivíduos da restinga, por hectare: ver dados da parcela PEIC))
          * As probabilidades de morrer são as mesmas para todos os indivíduos, e o número esperado de mortes por intervalo é um. Portanto temos para todos os indivíduos: 
                  * $$p_i \ = \ 1/J $$
                  * $$x_i0 \ = \ X/J $$, pois $$E[t_i] \ = \ J$$ : //inteiro//
  - Atribua a cada indivíduo sua espécie, e os valores de $p_i$ e $x_i$
  - Gere $x_i$ propágulos para cada indivíduo, mantendo a identidade da mãe de cada um.  
  - Sorteie os indivíduos que morrem: 
          * Gere $u_i$, números sorteados de um distribuição uniforme entre zero e um, sendo que $i$ vai de $1$ a  $J$ ;
          * Morrem todos os indivíduos em que $p_i>u_i$
  - Conte o número de mortos, $D$
  - Sorteie sem reposição $D$ propágulos do total de propágulos disponíveis, para ocupar o lugar dos mortos, .
  - Para a mãe de cada propágulo sorteado, sorteie com reposição um pai, entre os indivíduos da mesma espécie.
  - Calcule para cada combinação de mãe e pai a média do número de propágulos que eles produzem por intervalo:
          * $$\bar x_{ij} = \frac{x_i + x_j}{2}$$
  - Sorteie o valor de $x_i$ de cada filhote de uma normal com parâmetros:
          * $$\mu = \bar x_{ij}$$
          * $$\sigma = \bar x_{ij} * c_v$$  
  - Armazene para cada indivíduo: 
          * Identidade da espécie 
          * $$x_i$$
  - Armazene para cada ciclo: 
          * $$D$$        
  - Calcule para cada novo indivíduo:
          * $$p_i \ = \ x_i/X$$
  - Retorne ao passo 4

===== Códigos =====

==== Em R ====
=== Função principal ===
  * Código inicial do Alexandre, com as seguintes modificações feitas pelo Paulo Inácio:
      * A fertilidade dos filhos é tomada de uma normal discretizada truncada, com parâmetro $$mu$$ = média das fertilidades dos pais, e coeficiente de variação definido pelo usuário (ver discussão). Os limites do truncamento são 1 e X.
      * Pai sorteado entre os co-específicos (código original soretava entre todos os indivíduos).
      * Ao invés do usuário definir o total de indivíduos (J), agora há um argumento para definir o número de indivíduos por espécie (j). O total de indivíduos é calculado pela função como o produto $$J = jS$$, garantindo que J seja um múltiplo de S.
      * Opção de gravar os resultados a intervalos regulares de tempo, que é útil para não estourar memória nas simulações longas. O argumento ''ciclo'' controla o tempo total da simulação, e o argumento ''step'' o intervalo para gravação.
  * **Argumentos da função**:
      * ''S'' : Número de espécies inicial
      * ''j'' : número de individuos por espécie
      * ''X'' : esforço reprodutivo total((como a simulação é probabilística, este valor é na verdade uma esperança.))
      * ''cv'' : coeficiente de variação fenotípica da fertilidade individual 
      * ''ciclo'' : número total de ciclos que devem ocorrer. Em cada ciclo os adultos estão sujeitos a morrer, de acordo com a probabilidade que é fertilidades/esforço reprodutivo.
      * ''step'' : intervalos para gravação dos resulatdos. O dados da comunidade são gravados a cada ''step'' ciclos. Para gravar em cada ciclo, faça ''step=1''

  * **Saída da Função**: uma lista com so seguintes objetos:
    * ''tempo'': vetor com rótulos do tempo passado em ciclos a cada instante gravado 
    * ''sp.list'': matriz de J linhas por ciclos/step colunas, com a identidade de cada indivíduso a cada tempo gravado.
    * ''sementes'' : matriz de J linhas por ciclos/step colunas, com a fertilidade (n de propágulos produzidos/ciclo) de cada indivíduso a cada tempo gravado.
    * prob.morte: matriz de J linhas por ciclos/step colunas, com a probabilidade de morte de cada indivíduso a cada tempo gravado.
    * n.mortes: vetor com o número total de mortes ao longo de cada intervalo de ''step'' ciclos que é gravado.


<code>
simula.neutra.step=function(S= 100, j=10, X=1000, cv=0.1, ciclo=1e6, step=100){ 
  ## Tamanho da comunidade
  J <- S*j
  ##Verifica se X e multiplo de J
  if(abs(X/J - round(X/J)) > .Machine$double.eps^0.5)
    {stop("\n\tO potencial reprodutivo (X) precisa ser multiplo do tamanho da comunidade (J). Tente novamente!\n\n")} 
  ##Matrizes para guardar os resultados
  ## matriz da especie de cada individuo por ciclo
  ind.mat=matrix(nrow=J,ncol=1+ciclo/step) 
  ##Matriz de propagulos produzidos por individuo em cada ciclo
  prop.mat=matrix(nrow=J,ncol=1+ciclo/step)
  ## Matriz de probabilidade de morte de cada individuo, por ciclo
  dead.mat=matrix(nrow=J,ncol=1+ciclo/step)
  ##Vetor com numero de mortos por ciclo
  n.dead <- c()
  n.dead[1] <- 0
  ##CONDICOES INICIAIS##
  ##Deduzidas de acordo com o modelo de Hubbell:
  ## Todas as especies comecam com o mesmo numero de individuos (j=J/S)
  ind.mat[,1] <- rep(1:S,each=j)
  cod.sp <- ind.mat[,1]
  ##A probabilidade inicial de morte, tomada de uma geométrica,
  ##dado que o numero de mortes esperado por ciclo eh a mesma para todos (p=1/J)
  dead.mat[,1] <- 1/J
  p.death <- dead.mat[,1]
  ##O esforco reprodutivo inicial e deduzido da esperanca de tempo de vida da geometrica E[t]=J
  ## o que portanto resulta em X/J propagulos produzidos a cada ciclo, por todos
  prop.mat[,1] <- X/J
  n.propag <- prop.mat[,1]
  ##Aqui comecam as simulacoes
  for(i in 2:(1+ciclo/step)){
    n.mortes <- 0
    for(j in 1:step){
      ## Sorteio dos que morrerao
      morte=rbinom(J, 1, prob=p.death)
      ##Total de mortos, que e armazenado em n.dead
      D=sum(morte)
      n.mortes <- n.mortes+D
      ## Se ha mortos comeca aqui a rotina de substituicao dos valores
      if(D>0) 
        {
           ##vetor de propagulos: cada propagulo tem o codigo numerico do individuo
          seed.bank <- rep(1:J,n.propag)
          ##Indices dos individuos que morreram
          nascer= which(morte==1)
          ##Sorteio dos propagulos que irao repor os mortos
          mami=sample(seed.bank,D)
          ##Um vetor para armazenar o fenotipo do pai
          papi <- c()
          ##Um loop para sortear o pai entre os individuos da especie de cada
          ## propagulo-mae sorteado
          for(w in 1:D){
            papi[w] <- sample(n.propag[ cod.sp==cod.sp[mami[w]] ],1)
          }
          ##O valor esperado de propagulos dos filhotes eh
          ## a media do numero medio de propagulos produzidos pelo pai e pela mae
          medias.prop=(n.propag[mami]+papi)/2
          ##Codigos das especies dos mortos sao substituidos pelos nascidos
          ##dos propagulos sorteados
          cod.sp[nascer]<-cod.sp[mami]
          ## Numero medio de propagulos dos novos individuos nascidos e sorteada
          ## de uma normal discretizada e truncada entre 1 e J,
          ## com o coeficiente de variacao estabelecido
          n.propag[nascer] <- sapply(medias.prop,rnormt,cv=cv,min=1,max=X)
          ##A matriz de probabilidades de morrer eh atualizada para os novos individuos
          p.death[nascer] <- n.propag[nascer]/X
        }
    }
    ## A cada step ciclos os resultados sao gravados
    ind.mat[,i] <- cod.sp
    dead.mat[,i] <- p.death
    prop.mat[,i] <- n.propag
    n.dead[i] <- n.mortes
  }
  tempo <- seq(0,ciclo,by=step)
  colnames(ind.mat) <- tempo
  colnames(dead.mat) <- tempo
  colnames(prop.mat) <- tempo
  names(n.dead) <- tempo
  resulta=list(tempo=tempo,sp.list=ind.mat,sementes=prop.mat,prob.morte=dead.mat,n.mortes=n.dead)
  return(resulta)
}
</code>

=== Funções Acessórias ===

==Modelo Neutro de Hubbell==

Esta é a função ''simula.neutra.step'' editada para retirara o tradeoff reprodutivo e a variação neste caracter. O resulatdo é o modelo neutro do Hubbell, com a diferença que o número **esperado** de mortes a cada ciclo é um, mas há variação estocástica nisto. no modelo original, há sempre uma morte por ciclo.

<code>
simula.neutra.hub=function(S= 100, j=10, ciclo=1e4, step=100){ 
  ## Tamanho da comunidade
  J <- S*j
  ##Matrizes para guardar os resultados
  ## matriz da especie de cada individuo por ciclo
  ind.mat=matrix(nrow=J,ncol=1+ciclo/step) 
  ##Vetor com numero de mortos por ciclo
  n.dead <- c()
  n.dead[1] <- 0
  ##CONDICOES INICIAIS##
  ##Deduzidas de acordo com o modelo de Hubbell:
  ## Todas as especies comecam com o mesmo numero de individuos (j=J/S)
  ind.mat[,1] <- rep(1:S,each=j)
  cod.sp <- ind.mat[,1]
  ##A probabilidade inicial de morte, tomada de uma geométrica,
  ##dado que o numero de mortes esperado por ciclo eh a mesma para todos (p=1/J)
  p.death <- 1/J
  ##Aqui comecam as simulacoes
  for(i in 2:(1+ciclo/step)){
    n.mortes <- 0
    for(j in 1:step){
      ## Sorteio dos que morrerao
      morte=rbinom(J, 1, prob=p.death)
      ##Total de mortos, que e armazenado em n.dead
      D=sum(morte)
      n.mortes <- n.mortes+D
      ## Se ha mortos comeca aqui a rotina de substituicao dos valores
      if(D>0) 
        {
          ##Indices dos individuos que morreram
          nascer= which(morte==1)
          ##Quem substitui os mortos
          novos <- sample(1:J,D,replace=T)
          ##Substituindo
          cod.sp[nascer]<-cod.sp[novos]
        }
    }
    ## A cada step ciclos os resultados sao gravados
    ind.mat[,i] <- cod.sp
    n.dead[i] <- n.mortes
  }
  tempo <- seq(0,ciclo,by=step)
  colnames(ind.mat) <- tempo
  names(n.dead) <- tempo
  resulta=list(tempo=tempo,sp.list=ind.mat,n.mortes=n.dead)
  return(resulta)
}
</code>

==Normal Discretizada Truncada==
Sorteia n números inteiros de uma aproximação discreta da normal truncada em seus limites inferiores e superiores dados sua média e coeficiente de variação. 

  * **Pseudocódigo:**
     * Dados os limites de truncamento (no caso um mínimo de um filhote e máximo de X filhotes), divida este intervalo contínuo em trechos de tamanho um, cujo ponto médio são os valores inteiros (e.g. 0.5-1.5, 1.5-2.5 ...). 
     * A probabilidade associada a cada um destes intervalos é calculada para uma normal com os parâmetros especificados.
     * Use estas probabilidades para sortear com reposição os números inteiros que são os pontos médios de cada intervalo.

  * **Argumentos:**
    * ''mean'' : parâmetro $$mu$$ da normal, que corresponde à média
    * ''n'' : número de sorteios
    * ''cv'' : coeficiente de variação (desvio padrão / média)
    * ''min'' , ''max'' : limites mínimo e máximo de truncamneto

 
<code>
##Funcao de sorteio de uma normal discretizada e truncada
rnormt <- function(mean,n=1,cv,min,max){
  vals <- (min-0.5):(max+0.5)
  p <- pnorm(vals,mean=mean,sd=mean*cv)
  p2 <- diff(p)
  sample(min:max,n,prob=p2,replace=T)
}
</code>

==Riqueza de espécies==
Função simples que retorna o número de elementos únicos de um vetor alfanumérico. Aplicado ao vetor de coódigo de espécie de cada indivíduo retorna a riqueza.
<code>
rich <- function(x)length(unique(x))
</code>

==Rannk-abundance Plot==
Dado um vetor de abundãncia de espécies, faz um gráfico simples de Whitaker (log abundâncias x ranque).
<code>
w.plot <- function(x,...){
  plot(1:length(x),sort(x,decreasing=T),
       xlab="Rank order",ylab="Abundance",log="y",...)
}
</code>

==Box-plot e rank-abundance plot==
Função feita para explorar graficamente a distribuição da fertilidade. Plota box-plots da fertlidade dos indivíduso de cada espécie, sobre um gr[áfico de rank-abundância da fertilidade de cada espécie, num dado instante gravado da simulação (argumento ''tempo''). 
<code>
box.w <- function(tempo,cod.sp,n.propag){
  sp <- cod.sp[,tempo]
  prop <- n.propag[,tempo]
  ab <- tapply(sp,sp,length)
  par(mfrow=c(2,1))
  boxplot(prop~factor(sp,levels=names(ab)[order(ab,decreasing=T)]))
  w.plot(ab)
  par(mfrow=c(1,1))
  }
</code>

==Fertilidade Média por Espécie==
A simulação grava uma matriz com as fertilidades de cada indivíduo a cada ''step'' ciclos. Esta função resume estes dados, calculando uma estatística descritiva (o //default// é a média) da fertilidade de cada espécie, em cada instante de tempo.

  * **Argumentos:**
     * ''cod.sp'': matriz retronada pela simulação de identidade de cada indivíduo x tempo
     * ''n.propag'': matriz retronada pela simulação de fertilidade de cada indivíduo x tempo
     * ''fun'': função a aplicar a cada espécie por tempo.

  * **Saída**
     * Uma matriz das médias de fertilidade de cada espécie em cada tempo. **Atenção:** os  valores nesta matriz são não existentes (''NA'') a patir do tempo em que a espécie de sextingue.

<code>
fert.t <- function(cod.sp,n.propag,fun=mean){
  especie <- unique(cod.sp[,1])
  medias <- matrix(ncol=ncol(n.propag),nrow=length(especie))
  colnames(medias) <- colnames(n.propag)
  rownames(medias) <- especie
  for(i in 1:ncol(n.propag)){
    medias[,i] <- tapply(n.propag[,i],factor(cod.sp[,i],levels=especie),fun)
  }
  medias
}
</code>

===== Discussão =====
Novos tópicos com //"new blog entry"//. Para discutir um tópico, use o link //comments//.

{{blog>blog:neutral1?10}}