3 ## Probability Density Under Birth--Death Models
5 ## Copyright 2012 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 dyule <- function(x, lambda = 0.1, t = 1, log = FALSE)
12 tmp <- exp(-lambda * t)
13 res <- if (log) log(tmp) + (x - 1) * log(1 - tmp) else tmp * (1 - tmp)^(x - 1)
14 out.of.range <- x <= 0
15 if (any(out.of.range))
16 res[out.of.range] <- if (log) -Inf else 0
20 dbd <- function(x, lambda, mu, t, conditional = FALSE, log = FALSE)
22 if (length(lambda) > 1) {
24 warning("only the first value of 'lambda' was considered")
28 warning("only the first value of 'mu' was considered")
32 warning("only the first value of 't' was considered")
35 if (mu == 0) return(dyule(x, lambda, t, log))
37 ## for the unconditional case, we have to consider x=0 separately:
42 out.of.range <- x <= 0
45 res <- numeric(length(x))
47 ## the situation were speciation and extinction probabilities are equal:
52 res[] <- if (log) log(1 - eta) + (x - 1) * log(eta) else (1 - eta) * eta^(x - 1)
53 } else { # the unconditional case:
56 res[!zero] <- (1 - eta)^2 * eta^(x[!zero] - 1)
57 } else res[] <- (1 - eta)^2 * eta^(x - 1)
59 } else { # the general case with lambda != mu
61 ## this expression is common to the conditional and unconditional cases:
62 Ent <- exp((lambda - mu) * t)
66 res[] <- log(lambda - mu) - log(lambda * Ent - mu) +
67 (x - 1) * (log(lambda) + log(Ent - 1) - log(lambda * Ent - mu))
69 eta <- lambda * (Ent - 1)/(lambda * Ent - mu)
70 res[] <- (1 - eta) * eta^(x - 1)
72 } else { # finally, the unconditional case:
73 eta <- lambda * (Ent - 1)/(lambda * Ent - mu)
75 res[zero] <- eta * mu / lambda
76 res[!zero] <- (1 - mu * eta / lambda) * (1 - eta) * eta^(x[!zero] - 1)
77 } else res[] <- (1 - mu * eta / lambda) * (1 - eta) * eta^(x - 1)
81 if (any(out.of.range))
82 res[out.of.range] <- if (log) -Inf else 0
86 dbdTime <- function(x, birth, death, t, conditional = FALSE,
87 BIRTH = NULL, DEATH = NULL, fast = FALSE)
91 warning("only the first value of 't' was considered")
95 PrNt <- function(t, T, x) {
96 tmp <- exp(-RHO(t, T))
97 Wt <- tmp * (1 + INT(t))
98 (1/Wt)*(1 - 1/Wt)^(x - 1)
100 } else { # the unconditional case:
101 PrNt <- function(t, T, x)
103 tmp <- exp(-RHO(t, T))
104 Wt <- tmp * (1 + INT(t))
105 out <- numeric(length(x))
108 out[zero] <- 1 - tmp/Wt
109 out[!zero] <- (tmp/Wt^2)*(1 - 1/Wt)^(x[!zero] - 1)
110 } else out[] <- (tmp/Wt^2)*(1 - 1/Wt)^(x - 1)
114 case <- .getCase(birth, death, BIRTH, DEATH)
115 ff <- .getRHOetINT(birth, death, BIRTH, DEATH, case = case, fast = fast)
118 environment(RHO) <- environment(INT) <- environment()