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

1 2 |

`t` |
time |

`a0` |
total number of type 1 particles at |

`b0` |
total number of type 2 particles at |

`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 |

a matrix of the transition probabilities

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

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)
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.