I am totally stuck with this exercise on Markov chains. Maybe someone can help me :).
Red and green bacteria
A growth medium at time $t = 0$ has 500 red bacteria and 500 green bacteria. Each hour, each bacterium divides in two. A colour-blind predator eats exactly 1000 bacteria per hour.
- After a very long time, what is the probability distribution for the number $\alpha$ of red bacteria in the growth medium?
- Roughly how long will it take to reach this final state?
- Assume that the predator has a $1\%$ preference for green bacteria (implemented as you choose). Roughly how much will this change the final distribution?
Within the accuracy of this question, you may assume either that one bacterium reproduces and then one is eaten 1000 times per hour, or that at the end of each hour all the bacteria reproduce and then 1000 are consumed.
I want to get analytical solution. So, my attempt on former method,
where steps are hours is based on binomial distribution. I think that
probability for predator to consume number $n_r$ of red bacteria is:
$$
P_{\text{bin}} (n_r) = \frac{N!}{(N-n_r)!n_r!} (1-p)^{N-n_r} p^{n_r}
$$
where $N=1000$ and $p$ is probability for predator to eat red bacteria
which depends on current state.
If I say that current number of red bacteria is $\alpha$,
I can write
probability $p=\alpha/N$ and transition probability:
$$
P(\alpha+n|\alpha) = \frac{N!}{(N-\alpha+n)!\alpha-n!}(1-p)^{N-\alpha+n} p^{\alpha - n}
$$
However when I try to simulate this with simple python program, where
each step each bacteria duplicates itself and afterwards predator
randomly eats half of the population I get different results.
For example from simulation I get $P(501|500)=0.035$ and from analytical
formula $P(501|500)=0.025$. Maybe someone can explain me why I get this difference
and why binomial distribution is wrong choice.
After some benchmarking and inconsistencies I decided to switch to later method, where one bacterium reproduces and then one is eaten 1000 times per hour. At this case there are possible only three transitions from each state: \begin{eqnarray} P(\alpha + 1|\alpha) &=& P(\text{predator eats green})P(\text{red reproduces}) = \frac{N - \alpha}{N} \frac{\alpha}{N}\\ P(\alpha|\alpha) &=& P(\text{predator eats green})P(\text{green reproduces}) \\ && + P(\text{predator eats red})P(\text{red reproduces}) \\ &=& \frac{(N - \alpha)^2}{N^2} + \frac{\alpha^2}{N^2}\\ P(\alpha -1|\alpha) &=& P(\text{predator eats red})P(\text{green reproduces}) = \frac{\alpha}{N}\frac{N - \alpha}{N}\\ \end{eqnarray}
However further I don't know what to do next. Equilibrium distribution (state) is eigenvector $\vec\rho$ of transition matrix $$ P_{ij} = P(i|j) $$ with eigenvalue 1. How I can find this vector. Writing eigenvalue problem in index notation gives me: $$ \rho_\alpha = P(\alpha|\alpha-1)\rho_{\alpha - 1} + P(\alpha|\alpha)\rho_{\alpha} + P(\alpha|\alpha+1)\rho_{\alpha+1} $$ Or maybe I should use detailed balance condition: $$ P(\alpha|\beta)\rho_{\beta} = P(\beta|\alpha)\rho_{\alpha} $$ Sorry for my english :).
The idea here is to realize that eventually either all reds disappear, or all greens, whichever comes first. So the limiting distribution is extremal.
Assume the subset of 1000 bacteria the predator eats every time slot is independent and uniformly distributed over all options. So, after each doubling step, each one of the 2000 bacteria has a probability of $1/2$ of being eaten.
One interesting aspect of this is that the expected number of reds is always the same, regardless of the number of reds we start with. To see this, let $N(t)$ be the number of reds (for time slots $t \in \{0, 1, 2, \ldots\}$), and assume $N(0)=n_0 \in \{0, 1, \ldots, 1000\}$ is a given constant. Then for each step $t\geq 0$ we get: $$ E[N(t+1)|N(t)] = 2N(t)/2 = N(t) $$ So $E[N(t+1)]=E[N(t)]=n_0$ for all $t \in \{0, 1, 2, \ldots\}$.
This can be used to derive a simple expression for the probability that red eventually dies out, given a particular starting condition $n_0$.