from __future__ import division
import numpy as np
from mc_tools import MarkovChain
# Taken from www.math.wustl.edu/~feres/Math450Lect04.pdf
P = np.zeros((10, 10))
P[[0, 0], [0, 2]] = 1./2
P[1, 1], P[1, 6] = 1./3, 2./3
P[2, 0] = 1
P[3, 4] = 1
P[[4, 4, 4], [3, 4, 8]] = 1./3
P[5, 5] = 1
P[6, 6], P[6, 8] = 1./4, 3./4
P[[7, 7, 7, 7], [2, 3, 7, 9]] = 1./4
P[8, 1] = 1
P[[9, 9, 9], [1, 4, 9]] = 1./3
np.set_printoptions(precision=3)
print P
[[ 0.5 0. 0.5 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0.333 0. 0. 0. 0. 0.667 0. 0. 0. ] [ 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0.333 0.333 0. 0. 0. 0.333 0. ] [ 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. 0. 0.25 0. 0.75 0. ] [ 0. 0. 0.25 0.25 0. 0. 0. 0.25 0. 0.25 ] [ 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0.333 0. 0. 0.333 0. 0. 0. 0. 0.333]]
mc = MarkovChain(P)
print mc
<bound method MarkovChain.__repr__ of Markov chain with transition matrix P = [[ 0.5 0. 0.5 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0.333 0. 0. 0. 0. 0.667 0. 0. 0. ] [ 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0.333 0.333 0. 0. 0. 0.333 0. ] [ 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. 0. 0.25 0. 0.75 0. ] [ 0. 0. 0.25 0.25 0. 0. 0. 0.25 0. 0.25 ] [ 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0.333 0. 0. 0.333 0. 0. 0. 0. 0.333]]>
Let us perform a simulation with initial state 0
:
X0 = mc.simulate(init=0)
bins = list(range(11))
hist, bin_edges = np.histogram(X0, bins=bins)
x0 = hist/len(X0)
print x0
[ 0.665 0. 0.335 0. 0. 0. 0. 0. 0. 0. ]
Let us write a function:
def empirical_dist(init):
X = mc.simulate(init=init)
hist, bin_edges = np.histogram(X, bins=bins)
x= hist/len(X)
return x
The set of states [0, 2]
seems to be closed:
print empirical_dist(2)
[ 0.655 0. 0.345 0. 0. 0. 0. 0. 0. 0. ]
The set of states [1, 6, 8]
seems to be closed:
print empirical_dist(1)
[ 0. 0.389 0. 0. 0. 0. 0.348 0. 0.263 0. ]
print empirical_dist(6)
[ 0. 0.372 0. 0. 0. 0. 0.367 0. 0.261 0. ]
print empirical_dist(8)
[ 0. 0.39 0. 0. 0. 0. 0.349 0. 0.261 0. ]
3
and 4
seem to communicate with each other and eventually get absorbed into [1, 6, 8]
:
print empirical_dist(3)
[ 0. 0.394 0. 0.001 0.001 0. 0.343 0. 0.261 0. ]
print empirical_dist(4)
[ 0. 0.4 0. 0.001 0.002 0. 0.343 0. 0.254 0. ]
5
is an absorbind state:
print empirical_dist(5)
[ 0. 0. 0. 0. 0. 1. 0. 0. 0. 0.]
From 9
, the path seems to be absorbed into [1, 6, 8]
:
print empirical_dist(9)
[ 0. 0.385 0. 0. 0.002 0. 0.355 0. 0.255 0.003]
From 7
, the path gets absorbed into either [0, 2]
or [1, 6, 8]
, depending on the realization:
print empirical_dist(7)
[ 0. 0.378 0. 0.007 0.008 0. 0.354 0.001 0.252 0. ]
print empirical_dist(7)
[ 0. 0.377 0. 0.001 0.001 0. 0.355 0.002 0.264 0. ]
N=50
print np.mean([empirical_dist(7) for i in range(N)], axis=0)
[ 0.173 0.287 0.086 0.001 0.001 0. 0.256 0.001 0.193 0. ]
As expected, the Markov chain mc
is reducible:
mc.is_irreducible
False
Communication classes of mc
:
mc.num_communication_classes
6
mc.communication_classes
[array([0, 2]), array([1, 6, 8]), array([3, 4]), array([5]), array([9]), array([7])]
Recurrent classes of mc
:
mc.num_recurrent_classes
3
mc.recurrent_classes
[array([0, 2]), array([1, 6, 8]), array([5])]
Recurrent states of mc
:
recurrent_states = np.concatenate(mc.recurrent_classes)
print recurrent_states
[0 2 1 6 8 5]
Transient states of mc
:
transient_states = np.setdiff1d(np.arange(mc.n), recurrent_states)
print transient_states
[3 4 7 9]
Canonical form of P
:
permutation = np.concatenate((recurrent_states, transient_states))
print permutation
[0 2 1 6 8 5 3 4 7 9]
print P[permutation, :][:, permutation]
[[ 0.5 0.5 0. 0. 0. 0. 0. 0. 0. 0. ] [ 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0.333 0.667 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0.25 0.75 0. 0. 0. 0. 0. ] [ 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. ] [ 0. 0. 0. 0. 0.333 0. 0.333 0.333 0. 0. ] [ 0. 0.25 0. 0. 0. 0. 0.25 0. 0.25 0.25 ] [ 0. 0. 0.333 0. 0. 0. 0. 0.333 0. 0.333]]
Decompose P
into irreducible submatrices:
for i, recurrent_class in enumerate(mc.recurrent_classes):
print 'P{0} ='.format(i)
print P[recurrent_class, :][:, recurrent_class]
P0 = [[ 0.5 0.5] [ 1. 0. ]] P1 = [[ 0.333 0.667 0. ] [ 0. 0.25 0.75 ] [ 1. 0. 0. ]] P2 = [[ 1.]]
The Markov chain mc
above has three stationary distributions,
one for each of the three recurrent classes:
vecs = mc.stationary_distributions
print vecs
[[ 0.667 0. 0.333 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0.391 0. 0. 0. 0. 0.348 0. 0.261 0. ] [ 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. ]]
print np.dot(vecs, P)
[[ 0.667 0. 0.333 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0.391 0. 0. 0. 0. 0.348 0. 0.261 0. ] [ 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. ]]
Q = np.zeros((11, 11))
Q[0, [1, 5]] = 1./2
Q[1, [2, 10]] = 1./2
Q[2, 7] = 1.
Q[3, 4] = 1.
Q[4, 5] = 1.
Q[5, [2, 6]] = 1./2
Q[6, [7, 8]] = 1./2
Q[7, 4] = 1.
Q[8, 0] = 1.
Q[9, 4] = 1.
Q[10, [3, 8, 9]] = 1./3
print Q
[[ 0. 0.5 0. 0. 0. 0.5 0. 0. 0. 0. 0. ] [ 0. 0. 0.5 0. 0. 0. 0. 0. 0. 0. 0.5 ] [ 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. ] [ 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. ] [ 0. 0. 0.5 0. 0. 0. 0.5 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. 0. 0. 0.5 0.5 0. 0. ] [ 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. ] [ 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0.333 0. 0. 0. 0. 0.333 0.333 0. ]]
mc2 = MarkovChain(Q)
This Markov chain mc2
is irreducible:
mc2.is_irreducible
True
Whether it is aperiodic:
mc2.is_aperiodic
False
The priod mc2
:
mc2.period
4
Cyclic classes of mc2
:
cyclic_classes = mc2.cyclic_classes
print cyclic_classes
[array([0, 4]), array([1, 5]), array([ 2, 6, 10]), array([3, 7, 8, 9])]
Cyclic normal form of Q
:
permutation = np.concatenate(cyclic_classes)
print permutation
[ 0 4 1 5 2 6 10 3 7 8 9]
np.set_printoptions(precision=2)
print Q[permutation, :][:, permutation]
[[ 0. 0. 0.5 0.5 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0.5 0. 0.5 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0.5 0.5 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. ] [ 0. 0. 0. 0. 0. 0. 0. 0. 0.5 0.5 0. ] [ 0. 0. 0. 0. 0. 0. 0. 0.33 0. 0.33 0.33] [ 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]]
Powers of Q
:
R = Q[permutation, :][:, permutation]
R /= np.sum(R, axis=1, keepdims=True)
print R
[[ 0. 0. 0.5 0.5 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0.5 0. 0.5 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0.5 0.5 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. ] [ 0. 0. 0. 0. 0. 0. 0. 0. 0.5 0.5 0. ] [ 0. 0. 0. 0. 0. 0. 0. 0.33 0. 0.33 0.33] [ 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]]
print np.linalg.matrix_power(R, 5)
[[ 0. 0. 0.1 0.9 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0.12 0.88 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0.5 0.46 0.04 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0.5 0.44 0.06 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. 0. 0. 0. 0.75 0.25 0. ] [ 0. 0. 0. 0. 0. 0. 0. 0.04 0.69 0.23 0.04] [ 0. 0. 0. 0. 0. 0. 0. 0.03 0.71 0.24 0.03] [ 0.25 0.75 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0.25 0.75 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0.21 0.79 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0.25 0.75 0. 0. 0. 0. 0. 0. 0. 0. 0. ]]
S = np.linalg.matrix_power(R, 8)
for k in range(5):
S = S.dot(R)
print 'R^{0} ='.format(9 + k)
print S, '\n'
R^9 = [[ 0. 0. 0.12 0.88 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0.12 0.88 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0.5 0.44 0.06 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0.5 0.44 0.06 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. 0. 0. 0.02 0.72 0.24 0.02] [ 0. 0. 0. 0. 0. 0. 0. 0.02 0.72 0.24 0.02] [ 0. 0. 0. 0. 0. 0. 0. 0.02 0.72 0.24 0.02] [ 0.24 0.76 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0.24 0.76 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0.24 0.76 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0.24 0.76 0. 0. 0. 0. 0. 0. 0. 0. 0. ]] R^10 = [[ 0. 0. 0. 0. 0.5 0.44 0.06 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0.5 0.44 0.06 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. 0. 0. 0.02 0.72 0.24 0.02] [ 0. 0. 0. 0. 0. 0. 0. 0.02 0.72 0.24 0.02] [ 0.24 0.76 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0.24 0.76 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0.24 0.76 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0.12 0.88 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0.12 0.88 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0.12 0.88 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0.12 0.88 0. 0. 0. 0. 0. 0. 0. ]] R^11 = [[ 0. 0. 0. 0. 0. 0. 0. 0.02 0.72 0.24 0.02] [ 0. 0. 0. 0. 0. 0. 0. 0.02 0.72 0.24 0.02] [ 0.24 0.76 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0.24 0.76 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0.12 0.88 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0.12 0.88 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0.12 0.88 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0.5 0.44 0.06 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0.5 0.44 0.06 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0.5 0.44 0.06 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0.5 0.44 0.06 0. 0. 0. 0. ]] R^12 = [[ 0.24 0.76 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0.24 0.76 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0.12 0.88 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0.12 0.88 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0.5 0.44 0.06 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0.5 0.44 0.06 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0.5 0.44 0.06 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. 0. 0. 0.02 0.72 0.24 0.02] [ 0. 0. 0. 0. 0. 0. 0. 0.02 0.72 0.24 0.02] [ 0. 0. 0. 0. 0. 0. 0. 0.02 0.72 0.24 0.02] [ 0. 0. 0. 0. 0. 0. 0. 0.02 0.72 0.24 0.02]] R^13 = [[ 0. 0. 0.12 0.88 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0.12 0.88 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0.5 0.44 0.06 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0.5 0.44 0.06 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. 0. 0. 0.02 0.72 0.24 0.02] [ 0. 0. 0. 0. 0. 0. 0. 0.02 0.72 0.24 0.02] [ 0. 0. 0. 0. 0. 0. 0. 0.02 0.72 0.24 0.02] [ 0.24 0.76 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0.24 0.76 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0.24 0.76 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0.24 0.76 0. 0. 0. 0. 0. 0. 0. 0. 0. ]]
So far, periodicity is left undefined for reducible Markov chains:
mc3 = MarkovChain([[1, 0], [0, 1]])
mc3.is_irreducible
False
mc3.is_aperiodic
--------------------------------------------------------------------------- NotImplementedError Traceback (most recent call last) <ipython-input-49-e2bd2b21b5ad> in <module>() ----> 1 mc3.is_aperiodic /Users/oyama/Dropbox/Development/graph_tools/mc_tools.py in is_aperiodic(self) 133 if not self.is_irreducible: 134 raise NotImplementedError( --> 135 'Not defined for a reducible Markov chain' 136 ) 137 else: NotImplementedError: Not defined for a reducible Markov chain
mc3.period
--------------------------------------------------------------------------- NotImplementedError Traceback (most recent call last) <ipython-input-50-13254b9bcf14> in <module>() ----> 1 mc3.period /Users/oyama/Dropbox/Development/graph_tools/mc_tools.py in period(self) 142 if not self.is_irreducible: 143 raise NotImplementedError( --> 144 'Not defined for a reducible Markov chain' 145 ) 146 else: NotImplementedError: Not defined for a reducible Markov chain
mc3.cyclic_classes
--------------------------------------------------------------------------- NotImplementedError Traceback (most recent call last) <ipython-input-51-95e38e278a88> in <module>() ----> 1 mc3.cyclic_classes /Users/oyama/Dropbox/Development/graph_tools/mc_tools.py in cyclic_classes(self) 151 if not self.is_irreducible: 152 raise NotImplementedError( --> 153 'Not defined for a reducible Markov chain' 154 ) 155 else: NotImplementedError: Not defined for a reducible Markov chain