[BACK]Return to extinction_risk_1.R CVS log [TXT][DIR] Up to [local] / ext

Diff for /ext/extinction_risk_1.R between version 1.15 and 1.16

version 1.15, 2015/08/09 01:28:14 version 1.16, 2015/08/12 11:18:58
Line 87  ext1 <- function(dat, t = 100, ne = 1, alpha = 0.05, v
Line 87  ext1 <- function(dat, t = 100, ne = 1, alpha = 0.05, v
   s <- (1 / q) * sum((delta.log.n - mu * tau)^2 / tau) # ML estimate of variance    s <- (1 / q) * sum((delta.log.n - mu * tau)^2 / tau) # ML estimate of variance
   us <- q * s / (q - 1) # an unbiased estimate of variance    us <- q * s / (q - 1) # an unbiased estimate of variance
   xd <- log(ps[length(ps)] / ne)    xd <- log(ps[length(ps)] / ne)
     nn <- ps[-1]
     nn2 <- ps[-length(ps)]
     lnl <- - sum(log(nn * sqrt(2 * tau * pi))) - (q / 2) * log(s)
            - 1 / (2 * s) * sum((1 / tau) * (log(nn / nn2) - mu * tau)^2)
     AIC <- - 2 * lnl + 4 # 2 parameters, mu and s
   
   if (qq.plot == TRUE) {    if (qq.plot == TRUE) {
     qqnorm(delta.log.n)      qqnorm(delta.log.n)
     qqline(delta.log.n)      qqline(delta.log.n)
Line 191  ext1 <- function(dat, t = 100, ne = 1, alpha = 0.05, v
Line 196  ext1 <- function(dat, t = 100, ne = 1, alpha = 0.05, v
   
   if (verbose == TRUE) {    if (verbose == TRUE) {
     results <- list(ne = ne, t = t, verbose = verbose,      results <- list(ne = ne, t = t, verbose = verbose,
         AIC = AIC,
       n = q + 1,        n = q + 1,
       xd = xd,        xd = xd,
       Growth.rate = mu,        Growth.rate = mu,
Line 224  print.ext1 <- function(obj, digits = 5) {
Line 230  print.ext1 <- function(obj, digits = 5) {
       formatC(obj$Variance, digits = digits),        formatC(obj$Variance, digits = digits),
       formatC(obj$xd, digits = digits),        formatC(obj$xd, digits = digits),
       formatC(obj$n, digits = digits),        formatC(obj$n, digits = digits),
         formatC(obj$AIC, digits = digits),
       formatC(obj$Extinction.probability, digits = digits)),        formatC(obj$Extinction.probability, digits = digits)),
     c(paste("(",formatC(obj$lower.CL.mu, digits = digits),", ",      c(paste("(",formatC(obj$lower.CL.mu, digits = digits),", ",
       formatC(obj$upper.CL.mu, digits = digits),")", sep = ""),        formatC(obj$upper.CL.mu, digits = digits),")", sep = ""),
Line 231  print.ext1 <- function(obj, digits = 5) {
Line 238  print.ext1 <- function(obj, digits = 5) {
       formatC(obj$upper.CL.s, digits = digits),")", sep = ""),        formatC(obj$upper.CL.s, digits = digits),")", sep = ""),
       "-",        "-",
       "-",        "-",
         "-",
       paste("(",formatC(obj$lower.CL.P, digits = digits),", ",        paste("(",formatC(obj$lower.CL.P, digits = digits),", ",
       formatC(obj$upper.CL.P, digits = digits),")", sep = "")))        formatC(obj$upper.CL.P, digits = digits),")", sep = "")))
     dimnames(output.est) <- list(      dimnames(output.est) <- list(
Line 238  print.ext1 <- function(obj, digits = 5) {
Line 246  print.ext1 <- function(obj, digits = 5) {
         "Variance:",          "Variance:",
         "Current population size, xd = log(N(q) / ne):",          "Current population size, xd = log(N(q) / ne):",
         "Sample size, n = q + 1:",          "Sample size, n = q + 1:",
           "AIC:",
         paste("Probability of decline to", obj$ne, "within", obj$t, "years:")),          paste("Probability of decline to", obj$ne, "within", obj$t, "years:")),
       c("Estimate","95% confidence interval"))        c("Estimate","95% confidence interval"))
     print(output.est)      print(output.est)
Legend:
Removed from v.1.15  
changed lines
  Added in v.1.16