In what follows, the simulations find the number of bacteria at a given time, based on the starting number of bacteria and the rate at which each bacterium splits according to an exponential distribution and then the two results from the split do the same from that point forward. Note that if the rate for each bacterium is $\lambda$ then note the following points, which are used in the simulation.
- The expected number of bacteria after time $t$ is $e^{\lambda t}$ times the number at the start. If you want this to be $2^t$ then you need $\lambda=\log_e(2)$.
- With $n$ bacteria at any point in time, the time until the next split has an exponential distribution with rate $n \lambda$.
- Starting with $N_0=1$ bacterium, the number of bacteria at time $t$ is geometrically distributed on $1,2,3,4,\ldots$ with parameter $p=e^{-\lambda t}$, so $\mathbb P(N_t=k)=e^{-\lambda t}\left(1-e^{-\lambda t}\right)^{k-1}$. It has mean $e^{\lambda t}$ and variance $e^{\lambda t}\left(e^{\lambda t}-1\right)$.
- So starting with $N_0=n$ bacteria, the number of bacteria at time $t$ is negative binomially distributed on $n,n+1,n+2,\ldots$ with parameter $p=e^{-\lambda t}$ so $\mathbb P(N_t=k)={k-1 \choose n-1}e^{-n \lambda t}\left(1-e^{-\lambda t}\right)^{k-n}$. It has mean $n e^{\lambda t}$ and variance $ne^{\lambda t}\left(e^{\lambda t}-1\right)$.
The first simulation simply splits individual bacteria until the desired time point is reached for each of them:
numberbytimeA <- function(time, rate=1, start=1){
nextsplit <- rexp(start, rate)
while(min(nextsplit) <= time){
nextsplit <- c(nextsplit, nextsplit[nextsplit <= time])
nextsplit[nextsplit <= time] <- nextsplit[nextsplit <= time] +
rexp(sum(nextsplit <= time), rate)
}
return(length(nextsplit))
}
The next simulation add one bacterium at a time until the desired time point is reached:
numberbytimeB <- function(time, rate=1, start=1){
splittime <- 0
count <- start
while (splittime <= time){
splittime <- splittime + rexp(1, count*rate)
if (splittime > time){
return(count)
}
count <- count + 1
}
}
That can be quite slow, looping many times. An attempt to speed it up could be to generate more than enough splits to take up the time provided and then find how many there are when the time has been reached.
numberbytimeC <- function(time, rate=1, start=1,
maxend=start*exp(3+time*rate)){
splittimes <- cumsum(rexp(maxend-start+1,(start:maxend)*rate))
return(start + sum(splittimes <= time))
}
The fourth uses the geometric and negative binomial distributions to provide a massive speed-up when the time is long.
numberbytimeD <- function(time, rate=1, start=1){
return(start + rnbinom(1, start, exp(-rate*time)))
}
As far as I can tell by testing, all three produce the same results (allowing for the noise associated with simulation) but at different speeds.
So here is an actual simulation, over time $3$ with your rate $\log_e(2)$ starting with $1$ bacterium, and $10^5$ simulations. It uses the second function as the one I initially encoded, though the others would do as well but with different speeds.
set.seed(2022) # so you can reproduce my results
simB <- replicate(10^5, numberbytimeB(time=3, rate=log(2), start=1))
mean(simB) # expected 2^3 = 8
# 7.996
sd(simB) # expected sqrt(8*(8-1)) about 7.483
# 7.432131
plot(table(simB)/10^5) # empirical frequency of final number of bacteria
Note that when you start with $1$ bacterium, the most likely outcome is $1$ bacterium at the end, i.e. $0$ splits. As a geometric random variable, this is what we expect to happen, no matter how long the time period is or how fast the rate is; the probability of this outcome is $e^{-\lambda t}$ and tends towards $0$ as $t$ or $\lambda$ increase but other particular outcomes are always less likely.
N_{t}? – Henry Nov 10 '22 at 14:53