Lab 7-1: Markov Chains - Basic Examples#
(Steven Pestana, 2019. Derived from CEE599_MarkovChain_Lab.m, Jessica Lundquist, Oct 23, 2015)
(see Markov Chain Weather Example in Lecture notes, except note that we now have an initial condition, which would be like just saying today’s weather was dry)
import numpy as np
from scipy.linalg import eig
import scipy.stats as stats
from scipy import sparse
import matplotlib.pyplot as plt
%matplotlib inline
Example - City and Suburbs#
We are going to use a Markov Chain to look at the conditional probabilities of people moving between the city and the suburbs outside the city, at steps of 1 year.
Let’s define two states, in our model someone can either live in the city or the suburbs:
city = 0
suburbs = 1
Set up the Markov probability matrix, start with just an empty 2x2 array:
Pmarkov = np.zeros([2,2]) # create an empty 2-by-2 numpy array
Pmarkov
array([[0., 0.],
[0., 0.]])
Given someone lives in the city, they have 95% chance of staying and 5% chance of moving:
# Assign our given values to the table:
Pmarkov[0,0] = 0.95 # note that we use the array indices to describe probability of going from state 0 to state 0
Pmarkov[0,1] = 0.05 # probability of going from state 0 to state 1
Given someone lives in the suburbs, they have a 97% change of staying and a 3% chance of moving:
# Assign our given values to the table:
Pmarkov[1,0] = 0.03 # probability of going from state 1 to state 0
Pmarkov[1,1] = 0.97 # probability of going from state 1 to state 1
View the complete matrix:
print(Pmarkov)
[[0.95 0.05]
[0.03 0.97]]
Before going into matrix notation, let’s think this through. In the initial state, 60% of people live in the city & 40% live in the suburb.
percent_city0 = 0.6
percent_suburb0 = 0.4
Step 1st year using alebra
Take one step through the Markov chain. What fraction of people live in the city after one step?
This can be defined as the fraction of people who live in the city and stay in the city + the fraction of people who lived in the suburbs but move into the city.
\(p_{city,1}\) = 0.95 * 0.6 + 0.03 * 0.4
pcity1 = Pmarkov[0,0]*percent_city0 + Pmarkov[1,0]*percent_suburb0
print(pcity1)
0.582
What fraction of people live in the suburbs?
Similar to above, this is the fraction of people who lived in the city and moved to the suburbs plus the fraction of people who live in the suburbs and stay in the suburbs.
\(p_{suburb,1}\) = 0.05 * 0.6 + 0.97 * 0.4
psuburb1 = Pmarkov[0,1]*percent_city0 + Pmarkov[1,1]*percent_suburb0
print(psuburb1)
0.41800000000000004
Draw a picture on scratch paper of what you are doing in terms of matrix multiplication.
Step 1st year using matrix multiplication
Now use matrix algebra using numpy to get the same answer and check your work.
# initial state of populations for state 0 and state 1
percent0 = np.array([percent_city0, percent_suburb0])
# multiply initial states with Markov matrix
percent1 = np.dot(percent0, Pmarkov)
# check values of x1 here against what you got for pcity1 and pusuburb1
print(percent1)
[0.582 0.418]
Step 2nd year
Now we can make this the current state and take a second step with the matrix (using both one step and matrix notation to check our work):
# take a second step with the Markov matrix:
pcity2 = Pmarkov[0,0]*pcity1 + Pmarkov[1,0]*psuburb1
psuburb2=Pmarkov[0,1]*pcity1 + Pmarkov[1,1]*psuburb1
print(pcity2,psuburb2)
# take a second step with the Markov matrix (using np.dot)
percent2 = np.dot(percent1,Pmarkov)
print(percent2[0],percent2[1])
0.5654399999999999 0.43456000000000006
0.5654399999999999 0.43456000000000006
Good, we got the same answer (within rounding error), let’s just use matrix math for the rest of the lab.
Step 3rd year
We can also make multiple steps by raising the Markov matrix to a power:
# three steps all at once by taking our Pmarkov matrix to the 3rd power, then multipying by our original population percentage vector
percent3 = np.dot(percent0, np.linalg.matrix_power(Pmarkov,3))
print(percent3)
# we can also take a 3rd step by just multiplying through from the 2nd year
percent3_check = np.dot(percent2,Pmarkov)
print(percent3_check)
[0.5502048 0.4497952]
[0.5502048 0.4497952]
# 10 steps, then 100 steps:
percent10 = np.dot(percent0, np.linalg.matrix_power(Pmarkov,10))
percent100 = np.dot(percent0, np.linalg.matrix_power(Pmarkov,100))
print(percent10)
print(percent100)
[0.4727374 0.5272626]
[0.37505382 0.62494618]
Find the steady state probabilities using the eigenvector function
And now, use the eigenvector function in numpy to see if x100 is the same as the steady state matrix:
Here, we are using the concept that for a square matrix A, the column vector v is an eigenvector of A is theere is a scalar multiple \(\lambda\), such that multipliplying it times v is the same as applying the matrix multiplication of A on v.
\( \lambda v = Av \)
This fits the steady state matrix definition of
\( \pi = \pi P \)
for the case of \(\lambda\) = 1 and P = the transpose of A.
# So, we solve the eigenvalue problem of the transpose of the markov matrix A.
# See https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eig.html#scipy.linalg.eig
W, V = eig(Pmarkov.T)
V
array([[-0.70710678, -0.51449576],
[ 0.70710678, -0.85749293]])
Find the column of V associated with the eigenvalue of 1, but allow for some margin of round-off error
(from the function documentation: the column V[:,i] is the eigenvector corresponding to the eigenvalue W[i])
# First, we can inspect W by eye (since this is a small matrix)
W
array([0.92+0.j, 1. +0.j])
We want the second column (index 1)
V[:,1]
array([-0.51449576, -0.85749293])
But these need to be normalized first
p = V[:,1] / np.sum(V[:,1])
p
array([0.375, 0.625])
We can write a function to find this for us:
def eig_pi(X, err=1e-9):
'''Find the value of W closest to 1 (within some round-off error) and then look at the associated column of V.
Input: tuple of (W (eigenvalue matrix), V (eigenvectors)) that are output from scipy.linalg.eig()'''
[W, V] = X
# find the column where the eigenvalue, W, is within 1 +/- err
col1 = np.where((W>=1-err) & (W<=1+err))
# normalize the eigenvector, V, for that column, so that it sums to one
p = V[:,col1] / np.sum(V[:,col1])
return p
p = eig_pi(eig(Pmarkov.T))
print('Eigenvector Steady State Probabilities: {}'.format(p.T[0][0]))
# Yes, these are very close to the values we got for percent100
print('Markov Chain 100-step Probabilities: {}'.format(percent100))
Eigenvector Steady State Probabilities: [0.375 0.625]
Markov Chain 100-step Probabilities: [0.37505382 0.62494618]
Example - Create a Markov Matrix from data#
Now given a timeseries of four states, determine the Markov Matrix and use it to generate a timeseries
# first, read in data that I generated (sequence of states 1, 2, 3, and 4):
data = np.genfromtxt('../data/markov_random4.txt',dtype=int)
# This counts the transitions from each state to the next and marks that count
S = sparse.csr_matrix((np.ones_like(data[:-1]), (data[:-1], data[1:])), dtype=float)
# This converts those counts to matrix form
tm = S.todense()
print(tm)
[[ 0. 0. 0. 0. 0.]
[ 0. 16. 53. 22. 28.]
[ 0. 31. 34. 12. 56.]
[ 0. 13. 34. 48. 28.]
[ 0. 58. 12. 41. 13.]]
Take a look at the matrix above, what does this represent?
We had 16 counts when the series transitioned from 1 to 1, 53 counts when it transitioned from 1 to 2, etc. We want to transition these from counts to frequencies. To do this, we need to normalize the transition matrix to get probabilities
tm_norm = tm / tm.sum(axis=1)
print(tm_norm) # This is our normalized transition matrix.
[[ nan nan nan nan nan]
[0. 0.13445378 0.44537815 0.18487395 0.23529412]
[0. 0.23308271 0.2556391 0.09022556 0.42105263]
[0. 0.10569106 0.27642276 0.3902439 0.22764228]
[0. 0.46774194 0.09677419 0.33064516 0.10483871]]
/opt/hostedtoolcache/Python/3.7.17/x64/lib/python3.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in true_divide
"""Entry point for launching an IPython kernel.
# We take the above probabilities of transitions, and turn them into discrete CDF's.
# These will allow us to map random numbers generated from a uniform distribution into
# transitions that follow these probability rules.
tm_cdf = np.cumsum(tm_norm,1)
Now, we want a “random walk” for 2000 years that follows these transition probabilities.
n_years = 2000
q = np.random.uniform(0,1,n_years); # uniformly distributed random numbers n_years long
initialstate = 2; # give it an initial state, doesn't really matter which
Nrand = np.zeros_like(q) # initialize an array of the proper size, with the initial state
Nrand[0] = initialstate;
# Now, just like we did when we created monte carlo simulations from empirical CDFs,
# we use our uniform random numbers to look up the next state in the transition matrix
for i in range(1,n_years):
if q[i] <= tm_cdf[int(Nrand[i-1]),1]: #probability of transitioning from state i to 1
Nrand[i] = 1;
elif q[i] <= tm_cdf[int(Nrand[i-1]),2]: #transition from state i to 2
Nrand[i] = 2;
elif q[i] <= tm_cdf[int(Nrand[i-1]),3]: #transition from state i to 3
Nrand[i] = 3;
else:
Nrand[i] = 4;
# And how many times did state 3 appear 4 times?
Test3 = [Nrand[0:-3], Nrand[1:-2], Nrand[2:-1], Nrand[3:]] # stack our data 4 times, shifting it to the right by 1 each time
Test3 = np.stack(Test3, axis=1)
G2 = np.where((np.max(Test3, axis=1) == 3) & (np.min(Test3, axis=1) == 3))
# if both the maximum and the minimum are 3, then we have 4 threes in our sequence
frequencyof4threes = G2[0].size / Test3.shape[0]
print('Frequency of four 3s in a row = {}%'.format(100*np.round(frequencyof4threes,3)))
Frequency of four 3s in a row = 1.6%