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.

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

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.

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

M

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.

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

M

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.

• Start with a random state.

• Loop

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

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

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

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.

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

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\).

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

S.real

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

U.real

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]])

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

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.

from scipy.linalg import eig

S, U = eig(M.T)

S.real

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

U.real

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]])

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

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