I am trying to generate random deviates for the population size at time $t$ for a birth-death process with constant birth and death rates per individual and initial size $N_0 \gt 0$.
For the simple case of a pure birth process, this can be done drawing random deviates from a negative binomial (e.g., see https://math.stackexchange.com/questions/695192/elementary-proof-of-geometric-negative-binomial-distribution-in-birth-death-pr and references therein). Similarly, this can also be done for a birth-death model with immigration and initial population size of 0 (e.g., Ben Bolker's "Ecological models and data in R", p. 124, and Bailey's "The Elements of Stochastic Processes with Applications to the Natural Sciences", p. 99).
But I want to do it with arbitrary initial size, $N_0 \gt 0$. Based on the paper by Mather et al., 2012, "Fast stochastic algorithm for simulating evolutionary population dynamics", Bioinformatics (their algorithm 2) and by looking at the closed form solution for the birth-death process[1] I think this can be done in the following way:
First, draw a random number (call it $m$) from a binomial distribution of size $N_0$ with probability of success = 1 - probability of extinction, where probability of extinction = $(\mu e^{(\lambda-\mu)t} - 1)/(\lambda e^{(\lambda-\mu)t} - \mu)$, where $\lambda$ is the birth rate (per individual) and $\mu$ is the death rate (per individual). (This is the same as $\alpha$ in the expression in Bailey, p. 94).
If $m = 0$ end of story: population is extinct. Otherwise, draw a random number (call it $n$) from a negative binomial of size $m$ and some probability of success. The population size is $m + n$. The problem is that I am not quite sure what that probability of success is supposed to be in this negative binomial (I think it is not $e^{-(\lambda - \mu)t}$).
(I am aware I could use something like Gillespie's or Gibson and Bruck's Next Reaction Method, and related approaches. But I do not care about things that happen between the predefined $t$ interval; I just want the population size at $t$ in one single operation. Likewise, I want to avoid using an approximation via a sum of the intensities from the number of death and birth events from a Poisson; I'd like to get an exact simulation, which I think is possible).
[1] This is, e.g., equation 8.47 in p. 94 of Bailey, also reproduced in several other places, such as p. 566 of "Transition probabilities for general birth–death processes with applications in ecology, genetics,and evolution", by Crawford and Suchard, 2012, etc.
