Computing the population given a set of replication rates

So I'm looking to calculate the probability that a bacterial population dies out with the following conditions:

• Initial population is 1
• At each iteration the bacteria can die, do nothing, divide into 2 or divide into three. So each bacteria can change into either 0, 1, 2, 3 bacteria.
• Each event has a 25% likelihood.

So I know from simulation that the likelihood of the population dying out is ~36.3%. However, I'm looking for a method to calculate this analytically.

Some MATLAB code for the simulation:

ntries = 100000;
mnum = 10000;
res = ones(ntries, 1);

for i = 1:ntries
while res(i) > 0 && res(i) < mnum
res(i) = sum(randi(4,res(i), 1)-1);
end
end
mean(res==0)


Thanks,

Will

-
Are you sure it's not 41.4%? – joriki Apr 28 '11 at 20:24
Nope I reliably get 0.3632 although this assumes that after a population gets to 10,000 then it won't die out ... so that might be where I'm loosing 5%. – JudoWill Apr 28 '11 at 21:25
@JudoWill: Three people have independently arrived at 41.4%, so I daresay it's likely that the mistake is on your side :-) – joriki Apr 28 '11 at 21:28
@JudoWill: I've never used MATLAB, but I think the "res(i)" after "sum" shouldn't be there -- you want to sum up the random values; there's no reason to multiply them with the current number of bacteria. – joriki Apr 28 '11 at 21:36
yup, my bad ... I get ~41% now – JudoWill Apr 28 '11 at 21:37

Call $q_n$ the probability that the initial bacteria dies out at or before generation $n$ and condition by the size $i$ of the first generation. The event that every bacteria of this generation dies out at or before generation $n+1$ has probability $q_n^i$ and $i=0$, $1$, $2$ or $3$ with probability $1/4$, hence $q_{n+1}=g(q_n)$ with $$g(q)=q^0/4+q^1/4+q^2/4+q^3/4.$$ Since $q_0=0$, $q_n=g\circ g\circ\cdots\circ g(0)$ ($g$ composed $n$ times). One can check, for instance by considering the graph of $g$ on the interval $(0,1)$, that $q_n$ converges to the smallest nonnegative root $q$ of the equation $q=g(q)$.

Now $g(q)=q$ if and only if $0=q^3+q^2-3q+1=(q-1)(q^2+2q-1)$, hence $q=1$ or $(q+1)^2=2$, that is $q=-1\pm\sqrt{2}$. Finally the population dies out with probability $q=\sqrt2-1\approx41\%$.

See this chapter of lecture notes by Steve Lalley, section 1.2.

-
Either .41 or 41% :-) (Also $\approx$, not $=$) – joriki Apr 28 '11 at 20:26
@joriki: Yep. :-) Thanks. – Did Apr 28 '11 at 20:29

I find the probability to be $\sqrt{2}-1\approx 41.4\%$. If we let $p(n)$ be the probability that a population of $n$ dies out, we have $p(0)=1, p(n)=(p(n-1)+p(n)+p(n+1)+p(n+2))/4$. I first just did a numeric relaxation in Excel and recognized $p(1), p(2)=3-\sqrt{2}, p(3)=5\sqrt{2}-7.$ This was aided by the inverse symbolic calculator. The recurrence becomes $p(n+2)=-p(n+1)+3p(n)-p(n-1)$, whose characteristic polynomial has roots $1$, $\sqrt{2}-1$, and $-1-\sqrt{2}$. Finally we can check that $\sqrt{2}-1=(1+\sqrt{2}-1+3-2\sqrt{2}+5\sqrt{2}-7)/4$

-