dbd_prob: Transition probabilities of a death/birth-death process In MultiBD: Multivariate Birth-Death Processes

Description

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

Usage

 1 2 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

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 ## Not run: data(Eyam) loglik_sir <- function(param, data) { alpha <- exp(param) # Rates must be non-negative beta <- exp(param) # 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)

MultiBD documentation built on May 2, 2019, 11:50 a.m.