In [190]:

    

%pylab inline


    

    




Populating the interactive namespace from numpy and matplotlib

In [191]:

    

import numpy as np
from numpy import linalg as LA
import pylab as pl
from collections import defaultdict

s = defaultdict(dict) #states_transition
# sunny & cloudy example
# s[1][1] = 0.75
# s[1][2] = 0.25
# s[2][1] = 0.4
# s[2][2] = 0.6
s[1][1] = 0.7
s[1][2] = 0.3
s[1][3] = 0
s[2][1] = 0
s[2][2] = 0.85
s[2][3] = 0.15
s[3][1] = 0
s[3][2] = 0
s[3][3] = 1

initial_state = 1
n = 20


    

In [192]:

    

#create matrix array
n_states = len(s)

# this stuff is needed to sort matrix in form
# [ Q 0 ]
# [ R I ]
absorbing_states = [ i for i in xrange(1, n_states + 1) if s[i][i] == 1 ]
state_reindex = []

arr = []
initial_arr = []
states_to_enum = [ i for i in xrange(1, n_states + 1) if i not in absorbing_states ]
for to_state in states_to_enum:
    state_reindex.append(to_state)
    # append row
    arr.append([ s[from_state][to_state] for from_state in states_to_enum ])
    if to_state == initial_state:
        initial_arr.append([1])
    else:
        initial_arr.append([0])

        
a = np.matrix(arr)
print a

#create initial array
initial = np.matrix(initial_arr)


    

    




[[ 0.7   0.  ]
 [ 0.3   0.85]]

In [197]:

    

def noAbs():
    values = []
    values.append(initial)
    curent_matrix = 1
    for i in range(n):
        curent_matrix = curent_matrix * a
        values.append(curent_matrix * initial)

    for i in xrange(0, n_states):
        state_p = [ v.item(i) for v in values ]
        pl.bar(range(n + 1), state_p, label="state", width=3)
        for iv in range(len(state_p)):
            pl.text(iv, state_p[iv], "{0}".format(state_p[iv]))
        pl.title("state {0}".format(i+1))
        pl.show()

# no absorbing
if not absorbing_states:
    noAbs()
else:
    print a
    i_q = np.identity(len(states_to_enum)) - a
    print i_q
    N = LA.inv(i_q)
    print N


    

    




[[ 0.7   0.  ]
 [ 0.3   0.85]]
[[ 0.3   0.  ]
 [-0.3   0.15]]
[[ 3.33333333  0.        ]
 [ 6.66666667  6.66666667]]

In [193]:

    




    

In [193]:

    




    

In [193]: