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

## Description

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

## Usage

 ```1 2``` ```bbd_prob(t, a0, b0, lambda1, lambda2, mu2, gamma, A, 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` `lambda1` birth 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` upper bound for the total number of type 1 particles `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. 2015. "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``` ```## Not run: data(Eyam) # (R, I) in the SIR model forms a birth/birth-death process loglik_sir <- function(param, data) { alpha <- exp(param[1]) # Rates must be non-negative beta <- exp(param[2]) N <- data\$S[1] + data\$I[1] + data\$R[1] # Set-up SIR model with (R, I) brates1 <- function(a, b) { 0 } brates2 <- function(a, b) { beta * max(N - a - b, 0) * b } drates2 <- function(a, b) { 0 } trans21 <- function(a, b) { alpha * b } sum(sapply(1:(nrow(data) - 1), # Sum across all time steps k function(k) { log( bbd_prob( # Compute the transition probability matrix t = data\$time[k + 1] - data\$time[k], # Time increment a0 = data\$R[k], b0 = data\$I[k], # From: R(t_k), I(t_k) brates1, brates2, drates2, trans21, A = data\$R[k + 1], B = data\$R[k + 1] + data\$I[k] - data\$R[k], computeMode = 4, nblocks = 80 # Compute using 4 threads )[data\$R[k + 1] - data\$R[k] + 1, data\$I[k + 1] + 1] # To: R(t_(k+1)), I(t_(k+1)) ) })) } loglik_sir(log(c(3.204, 0.019)), Eyam) # Evaluate at mode ## End(Not run) ```

MultiBD documentation built on May 30, 2017, 3:08 a.m.