fix in birthdeath()

git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@186 6e262413-ae40-0410-9e79-b911bd7a66b7

diff --git a/DESCRIPTION b/DESCRIPTION
index a9eed554c5d4e05560fb6eb7b4dee7644ed14bce..2877fec05708a6d583fa72b7f0d7eeca190b268e 100644 (file)
--- 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 18485974832ad5157a9c413d79334e09afdce807..3660ac41656c766a522d731a10cfc489ae01c384 100644 (file)
--- 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 3aee6313863e48bf3a9b6d0bc69f2fd2fa092a05..7d53fd29218613b8a3047a8182ceee2c8556b46b 100644 (file)
--- 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 ce6863ecdde096270d3e6b73fdf88f58d8a7d76b..33f846993f5586383387a6a0cad685dad75d29b9 100644 (file)
--- 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}