Maximum size of a subcritical birth-death process Beginning with a population of $n_0$ individuals, let each individual have a probability $p$ to survive until it replicates into two independent and identical individuals, where $p<\frac12$.
It follows that the population goes extinct in the long-run with probability 1, and the expected number of replications before extinction is $n_0p/(1-2p)$.
The population must reach some maximum size $N$ before going extinct, where $N\geq1$. What is the expected value of this maximum size $N$?

I have written this $\texttt{R}$ code to compute a point estimate for the mean and variance of $N$ given fixed values for $p$ and $n_0$:
rm(list=ls())

N <- 10000
maxpop <- numeric(N)
p <- 0.3
n_0 <- 2

for(n in 1:N) {
  pop <- n_0
  maxpop[n] <- n_0
  while(pop > 0) {
    if(runif(1)>p) {
      pop <- pop - 1
    }
    else {
      pop <- pop + 1
      maxpop[n] <- max(maxpop[n], pop+1)
    }
  }
}

print(mean(maxpop))
print(var(maxpop))

Perhaps it can be useful.
