## Nee 1994 equations.
make.cache.bd.nee <- function(tree=NULL, sampling.f=NULL,
unresolved=NULL, times=NULL) {
if ( !is.null(sampling.f) && !is.null(unresolved) )
stop("Cannot specify both sampling.f and unresolved")
sampling.f <- check.sampling.f(sampling.f, 1)
if ( !is.null(tree) && !is.null(unresolved) )
unresolved <- check.unresolved.bd(tree, unresolved)
if ( is.null(times) ) { # proper tree
tree <- check.tree(tree)
times <- branching.times(tree)
} else { ## Special case...
if ( !is.null(tree) )
stop("Can only supply times if tree=NULL")
}
times <- as.numeric(sort(times, decreasing=TRUE))
x <- c(NA, times)
N <- length(x)
n.node <- length(times) # == tree$Nnode
const <- lfactorial(N - 1)
list(N=N, x=x, f=sampling.f, unresolved=unresolved, n.node=n.node,
const=const)
}
make.all_branches.bd.nee <- function(cache, control) {
N <- cache$N
x <- cache$x
f <- cache$f
unresolved <- cache$unresolved
sum.x.3.N <- sum(x[3:N])
function(pars, intermediates, preset=NULL) { # preset ignored
if ( pars[2] == pars[1] )
pars[1] <- pars[1] + 1e-12
r <- pars[[1]] - pars[[2]]
a <- pars[[2]] / pars[[1]]
if ( f < 1 )
loglik <-
(N - 2) * log(f*abs(r)) + N * log(abs(1 - a)) + r * sum.x.3.N -
2*sum(log(abs(f*exp(r * x[2:N]) - a + 1 - f)))
else
loglik <-
(N - 2) * log( abs(r)) + N * log(abs(1 - a)) + r * sum.x.3.N -
2*sum(log(abs( exp(r * x[2:N]) - a )))
if ( !is.null(unresolved) ) {
ert <- exp(r * unresolved$t)
loglik <- loglik +
sum((unresolved$n-1) * (log(abs(ert - 1)) - log(abs(ert - a))))
}
## At some point, this should be e.root only. See derivations.pdf
e <- log(f * f * r * (1 - a)) -
2*log(abs(exp(-r * x[2])*(a - 1 + f) - f))
list(vals=c(loglik, e))
}
}
rootfunc.bd.nee <- function(res, pars, condition.surv,
intermediates, const) {
if ( intermediates )
stop('intermediates cannot be produced -- use method="ode"')
vals <- res$vals
loglik <- vals[[1]] + const
if ( !condition.surv ) # notice this is opposite to usual!
loglik <- loglik + vals[[2]]
loglik
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.