]> git.donarmstrong.com Git - ape.git/blob - R/yule.time.R
some small changes
[ape.git] / R / yule.time.R
1 ## yule.time.R (2009-02-20)
2
3 ##    Fits the Time-Dependent Yule Model
4
5 ## Copyright 2009 Emmanuel Paradis
6
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
9
10 yule.time <- function(phy, birth, BIRTH = NULL, root.time = 0,
11                       opti = "nlm", start = 0.01)
12 {
13     opti <- pmatch(opti, c("nlm", "nlminb", "optim"))
14     if (is.na(opti)) stop("ambiguous argument 'opti'")
15     LAMBDA <- function() x
16     body(LAMBDA) <- body(birth)
17     formals(LAMBDA) <- alist(t=)
18     BT <- branching.times(phy)
19     T <- BT[1]
20     x <- BT[1] - BT + root.time
21     m <- phy$Nnode
22
23     paranam <- c(names(formals(birth)))
24     np <- length(paranam)
25     start <- rep(start, length.out = np)
26
27     ## Foo is always vectorized
28     if (is.null(BIRTH)) {
29         Foo <- function(x) {
30             n <- length(x)
31             res <- numeric(n)
32             for (i in 1:n)
33                 res[i] <- integrate(LAMBDA, x[i], T)$value
34             res
35         }
36     } else {
37         environment(BIRTH) <- environment()
38         Foo <- function(x) BIRTH(T) - BIRTH(x)
39     }
40
41     half.dev <- function(p) {
42         for (i in 1:np)
43             assign(paranam[i], p[i], pos = sys.frame(1))
44         root.term <-
45             if (is.null(BIRTH)) integrate(LAMBDA, x[1], T)$value
46             else BIRTH(T) - BIRTH(x[1])
47         sum(Foo(x)) + root.term - sum(log(LAMBDA(x[2:m])))
48     }
49
50     switch(opti,
51         {
52             out <- nlm(half.dev, start, hessian = TRUE)
53             est <- out$estimate
54             se <- sqrt(diag(solve(out$hessian)))
55             loglik <- lfactorial(m) - out$minimum
56         },{
57             out <- nlminb(start, half.dev)
58             est <- out$par
59             se <- NULL
60             loglik <- lfactorial(m) - out$objective
61         },{
62             out <- optim(start, half.dev, hessian = TRUE,
63                          control = list(maxit = 1000), method = "BFGS")
64             est <- out$par
65             se <- sqrt(diag(solve(out$hessian)))
66             loglik <- lfactorial(m) - out$value
67         })
68     names(est) <- paranam
69     if (!is.null(se)) names(se) <- paranam
70     structure(list(estimate = est, se = se, loglik = loglik), class = "yule")
71 }