===================================================================
RCS file: /cvs/ext/extinction_risk_1.R,v
retrieving revision 1.15
retrieving revision 1.16
diff -u -p -r1.15 -r1.16
--- ext/extinction_risk_1.R	2015/08/09 01:28:14	1.15
+++ ext/extinction_risk_1.R	2015/08/12 11:18:58	1.16
@@ -1,5 +1,5 @@
 # extinction_risk_1.R
-# $Id: extinction_risk_1.R,v 1.15 2015/08/09 01:28:14 hako Exp $
+# $Id: extinction_risk_1.R,v 1.16 2015/08/12 11:18:58 hako Exp $
 #
 # Author: Hiroshi Hakoyama <hako@affrc.go.jp>
 # Copyright (c) 2013-2015 Hiroshi Hakoyama <hako@affrc.go.jp>, All rights reserved.
@@ -87,7 +87,12 @@ 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
   us <- q * s / (q - 1) # an unbiased estimate of variance
   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) {
     qqnorm(delta.log.n)
     qqline(delta.log.n)
@@ -191,6 +196,7 @@ ext1 <- function(dat, t = 100, ne = 1, alpha = 0.05, v
 
   if (verbose == TRUE) {
     results <- list(ne = ne, t = t, verbose = verbose,
+      AIC = AIC,
       n = q + 1,
       xd = xd,
       Growth.rate = mu,
@@ -224,6 +230,7 @@ print.ext1 <- function(obj, digits = 5) {
       formatC(obj$Variance, digits = digits), 
       formatC(obj$xd, digits = digits),
       formatC(obj$n, digits = digits),
+      formatC(obj$AIC, digits = digits),
       formatC(obj$Extinction.probability, digits = digits)), 
     c(paste("(",formatC(obj$lower.CL.mu, digits = digits),", ", 
       formatC(obj$upper.CL.mu, digits = digits),")", sep = ""), 
@@ -231,6 +238,7 @@ print.ext1 <- function(obj, digits = 5) {
       formatC(obj$upper.CL.s, digits = digits),")", sep = ""), 
       "-",
       "-",
+      "-",
       paste("(",formatC(obj$lower.CL.P, digits = digits),", ", 
       formatC(obj$upper.CL.P, digits = digits),")", sep = "")))
     dimnames(output.est) <- list(
@@ -238,6 +246,7 @@ print.ext1 <- function(obj, digits = 5) {
         "Variance:",
         "Current population size, xd = log(N(q) / ne):",
         "Sample size, n = q + 1:",
+        "AIC:",
         paste("Probability of decline to", obj$ne, "within", obj$t, "years:")),
       c("Estimate","95% confidence interval"))
     print(output.est)