11. Markov Chain, Stationary Distribution

When we have a matrix that represents transition probabilities or a Markov chain, it is often of interest to find the marginal probabilities of the states. Such probabilities are called the stationary distribution which represents the probabilities of the states in the long run. There are many ways to compute the stationary distribution:

• sampling,

• Power Method or

• matrix decomposition.

The following matrix represents 5 states. We will normalize them so that the matrix represents probabilities and each row sums to 1.

[1]:

import numpy as np

M = np.array([
[0, 1, 1, 1, 1],
[0, 0, 0, 0, 0],
[0, 1, 0, 0, 1],
[0, 1, 1, 0, 1],
[0, 1, 0, 0, 0]
])
M

[1]:

array([[0, 1, 1, 1, 1],
[0, 0, 0, 0, 0],
[0, 1, 0, 0, 1],
[0, 1, 1, 0, 1],
[0, 1, 0, 0, 0]])


Here we detect any row with all zeroes and set all the elements to 1.

[2]:

for r, s in enumerate(M.sum(axis=1)):
if s == 0:
M[r, :] = 1

[3]:

M

[3]:

array([[0, 1, 1, 1, 1],
[1, 1, 1, 1, 1],
[0, 1, 0, 0, 1],
[0, 1, 1, 0, 1],
[0, 1, 0, 0, 0]])


Finally, we normalize the rows to probabilities in such a way that they sum to 1.

[4]:

M = np.array([M[r,:] * s for r, s in enumerate(1 / M.sum(axis=1))])

[5]:

M

[5]:

array([[0.        , 0.25      , 0.25      , 0.25      , 0.25      ],
[0.2       , 0.2       , 0.2       , 0.2       , 0.2       ],
[0.        , 0.5       , 0.        , 0.        , 0.5       ],
[0.        , 0.33333333, 0.33333333, 0.        , 0.33333333],
[0.        , 1.        , 0.        , 0.        , 0.        ]])


11.1. Sampling

Sampling is easy.

• Loop

• Increment the state (signal that we have seen the state)

• Sample uniformly $$p$$ from $$[0, 1)$$ and pick the next state

[6]:

from random import choice
import bisect

np.random.seed(37)

indexes = {r: [c for c in range(M.shape[1]) if M[r][c] > 0] for r in range(M.shape[0])}
cumsum = {r: M[r][c].cumsum() for r, c in indexes.items()}

c = np.zeros(M.shape[0]) + 1e-100
r = choice(range(M.shape[0]))

for it in range(30_000):
prev = c / c.sum()
c[r] += 1
curr = c / c.sum()
d = np.linalg.norm(prev - curr, 1)

if d < 0.001:
print(f'num of iterations: {it}')
break

p = np.random.random()
i = bisect.bisect_left(cumsum[r], p)
r = indexes[r][i]

c = c / c.sum()
c

num of iterations: 1151

[6]:

array([0.09548611, 0.42534722, 0.14670139, 0.11024306, 0.22222222])


11.2. Power Method

Using power iteration, we can also find the stationary distribution.

[7]:

X = M.dot(M)
for it in range(10):
Y = X.dot(X)

x = np.diag(X)
y = np.diag(Y)
d = np.linalg.norm(x - y, 1)

if d < 0.001:
print(f'num of iterations: {it}')
break

X = Y

np.diag(X)

num of iterations: 3

[7]:

array([0.08759149, 0.43795601, 0.1459856 , 0.10948932, 0.2189781 ])


11.3. Numpy, eig

We can also use eigen decomposition from Numpy. The stationary probability will be the normalized eigenvector associated with the eigenvalue closest to 1. Note that we have to transpose $$M$$.

[8]:

S, U = np.linalg.eig(M.T)

[9]:

S.real

[9]:

array([ 1.        , -0.19142147, -0.19142147, -0.20857853, -0.20857853])

[10]:

U.real

[10]:

array([[-0.16531686, -0.16178269, -0.16178269, -0.6152745 , -0.6152745 ],
[-0.8265843 ,  0.7453579 ,  0.7453579 ,  0.6968425 ,  0.6968425 ],
[-0.2755281 , -0.2841635 , -0.2841635 ,  0.09145666,  0.09145666],
[-0.20664607, -0.28528506, -0.28528506,  0.0053999 ,  0.0053999 ],
[-0.41329215, -0.01412664, -0.01412664, -0.17842455, -0.17842455]])

[11]:

(U[:,np.isclose(S, 1)][:,0] / U[:,np.isclose(S, 1)][:,0].sum()).real

[11]:

array([0.08759124, 0.4379562 , 0.1459854 , 0.10948905, 0.2189781 ])


11.4. Scipy, eig

Eigen decomposition from Scipy works the same way as Numpy.

[12]:

from scipy.linalg import eig

S, U = eig(M.T)

[13]:

S.real

[13]:

array([ 1.        , -0.19142147, -0.19142147, -0.20857853, -0.20857853])

[14]:

U.real

[14]:

array([[-0.16531686, -0.16178269, -0.16178269, -0.6152745 , -0.6152745 ],
[-0.8265843 ,  0.7453579 ,  0.7453579 ,  0.6968425 ,  0.6968425 ],
[-0.2755281 , -0.2841635 , -0.2841635 ,  0.09145666,  0.09145666],
[-0.20664607, -0.28528506, -0.28528506,  0.0053999 ,  0.0053999 ],
[-0.41329215, -0.01412664, -0.01412664, -0.17842455, -0.17842455]])

[15]:

(U[:,np.isclose(S, 1)][:,0] / U[:,np.isclose(S, 1)][:,0].sum()).real

[15]:

array([0.08759124, 0.4379562 , 0.1459854 , 0.10948905, 0.2189781 ])