Python code explanation for stationary distribution of a Markov chain

Question

I have got this code:

import numpy as np
from scipy.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]

But I can't understand the line:

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

Can anyone explain the part: U[:, np.where(np.abs(S - 1.) < 1e-8)[0][0]].flat ?

I know that the routine returns S: eigenvalue, U : eigenvector. I need to find the eigenvector corresponding to the eigenvalue 1. I have wrote the code below:

for i in range(len(S)):
    if S[i] == 1.0:
        j = i

 matrix = np.array(U[:, j].flat)

I am getting output:

:  [ 0.6144763   0.57607153  0.53766676  0.03840477]

but it does not give the same output. why?!

Answer 1

How to find a stationary distribution.

Ok, I came to this post looking to see if there was a built-in method to find the stationary distribution. It looks like there's not. So, for anyone coming in from Google, this is how I would find the stationary distribution in this circumstance:

import numpy as np

#note: the matrix is row stochastic.
#A markov chain transition will correspond to left multiplying by a row vector.
Q = np.array([
    [.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]])

#We have to transpose so that Markov transitions correspond to right multiplying by a column vector.  np.linalg.eig finds right eigenvectors.
evals, evecs = np.linalg.eig(Q.T)
evec1 = evecs[:,np.isclose(evals, 1)]

#Since np.isclose will return an array, we've indexed with an array
#so we still have our 2nd axis.  Get rid of it, since it's only size 1.
evec1 = evec1[:,0]

stationary = evec1 / evec1.sum()

#eigs finds complex eigenvalues and eigenvectors, so you'll want the real part.
stationary = stationary.real

What that one weird line is doing.

Let's break that line into parts:

#Find the eigenvalues that are really close to 1.
eval_close_to_1 = np.abs(S-1.) < 1e-8

#Find the indices of the eigenvalues that are close to 1.
indices = np.where(eval_close_to_1)

#np.where acts weirdly.  In this case it returns a 1-tuple with an array of size 1 in it.
the_array = indices[0]
index = the_array[0]

#Now we have the index of the eigenvector with eigenvalue 1.
stationary = U[:, index]

#For some really weird reason, the person that wrote the code
#also does this step, which is completely redundant.
#It just flattens the array, but the array is already 1-d.
stationary = np.array(stationary.flat)

If you compress all these lines of code into one line you get stationary = np.array(U[:, np.where(np.abs(S-1.)<1e-8)[0][0]].flat)

If you remove the redundant stuff you get stationary = U[:, np.where(np.abs(S - 1.) < 1e-8)[0][0]]

Why your code gives a different stationary vector.

As @Forzaa pointed out, your vector cannot represent a vector of probabilities because it does not sum to 1. If you divide it by its sum, you'll get the vector the original code snippet has.

Just add this line:

stationary = matrix/matrix.sum()

Your stationary distribution will then match.