Computing a stationary distribution of a Markov chain¶
Let stochastic matrix $ P $ denote the transition probability matrix of a Markov chain which takes values from $ S=\{0,\dots,n-1\} $
Assume that $ P $ is aperiodic and irreducible, so there exists a unique stationary distribution $ \psi^\star $ such that
$$ \psi^\star = \psi^\star \cdot P $$Our task is to compute $ \psi^\star $
Method 1: successive approximations¶
Due to convergence results, we can use the following algorithm (see previous video on Markov chains)
- Start from arbitrary distribution $ \psi_0 $
- Compute the updated distribution $ \psi_t = \psi_{t-1} P $ until $ \psi_t $ and $ \psi_{t-1} $ are indistinguishable
In [1]:
import numpy as np
P = np.array([[.5,.3,.2],[.5,.4,.1],[.1,.1,.8]])
ψ = np.array([1,0,0])
In [2]:
def stationary_sa(P,psi0,tol=1e-6,maxiter=100,callback=None):
'''Computes stationary distribution for the Markov chain given by transition probability
matrix P, with given maximum number of iterations, and convergence tolerance.
callback function is called at each iteration if given.
Method: successive approximations
'''
pass
stationary_sa(P,ψ)
In [3]:
def stationary_sa(P,psi0=[None,],tol=1e-6,maxiter=100,callback=None):
'''Computes stationary distribution for the Markov chain given by transition probability
matrix P, with given maximum number of iterations, and convergence tolerance.
callback function is called at each iteration if given.
Method: successive approximations
'''
if psi0[0]==None:
# degenrate initial distribution
psi0 = [0,]*P.shape[0]
psi0[0]=1.0
P,psi0 = np.asarray(P),np.asarray(psi0) # convert to np arrays (in case lists were passed)
assert np.all(np.abs(P.sum(axis=1)-1)<1e-10), 'Passed P is not stochastic matrix'
assert np.abs(psi0.sum()-1)<1e-10, 'Passed probabilities do not sum up to one'
for i in range(maxiter): # main loop
psi1 = psi0 @ P # update approximation of psi^star
err = np.amax(abs(psi0-psi1)) # error is the max absolute deviation element-wise
if callback != None: callback(err=err,iter=i,psi=psi1)
if err<tol:
break # break out if converged
psi0 = psi1 # get ready to the next iteration
else:
raise RuntimeError('Failed to converge in %d iterations'%maxiter)
return psi1
def printme(**kwargs):
print('iter %d, psi = %r'%(kwargs['iter'],kwargs['psi']))
stationary_sa(P,callback=printme)
iter 0, psi = array([0.5, 0.3, 0.2]) iter 1, psi = array([0.42, 0.29, 0.29]) iter 2, psi = array([0.384, 0.271, 0.345]) iter 3, psi = array([0.362 , 0.2581, 0.3799]) iter 4, psi = array([0.34804, 0.24983, 0.40213]) iter 5, psi = array([0.339148, 0.244557, 0.416295]) iter 6, psi = array([0.333482 , 0.2411967, 0.4253213]) iter 7, psi = array([0.32987148, 0.23905541, 0.43107311]) iter 8, psi = array([0.32757076, 0.23769092, 0.43473833]) iter 9, psi = array([0.32610467, 0.23682143, 0.4370739 ]) iter 10, psi = array([0.32517044, 0.23626736, 0.4385622 ]) iter 11, psi = array([0.32457512, 0.2359143 , 0.43951058]) iter 12, psi = array([0.32419577, 0.23568931, 0.44011492]) iter 13, psi = array([0.32395403, 0.23554595, 0.44050002]) iter 14, psi = array([0.32379999, 0.23545459, 0.44074542]) iter 15, psi = array([0.32370183, 0.23539638, 0.44090179]) iter 16, psi = array([0.32363928, 0.23535928, 0.44100144]) iter 17, psi = array([0.32359943, 0.23533564, 0.44106493]) iter 18, psi = array([0.32357403, 0.23532058, 0.4411054 ]) iter 19, psi = array([0.32355784, 0.23531098, 0.44113118]) iter 20, psi = array([0.32354753, 0.23530486, 0.44114761]) iter 21, psi = array([0.32354096, 0.23530096, 0.44115808]) iter 22, psi = array([0.32353677, 0.23529848, 0.44116475]) iter 23, psi = array([0.3235341, 0.2352969, 0.441169 ]) iter 24, psi = array([0.3235324 , 0.23529589, 0.44117171]) iter 25, psi = array([0.32353132, 0.23529525, 0.44117344]) iter 26, psi = array([0.32353062, 0.23529484, 0.44117454]) iter 27, psi = array([0.32353018, 0.23529458, 0.44117524])
Out[3]:
array([0.32353018, 0.23529458, 0.44117524])
Method 2: direct solution with linear solver¶
$ \psi^\star \in \mathbb{R}^n $ is a row vector of probabilities, $ \sum \psi^\star_i = 1 $
$$ \psi^\star = \psi^\star \cdot P $$Write a linear system of equations:
$$ \psi^\star (I - P) = 0, \quad \psi^\star $$$$ (I - P') \tilde{\psi}^\star = 0, \quad \tilde{\psi}^\star \text{ is a column vector!} $$We can use the second form with the standard np.linalg.solver
In [4]:
# have to be careful with the zero solution!
m = P.shape[0]
A = np.eye(m) - P.T
print(np.linalg.solve(A,np.zeros(m)))
print(np.linalg.matrix_rank(A))
[ 0. -0. -0.] 2
Include the linear constraint into the system itself¶
Let $ e $ be the column vector of appropriate length with all elements equal 1
The system
$$ \left( \begin{array}{cc} I - P' & e \\ e' & 1 \end{array} \right) \cdot \left( \begin{array}{c} \psi \\ \epsilon \end{array} \right) = \left( \begin{array}{c} e \\ 2 \end{array} \right) $$has full rank and therefore can be solved directly
In [5]:
def stationary_linalg(P,psi0=None):
'''Computes stationary distribution for the Markov chain given by transition probability
matrix P. Method: linear solver.
'''
pass
In [6]:
def stationary_linalg(P,psi0=[None,]):
'''Computes stationary distribution for the Markov chain given by transition probability
matrix P. Method: linear solver.
'''
if psi0[0]==None:
# degenrate initial distribution
psi0 = [0,]*P.shape[0]
psi0[0]=1.0
P,psi0 = np.asarray(P),np.asarray(psi0) # convert to np arrays (in case lists were passed)
assert np.all(np.abs(P.sum(axis=1)-1)<1e-10), 'Passed P is not stochastic matrix'
assert np.abs(psi0.sum()-1)<1e-10, 'Passed probabilities do not sum up to one'
m = P.shape[0] # dimension of the problem
A = np.ones((m+1,m+1)) # square matrix of ones
A[:-1,:-1] = np.eye(m)-P.T
b = np.ones(m+1)
b[-1]=2
res = np.linalg.solve(A,b)
return res[:-1]
print(stationary_linalg(P))
[0.32352941 0.23529412 0.44117647]
Further learning resources¶
- Post on Brilliant.org https://brilliant.org/wiki/stationary-distributions/