There's only one Mr Tripletoddletrouble. In fact, unless he has
a son to pass his surname down to, he'll be the last
Mr Tripletoddletrouble. Social mores of his time and place sadly
disallow even such an exquisite surname to survive by passing through the
female line.
Mr Tripletoddletrouble has a rare and mathematically convenient
genetic condition, which any future generations of
Tripletoddletroubles will inherit: if he fathers any sons at all,
and there's a 50% chance he will, it will be one set of triplets.
So at every step of the family tree we see, equally likely, either three sons or none.
The bad news for
onomatologists is
this leaves a 50% chance of a wonderful name becoming extinct among the
next generation. The good news is that every Mr
Tripletoddletrouble has, on average, 1.5 sons — since this is
safely above one, the expected population of Tripletoddletroubles enjoys
exponential growth, and there is a positive probability their surname
will survive forever.
What's the probability that the Tripletoddletrouble surname will, eventually, become extinct?
Here's a quick R simulation.
set.seed(123)
nsims <- 1e5
ngens <- 20
simulate_extinction <- function(ngens) {
nsurvivors <- c(1, rep(NA, ngens - 1))
for (gen in seq_len(ngens - 1)) {
nsurvivors[gen + 1] <- 3 * rbinom(1, nsurvivors[gen], 0.5)
}
extinct <- (!is.na(nsurvivors) & nsurvivors == 0) # rbinom gives NA if population huge
return(extinct)
}
pextinct <- rowMeans(replicate(n = nsims, simulate_extinction(ngens)))
plot(pextinct, xlab = "Generation number", ylab = "Probability of extinction")
abline(h = (sqrt(5) - 1)/2, col = "red")
sprintf("Estimated probability of extinction = %f", pextinct[ngens])
The red line in the plot is at
$$\varphi - 1 = \varphi^{-1} = \frac{\sqrt{5} - 1}{2} \approx 0.618034 $$
[1] "Estimated probability of extinction = 0.618150"
This is a question about branching processes. Indeed, one of the earliest investigations into their stochastic behaviour originated in Victorian concerns about the extinction of unusual surnames. The resulting Galton-Watson process paper is available online:
"Either 3 or 0 offspring, equally likely" is arguably the simplest branching process with non-trivial probability of extinction. We need at least two outcomes if chance is to play a role, including zero offspring for extinction to be possible. "Either 1 or 0 offspring" is clearly doomed: with no branches budding off, the family line becomes extinct the first time there are no sons. "Either 2 or 0 offspring, equally likely" gives a mean of exactly one offspring to replace each individual. When fate balances on this knife-edge, it turns out extinction is certain in the long run, even if the family tree successfully buds a few times. We could tweak the offspring distribution to produce a range of desired extinction probabilities, but only by introducing unequal probabilities or more than two outcomes. This set-up doesn't feel artificially "tuned" to shoe-horn in the golden ratio.
Let's find the probability of ultimate extinction, $\theta$, algebraically. Intuitively, this probability splits into two parts: either the original Mr Tripletoddletrouble has no sons and his line becomes extinct immediately, or he successfully has three sons, but each of their three lines eventually becomes extinct. Since a son is in the same position as the original Mr Tripletoddletrouble, their lines also each have extinction probability $\theta$. Since we are concerned only with direct male descendants, each line's fate is independent of the others. Given that there are three sons, the probability the surname becomes extinct is therefore $\theta^3$.
From the tree diagram, we see the probability of extinction $\theta$ must obey the equation
$$\theta = \frac{1}{2} + \frac{1}{2}\theta^3 \tag{1}$$
which we can solve (and will, shortly). First let's tie this into some wider theory of branching processes. The number of offspring of any individual is a random variable with probability distribution $p_0 = p_3 = 0.5$ and $p_n = 0$ otherwise, so its probability generating function is:
$$\Pi(s) = \sum_n p_n s^n = \frac{1}{2} + \frac{1}{2} s^3 $$
Looks familiar? No coincidence. More later. By evaluating the derivative of the pgf at $s=1$, we get the mean number of offspring. This number $R_0 = \Pi'(1)$ is important in population ecology and human demography, where it's called the net reproduction rate (it's usually defined as the mean number of daughters produced by each female, rather than sons by each male — maternity is easier to track than paternity, and in many species females can reproduce by parthenogenesis), while in epidemiology it's the basic reproduction number (mean number of infections directly generated by one infected individual, in a fully susceptible population). If $\Pi'(1) \leq 1$ then ultimate extinction is certain. If $\Pi'(1) \gt 1$ the probability of extinction is below one. We have
$$\Pi'(s) = \frac{3s^2}{2} \implies \Pi'(1) = 1.5 > 1$$
so the surname has positive probability of survival. How many Tripletoddletroubles survive in each generation? Take one individual as "generation zero", and let $Z_n$ be the number of descendants after $n$ generations. $Z_n$ is a random variable whose probability distribution can be read off from the coefficients of its pgf $\Pi_{n}(s)$, which we find by iteratively applying $\Pi$, the offspring pgf, $n$ times:
$$\Pi_{n}(s) = \Pi(\dots\Pi(\Pi(s))\dots) $$
Why? $Z_{n}$ is the sum of the offspring of the $Z_{n-1}$ survivors in the previous generation. The numbers of offspring from each survivor are independent, identically distributed (iid) random variables with pgf $\Pi$, and the number of them we are adding up has pgf $\Pi_{n-1}$, so by the rule for the pgf of the sum of a random number of iid variables (proof in this answer), $Z_n$ has pgf $\Pi_{n}(s) = \Pi_{n-1}(\Pi(s))$. For example, after two generations
$$\Pi_2(s) = \Pi(\Pi(s)) = \frac{1}{2} + \frac{1}{2} \left(\frac{1}{2} + \frac{s^3}{2} \right)^3 = \frac{9}{16} + \frac{3s^3}{16} + \frac{3s^6}{16} + \frac{s^9}{16} $$
so there's a $\frac{1}{16}$ chance of nine descendants but $\Pi_2(0) = \frac{9}{16}$ chance that extinction has already occurred. $\mathbb{E}(Z_2)$, the expected number of descendants after two generations, is found by $\Pi'_2(1) = 2.25$. It's no coincidence this equals $1.5^2$.
The mean and variance of the number of offspring from a single individual are $\mu = \Pi'(1)$ and $\sigma^2 = \Pi''(1) + \mu - \mu^2$. You can prove by induction that $\mathbb{E}(Z_n) = \mu^n$. Now it's obvious why ultimate extinction is certain when $\mu < 1$. With $\mu = 1.5$ we see exponential growth on average, despite our high chance of early extinction. Essentially, chains of surname transmission tend to either fizzle out or blow up, and $\mu = 1.5$ guarantees enough chance of blowing up that extinction is not inevitable. Good news for the Tripletoddletroubles; bad news if we switch context from surnames to infectious diseases with $R_0 > 1$. The way chains of infection can randomly either "go big or go home", rather than follow a deterministic rule like "each case infects exactly two susceptibles", relates to the epidemiological idea of overdispersion due to clustering or super-spreading events. The variance of the number of descendants after $n$ generations can be considerable, as $Z_n$ might be enormous or zero. Again by induction, we find:
$$\operatorname{Var}(Z_n) =
\begin{cases}
\frac{\mu^{n-1} \sigma^2 \left(\mu^n - 1\right)}{\mu - 1}, & \mu \neq 1 \\[2ex]
n \sigma^2, & \mu = 1
\end{cases}$$
In general, the probability of ultimate extinction is the smallest positive solution, $\theta^{*}$, of the equation $\theta = \Pi(\theta)$. That's exactly equation $(1)$ we derived above! But how did we know which solution to take? The probability of extinction by generation $n$ is $\Pi_n(0)$, since that's the constant or $s^0$ term of the pgf of $Z_n$, hence represents $\Pr(Z_n = 0)$. The probability of ultimate extinction must be $\lim_{n \to \infty} \Pi_n(0)$ which we can find using a cobweb plot of $y=\Pi(x)$ and $y=x$ for $0 \le x \le 1$. Since $\Pi(0) = p_0$, the probability an individual has no offspring, we can assume the y-intercept is between $0 \lt \Pi(0) \le 1$ (if $p_0 = 0$ then extinction is clearly impossible). So $y = \Pi(x)$ starts above $y=x$, and the first time it intersects $y=x$ must be from above. Since $\Pi(x)$ and its derivatives have only non-negative coefficients, its graph is increasing and convex on $0 \le x \le 1$. This means it can intersect $y=x$ at most twice in this interval: once from above, then again from below. $\Pi(1) = \sum p_n = 1$ so the graphs certainly intersect at $(1,1)$.
This intersection's nature depends on the slope $\Pi'(1)$, which represents the mean number of offspring $\mu$ (biologically, $R_0$). If $\Pi'(1) > 1$ it must be steeper than $y=x$ so $y = \Pi(x)$ is hitting the line from below, in which case there must have been an earlier intersection in $0 \lt x \lt 1$. If $\Pi'(1) < 1$ it's shallower so hitting from above, and there's no earlier root. If $\Pi'(1) = 1$ the two curves just touch at $(1,1)$, but $y = \Pi(x)$ must have been shallower before (its average slope over $0 \le x \le 1$ is $1 - p_0$ so below one), hence approaches the line from above and there can be no earlier root. This is why if $\mu=1$ but $p_0 > 0$, ultimate extinction has probability one.
To find $\lim_{n \to \infty} \Pi_n(0)$ graphically, read off horizontally from the y-intercept at $y = \Pi(0)$ to the $y=x$ line, where now $x = \Pi(0)$. Then read off vertically to the $y = \Pi(x)$ graph, where now $y = \Pi(\Pi(x)) = \Pi_2(x)$. Read off horizontally to the line so $x = \Pi_2(x)$. Read off vertically to the curve so $y = \Pi(\Pi_2(x)) = \Pi_3(x)$. Note that all horizontal readings are rightwards and vertical readings are upwards, since $y = \Pi(x)$ is increasing so each vertical positions is above the previous one. This procedure must converge to the first (i.e. smallest positive $x$) intersection point $x = \Pi(x)$, where $y = \Pi(x)$ hits $y = x$ from above. We illustrate the three cases $\mu = 1.5$ ($p_0 = p_3 = \frac{1}{2}$), $\mu = 1$ ($p_0 = \frac{2}{3}, p_3 = \frac{1}{3}$) and $\mu = 0.5$ ($p_0 = \frac{5}{6}, p_3 = \frac{1}{6}$). The dotted blue line is the tangent to $y = \Pi(x)$ at $(1, 1)$, and shows the role of its slope $\Pi'(1) = \mu$ in determining whether there was an earlier intersection.
We need the smallest positive solution $\theta*$ of $(1)$. Moving $\theta$ to the right-hand side and doubling to clear out the fractions, we obtain:
$$0 = \theta^3 - 2 \theta + 1 = (\theta - 1)(\theta^2 + \theta - 1)$$
The solutions are $-\varphi < \varphi^{-1} < 1$ so the smallest positive solution is $\theta^{*} = \varphi^{-1}$.
Time to reveal the "fiddle". This link to the golden ratio isn't a result
I recall seeing before,
but I reverse-engineered it by thinking about the required factorisation of
the final equation. Since
$\Pi(1) = \sum p_n = 1$, we always have $\theta = 1$ as a root of
$\theta = \Pi(\theta)$, so $(\theta - 1)$ must appear as a factor once we
set one side to zero. I also knew what quadratic I wanted to see.
After that I worked back to try to form a valid pgf. Negative
coefficients are disallowed; positive coefficients just
needed normalising so they sum to unity. I hoped
that the resulting probability distribution for the offspring would be
a "nice" one — which I think it is!
R code for cobweb plot
ngens <- 100
par(mfrow=c(1, 3), pty = "s", xaxs = "i", yaxs = "i")
for(p0 in c(1/2, 2/3, 5/6)) {
pgf <- function(x) {p0 + (1-p0)*x^3}
mu <- p0*0 + (1-p0)*3
plot(pgf, xlim = c(0,1), ylim = c(0,1), xlab = "", ylab = "",
main = paste0("Mean offspring = ", mu))
segments(0, 0, 1, 1)
abline(1 - mu, mu, col = "blue", lty = "dotted")
pextinct <- c(0, rep(NA, ngens))
for (n in seq_len(ngens)) {
pextinct[n + 1] <- pgf(pextinct[n])
segments(pextinct[n], pextinct[n], pextinct[n], pextinct[n + 1], col = "red")
segments(pextinct[n], pextinct[n + 1], pextinct[n + 1], pextinct[n + 1], col = "red")
}
print(pextinct)
}
