===================================================================
RCS file: /cvs/ext/extinction_risk_1.R,v
retrieving revision 1.20
retrieving revision 1.21
diff -u -p -r1.20 -r1.21
--- ext/extinction_risk_1.R	2015/08/15 06:50:20	1.20
+++ ext/extinction_risk_1.R	2015/08/16 13:59:30	1.21
@@ -1,5 +1,5 @@
 # extinction_risk_1.R
-# $Id: extinction_risk_1.R,v 1.20 2015/08/15 06:50:20 hako Exp $
+# $Id: extinction_risk_1.R,v 1.21 2015/08/16 13:59:30 hako Exp $
 #
 # Author: Hiroshi Hakoyama <hako@affrc.go.jp>
 # Copyright (c) 2013-2015 Hiroshi Hakoyama <hako@affrc.go.jp>, All rights reserved.
@@ -206,25 +206,38 @@ RejectionRate <- function(mu.obs, xd, s.obs, tt, tq, q
 AsymptoticallyExactConfidenceInterval <- function(mu.obs, xd, s.obs, tt, tq, qq, alpha, n.sim, alpha.sim, ci.function = ConfidenceInterval) {
   P.obs <- Pr(W(mu.obs, xd, s.obs, tt), Z(mu.obs, xd, s.obs, tt))
   if (P.obs == 0 || P.obs == 1) {
-    return(c(ci.function(mu.obs, xd, s.obs, tt, tq, qq, alpha), NaN, NaN))
+    return(ci.function(mu.obs, xd, s.obs, tt, tq, qq, alpha))
   }
   res <- RejectionRate(mu.obs, xd, s.obs, tt, tq, qq, alpha, alpha, n.sim, ci.function)
-  if (is.nan(res[[1]])) return(c(NaN, NaN, alpha, 0))
   gg <- res[[1]]
   pp <- res[[2]]
-  ai <- alpha
-  i <- 0
+#    cat(alpha, gg, "\n")
+  if (gg > alpha) {
+    ai.lower <- 0
+    ai.upper <- alpha
+    ai <- alpha / 2
+  } else {
+    ai.lower <- alpha
+    ai.upper <- 1
+    ai <- (1 + alpha) / 2
+  }
   while (pp < alpha.sim) {
-    i <- i + 1
     res <- RejectionRate(mu.obs, xd, s.obs, tt, tq, qq, ai, alpha, n.sim, ci.function)
-    if (is.nan(res[[1]])) return(c(NaN, NaN, ai, i))
     gg <- res[[1]]
     pp <- res[[2]]
-    ai <- ai * alpha / gg
-    if (ai < 10^(-10) || ai > 1) return(c(NaN, NaN, ai, i))
-    if (i > 20) return(c(NaN, NaN, ai, i))
+
+    if (gg > alpha) {
+      ai.lower <- ai.lower
+      ai.upper <- ai
+      ai <- (ai.lower + ai.upper) / 2
+    } else {
+      ai.lower <- ai
+      ai.upper <- ai.upper
+      ai <- (ai.lower + ai.upper) / 2
+    }
+#  cat(ai, gg, "\n")
   }
-  c(res[[3]], res[[4]], ai, i)
+  c(res[[3]], res[[4]], ai, gg, pp)
 }
 #-------------------------------------------------------------------
 print.ext1 <- function(obj, digits = 5) {
@@ -266,8 +279,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