adding bd.time

[ape.git] / R / CDF.birth.death.R

diff --git a/R/CDF.birth.death.R b/R/CDF.birth.death.R
index 6f58fa46ca43dba9dfb29e12a92e8223c6f5d5c3..3608a35f1d08b3bb1a63232cf54b25de1ef464bc 100644 (file)
--- a/R/CDF.birth.death.R
+++ b/R/CDF.birth.death.R
@@ -1,4 +1,4 @@
-## CDF.birth.death.R (2010-08-09)
+## CDF.birth.death.R (2010-09-27)
 
 ##    Functions to simulate and fit
 ##  Time-Dependent Birth-Death Models
@@ -8,9 +8,9 @@
 ## This file is part of the R-package `ape'.
 ## See the file ../COPYING for licensing issues.
 
+integrateTrapeze <- function(FUN, from, to, nint = 10)
 ## compute an integral with a simple trapeze method
 ## (apparently, Vectorize doesn't give faster calculation)
-integrateTrapeze <- function(FUN, from, to, nint = 10)
 {
     x <- seq(from = from, to = to, length.out = nint + 1)
     ## reorganized to minimize the calls to FUN:
@@ -41,24 +41,6 @@ integrateTrapeze <- function(FUN, from, to, nint = 10)
     }
 }
 
-###.getRHO <- function(birth, death, BIRTH = NULL, DEATH = NULL, case)
-###{
-###    ## build the RHO(), \rho(t), function
-###    switch (case,
-###            function(t1, t2) (t2 - t1)*(death - birth), # case 1
-###            function(t1, t2) # case 2
-###                integrate(function(t) death(t) - birth(t), t1, t2)$value,
-###            function(t1, t2)
-###                DEATH(t2) - BIRTH(t2) - DEATH(t1) + BIRTH(t1), # case 3
-###            function(t1, t2)
-###                death - integrate(birth, t1, t2)$value, # case 4
-###            function(t1, t2)
-###                integrate(death, t1, t2)$value - birth, # case 5
-###            function(t1, t2) death - BIRTH(t2) + BIRTH(t1), # case 6
-###            function(t1, t2) DEATH(t2) - DEATH(t1) - birth # case 7
-###        )
-###}
-
 .getRHOetINT <- function(birth, death, BIRTH = NULL, DEATH = NULL, case, fast)
 {
     ## build the RHO(), \rho(t), and INT(), I(t), functions
@@ -392,3 +374,110 @@ rbdtree <-
     foo(1L, 0)
     .makePhylo(edge[1:j, ], edge.length[1:j], i)
 }
+
+bd.time <- function(phy, birth, death, BIRTH = NULL, DEATH = NULL,
+                    ip, lower, upper, fast = FALSE, boot = 0, trace = 0)
+{
+    guess.bounds <- if (missing(lower)) TRUE else FALSE
+    BIG <- 1e10
+    PrNt <- function(t, T, x)
+    {
+        tmp <- exp(-RHO(t, T))
+        Wt <- tmp * (1 + INT(t))
+        out <- numeric(length(x))
+        zero <- x == 0
+        if (length(zero)) {
+            out[zero] <- 1 - tmp/Wt
+            out[!zero] <- (tmp/Wt^2)*(1 - 1/Wt)^(x[!zero] - 1)
+        } else out[] <- (tmp/Wt^2)*(1 - 1/Wt)^(x - 1)
+        out
+    }
+
+    case <- .getCase(birth, death, BIRTH, DEATH)
+
+    if (is.function(birth)) {
+        paranam <- names(formals(birth))
+        if (guess.bounds) {
+            upper <- rep(BIG, length(paranam))
+            lower <- -upper
+        }
+        formals(birth) <- alist(t=)
+        environment(birth) <- environment()
+        if (!is.null(BIRTH)) environment(BIRTH) <- environment()
+    } else {
+        paranam <- "birth"
+        if (guess.bounds) {
+            upper <- 1
+            lower <- 0
+        }
+    }
+
+    if (is.function(death)) {
+        tmp <- names(formals(death))
+        np2 <- length(tmp)
+        if (guess.bounds) {
+            upper <- c(upper, rep(BIG, np2))
+            lower <- c(lower, rep(-BIG, np2))
+        }
+        paranam <- c(paranam, tmp)
+        formals(death) <- alist(t=)
+        environment(death) <- environment()
+        if (!is.null(DEATH)) environment(DEATH) <- environment()
+    } else {
+        paranam <- c(paranam, "death")
+        if (guess.bounds) {
+            upper <- c(upper, .1)
+            lower <- c(lower, 0)
+        }
+    }
+
+    np <- length(paranam)
+
+    ff <- .getRHOetINT(birth, death, BIRTH, DEATH, case = case, fast = fast)
+    RHO <- ff[[1]]
+    INT <- ff[[2]]
+    environment(RHO) <- environment(INT) <- environment()
+
+    x <- branching.times(phy)
+    n <- length(x)
+    Tmax <- x[1]
+    x <- Tmax - x # change the time scale so the root is t=0
+    x <- sort(x)
+
+    foo <- function(para) {
+        for (i in 1:np)
+            assign(paranam[i], para[i], pos = sys.frame(1))
+        p <- CDF.birth.death(birth, death, BIRTH, DEATH, Tmax = Tmax,
+                             x = x, case = case, fast = fast)
+        ## w is the probability of the observed tree size (= number of tips)
+        w <- PrNt(0, Tmax, Ntip(phy))
+        ## p is the expected CDF of branching times
+        ## ecdf(x)(x) is the observed CDF
+        sum((1:n/n - p)^2)/w # faster than sum((ecdf(x)(x) - p)^2)/w
+    }
+
+    if (missing(ip)) ip <- (upper - lower)/2
+
+    out <- nlminb(ip, foo, control = list(trace = trace, eval.max = 500),
+                  upper = upper, lower = lower)
+
+    names(out$par) <- paranam
+    names(out)[2] <- "SS"
+
+    if (boot) { # nonparametric version
+        PAR <- matrix(NA, boot, np)
+        i <- 1L
+        while (i <= boot) {
+            cat("i =", i, "\n")
+            x <- sort(sample(x, replace = TRUE))
+            o <- try(nlminb(ip, foo, control = list(trace = 0, eval.max = 500),
+                            upper = upper, lower = lower))
+            if (class(o) == "list") {
+                PAR[i, ] <- o$par
+                i <- i + 1L
+            }
+        }
+        out$boot <- PAR
+    }
+    out
+}