X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=R%2Fbirthdeath.R;h=3aee6313863e48bf3a9b6d0bc69f2fd2fa092a05;hb=f71a875dd49951d5c72f55b8578b52cdaf23f956;hp=82decda388a03454f9b988e1f5dd9f39de834846;hpb=4ceef408de61dc86f0a93b0396aecc6e30cc0d70;p=ape.git

diff --git a/R/birthdeath.R b/R/birthdeath.R
index 82decda..3aee631 100644
--- a/R/birthdeath.R
+++ b/R/birthdeath.R
@@ -1,4 +1,4 @@
-## birthdeath.R (2009-05-10)
+## birthdeath.R (2010-11-04)
 
 ##   Estimation of Speciation and Extinction Rates
 ##             with Birth-Death Models
@@ -6,7 +6,7 @@
 ## birthdeath: standard model
 ## bd.ext: extended version
 
-## Copyright 2002-2009 Emmanuel Paradis
+## Copyright 2002-2010 Emmanuel Paradis
 
 ## This file is part of the R-package `ape'.
 ## See the file ../COPYING for licensing issues.
@@ -17,6 +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)
         -2 * (lfactorial(N - 1)
               + (N - 2) * log(r)
               + r * sum(x[3:N])
@@ -78,39 +79,47 @@ print.birthdeath <- function(x, ...)
     cat("      b - d: [", x$CI[2, 1], ", ", x$CI[2, 2], "]", "\n\n", sep = "")
 }
 
-bd.ext <- function(phy, S)
+bd.ext <- function(phy, S, conditional = TRUE)
 {
-    if (!inherits(phy, "phylo")) stop('object "phy" is not of class "phylo"')
+    if (!inherits(phy, "phylo"))
+        stop('object "phy" is not of class "phylo"')
     if (!is.null(names(S))) {
         if (all(names(S) %in% phy$tip.label)) S <- S[phy$tip.label]
-        else warning('the names of argument "S" and the names of the tip labels
-did not match: the former were ignored in the analysis.')
+        else warning('the names of argument "S" and the tip labels
+did not match: the former were ignored.')
     }
     N <- length(S)
     x <- branching.times(phy)
     x <- c(x[1], x)
     trm.br <- phy$edge.length[phy$edge[, 2] <= N]
-    dev <- function(a, r)
-    {
-        -2 * (lfactorial(N - 1)
-              + (N - 2) * log(r)
-              + (3 * N) * log(1 - a)
-              + 2 * r * sum(x[2:N])
-              - 2 * sum(log(exp(r * x[2:N]) - a))
-              + r * sum(trm.br)
-              + sum((S - 1) * log(exp(r * trm.br) - 1))
-              - sum((S + 1) * log(exp(r * trm.br) - a)))
-    }
-    out <- nlm(function(p) dev(p[1], p[2]), c(0, 0.2), hessian = TRUE)
-    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))))
-    }
-    else {
-        para <- out$estimate
-        se <- sqrt(diag(solve(out$hessian)))
+    if (conditional) {
+        dev <- function(a, r) {
+            if (a >= 1 || a < 0 || r <= 0) return(1e50)
+            ert <- exp(r * trm.br)
+            zeta <- (ert - 1)/(ert - a)
+            -2 * (lfactorial(N - 1)
+                  + (N - 2) * log(r)
+                  + N * log(1 - a)
+                  + 2 * r * sum(x[2:N])
+                  - 2 * sum(log(exp(r * x[2:N]) - a))
+                  + sum(log(1 - zeta) + (S - 1)*log(zeta)))
+        }
+    } else {
+        dev <- function(a, r) {
+            if (a >= 1 || a < 0 || r <= 0) return(1e50)
+            -2 * (lfactorial(N - 1)
+                  + (N - 2) * log(r)
+                  + (3 * N) * log(1 - a)
+                  + 2 * r * sum(x[2:N])
+                  - 2 * sum(log(exp(r * x[2:N]) - a))
+                  + r * sum(trm.br)
+                  + sum((S - 1) * log(exp(r * trm.br) - 1))
+                  - sum((S + 1) * log(exp(r * trm.br) - a)))
+        }
     }
+    out <- nlm(function(p) dev(p[1], p[2]), c(0.1, 0.2), hessian = TRUE)
+    para <- out$estimate
+    se <- sqrt(diag(solve(out$hessian)))
     Dev <- out$minimum
     cat("\nExtended Version of the Birth-Death Models to\n")
     cat("    Estimate Speciation and Extinction Rates\n\n")