Chapter 9 Branching process
The branching process is a model of population growth where individuals independently reproduce and then are replaced by the next generation. The stochastic process is \((X_n)_{n\in\mathbb N_0}\) where \(X_n\) is the population size at generation \(n\). We will usually take \(X_0=1\), meaning that the population starts with one original ancestor.
For each individual in the population, independently of all other individuals, the number of offspring is drawn from the same offspring distribution \[ \mathbb P(Y=k)=p_k,\quad k=0,1,2,\ldots. \] Here \(Y\) is a generic random variable representing the number of children left by one individual. The numbers \(p_0,p_1,p_2,\ldots\) must satisfy \[ \sum_{k=0}^\infty p_k=1. \] The case \(p_0>0\) is important because it means there is a positive probability that a single individual has no children, so extinction is possible. If \(p_0=0\), extinction is impossible once the population is positive.
If \(X_n=k\), then the next generation is the sum of the offspring counts from those \(k\) individuals: \[ X_{n+1}=\sum_{j=1}^{X_n}Y_j^{(n)}, \] where the variables \(Y_j^{(n)}\) are independent copies of \(Y\). If \(X_n=0\), then the sum is empty and \(X_{n+1}=0\). Thus \(0\) is an absorbing state.
This process is time-homogeneous because the same reproduction rule is used at every generation. In Markov chain notation, \[ \mathbb P(X_{n+1}=j\mid X_n=i) \] does not depend on \(n\).
9.1 Mean offspring number
The most important single number associated with the offspring distribution is its mean \[ \mu=\mathbb E(Y)=\sum_{k=0}^\infty k p_k. \] We can derive the expected size of the population at time \(n\) by conditioning on \(X_n\). If \(X_n=k\), then \[ \mathbb E(X_{n+1}\mid X_n=k) =\mathbb E\left(\sum_{j=1}^k Y_j^{(n)}\right) =k\mu. \] Therefore, \[ \mathbb E(X_{n+1}\mid X_n)=\mu X_n. \] Taking expectations on both sides gives \[ \mathbb E(X_{n+1})=\mu\mathbb E(X_n). \] Since \(X_0=1\), we obtain \[ \mathbb E(X_n)=\mu^n. \]
Expected population size. For a branching process with \(X_0=1\) and mean offspring number \(\mu\), \[ \mathbb E(X_n)=\mu^n. \] More generally, if \(X_0=m\), then \(\mathbb E(X_n)=m\mu^n\).
This formula is useful, but it does not by itself tell us whether extinction occurs. When \(\mu>1\), the expected population size grows without bound, but there can still be a positive probability of extinction. When \(\mu\leq 1\), the population goes extinct with probability one under the usual nontrivial assumptions.
9.2 Extinction and survival
Let \[ \pi_0=\mathbb P(X_n=0\text{ eventually}\mid X_0=1) \] be the probability that the population eventually goes extinct when starting from one individual. The event of survival forever is \[ \{X_n>0\text{ for all }n\geq0\}, \] and its probability is \(1-\pi_0\).
If \(\mu<1\), then extinction follows quickly from Markov’s inequality: \[ \mathbb P(X_n>0)=\mathbb P(X_n\geq1)\leq \mathbb E(X_n)=\mu^n. \] As \(n\to\infty\), \(\mu^n\to0\), so \[ \mathbb P(X_n>0)\to0. \] This implies extinction with probability one.
The critical case \(\mu=1\) is more delicate. Except for the trivial case \(p_1=1\), where \(X_n=1\) forever, the branching process also goes extinct with probability one.
Extinction criterion. For a nontrivial branching process with \(X_0=1\), extinction occurs with probability one if \(\mu\leq1\). If \(\mu>1\), then the process survives forever with positive probability, so \(\pi_0<1\).
9.3 Probability generating functions
A useful tool for discrete nonnegative random variables is the probability generating function. If \(Y\) has offspring distribution \(p_0,p_1,p_2,\ldots\), define \[ p(s)=\mathbb E(s^Y)=\sum_{k=0}^\infty p_k s^k, \] for \(0\leq s\leq1\). We use the convention \(0^0=1\), so \(p(0)=p_0\). Also, \[ p(1)=\sum_{k=0}^\infty p_k=1, \] and \[ p'(1)=\sum_{k=1}^\infty k p_k=\mu. \]
The probability generating function is increasing on \([0,1]\), and it is convex since \[ p''(s)=\sum_{k=2}^\infty k(k-1)p_k s^{k-2}\geq0. \] This convexity is what gives the useful graphical picture for extinction probabilities.
Probability generating function. The probability generating function for the offspring distribution is \[ p(s)=\sum_{k=0}^\infty p_k s^k=\mathbb E(s^Y). \] It satisfies \(p(0)=p_0\), \(p(1)=1\), and \(p'(1)=\mu\).
9.4 The equation for the extinction probability
We can derive an equation for \(\pi_0\) by conditioning on the first generation. If \(X_1=k\), then the population consists of \(k\) independent family lines. Each family line eventually goes extinct with probability \(\pi_0\). Thus, given \(X_1=k\), the whole population eventually goes extinct with probability \(\pi_0^k\). Therefore, \[ \pi_0 =\sum_{k=0}^\infty \mathbb P(\text{extinction}\mid X_1=k)p_k =\sum_{k=0}^\infty \pi_0^k p_k =p(\pi_0). \] So the extinction probability is a solution of \[ s=p(s). \] Since \(p(1)=1\), the equation always has the solution \(s=1\). The question is whether there is a smaller solution in \([0,1)\).
Extinction probability equation. The extinction probability \(\pi_0\) is the smallest solution in \([0,1]\) of \[ s=p(s). \] If \(\mu\leq1\), this smallest solution is \(\pi_0=1\). If \(\mu>1\), then there is a solution \(\pi_0\in[0,1)\), and this is the extinction probability.
The graphical interpretation is to plot \(y=p(s)\) and \(y=s\) on the same axes. The extinction probability is the first intersection point in \([0,1]\). Since \(p(s)\) is convex and \(p(1)=1\), the slope \(p'(1)=\mu\) determines whether a smaller intersection can occur. When \(\mu\leq1\), the only intersection in \([0,1]\) is \(s=1\). When \(\mu>1\), the graph of \(p(s)\) crosses the line \(y=s\) at a smaller value, and that smaller value is \(\pi_0\).
9.5 Finite-time extinction probabilities
Let \[ \theta_n=\mathbb P(X_n=0\mid X_0=1). \] These probabilities give the chance that the population has gone extinct by generation \(n\). We have \[ \theta_0=0, \] because the initial population is one. Also, \[ \theta_1=p_0=p(0). \] The same conditioning argument used above gives the recursion \[ \theta_{n+1}=p(\theta_n). \] Thus \[ \theta_n=(\underbrace{p\circ p\circ\cdots\circ p}_{n\text{ compositions}})(0). \] The sequence \(\theta_n\) is increasing and bounded above by \(1\), so it converges. Its limit is the extinction probability: \[ \lim_{n\to\infty}\theta_n=\pi_0. \]
9.6 Generating function of \(X_n\)
Let \[ p_n(s)=\mathbb E(s^{X_n}) \] be the probability generating function of \(X_n\). Since \(X_0=1\), we have \[ p_0(s)=s. \] Conditioning on \(X_n\) gives \[ \begin{aligned} p_{n+1}(s) &=\mathbb E(s^{X_{n+1}})\\ &=\mathbb E\left[\mathbb E(s^{X_{n+1}}\mid X_n)\right]. \end{aligned} \] If \(X_n=k\), then \(X_{n+1}\) is the sum of \(k\) independent offspring variables, so \[ \mathbb E(s^{X_{n+1}}\mid X_n=k)=\big(p(s)\big)^k. \] Therefore, \[ p_{n+1}(s)=\mathbb E\left[(p(s))^{X_n}\right]=p_n(p(s)). \] Inductively, \[ p_n(s)=p(p(\cdots p(s)\cdots)), \] with \(n\) total compositions of \(p\).
In particular, \[ \mathbb P(X_n=0)=p_n(0), \] which agrees with the recursion \(\theta_{n+1}=p(\theta_n)\).
9.7 Examples
Example. Suppose \[ p_0=\frac14,\qquad p_2=\frac34, \] and \(p_k=0\) otherwise. Then \[ \mu=0\cdot\frac14+2\cdot\frac34=\frac32>1, \] so survival has positive probability.
The generating function is \[ p(s)=\frac14+\frac34s^2. \] The extinction probability solves \[ s=\frac14+\frac34s^2. \] Equivalently, \[ 3s^2-4s+1=0=(3s-1)(s-1). \] Thus the possible solutions in \([0,1]\) are \(s=1/3\) and \(s=1\). The extinction probability is the smaller one: \[ \pi_0=\frac13. \] The survival probability is \[ 1-\pi_0=\frac23. \]
Example. Suppose \[ p_0=p_1=p_2=p_3=\frac14. \] Then \[ \mu=0\cdot\frac14+1\cdot\frac14+2\cdot\frac14+3\cdot\frac14=\frac32>1, \] so survival is possible. The generating function is \[ p(s)=\frac14(1+s+s^2+s^3). \] The equation \(s=p(s)\) becomes \[ s=\frac14(1+s+s^2+s^3), \] or \[ s^3+s^2-3s+1=0. \] Factoring gives \[ s^3+s^2-3s+1=(s-1)(s^2+2s-1). \] The roots are \[ s=1,\qquad s=-1+\sqrt2,\qquad s=-1-\sqrt2. \] Only \(1\) and \(\sqrt2-1\) are in \([0,1]\), so \[ \pi_0=\sqrt2-1. \] The survival probability is \[ 1-\pi_0=2-\sqrt2. \]
Example. Suppose
\[
p_0=0.2,\qquad p_1=0.1,\qquad p_2=0.4,\qquad p_3=0.3.
\]
Then
\[
p(s)=0.2+0.1s+0.4s^2+0.3s^3.
\]
To find the extinction probability, we solve \(s=p(s)\), or
\[
0=0.2-0.9s+0.4s^2+0.3s^3.
\]
In R, we can solve this polynomial using polyroot().
## [1] 0.257334-7.888609e-31i 1.000000+9.525013e-17i -2.590667-9.525013e-17iThe mean number of offspring is \[ \mu=0(0.2)+1(0.1)+2(0.4)+3(0.3)=1.8>1. \] Thus the extinction probability is the root in \([0,1)\), approximately \[ \pi_0\approx0.257334. \]
9.8 Simulating a branching process in R
The following code simulates a branching process with a finite offspring distribution. The vector p stores the probabilities \(p_0,p_1,\ldots,p_m\).
p <- c(0.25, 0, 0.75) # p0, p1, p2
maxgen <- 10
x <- numeric(maxgen + 1)
x[1] <- 1
for (n in 1:maxgen){
if (x[n] == 0){
x[n+1] <- 0
} else {
offspring <- sample(0:(length(p)-1), size=x[n], replace=TRUE, prob=p)
x[n+1] <- sum(offspring)
}
}
plot(0:maxgen, x, type="b", pch=20, xlab="generation", ylab="population size")
To estimate the extinction probability by simulation, we can run many independent copies and record the proportion that hit zero before a large final generation.
sim_branch <- function(p, maxgen = 10, max_pop = 1e6){
x <- 1
for (n in 1:maxgen){
if (x == 0) return(0)
# If the population has become very large, treat it as surviving
# for purposes of estimating extinction probability.
if (x > max_pop) return(Inf)
offspring <- sample(0:(length(p)-1), size = x, replace = TRUE, prob = p)
x <- sum(offspring)
}
x
}
p <- c(0.25, 0, 0.75)
nsim <- 1000
final_sizes <- replicate(nsim, sim_branch(p, maxgen = 100))
mean(final_sizes == 0)## [1] 0.341For this example the theoretical extinction probability is \(1/3\), and the simulation estimate should be reasonably close when nsim is large.