Every Friday, FiveThirtyEight.com posts a weekly Riddler challenge. On January 27th, 2017, the following Riddler Classic was posted:

At the beginning of time, there is a single microorganism. Each day, every member of this species either splits into two copies of itself or dies. If the probability of multiplication is p, what are the chances that this species goes extinct?

This is a textbook branching process. That is, one unit branches into a random number of new units, each of which are indistinguishable and independent from the parent unit, and thus, the branching process renews.

Notice that in a usual branching process, organisms’ subsequent generations are considered all of the offspring it will ever have before dying, where in our process, the organism splits, and hence one of the “child” offspring is in fact the parent. This technicality won’t affect our logic in what follows.

Let \(u_n\) be the probability that the population goes extinct at generation \(n\). Suppose any individual microorganism in generation \(n\) has \(k\) offpsring (which will happen with probability \(p_k\)), thus creating \(k\) new identical branching process subpopulations. The probability that all of these new subpopulations go extinct is thus \((u_{n-1})^k\). We now have the following property due to the Law of Total Probability:

\[u_n = \sum_{k = 0}^{\infty} p_k (u_{n-1})^k\]

Notice the similarity to a distribution’s probability generation function (pgf):

\[\phi(s) = E[s^X] = \int s^X dP\]

If we let \(n \rightarrow \infty\) (let’s eschew discussing convergence conditions here), we can thus find the extinction probability until the end of time by solving for the smallest \(u^*\) satisfying

\[u^* = \phi(u^*)\] The Riddler states that a single microorganism multiplies into two organisms with probability \(p\), and thus dies without splitting with probability \(1 - p\). So, \(p_2 = p\), \(p_0 = 1 - p\), and \(p_k = 0\) for any \(k \ne 0, 2\). Thus:

\[\phi(u^*) = p(u^*)^2 + (1 - p) = u^*\] And thus we get

\[u^* = \min\left\{1, \frac{1 - p}{p}\right\}\]

As proof, I also simulated this situation for various values of \(p\) in R.

library(dplyr)
library(ggplot2)

runSim <- function(i, p) {
  cells <- 1
  while (TRUE) {
    cells <- sum(2*rbinom(cells, 1, p))
    if (cells == 0) return(TRUE)
    if (cells > 1000) return(FALSE)
  }
}

p_index <- seq(0.001, 1, 0.001)

runAggSim <- function(p) mean(sapply(1:1e3, function(i) runSim(i, p)))

probs <- sapply(p_index, runAggSim)

truep <- sapply(p_index, function(p) min(1, (1 - p)/p))

df <- data.frame(p_index, probs, truep)

ggplot(df) +
  theme_minimal() +
  geom_line(aes(p_index, truep), color = 'red', size = 3, alpha = 1/2) +
  geom_line(aes(p_index, probs), size = 1/4) +
  xlab('p') + ylab('Probability of Extinction')