# Markov Chain for combined events?

I was wondering how to apply Markov steady-state chains to a simple case in genetics, eye colour.

Assuming the probabilities of being brown, green or blue eyed, given your parents were a certain colour combination are:

Br + Br - Br 75%  , Gr 18.75% , Bl 6.25%
Gr + Br - Br 50%  , Gr 37.5%  , Bl 12.5%
Bl + Br - Br 50%  , Gr 0%     , Bl 50%
Gr + Gr - Br 0.5% , Gr 74.75% , Bl 24.75%
Gr + Bl - Br 0%   , Gr 50%    , Bl 50%
Bl + Bl - Br 0%   , Gr 1%     , Bl 99%


Corresponding matrix:

{{3/4, 1/2, 1/2, 0.005, 0, 0}, {0.1875, 0.375, 0, 0.7475, 1/2, 1/
100}, {0.0625, 0.125, 1/2, 0.2475, 1/2, 99/100}}


I want to find the steady state for the eye colour of the population.

• If you've made any attempt at solving this yourself it would be a good idea to share that attempt. Providing data that's been formatted for Mathematica is also appreciated. Aug 28, 2013 at 1:43
• This question requires an initial gentype frequency. some assumption or declaration of probabilities of genotypes mating and then use of Mendelian genetic rules to develop transition matrix for genotypes (as states). The eye colors can be inferred from the genotype-phenotype relationships after determining stationary(steady state) distribution for given initial state. The matrix you have, though rows add to 1, is not a transition matrix, eg entry (1,1) is not transition from state 1 to state 1 etc. Aug 28, 2013 at 4:20
• I haven't made any attempts that led anywhere. I saw the table of values on twitter and immediately wondered what the steady-state would be. I was willing to assume that any one person can marry any other person on Earth to remove any regional differences there may be, but I was hoping that it was possible to somehow use the steady-state as the probability that person A is Br/Gr/Bl and same for person B and solve? I haven't studied Markov chains in that much detail and it seems like it is a lot more complicated than I had hoped. Aug 28, 2013 at 8:14
• I thought about setting: P(Brown)=0.75*P(Brown)*P(Brown)*... etc. It would be assuming that there is equal likelihood of a person with brown eyes having a child with a person with green, blue or brown eyes (so there is no bias from a particular eye colour to a particular eye colour) and that eye colour is spread the same universally. It would seem to lead to a solution but I don't have a maths package and there are too many terms for Wolfram Alpha. Sep 9, 2013 at 17:03

@YayaManisya's answer can be acquired by numerical fixed point calculation:

FixedPoint[
Apply[# / Total@# & @*
Function[{br, gr, bl},
{br br, 2 gr br, 2 bl br, gr gr, 2 gr bl, bl bl}.
{{.75, .1875, .0625}, {.5, .375, .125}, {.5, 0, .5},
{.005, .7475, .2475}, {0, .5, .5}, {0, .01, .99}}]],
{1, 1, 1} / 3]


{0.0218527, 0.154522, 0.823625}

Here we repeatedly apply a function which computes a vector of pairwise probabilities for parents, and multiply that with the corresponding offspring probability matrix, and in the end normalizes the sum of probabilities to 1, mostly to avoid numerical errors destroying the result. FixedPoint nests these calls until the result doesn't change (much).

This result doesn't seem to vary in this case depending on the starting population.

Change in population eye color can be seen here for one starting population:

NestList[
Apply[# / Total@# & @*
Function[{br, gr, bl},
{br br, 2 gr br, 2 bl br, gr gr, 2 gr bl, bl bl}.
{{.75, .1875, .0625}, {.5, .375, .125}, {.5, 0, .5},
{.005, .7475, .2475}, {0, .5, .5}, {0, .01, .99}}]],
{0.9, 0.075, 0.025}, 100] // Transpose //
ListLinePlot[#, PlotStyle -> {Brown, Green, Blue}] &


Also, the solution can be found symbolically (especially now when we are convinced there's just one steady state):

Solve[0 <= br <= 1 && 0 <= gr <= 1 && 0 <= bl <= 1 &&
br + gr + bl == 1 &&
{br, gr, bl} == {br br, 2 gr br, 2 bl br, gr gr, 2 gr bl, bl bl}.
Rationalize@{{.75, .1875, .0625}, {.5, .375, .125}, {.5, 0, .5},
{.005, .7475, .2475}, {0, .5, .5}, {0, .01, .99}},
{br, gr, bl}] // FullSimplify // ToRadicals


$$\text{br}\to \frac{2 \left(-342-46790 \sqrt{2}+5 \sqrt{4438369759-409761860 \sqrt{2}}\right)}{22339241},\\\text{ gr}\to \frac{10 \left(-93580-342 \sqrt{2}+5 \sqrt{2 \left(4438369759-409761860 \sqrt{2}\right)}\right)}{22339241} ,\\\text{bl}\to \frac{46551450+194000 \sqrt{2}-22339241 \sqrt{\frac{87264648203600}{499041 688456081}+\frac{9394337092000 \sqrt{2}}{499041688456081}}}{44678 482}$$

These values correspond to the numerical values achieved with FixedPoint.

Here's an easy and probably useless implementation. Imagine nine states in a Markov chain. Each state is an ordered pair, drawn from $(Br,\ Gr,\ Bl)$ that describes the eye color of a person and their mate. If you assume that a person chooses a mate with one of the three eye colors with equal probability, then I think the following MarkovProcess tracks a single path through an infinite family tree (or the 100 generations I simulate below):

proc = DiscreteMarkovProcess[Table[1/9, {9}],
Flatten@Table[Table[p, {3}], {p, #}]/3 & /@ {
{0.75, 0.1875, 0.0625},
{0.5, 0.375, 0.125},
{0.5, 0, 0.5},
{0.5, 0.375, 0.125},
{0.005, 0.7475, 0.2475},
{0, 0.5, 0.5},
{0.5, 0, 0.5},
{0, 0.5, 0.5},
{0, 0.01, 0.99}}];
ListPlot[RandomFunction[proc, {0, 100}],
Filling -> Axis, TicksStyle -> 16, ImageSize -> 500,
Ticks -> {Automatic, Flatten[Table[{3 (i - 1) + j,
StringJoin[{"(", {"Br", "Gr", "Bl"}[[i]], ",", {"Br", "Gr", "Bl"}[[j]], ")"}]},
{i, 3}, {j, 3}], 1]}
]


Hardly useful! But maybe interesting as a simple example of a Markov Chain simulation? Of course, there's no steady state here....

This is a simulation and perhaps will motivate better answers or can be adapted. In the following:

• the starting population is (1/3,1/3,1/3) blue:brown:green

• each generation produces one child from one coupling with probability of color as per OP

• no deaths

Again this is just to a toy example to illustrate evolution over time of a population with the above assumptions. The gif is a downsampling of vis/@sim[100,1000]:

mat = Transpose[{{3/4, 1/2, 1/2, 0.005, 0, 0}, {0.1875, 0.375, 0,
0.7475, 1/2, 1/100}, {0.0625, 0.125, 1/2, 0.2475, 1/2, 99/100}}];
mnpar = {{"Blue", "Blue"} -> mat[[6]], {"Blue", "Brown"} ->
mat[[3]], {"Blue", "Green"} -> mat[[5]], {"Brown", "Brown"} ->
mat[[1]], {"Brown", "Green"} -> mat[[2]], {"Green", "Green"} ->
mat[[4]]};
fun[1, 0, 0] := "Blue";
fun[0, 1, 0] := "Brown";
fun[0, 0, 1] := "Green";
sim[n_, num_] :=
Module[{rv =
RandomVariate[MultinomialDistribution[1, {1/3, 1/3, 1/3}], n]},
Map[fun,
NestList[{RandomVariate[
MultinomialDistribution[
1, (Sort[fun @@@ RandomSample[#, 2]] /. mnpar)]]}~Join~# &,
rv, num], {2}]
]
vis[lst_] := Module[{col = {"Blue", "Brown", "Green"}, tly},
tly = {#, Count[lst, #]} & /@ col;
BarChart[#2/Total[tly[[All, 2]]] & @@@ tly, ChartLabels -> col,
PlotLabel -> Row[{"Population:", Total[tly[[All, 2]]]}],
LabelingFunction -> Above, PlotRange -> {0, 0.5}, Frame -> True,
ChartStyle -> {Blue, Brown, Green}]]


Actually you just need to assign unknowns to the steady state probabilities and sum them up accordingly.

Let $P(W)$, $P(G)$ and $P(B)$ represent the probabilities of being brown, green and blue-eyed, and thus

$P(W) + P(G) + P(B) = 1$ ,

assuming that you get no other eye colours in children due to illness etc.

Then you will get

$\begin{eqnarray}P(W) = &&0.75[P(W)]^2 + [P(W)][P(G)] \\ &&+[P(W)][P(B)] + 0.005[P(G)]^2\end{eqnarray}$

$\begin{eqnarray}P(G) = &&0.1875[P(W)]^2 + 0.75[P(W)][P(G)] + 0.7475[P(G)]^2 \\ &&+ [P(G)][P(B)] + 0.01[P(B)]^2\end{eqnarray}$

$\begin{eqnarray}P(B) = &&0.0625[P(W)]^2 + 0.25[P(W)][P(G)] + 0.2475[P(G)]^2\\ &&+ [P(W)][P(B)] + [P(G)][P(B)] + 0.99[P(B)]^2\end{eqnarray}$

From here on you can use normal calculation method to find the answer. Simplify the equations first, and after that you can use any calculator that enable square root calculation to find the final numerical answers.

I get

$\begin{eqnarray} P(W) &&= 1 / (1 + 2.5 * \sqrt{199 + 20 * \sqrt{2}} + 5 * \sqrt{2})\\ &&= 0.0218527 \\ &&= 2.2\% \text{ of the population,}\end{eqnarray}$

$\begin{eqnarray}P(G) &&= P(W) * 5 * \sqrt{2}\\ &&= 0.1545219 \\ &&= 15.5\% \text{ of the population,}\end{eqnarray}$

$\begin{eqnarray}P(B) &&= 1 - P(W) - P(G) \\ &&= 82.3\% \text{ of the population.}\end{eqnarray}$

• This site is about the Mathematica software. Can you complement you answer with valid MMA code to carry out the calculations? Sep 13, 2018 at 17:42
• Welcome to Mathematica.SE! I have edited your answer using MathJax to make it more readable. Please check, if I introduced any errors.
– Johu
Sep 13, 2018 at 17:55