R/calc_expected_sptree_values.R

In wadedismukes/treeduckenTools: treeduckenTools

#' Calculate expected tree length
#'
#' @details Calculates expected tree length if death_rate is non-zero given birth rate,
#' death rate, and number of tips using the following equation from Gernhard 2008.
#' @param birth_rate species birth rate (real number)
#' @param death_rate species death rate (real number)
#' @param ntips number of tips in the present (integer)
#' @param k speciation event
#' @return Numeric value of tree length given the input parameters
#' @examples
#' calc_expected_bdtree_length(0.05, 0.01, 10)
#' calc_expected_pbtree_length(0.1, 10)
calc_expected_bdtree_length <- function(birth_rate, death_rate, ntips, k = 1) {
    if (death_rate > 1e-04 && death_rate != birth_rate) {
        rho <- death_rate / birth_rate
        outer_prod <- ( ( k + 1) / birth_rate) * choose( ntips, k + 1) * ( -1) ^ k
        outer_sum <- 0
        for (i in 0:(ntips - k - 1)) {
            sum_terms <- choose( ntips - k - 1, i) *
                ( ( k + i + 1) * rho) ^ ( -1) * ( 1 / rho - 1) ^ ( k + i)
            inner_sum <- 0
            for (j in 1:(k + i)) {
                inner_sum_terms <- choose( k + i, j) * ( 1 / j) * ( -1) ^ j *
                    ( 1 - ( 1 / ( 1 - rho)) ^ j)
                inner_sum <- inner_sum + inner_sum_terms
            }
            sum_terms <- sum_terms * ( log( 1 / (1 - rho)) - inner_sum)
            outer_sum <- outer_sum + sum_terms
        }
        exp_tree_length <- outer_sum * outer_prod
    } else if (death_rate > 1e-04 && death_rate == birth_rate) {
        exp_tree_length <- ( ntips - k) / ( birth_rate * k)
    } else {
        exp_tree_length <- calc_expected_pbtree_length(birth_rate, ntips)
    }
    exp_tree_length
}

calc_expected_pbtree_length <- function(birth_rate, ntips, k = 1) {
    exp_tree_length <- 0
    for (i in (k + 1):ntips) {
        exp_tree_length <- exp_tree_length + (birth_rate * i) ^ (-1)
    }
    exp_tree_length
}

calc_expected_time_spec_event <- function(birth_rate, death_rate, ntips, k) {
    return(calc_expected_bdtree_length(birth_rate, death_rate, ntips, k))
}