Extinction Probability Part I: Theory

Now that we have a handle on the PGF, we can try to use it to solve our original question: what is the probability of extinction for our branching process?

Let’s start by finding the PGF of a random sum of random variables. As you might be able to guess, this will be useful for branching processes, where the total reproduction of any single generation is a random sum of random variables, (it’s a ‘random sum’ because it depends on how many members of that generation there are). Let \(Z=X_1+X_2+\ldots +X_N\), where the \(X\) terms are i.i.d. random variables with support 0, 1, 2,… and \(N\) is a random variable with support 0,1,2,…

What is the PGF of \(Z\)? Let’s begin with the definition: \[ G_Z(s) = E(s^Z) = \sum_{z=0}^\infty P(Z = z)s^z. \] Is there anything we wish we knew? Well, to start, it would be helpful to know \(N\), since then we would at least know how many \(X\) terms go into \(Z\). Let’s try, then, conditioning on \(N\), (being sure to sum over all values of \(N\) ): \[ =\sum_{n=0}^\infty \sum_{z=0}^\infty P(Z=z\ | \ N=n)P(N=n)s^z. \]

Since \(P(N=n)\) has nothing to do with \(Z\), we can factor that out of the first sum: \[ =\sum_{n=0}^\infty P(N=n)\sum_{z=0}^\infty P(Z=z\ | \ N=n)s^z. \]

Hold on a moment … that \(\sum_{z=0}^\infty P(Z=z\ |\ N=n)s^z\) term is just the PGF of \(Z\), conditioning on \(N=n\). Further, because \(N=n\) is fixed (we are conditioning on it), \(Z\) is just the sum of \(n\) total \(X\) random variables. So, we are looking for the PGF of a random variable \(Z\) that is the sum of independent random variables \(X_1+X_2+\ldots +X_n\). As we have seen above, the PGF of \(Z\) is just, then, the product of the PGFs of the \(X\) terms. Since the \(X\) terms have identical distributions (they are iid.), all of their PGFs are the same, and we can thus just write: \[ = \sum_{n=0}^\infty P(N = n)(G_X (s))^n. \] Now - and here’s where it gets really funky - doesn’t this equation just look like the PGF of \(N\), except we have the term \(\left( G_X(s)\right)\) where \(s\) would usually be? That means, dear reader, that this comes out to: \[ G_N\left( G_X(s)\right). \] Incredible! The PGF of \(Z\) is just the PGF of \(X\) fed into the PGF of \(N\). Formally, repeating the prompt from above, if \(Z=X_1+X_2+\ldots +X_N\), where the \(X\) terms are i.i.d. random variables with support 0, 1, 2, … and \(N\) is a random variable with support 0,1,2, … , where the terms are i.i.d. random variables with support and \(N\) is a random variable with support 0, 1, 2,…, then we have: \[ G_Z(s)=G_N\left( G_X(s)\right). \] Exercise 8 Let \(X_1, X_2, \ldots\) be a sequence of i.i.d. Bernoulli random variables with parameter \(p\). Let \(N\) be a Poisson random variable with parameter \(\lambda\), which is independent of the \(X_i\).

  1. Find the PGF of \(Z=\sum_{i=1}^N X_i\).

  2. Use (a) to identify the distribution of \(Z\).

Let’s now see how this property can be useful with extinction probabilities for branching processes. Let’s define our branching process as follows: \[ Z_t = X_{1,t−1} + X_{2,t−1}+. . . +X_{Z_{t−1},t−1}, \]

where \(Z_t\) is the size of the population at time \(t\) and \(X_{i,t-1}\) gives the number of offspring produced by the \(i^{th}\) member of the population at time \(t-1\) (naturally, the total sum of these offspring constitute the population at time \(t\)). There are \(Z_{t-1}\) of these \(X\) terms because, well, there are \(Z_{t-1}\) cells in the population at time \(t-1\)!

Since we have that \(Z_t\) is the sum of a random number \(Z_{t-1}\) of random variables, we have, using our result from above, that the PGF of \(Z_t\) is given by: \[ G_{Z_t}(s)=G_{Z_{t-1}}\left( G_{Z_0}(s)\right). \] We write \(G_{Z_0}(s)\) here, or the generating function of the first cell in the population, because that is the generating function of all of the \(X\) terms; each \(X\) term generates offspring independently and with the same distribution as the first cell!

From here, we have some interesting iterations. We can use the above formula to find the PGF of \(Z_{t-1}\) in the same fashion: \[ G_{Z_{t−1}} (s) = G_{Z_{t−2}}\left(Π_{Z_0} (s)\right), \] which, if we plug in to our expression for \(G_{Z_t}(s)\) above, yields: \[ G_{Z_{t}} (s) = G_{Z_{t−2}}\left(G_{Z_0} (G_{Z_0}(s))\right). \]

If we continue in this way, we will get: \[ G_{Z_t}(s) = G_{Z_0}\left(G_{Z_0}( \ldots G_{Z_0}(s)\ldots)\right), \] where there are \(t\) of the \(G_{Z_0}\) terms.

It’s certainly complicated up to this point; it might help to review the following Video 1 recap before pressing on.

Excellent! We now have the PGF of the branching process \(Z_t\). Let’s now turn towards using this information to find \(p_e\), or the probability that the population goes extinct (with the help of material from this source).

We will need two key facts to find \(p_e\):

  1. Recall that \(G_{Z_n}(0)=P(Z_n=0)\) (from above; plugging into the PGF of \(Z\) gives \(P(Z=0)\). We also know that: \[ p_e=\lim_{n\to \infty} P(Z_n=0), \]

because the right hand side of this equation is essentially giving the probability of the process going to zero as \(n\) goes to \(\infty\) (lots of time passes!) which is equivalent to the extinction probability \(p_e\). Plugging in \(G_{Z_n}(0)=P(Z_n=0)\) we have: \[ p_e=\lim_{n\to \infty} G_{Z_n}(0). \]

  1. Recall what we showed above: \[ G_{Z_n}(s) = G_{Z_0}\left(G_{Z_0}( \ldots G_{Z_0}(s)\ldots)\right). \] This fact implies that: \[ G_{Z_{n+1}}(s) = G_{Z_0}\left(G_{Z_n}(0)\right). \] That is, we just add another \(G_{Z_0}\) to go from \(G_{Z_n}(0)\) to \(G_{Z_{n+1}}(0)\).

With these two facts, we can begin a “chain of equality” that will lead to an equation for \(p_e\). Thanks to fact (1.) we have \[ p_e=\lim_{n\to \infty} G_{Z_n}(0)=\lim_{n\to \infty} G_{Z_{n+1}}(0), \]

where the second equality means that adding 1 in the subscript doesn’t change anything here ( \(n\) and \(n+1\) look much the same in the limit as \(n\) goes to \(\infty\)).

Next, using fact (2.), we have \[ \lim_{n\to \infty} G_{Z_{n+1}}(0)=\lim_{n\to \infty} G_{Z_0}(G_{Z_{n}}(0)). \] Now, \[ \lim_{n\to \infty}G_{Z_0}(G_{Z_n}(0))=G_{Z_0}(\lim_{n\to \infty}G_{Z_n}(0)). \] We can do this because \(G_{Z_0}\), in the processes we are working with, is a continuous function. Finally, plugging in the result from (1.), we have: \[ G_{Z_0}\left( \lim_{n\to \infty} G_{Z_n}(0)\right)=G_{Z_n}(p_e), \] and thus, the two ends of our chain of equality are: \[ p_e=G_{Z_0}(p_e). \] Wow! This might not seem like much, but it means that we can set \(p_e\) equal to the PGF of the offspring distribution, \(G_{Z_0}\), with \(p_e\) as the input, and solve the resulting equatiion for \(p_e\). That is, we have an easy way to find our extinction probability by using the PGF!

We’ll do some examples shortly to let the power of this result sink in, but first we need a bit more housekeeping. Specifically, we can assert that \(p_e\), which is a probability and thus is between 0 and 1, is the smallest solution between 0 and 1 to the equation: \[ w=G_{Z_0}(w). \] That is, there is no other valid probability that is the solution to the above and is smaller than \(p_e\).

We can prove this with some clever math. Imagine some value $k4 (a value between 0 and 1 ) that was also a solution, so we had: \[ k=G_{Z-0}(k). \]

We know that, from our original definitions, \(G_{Z_0}\) is an increasing function, so since \(0\le k\) we must have: \[ G_{Z_0}(0)\le G_{Z_0}(k)=k. \] We can continue in this way, addinig \(G_{Z_0}\) terms to both sides. The \(k\) doesn’t change, but we start to get \[ G_{Z_0}\left(G_{Z_0} (. . . G_{Z_0} (0). . .)\right) ≤ k \]

Plugging \(G_{Z_n}(0)=G_{Z_0}\left(G_{Z_0} (. . . G_{Z_0} (0). . .)\right) ≤ k\) from above, we get: \[ G_{Z_n}(0)\le k, \] and letting \(n\) go to infinity, we get: \[ \lim_{n\to \infty} G_{Z_0}(0)\le k. \]

Of course, remember from fact (1.) above that we have \(\lim_{n\to \infty} G_{Z_0}(0)=p_e\), so we can just write: \[ p_e\le k, \]

Which proves that \(p_e\) is smaller than any other solution to: \[ w=G_{Z_0}(w). \]

This will be very useful if we find many solutions to this problem; just take the smallest one to get the correct probability of extinction!

Exercise 9 Given a branching process with the following offspring distributions, determine the extinction probability \(p_e\). Suppose \(Z_0=1\).

  1. \(p_0=0.25\), \(p_1=0.4\), \(p_2=0.35\).

  2. \(p_0=0.5\), \(p_1=0.1\), \(p_3=0.4\).

Exercise 10 Consider the branching process with offspring distribution as in Exercise 9 and suppose \(Z_0=1\). What is the probability that the population is extinct in the second generation \((Z_2=0)\), given taht it did not die out in the first generation \((Z_1>0)\)?

Exercise 11 A branching process has offspring distribution \((p_0, p_1, p_3)=(0.25, 0.25, 0.5)\). Find the following:

  1. The offspring PGF \(G_{Z_0}(s)\).

  2. The extinction probability \(p_e\).

  3. \(G_{Z_2}(s)\).

  4. \(P(Z_2=0)\).

Extinction Probability Part II: Examples

Let’s do some examples of finding the extinction probability by using the equation \(p_e=G_{Z_0}(p_e)\) to introduce some concreteness. Specifically, let’s go back to our process from Lecture 4 before:

  1. \(p_0 = 1/4\), \(p_1 = 1/2\), and \(p_2 = 1/4\) (and \(p_i = 0\) for \(i \ge 3\)).

Let’s calculate the PGF for this process: \[ E\left(s^{Z_0}\right)=G_{Z_0}(s)=\sum_{z=0}^2 P(Z_0=z)s^z = \frac{1}{4}+\frac{s}{2}+\frac{s^2}{4}. \]

We then use our result from earlier to solve for \(p_e\): \[ p_e=G_{Z_0}(p_e)=\frac{1}{4}+\frac{s}{2}+\frac{s^2}{4}. \] We can use R to code the quadratic equation and find the roots of the above:

FUN_quad <- function(a, b, c) {
 pos_root <- (-b + sqrt((b ^ 2) - 4 * a * c)) / (2 * a)
 neg_root <- (-b - sqrt((b ^ 2) - 4 * a * c)) / (2 * a)
 print(pos_root); print(neg_root)
}
FUN_quad(1 / 4, -1 / 2, 1 / 4)
## [1] 1
## [1] 1

We know that \(p_e\) is the smaller of these two solutions, and the solutions are the same, so \(p_e=1\). That is, the population will go extinct with certainty! Let’s analyze a branching process with a more general structure.

  1. Xou may have noticed that the offspring distribution of the process in the previous example was \(Bin(2,1/2)\). Can we find a general solution for a process with offspring distribution \(Bin(2,p)\)? That is, what happens to the probability of extinction as \(p\), the probability of success (in each of the two independent trials for generating offspring), changes.

Recall that we found the PGF of a \(Bin(n,p)\) random variable earlier, so we can use that result here: \[ G_{Z_0}(s)=(1-p+p\cdot s)^2. \]

Using our result from earlier to find \(p_e\) and letting \(q=1-p\), we can write: \[ G_{Z_0}(p_e)=(q+p\cdot p_e)^2=p_e. \] The two solutonns are 1 and \(\frac{q^2}{p^2}\). Remember that \(p_e\), the probability of extinction, is the minimum of our solutions. Therefore, we have \[ p_e=\min\left\{1, \ \frac{q^2}{p^2}\right\}. \] Xou may check out the video for some details that might help build understanding.

Let’s see if we can get some intuition around this. If \(p\) is very large, then get \(\frac{q^2}{p^2}\) gets small, which means \(p_e\) gets small; this makes sense, because \(p\) is the probability of generating offspring (on each of the 2 independent trials) and thus a higher value of \(p\) means a lower extinction probability \(p_e\) (and vice versa for \(q\); when \(p\) is large, \(p_e\) is large). It also makes sense that \(p_e\) maxes out at 1; it’s a probability, it can’t go any higher than 1.

We can easily visualize this result in R. Unsurprisingly, a larger \(p\) reduces the probability of extinction; interestingly, if \(p\le 0.5\), then we have \(p_e=1\). That is, if \(p\) is 0.5 or less, the population will go extinct. This is illustrated next.

library("xts")
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
#library("roll")
library("ggplot2")
library("data.table")
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:xts':
## 
##     first, last
p <- seq(from = .01, to = 1, by = .01)
data <- data.table(p = p,
 p_e = (1 - p) ^ 2 / p ^ 2)
data$p_e[data$p_e > 1] <- 1
ggplot(data, aes(x = p, y = p_e)) +
 geom_point() +
 theme_bw() +
 xlab("p") + ylab("Extinction Probability") +
 ggtitle("Extinction Probability for Bin(2, p) offspring distribution")