Archive for Galton-Watson extinction

extinction minus one

Posted in Books, Kids, pictures, R, Statistics, University life with tags , , , , , , , , , , , , , , , on March 14, 2022 by xi'an

The riddle from The Riddler of 19 Feb. is about the Bernoulli Galton-Watson process, where each individual in the population has one or zero descendant with equal probabilities: Starting with a large population os size N, what is the probability that the size of the population on the brink of extinction is equal to one? While it is easy to show that the probability the n-th generation is extinct is

\mathbb{P}(S_n=0) = 1 - \frac{1}{2^{nN}}

I could not find a way to express the probability to hit one and resorted to brute force simulation, easily coded

for(t in 1:(T<-1e8)){N=Z=1e4 
while(Z>1)Z=rbinom(1,Z,.5)
F=F+Z}
F/T

which produces an approximate probability of 0.7213 or 0.714. The impact of N is quickly vanishing, as expected when the probability to reach 1 in one generation is negligible…

However, when returning to Dauphine after a two-week absence, I presented the problem with my probabilist neighbour François Simenhaus, who immediately pointed out that this probability was more simply seen as the probability that the maximum of N independent geometric rv’s was achieved by a single one among the N. Searching later a reference for that probability, I came across the 1990 paper of Bruss and O’Cinneide, which shows that the probability of uniqueness of the maximum does not converge as N goes to infinity, but rather fluctuates around 0.72135 with logarithmic periodicity. It is only when N=2^n that the sequence converges to 0.721521… This probability actually writes down in closed form as

N\sum_{i=1}^\infty 2^{-i-1}(1-2^{-i})^{N-1}

(which is obvious in retrospect!, albeit containing a typo in the original paper which is missing a ½ factor in equation (17)) and its asymptotic behaviour is not obvious either, as noted by the authors.

On the historical side, and in accordance with Stiegler’s law, the Galton-Watson process should have been called the Bienaymé process! (Bienaymé was a student of Laplace, who successively lost positions for his political idea, before eventually joining Académie des Sciences, and later founding the Société Mathématique de France.)

another Galton-Watson riddle

Posted in Statistics with tags , , , on February 3, 2017 by xi'an

The riddle on the Riddler this week is definitely a classic, since it rephrases the standard Galton-Watson branching process (which should have been called Bienaymé‘s process, as he established the relation before Watson, while the jack-of-all-trades Francis Galton only posed the question):

At the beginning, 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?

As is easily seen from the moment generating function, the species goes instinct if p≤½. Actually, I always found it intriguing [intuitively] that the value ½ is included in the exclusion range!

a Galton-Watson riddle

Posted in R, Travel with tags , , , on December 30, 2016 by xi'an

riddleThe Riddler of this week has an extinction riddle which summarises as follows:

One observes a population of N individuals, each with a probability of 10⁻⁴ to kill the observer each day. From one day to the next, the population decreases by one individual with probability

K√N 10⁻⁴

What is the value of K that leaves the observer alive with probability ½?

Given the sequence of population sizes N,N¹,N²,…, the probability to remain alive is

(1-10^{-4})^{N+N^1+\ldots}

where the sum stops with the (sure) extinction of the population. Which is the moment generating function of the sum. At x=1-10⁻⁴. Hence the problem relates to a Galton-Watson extinction problem. However, given the nature of the extinction process I do not see a way to determine the distribution of the sum, except by simulation. Which returns K=27 for the specific value of N=9.

N=9
K=3*N
M=10^4
vals=rep(0,M)
targ=0
ite=1
while (abs(targ-.5)>.01){

for (t in 1:M){
  gen=vals[t]=N
  while (gen>0){
   gen=gen-(runif(1)<K*sqrt(gen)/10^4)
   vals[t]=vals[t]+gen}
  }
  targ=mean(exp(vals*log(.9999)))
print(c(as.integer(ite),K,targ))
  if (targ<.5){ K=K*ite/(1+ite)}else{
     K=K/(ite/(1+ite))}
  ite=ite+1}

The solution proposed on The Riddler is more elegant in that the fixed point equation is

\prod_{J=1}^9 \frac{K \cdot \sqrt{J}}{K \cdot \sqrt{J} + J}=\frac{1}{2}

with a solution around K=27.