msuchard / multibd Goto Github PK

Hi, I'm not 100% clear on the form of the returned matrix - am I right in thinking that after running

mat <- dbd_prob( ... a0 = 2, b0 = 5 ...)

(with all the other parameters)

Then mat[1,7] will be the probability of going from state (2,5) to state (0,6) (in the given time and with the given parameters)? And mat[3,9] would be the probability of going to state (2,8)?

I'm a bit confused because in the 'short example' in the readme, the indexing is done like
[1, data$I[k + 1] + 1] # To: S(t_(k+1)), I(t_(k+1))

and I don't completely understand what's going on here, particularly why the first index is 1.

If I understood the model correctly, the following rates in a birth/birth-death process should correspond to a simple (linear) birth-and-death process:

lambda_1 = 0
lambda_2 = b * lambda
mu_2 = b * mu
trans_12 = 0

but when I try to sum all the probabilities from 0 to a large number, I get a value very far from 1:

library(MultiBD)

lambda <- 1.0
mu <- 2.0

brates1 <- function(a, b) { 0 }
brates2 <- function(a, b) { b  * lambda }
drates2 <- function(a, b) { b * mu }
trans21 <- function(a, b) { 0 }

p <- bbd_prob(t=2.0, a0=0, b0=100, brates1, brates2, drates2, trans21, A=0, B=1000)

sprintf("%.20f", sum(p))

Note that most of the values are zero, so increasing the upper bound does not improve the accuracy.
I don't understand if I misinterpreted the model or there is a bug in the code.

If it is not a bug, what rates should I define to obtain a simple death-and-birth process?