## Wednesday, September 5, 2012

### The Stationary Distribution: A third way

Ian Carroll helpfully reminded me that I had neglected the easiest (and most computationally efficient) way of calculating the stationary distribution of the Markov model in the previous post, and I would be remiss if I didn't demonstrate it and the logic behind its use. This time we take advantage of the definition of an eigenvector. Recall that our goal is to find the probability distribution of states that, when left-multiplied by the transition probability matrix, results in the same distribution. In math-speak, this is the following: $${\bf v^TP = v^T}$$

Now, with that in mind, recall (or learn for the first time if you're one of today's lucky 10,000) that the eigenvector ${\bf u}$ and corresponding eigenvalue $\lambda$ for some square matrix ${\bf A}$ are defined as follows: $${\bf Au} = \lambda {\bf u}$$

More specifically, these are called the right eigenvectors/eigenvalues, since they are defined when the vector is multiplied on the right. We can also define the left eigenvectors/eigenvalues as the ${\bf u}$, $\lambda$ combos that satisfy the following: $${\bf u^TA} = \lambda {\bf u^T}$$ Now, if we plug in our stationary distribution vector ${\bf v}$ for ${\bf u}$, our transition matrix ${\bf P}$ for ${\bf A}$, and let $\lambda = 1$, we get the familiar ${\bf v^TP = v^T}$. Also, finding the left eigenvalues and eigenvectors for a matrix is the same as finding the right-side versions for the transpose of that matrix. This means that, computationally, we can take the eigen decomposition of ${\bf P^T}$, find the column in the ${\bf U}$ matrix corresponding to an eigenvalue of 1, normalize that vector (since it can be scaled by any arbitrary constant), and that will give the stationary distribution. I demonstrate this below:

import numpy as np
from numpy.linalg import eig

transition_mat = np.matrix([
[.95,.05,0.,0.],\
[0,0.9,0.09,0.01],\
[0,0.05,0.9,0.05],\
[0.8,0,0.05,0.15]])

S,U = eig(transition_mat.T)
stationary = np.array(U[:,np.where(np.abs(S-1.) < 1e-8)[0][0]].flat)
stationary = stationary / np.sum(stationary)

>>> print stationary
[ 0.34782609  0.32608696  0.30434783  0.02173913]


We get the same stationary distribution as before without having to perform and extra exponentiation or matrix multiplication! The glories of the eigen decomposition never cease to amaze!

#### 1 comment:

1. Oh thank you so much, very useful !