simula.neutra.step: simula.neutra.step

In adalardo/neutraPack: Neutral model simulation

Usage

simula.neutra.step(S = 100, j = 10, X = 1000, dp = 0.1, ciclo = 1e+06, step = 100)

Arguments

S
j
X
dp
ciclo
step

Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (S = 100, j = 10, X = 1000, dp = 0.1, ciclo = 1e+06, 
    step = 100) 
{
    t0 = proc.time()[[3]]
    J <- S * 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")
    }
    ind.mat = matrix(nrow = J, ncol = 1 + ciclo/step)
    prop.mat = matrix(nrow = J, ncol = 1 + ciclo/step)
    dead.mat = matrix(nrow = J, ncol = 1 + ciclo/step)
    n.dead <- c()
    n.dead[1] <- 0
    ind.mat[, 1] <- rep(1:S, each = j)
    cod.sp <- ind.mat[, 1]
    dead.mat[, 1] <- 1/J
    p.death <- dead.mat[, 1]
    prop.mat[, 1] <- X/J
    n.propag <- prop.mat[, 1]
    for (i in 2:(1 + ciclo/step)) {
        n.mortes <- 0
        for (a in 1:step) {
            morte = rbinom(J, 1, prob = p.death)
            D = sum(morte)
            n.mortes <- n.mortes + D
            if (D > 0) {
                seed.bank <- rep(1:J, n.propag)
                nascer = which(morte == 1)
                mami = sample(seed.bank, D)
                papi <- c()
                for (w in 1:D) {
                  papi[w] <- sample(n.propag[cod.sp == cod.sp[mami[w]]], 
                    1)
                }
                medias.prop = (n.propag[mami] + papi)/2
                cod.sp[nascer] <- cod.sp[mami]
                n.propag[nascer] <- sapply(medias.prop, rnormt, 
                  dp = dp, min = 1, max = X)
                p.death[nascer] <- n.propag[nascer]/X
            }
        }
        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)
    t1 = proc.time()[[3]]
    cat("\n\t tempo de processamento: ", round((t1 - t0)/60, 
        2), "\n")
    attributes(resulta)$start = list(especies = S, individuos = j, 
        nprop = X, sd = dp, ciclos = ciclo, passos = step)
    return(resulta)
  }

adalardo/neutraPack documentation built on May 10, 2019, 5:13 a.m.