-## CDF.birth.death.R (2010-08-09)
+## CDF.birth.death.R (2010-09-27)
## Functions to simulate and fit
## Time-Dependent Birth-Death Models
## 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:
}
}
-###.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
.makePhylo(edge[1:j, ], edge.length[1:j], i)
}
-drop.fossil <- function(phy, tol = 0)
+drop.fossil <- function(phy, tol = 1e-8)
{
n <- Ntip(phy)
x <- dist.nodes(phy)[n + 1, ][1:n]
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("\rDoing bootstrap no.", 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
+}