2
$\begingroup$

I have a second order Markov chain with 4 states {A,T,C,G} (the 4 DNA nucleotides).

the transition matrix looks like this:

    A    T    C    G
AA[0.1, 0.6, 0.2, 0.1]
AT[0.3, 0.1, 0.5, 0.1]
AC[0.5, 0.3,  0,  0.2]
AG[..., ..., ..., ...]
TA[..., ..., ..., ...]
TT[..., ..., ..., ...]
TC[..., ..., ..., ...]
TG[..., ..., ..., ...]
CA[..., ..., ..., ...]
CT[..., ..., ..., ...]
CG[..., ..., ..., ...]
GA[..., ..., ..., ...]
GT[..., ..., ..., ...]
GC[..., ..., ..., ...]
GG[..., ..., ..., ...]

I wanted to calculate the stationary probability vector for the 4 states to which this matrix converges. The Markov chain is regular.

In case of first order Markov chains this is easily done by calculating the limit of $P^n$ with $n\rightarrow \infty$.

I do not know how to approach the problem in case of second order Markov chains.

Also, having a limited dataset from which to determine the transition matrix, can I consider the stationary distribution of the 4 nucleotides as being the theoretical distribution I would have if I had a much larger pool from which to draw (with the same transition matrix)?

In other words, can I consider the stationary distribution like an estimation of the theoretical nucleotide frequency given the transition matrix obtained from limited data?

$\endgroup$
4
$\begingroup$

A second order Markov chain is a random process $(X_n)_n$ on an alphabet $A$, whose distribution is specified by its transition probabilities $Q(x\mid y,z)=P[X_n=x\mid X_{n-1}=y,X_{n-2}=z]$, for every $(x,y,z)$ in $A\times A\times A$ (and by an initial distribution on $A\times A$). A stationary distribution of $(X_n)$ is a probability measure $\pi$ on $A\times A$ such that, if $\pi(x,y)=P[X_n=x,X_{n-1}=y]$ for every $(x,y)$ in $A\times A$ and some $n$, then $\pi(x,y)=P[X_{n+1}=x,X_{n}=y]$ for every $(x,y)$ in $A\times A$.

Thus, one asks that, for every $(x,y)$ in $A\times A$, $$ \pi(x,y)=\sum_{z\in A}Q(x\mid y,z)\pi(y,z). $$ As in the first order case, this linear system, together with the normalizing condition $$ \sum_{(x,y)\in A\times A}\pi(x,y)=1, $$ fully determines $\pi$ as soon as $(X_n)_n$ is irreducible. A new feature, absent of the first order case, is that every stationary distribution $\pi$ has identical marginals, that is, for every $x$ in $A$, $$ \varrho(x)=\sum_{y\in A}\pi(x,y)=\sum_{y\in A}\pi(y,x). $$ Finally, the MLE of $\pi$ based on $(X_k)_{0\leqslant k\leqslant n}$ is $\hat\pi_n$ defined by $$ \hat\pi_n(x,y)=\frac1n\sum_{k=1}^n\mathbf 1_{X_k=x,X_{k-1}=y}. $$ The MLE is consistent, that is, $\hat\pi_n(x,y)\to\pi(x,y)$ almost surely, for every $(x,y)$ in $A\times A$, when $n\to\infty$. In particular, the frequency of $x$ in $A$ stabilizes, since $$ \frac1n\sum_{k=1}^n\mathbf 1_{X_k=x}=\sum_{y\in A}\hat\pi_n(x,y)\to\varrho(x). $$

$\endgroup$
29
  • $\begingroup$ Would it be correct to think of the stationary distribution pi as the theoretical frequency distribution, of every couple (x,y) in AxA, in an infinitely long chain generated from the starting transition matrix? Also, would there be a way to extrapolate the frequencies of x in A instead of (x,y) in AxA? $\endgroup$ May 6 '13 at 13:29
  • $\begingroup$ Would it be correct to... Yes. // extrapolate the frequencies of x in A... See Edit. $\endgroup$
    – Did
    May 6 '13 at 15:04
  • $\begingroup$ i apologize but i am not a mathematician and i have trouble understanding the 2 last formulas, in particular i am not familiar with the rho notation. if you could elaborate a little bit it would be awesome. Also, do you think that extrapolating transition matrix from limited data and using it to calculate the limit frequencies of each state is more effective than directly using the frequencies calculated from the limited data? $\endgroup$ May 6 '13 at 15:24
  • 1
  • 1
    $\begingroup$ @Bernhard Of course they do not (this was not clear to me from your previous comment that this is what you were after) since every second order Markov chain $(X_n)$ is associated to a first order Markov chain $(Y_n)$ by defining $Y_n=(X_n,X_{n+1})$. Thus, take the result you are interested in in one of the canonical sources you cite, for $(Y_n)$, and trivially deduce its version for $(X_n)$. For example, $(Y_n)$ has a stationary distribution $\pi$, hence $(X_n)$ has a stationary distribution $\pi'$, which one gets by projecting $\pi$, that is, $$\pi'(x)=\sum_y\pi(x,y)=\sum_y\pi(y,x).$$ $\endgroup$
    – Did
    Aug 4 '16 at 9:13
-1
$\begingroup$

Since your second-order Markov chain is regular, you can still compute $\lim_{n\to \infty} Q^n$, where

$$ Q^n(x|y,z) = \mathbb{P}(X_{n+1}=x|X_1=y,X_0=z). $$

In the limit, you will have $\lim_{n\to \infty}Q^n(x|y,z)=\varrho(x)$ for every pair $(y,z)$. More information can be found in Chapter 7 of this book (especially equation 7.1.3).

Note however, that regularity or irreducibility of $Q$ do not imply that a unique stationary distribution $\pi(x,y)$ exists:

Consider a second-order Markov chain on $\{1,2,3,4\}$. Consider further, that there are two possible classes of cycles this Markov chain may go through: 1-2-3-4-1 and 1-2-3-1 (to break periodicity), or 1-4-3-2-1 and 1-3-2-1. From all pairs of states, the Markov chain moves to any of these two cycles and remains in them. One can show that this second-order Markov chain is regular and that $Q^n$ converges to a matrix with identical rows with entries $\varrho(x)$. But there is no unique invariant distribution $\pi(x,y)$, since there are two recurrent classes of tuples $(x,y)$ of states.

More generally, as mentioned in this book on page 173, if $x$ is a recurrent state for $(X_n)$, it need not be the case that a state $(x,y)$ for $(X_n,X_{n+1})$ is recurrent.

$\endgroup$
3
  • $\begingroup$ The transition Q given in the example is obviously reducible, two communication classes being {12,23,34,41,31} and {14,43,32,21,13} (the status of states 24 and 42 being unclear). $\endgroup$
    – Did
    Oct 20 '16 at 8:49
  • $\begingroup$ A second-order Markov chain is irreducible if, for every $z,y,x$ there is a $n=n(z,y,x)$ such that $Q^n(x|y,z)>0$. Irreducibility is defined on the states of $X$, i.e., on $A$, and not on tuples of states (on $A\times A$). The example I give is clearly irreducible. $\endgroup$
    – Bernhard
    Oct 20 '16 at 9:15
  • $\begingroup$ See comment on my answer. (Once again, and for the last time, next time you want to explain your own theories à propos a question, please do so by first posting an answer, not by polluting the comment thread of another answer. Thanks in advance.) $\endgroup$
    – Did
    Oct 20 '16 at 10:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.