R/dbd.R

   1 ## dbd.R (2012-03-19)
   2
   3 ##   Probability Density Under Birth--Death Models
   4
   5 ## Copyright 2012 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 dyule <- function(x, lambda = 0.1, t = 1, log = FALSE)
  11 {
  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
  17     res
  18 }
  19
  20 dbd <- function(x, lambda, mu, t, conditional = FALSE, log = FALSE)
  21 {
  22     if (length(lambda) > 1) {
  23         lambda <- lambda[1]
  24         warning("only the first value of 'lambda' was considered")
  25     }
  26     if (length(mu) > 1) {
  27         mu <- mu[1]
  28         warning("only the first value of 'mu' was considered")
  29     }
  30     if (length(t) > 1) {
  31         t <- t[1]
  32         warning("only the first value of 't' was considered")
  33     }
  34
  35     if (mu == 0) return(dyule(x, lambda, t, log))
  36
  37     ## for the unconditional case, we have to consider x=0 separately:
  38     if (!conditional) {
  39         zero <- x == 0
  40         out.of.range <- x < 0
  41     } else {
  42         out.of.range <- x <= 0
  43     }
  44
  45     res <- numeric(length(x))
  46
  47     ## the situation were speciation and extinction probabilities are equal:
  48     if (lambda == mu) {
  49         tmp <- lambda * t
  50         eta <- tmp/(1 + tmp)
  51         if (conditional) {
  52             res[] <- if (log) log(1 - eta) + (x - 1) * log(eta) else (1 - eta) * eta^(x - 1)
  53         } else { # the unconditional case:
  54             if (length(zero)) {
  55                 res[zero] <- eta
  56                 res[!zero] <- (1 - eta)^2 * eta^(x[!zero] - 1)
  57             } else res[] <- (1 - eta)^2 * eta^(x - 1)
  58         }
  59     } else { # the general case with lambda != mu
  60
  61         ## this expression is common to the conditional and unconditional cases:
  62         Ent <- exp((lambda - mu) * t)
  63
  64         if (conditional) {
  65             if (log) {
  66                 res[] <- log(lambda - mu) - log(lambda * Ent - mu) +
  67                     (x - 1) * (log(lambda) + log(Ent - 1) - log(lambda * Ent - mu))
  68             } else {
  69                 eta <- lambda * (Ent - 1)/(lambda * Ent - mu)
  70                 res[] <- (1 - eta) * eta^(x - 1)
  71             }
  72         } else { # finally, the unconditional case:
  73             eta <- lambda * (Ent - 1)/(lambda * Ent - mu)
  74             if (length(zero)) {
  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)
  78         }
  79     }
  80
  81     if (any(out.of.range))
  82         res[out.of.range] <- if (log) -Inf else 0
  83     res
  84 }
  85
  86 dbdTime <- function(x, birth, death, t, conditional = FALSE,
  87                     BIRTH = NULL, DEATH = NULL, fast = FALSE)
  88 {
  89     if (length(t) > 1) {
  90         t <- t[1]
  91         warning("only the first value of 't' was considered")
  92     }
  93
  94     if (conditional) {
  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)
  99         }
 100     } else { # the unconditional case:
 101         PrNt <- function(t, T, x)
 102         {
 103             tmp <- exp(-RHO(t, T))
 104             Wt <- tmp * (1 + INT(t))
 105             out <- numeric(length(x))
 106             zero <- x == 0
 107             if (length(zero)) {
 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)
 111             out
 112         }
 113     }
 114     case <- .getCase(birth, death, BIRTH, DEATH)
 115     ff <- .getRHOetINT(birth, death, BIRTH, DEATH, case = case, fast = fast)
 116     RHO <- ff[[1]]
 117     INT <- ff[[2]]
 118     environment(RHO) <- environment(INT) <- environment()
 119     Tmax <- t
 120     PrNt(0, t, x)
 121 }