From 477a8f1b7e5841202ef29d3d8af3c93acd35c043 Mon Sep 17 00:00:00 2001
From: paradis <paradis@6e262413-ae40-0410-9e79-b911bd7a66b7>
Date: Fri, 20 Apr 2012 04:11:23 +0000
Subject: [PATCH] fix in birthdeath()

git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@186 6e262413-ae40-0410-9e79-b911bd7a66b7
---
 DESCRIPTION        |  4 +--
 NEWS               | 10 ++++++++
 R/birthdeath.R     | 63 ++++++++++++++++++++++++++++------------------
 man/subtreeplot.Rd |  8 +++---
 4 files changed, 55 insertions(+), 30 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index a9eed55..2877fec 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: ape
-Version: 3.0-2
-Date: 2012-04-05
+Version: 3.0-3
+Date: 2012-04-20
 Title: Analyses of Phylogenetics and Evolution
 Author: Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong, Richard Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel, Christoph Heibl, Daniel Lawson, Vincent Lefort, Pierre Legendre, Jim Lemon, Yvonnick Noel, Johan Nylander, Rainer Opgen-Rhein, Andrei-Alin Popescu, Klaus Schliep, Korbinian Strimmer, Damien de Vienne
 Maintainer: Emmanuel Paradis <Emmanuel.Paradis@ird.fr>
diff --git a/NEWS b/NEWS
index 1848597..3660ac4 100644
--- a/NEWS
+++ b/NEWS
@@ -1,6 +1,16 @@
 		CHANGES IN APE VERSION 3.0-2
 
 
+BUG FIXES
+
+    o birthdeath() now catches errors and warnings much better so that
+      a result is returned in most cases.
+
+
+
+		CHANGES IN APE VERSION 3.0-2
+
+
 NEW FEATURES
 
     o The new function alex (alignment explorator) zooms in a DNA
diff --git a/R/birthdeath.R b/R/birthdeath.R
index 3aee631..7d53fd2 100644
--- a/R/birthdeath.R
+++ b/R/birthdeath.R
@@ -1,4 +1,4 @@
-## birthdeath.R (2010-11-04)
+## birthdeath.R (2012-04-20)
 
 ##   Estimation of Speciation and Extinction Rates
 ##             with Birth-Death Models
@@ -6,7 +6,7 @@
 ## birthdeath: standard model
 ## bd.ext: extended version
 
-## Copyright 2002-2010 Emmanuel Paradis
+## Copyright 2002-2012 Emmanuel Paradis
 
 ## This file is part of the R-package `ape'.
 ## See the file ../COPYING for licensing issues.
@@ -17,7 +17,7 @@ birthdeath <- function(phy)
     N <- length(phy$tip.label)
     x <- c(NA, branching.times(phy))
     dev <- function(a, r) {
-        if (r <= 0) return(1e50)
+        if (r < 0 || a > 1) return(1e100)
         -2 * (lfactorial(N - 1)
               + (N - 2) * log(r)
               + r * sum(x[3:N])
@@ -28,32 +28,47 @@ birthdeath <- function(phy)
     if (out$estimate[1] < 0) {
         out <- nlm(function(p) dev(0, p), 0.2, hessian = TRUE)
         para <- c(0, out$estimate)
-        se <- c(0, sqrt(diag(solve(out$hessian))))
+        inv.hessian <- try(solve(out$hessian))
+        se <-
+            if (class(inv.hessian) == "try-error") NA
+            else sqrt(diag(inv.hessian))
+        se <- c(0, se)
     }
     else {
         para <- out$estimate
-        se <- sqrt(diag(solve(out$hessian)))
+        inv.hessian <- try(solve(out$hessian))
+        se <-
+            if (class(inv.hessian) == "try-error") c(NA, NA)
+            else sqrt(diag(inv.hessian))
     }
     Dev <- out$minimum
-    ## compute the 95 % profile likelihood CIs
-    ## (not very clean... but seems to work -- EP 2003-03-29)
-    CI <- matrix(NA, 2, 2)
-    foo <- function(p) dev(p, para[2]) - 3.84 - Dev
-    inc <- 1e-2
-    lo <- para[1] - inc
-    up <- para[1] + inc
-    while (foo(lo) < 0) lo <- lo - inc
-    while (foo(up) < 0) up <- up + inc
-    CI[1, 1] <- uniroot(foo, lower = lo, upper = para[1])$root
-    if (CI[1, 1] < 0) CI[1, 1] <- 0
-    CI[1, 2] <- uniroot(foo, lower = para[1], upper = up)$root
-    foo <- function(p) dev(para[1], p) - 3.84 - Dev
-    lo <- para[2] - inc
-    up <- para[2] + inc
-    while (foo(lo) < 0) lo <- lo - inc
-    while (foo(up) < 0) up <- up + inc
-    CI[2, 1] <- uniroot(foo, lower = lo, upper = para[2])$root
-    CI[2, 2] <- uniroot(foo, lower = para[2], upper = up)$root
+
+    ## 95% profile likelihood CIs
+
+    ## which: index of the parameter (1 or 2)
+    ## s: sign of the increment (-1 or +1)
+    foo <- function(which, s) {
+        i <- 0.1
+
+        if (which == 1) {
+            p <- para[1] + s * i
+            bar <- function() dev(p, para[2])
+        } else { # which == 2
+            p <- para[2] + s * i
+            bar <- function() dev(para[1], p)
+        }
+
+        while (i > 1e-9) {
+            while (bar() < Dev + 3.84) p <- p + s * i
+            p <- p - s * i
+            i <- i / 10
+        }
+        p
+    }
+
+    CI <- mapply(foo, c(1, 2, 1, 2), c(-1, -1, 1, 1))
+    dim(CI) <- c(2, 2)
+
     names(para) <- names(se) <- rownames(CI) <- c("d/b", "b-d")
     colnames(CI) <- c("lo", "up")
     obj <- list(tree = deparse(substitute(phy)), N = N,
diff --git a/man/subtreeplot.Rd b/man/subtreeplot.Rd
index ce6863e..33f8469 100644
--- a/man/subtreeplot.Rd
+++ b/man/subtreeplot.Rd
@@ -30,12 +30,12 @@ subtreeplot(x, wait=FALSE, ...)
 \examples{
 \dontrun{
 #example 1: simple
-tree1<-rtree(50) #random tree with 50 leaves
-tree2<-subtreeplot(tree1, wait=TRUE) # on exit, tree2 will be a subtree of tree1.
+tree1 <- rtree(50)
+tree2 <- subtreeplot(tree1, wait = TRUE) # on exit, tree2 will be a subtree of tree1
 
 #example 2: more than one zoom
-tree1<-rtree(60)
-tree2<-subtreeplot(subtreeplot(subtreeplot(tree1))) #allows three succssive zooms.
+tree1 <- rtree(60)
+tree2 <- subtreeplot(subtreeplot(subtreeplot(tree1))) # allow three succssive zooms
 }
 }
 \keyword{hplot}
-- 
2.39.2