Exponential and gamma distributions I am trying to make a simulation of bacterium splits is exponentially distributed in R.
I am assuming that initially, bacterias have an exponential distribution. Let's say:
N_{0} = 1 
N_{1} = 2^r with r ~ exp(\lambda)
N_{2} = N_{1} * 2^r
...
N_{t} = N_{t-1} *...* N_{1} * 2^r

So, what I would like to do is a simulation where I am able to see this. In total I would need to do two simulations, the first part showing the exponential simulation of all the
N_{t}, and then a second simulation showing log2(N_{t}).
For the first part, I tried, first simulating the exponential distribution of all the N_{t}:
X <- data.frame(Cell=paste("C", 1:1000, sep=""), DT=1)
X$N <- 1
TOT <- sum(X$N)
cell.max <- 10e10
T <- 1000
t <- 1
while (TOT <= cell.max) {
    X[,t+3]<- X[,t+2] * (2 ^ (rexp(nrow(X), X$DT)))
    TOT <- sum(apply(X[,-c(1:2), drop = F], 1, sum))
    t <- t+1
} 

plot(density((t/(log2(X$N)))))

But I have errors in the plotting density, and I am not sure If I am correctly doing this simulation.
For the second simulation, I am completely lost on how to do it.
Could I have some hints/help/suggestions to do both simulations, please? All are welcome.

Relationship between sum and maximum of exponential distributed variables.
If the time untill a division for a single bacteria is exponentially distributed with rate $\lambda$ then the time untill a division for $n$ bacteria has a rate $n\lambda$ and the waiting time for the next division occuring decreases when there are more bacteria. This makes that the distribution of the waiting time is much like a sum of exponential distributions with variable rate, which is equivalent to the maximum of exponential distributions
$$max(s_1, s_2, \dots , s_n) \qquad \sim \qquad \sum_{k=1}^n t_k \\ \text{where $s_i \sim Exp(\lambda)$ and $t_k \sim Exp(\lambda(n+1-k))$}$$
(see question Proof that the limit of a series of exponential distributed random variables follows a Gumbel distribution from which is the following image)


The image gives an example for 9 bins. It illustrates the bins as a 2D plane (filled with points according to a Poisson process), in which the bins are horizontal strips in this 2D plane ... The horizontal direction could be viewed as the 'time'. We marked the balls red when they are the first to fill the bin.


The markings $s_i$ and $t_i$ in the image relate to variables ... They are used to express the waiting time to fill all the bins with at least 1 ball.

Time until n divisions
The time until n divisions is similar to the maximum of $n$ exponentially distributed variables. This follows approximately a Gumbel distribution.
The simulation below demonstrates this:

set.seed(1)

# This function simulates the cell division
# by tracking for each cell the time when it is gonna split. 
# The algorithm starts with a vector of times untill the next split
# and keeps adding times for new cells after splitting

sim = function(lambda = 1, N_end = 1000) {

  ### the initial vector of length one when n= 1
  t_vector = c(rexp(1,lambda))

  ### start settings 
  t = 0
  n = 1

  ### loop untill growth has reached N_end
  while (n < N_end) {
     # find the next cell to divide (the cell with lowest t)
     next_d = which.min(t_vector) 

     # advance the time
     t = t_vector[next_d]

     # compute new time (for next cell division) for the dividing cell
     t_vector[next_d] = t + rexp(1,lambda)

     # compute time (for next cell division) for the newly created cell
     # this time is added to the vector of times
     t_vector = c(t_vector, t + rexp(1,lambda))

     # increase the counter for the cell number
     n = n + 1
  }
  return(t)
}

### make simulations
t_vec = replicate(5*10^3,sim(N_end = 100))

### plot simulations as histogram
hist(t_vec, breaks = seq(0, ceiling(max(t_vec)), 0.2), freq = 0, main = "histogram of simulated data, with Gumbel density fit")

### compute and plot gumbel density fit
pars = evir::gumbel(t_vec)$par.ests
ts = seq(0,20,0.1)
z = (ts-pars[2])/pars[1]
lines(ts, 1/pars[1] * exp(-z-exp(-z)) )

Number of divisions in time t
There is a correspondence between waiting time $t$ until $n$ events, and number of events $n$ per time $t$. (example: https://stats.stackexchange.com/a/450135)
The number of divisions $n$ within time t is equivalent to generating the image above, and keep adding bins untill one of the bins has it's first ball for a time above t. The probability for this is $exp(-\lambda t)$ for each bin. The total number of bins/bacteria untill you get a case with $s_i > t$ is geometric distributed (the geometric distribution describes the probability distribution of the number $n$ of Bernoulli succesvol trials needed to get one failure). With $p = e^{-\lambda t}$ the probability for a failure.
$$P(N = n) = (1-p)^np = (1-e^{-\lambda t})^n e^{-\lambda t}$$
This is the same result as explained in the link from Henry's answer, but argued via a different path.
