Nothing
## yule.time.R (2009-02-20)
## Fits the Time-Dependent Yule Model
## Copyright 2009 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
yule.time <- function(phy, birth, BIRTH = NULL, root.time = 0,
opti = "nlm", start = 0.01)
{
opti <- pmatch(opti, c("nlm", "nlminb", "optim"))
if (is.na(opti)) stop("ambiguous argument 'opti'")
LAMBDA <- function() x
body(LAMBDA) <- body(birth)
formals(LAMBDA) <- alist(t=)
BT <- branching.times(phy)
T <- BT[1]
x <- BT[1] - BT + root.time
m <- phy$Nnode
paranam <- c(names(formals(birth)))
np <- length(paranam)
start <- rep(start, length.out = np)
## Foo is always vectorized
if (is.null(BIRTH)) {
Foo <- function(x) {
n <- length(x)
res <- numeric(n)
for (i in 1:n)
res[i] <- integrate(LAMBDA, x[i], T)$value
res
}
} else {
environment(BIRTH) <- environment()
Foo <- function(x) BIRTH(T) - BIRTH(x)
}
half.dev <- function(p) {
for (i in 1:np)
assign(paranam[i], p[i], pos = sys.frame(1))
root.term <-
if (is.null(BIRTH)) integrate(LAMBDA, x[1], T)$value
else BIRTH(T) - BIRTH(x[1])
sum(Foo(x)) + root.term - sum(log(LAMBDA(x[2:m])))
}
switch(opti,
{
out <- nlm(half.dev, start, hessian = TRUE)
est <- out$estimate
se <- sqrt(diag(solve(out$hessian)))
loglik <- lfactorial(m) - out$minimum
},{
out <- nlminb(start, half.dev)
est <- out$par
se <- NULL
loglik <- lfactorial(m) - out$objective
},{
out <- optim(start, half.dev, hessian = TRUE,
control = list(maxit = 1000), method = "BFGS")
est <- out$par
se <- sqrt(diag(solve(out$hessian)))
loglik <- lfactorial(m) - out$value
})
names(est) <- paranam
if (!is.null(se)) names(se) <- paranam
structure(list(estimate = est, se = se, loglik = loglik), class = "yule")
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.