===================================================================
RCS file: /cvs/ext/extinction_risk_1.R,v
retrieving revision 1.16
retrieving revision 1.17
diff -u -p -r1.16 -r1.17
--- ext/extinction_risk_1.R	2015/08/12 11:18:58	1.16
+++ ext/extinction_risk_1.R	2015/08/12 15:29:53	1.17
@@ -1,5 +1,5 @@
 # extinction_risk_1.R
-# $Id: extinction_risk_1.R,v 1.16 2015/08/12 11:18:58 hako Exp $
+# $Id: extinction_risk_1.R,v 1.17 2015/08/12 15:29:53 hako Exp $
 #
 # Author: Hiroshi Hakoyama <hako@affrc.go.jp>
 # Copyright (c) 2013-2015 Hiroshi Hakoyama <hako@affrc.go.jp>, All rights reserved.
@@ -87,12 +87,9 @@ 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
-  
+  l <- - (q / 2) * log(2 * pi * s) - sum((delta.log.n - mu)^2) / (2 * s) # ML for the distribution of x = log.n (not for n).
+  AIC <- - 2 * l + 4 # AIC for x
+
   if (qq.plot == TRUE) {
     qqnorm(delta.log.n)
     qqline(delta.log.n)
@@ -262,8 +259,8 @@ print.ext1 <- function(obj, digits = 5) {
 #
 # Examples
 # Yellowstone grizzly bears (from Dennis et al., 1991)
-  # dat <- data.frame(Year = c(1959, 1960, 1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968, 1969, 1970, 1971, 1972, 1973, 1974, 1975, 1976, 1977, 1978, 1979, 1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987), 
- # Population = c(44, 47, 46, 44, 46, 45, 46, 40, 39, 39, 42, 44, 41, 40, 33, 36, 34, 39, 35, 34, 38, 36, 37, 41, 39, 51, 47, 57, 47))
+  dat <- data.frame(Year = c(1959, 1960, 1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968, 1969, 1970, 1971, 1972, 1973, 1974, 1975, 1976, 1977, 1978, 1979, 1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987), 
+ Population = c(44, 47, 46, 44, 46, 45, 46, 40, 39, 39, 42, 44, 41, 40, 33, 36, 34, 39, 35, 34, 38, 36, 37, 41, 39, 51, 47, 57, 47))
 # The probability of extinction (of decline to population size 1) within 100 years
 # ext1(dat, t = 100)
 # The probability of decline to 10 individuals within 100 years