MultiBD (version 0.2.0)

dbd_prob: Transition probabilities of a death/birth-death process

Description

Computes the transition pobabilities of a death/birth-death process using the continued fraction representation of its Laplace transform

Usage

dbd_prob(t, a0, b0, mu1, lambda2, mu2, gamma, a = 0, B, nblocks = 256, tol = 1e-12, computeMode = 0, nThreads = 4, maxdepth = 400)

Arguments

t

time

a0

total number of type 1 particles at t = 0

b0

total number of type 2 particles at t = 0

mu1

death rate of type 1 particles (a two variables function)

lambda2

birth rate of type 2 particles (a two variables function)

mu2

death rate function of type 2 particles (a two variables function)

gamma

transition rate from type 2 particles to type 1 particles (a two variables function)

a

lower bound for the total number of type 1 particles (default a = 0)

B

upper bound for the total number of type 2 particles

nblocks

number of blocks

tol

tolerance

computeMode

computation mode

nThreads

number of threads

maxdepth

maximum number of iterations for Lentz algorithm

Value

a matrix of the transition probabilities

References

Ho LST et al. 2016. "Birth(death)/birth-death processes and their computable transition probabilities with statistical applications". In review.

Examples

## Not run: 
# data(Eyam)
# 
# loglik_sir <- function(param, data) {
#   alpha <- exp(param[1]) # Rates must be non-negative
#   beta  <- exp(param[2])
# 
#   # Set-up SIR model
#   drates1 <- function(a, b) { 0 }
#   brates2 <- function(a, b) { 0 }
#   drates2 <- function(a, b) { alpha * b     }
#   trans12 <- function(a, b) { beta  * a * b }
# 
#   sum(sapply(1:(nrow(data) - 1), # Sum across all time steps k
#              function(k) {
#                log(
#                  dbd_prob(  # Compute the transition probability matrix
#                    t  = data$time[k + 1] - data$time[k], # Time increment
#                    a0 = data$S[k], b0 = data$I[k],       # From: S(t_k), I(t_k)
#                    drates1, brates2, drates2, trans12,
#                    a = data$S[k + 1], B = data$S[k] + data$I[k] - data$S[k + 1],
#                    computeMode = 4, nblocks = 80         # Compute using 4 threads
#                  )[1, data$I[k + 1] + 1]                 # To: S(t_(k+1)), I(t_(k+1))
#                )
#              }))
#   }
# 
#   loglik_sir(log(c(3.204, 0.019)), Eyam) # Evaluate at mode
# ## End(Not run)

# Birth-death-shift model for transposable elements

lam = 0.0188; mu = 0.0147; v = 0.00268; # birth, death, shift rates

drates1 <- function(a, b) { mu * a }
brates2 <- function(a, b) { lam * (a + b) }
drates2 <- function(a, b) { mu * b }
trans12 <- function(a, b) { v * a }

# Get transition probabilities
p <- dbd_prob(t = 1, a0 = 10, b0 = 0,
              drates1, brates2, drates2, trans12,
              a = 0, B = 50)