-
Notifications
You must be signed in to change notification settings - Fork 1
Markov Chains
The Markov Chain class in the stochastic processes module supports the analysis of discrete time markov chains. All of the Markov Chain procedures work with both rational and floating point numbers. Some procedures required the user to explicity specify that rational output is required. Matrix operations on rational numbers may be slower, and are more memory intensive. However, no accuracy is lost due to machine precision. The default output for Markov Chain procedures is a numpy array. The vector_display and matrix_display functions allow for the output to be displayed with the names of the states in the state space.
Markov Chains are implemented as a class object. A transition probability matrix is required to initiate the Markov Chain. Optionally, an initial state for the Markov Chain can be specified and the states in the state space can be name. Procedures that deal with conditional probability usually require the initial conditions for the Markov Chain. Others do not. If no names are chosen for the states in the state space, they are name 0,1,...,n-1 by default.
MarkovChain(P,[init],[states])In [3]: Y = MarkovChain(P=[[.97,.03],[.04,.96]],
init=[.7,.3],
states=['red','blue'])
In [4]: Y.display()
The transition probability matrix:
red blue
red 0.97 0.03
blue 0.04 0.96
----------------------------------------
The initial system state:
Prob
red 0.7
blue 0.3The absorption probability procedure computes the probability of absorption into the specified state, given any other starting state for the Markov Chain. The output is a vector where v[i] is the probability of absorption, given that i is the starting state. The state specified for this procedure must be an absorbing state.
X.absorbing(state)In [8]: P = np.array([
...: [1,0,0,0],
...: [Rational(3,10),Rational(4,10),Rational(1,10),Rational(2,10)],
...: [0,0,1,0],
...: [Rational(2,10),Rational(3,10),Rational(4,10),Rational(1,10)]
...: ])
In [9]: X = MarkovChain(P, states = ['red','blue','black','green'])
In [10]: abs_red = X.absorption_prob('red')
In [11]: vector_display(abs_red,X.state_space)
Out[11]:
Prob
red 1
blue 31/48
black 0
green 7/16This procedure computes the expected number of steps to absorption for the Markov Chain, given any starting point for the Markov Chain. This procedure requires the Markov Chain the have at least one absorbing state. The output is a vector where v[i] is the expected number of steps to absorption, given that i is the starting state
X.absorption_steps()In [12]: abs_steps = X.absorption_steps()
In [13]: vector_display(abs_steps, X.state_space)
Out[13]:
Prob
red 0
blue 55/24
black 0
green 15/8This procedures identifies equivalence classes in the Markov Chain and labels them as transient or recurrent. The method for classifying states is described in Weiss, B. 'A non-recursive algorithm for classifying the states of a finite Markov chain'. 1987. European Journal of Operations Research. The output is a dictionary of of labeled equivalence classes.
X.classify_states()In [14]: X = MarkovChain(P =[
....: [0, 0.1, 0.5, 0, 0.2, 0.2, 0, 0],
....: [0.1, 0.2, 0, 0, 0.6, 0, 0, 0.1],
....: [0, 0, 1, 0, 0, 0, 0, 0],
....: [0, 0, 0, 0, 1, 0, 0, 0],
....: [0, 0, 0, 1, 0, 0, 0, 0],
....: [0, 0, 0, 0, 0, 0.2, 0, 0.8],
....: [0, 0.3, 0, 0, 0, 0.4, 0.1, 0.2],
....: [0, 0, 0, 0, 0, 0.6, 0, 0.4]
....: ],
....: states = ['blue','green','purple','red',
....: 'yellow', 'black', 'orange', 'white'])
In [15]: X.classify_states()
Out[15]:
{'Recurrent 1': ['red', 'yellow'],
'Recurrent 2': ['black', 'white'],
'Recurrent 3': ['purple'],
'Transient': ['blue', 'green', 'orange']}This procedure computes the long run fraction of time spent in each state, given that the Markov Chain starts at any other state. These probabilities can be computed for both reducible and irreducible Markov Chains. The user must specify the method as rational if rational output is desired. The floating point method is less computationally expensive and uses less memory.
X.long_run_probs()In [16]: P = np.array([
....: [Rational(2,10),0,0,Rational(8,10)],
....: [0,1,0,0],
....: [Rational(1,10),Rational(4,10),Rational(2,10),Rational(3,10)],
....: [Rational(7,10),0,0,Rational(3,10)]
....: ])
In [17]: X = MarkovChain(P, states = ['red','blue','black','green'])
In [18]: Pi = X.long_run_probs(method = 'rational')
In [19]: matrix_display(Pi,X.state_space)
Out[19]:
red blue black green
red 7/15 0 0 8/15
blue 0 1 0 0
black 7/30 1/2 0 4/15
green 7/15 0 0 8/15This computes the probability of reaching a sequence of states, given the specified conditions. The desired states are entered as a list of tuples, consisting of two elements. The first element is the time period in which the state is realized, and the second element is the name of the state. A second list of sequences can be entered to compute a conditional probability. If the probability is not conditional, the initial conditions must be specified for the Markov Chain.
X.probability(states,[given])In [23]: Y = MarkovChain(P=
....: [[Rational(97,100), Rational(3,100)],
....: [Rational(4,100), Rational(96,100)]],
....: init = [Rational(7,10), Rational(3,10)],
....: states = ['red','blue'])
In [24]: Y.probability([(2,'blue'),(1,'red'),(0,'red')])
Out[24]:
2037
──────
100000
In [25]: Y.probability([(2,'blue')],given=[(1,'blue')])
Out[25]:
24
──
25This procedure computes a matrix B, where B[i,j] = True if state j is reachable from state i, and is false otherwise. The method for computing this matrix is described in Weiss, B. 'A non-recursive algorithm for classifying the states of a finite Markov chain'. 1987. European Journal of Operations Research. The output is a dictionary of of labeled equivalence classes.
X.reachability()In [26]: X = MarkovChain(P =[
....: [0, 0.1, 0.5, 0, 0.2, 0.2, 0, 0],
....: [0.1, 0.2, 0, 0, 0.6, 0, 0, 0.1],
....: [0, 0, 1, 0, 0, 0, 0, 0],
....: [0, 0, 0, 0, 1, 0, 0, 0],
....: [0, 0, 0, 1, 0, 0, 0, 0],
....: [0, 0, 0, 0, 0, 0.2, 0, 0.8],
....: [0, 0.3, 0, 0, 0, 0.4, 0.1, 0.2],
....: [0, 0, 0, 0, 0, 0.6, 0, 0.4]
....: ])
In [27]: B = X.reachability()
In [28]: matrix_display(B,X.state_space)
Out[28]:
0 1 2 3 4 5 6 7
0 True True True True True True False True
1 True True True True True True False True
2 False False True False False False False False
3 False False False True True False False False
4 False False False True True False False False
5 False False False False False True False True
6 True True True True True True True True
7 False False False False False True False TrueThe steady state procedure computes the long run probability that the Markov Chain is in each state. This procedure only works for irreducible Markov Chains. For reducible Markov Chains, use the long run probability procedure. The user must specify the method as rational if rational output is desired. The floating point method is less computationally expensive and uses less memory.
X.steady_state([method])In [29]: X = MarkovChain([[Rational(3,10),Rational(3,10),Rational(4,10),0,0],
....: [0,Rational(3,10),Rational(3,10),Rational(4,10),0],
....: [0,0,Rational(3,10),Rational(3,10),Rational(4,10)],
....: [Rational(3,10),Rational(3,10),Rational(4,10),0,0],
....: [Rational(3,10),Rational(3,10),Rational(4,10),0,0]],
....: states = ['blue','green','black', 'yellow','orange'])
In [30]: steady = X.steady_state(method='rational')
In [31]: vector_display(steady,X.state_space)
Out[31]:
Prob
blue 147/1070
green 21/107
black 37/107
yellow 39/214
orange 74/535This procedure computes the transition probability matrix for n steps of the Markov Chain. The user must specify the method as rational if rational output is desired. The floating point method is less computationally expensive and uses less memory.
X.trans_mat(n,[method])In [34]: X = MarkovChain([[.3,.3,.4,0,0],
[0,.3,.3,.4,0],
[0,0,.3,.3,.4],
[.3,.3,.4,0,0],
[.3,.3,.4,0,0]])
In [35]: X3 = X.trans_mat(n=3)
In [36]: matrix_display(X3,X.state_space)
Out[36]:
0 1 2 3 4
0 0.147 0.201 0.349 0.171 0.132
1 0.135 0.198 0.345 0.186 0.136
2 0.126 0.189 0.342 0.195 0.148
3 0.147 0.201 0.349 0.171 0.132
4 0.147 0.201 0.349 0.171 0.132