Posted March 10, 2013 by Hunter

# Python Script for Markov Chains Consider the following problem in order to understand why we would need to create a Python Script for Markov Chains:
A certain attire-fickle girl changes her outfit every hour. Her opinion for what outfit she should wear is always dependant only upon the outfit she is currently wearing. When she is wearing a blue-based outfit, she has a 30% likelihood of remaining in a blue-based outfit, 20% likelihood of changing into a green outfit, and a 50% likelihood of changing into a red outfit. When she is in a green outfit, she has a 10% chance of changing into a blue outfit, a 10% chance of not changing, and an 80% chance of changing into a red outfit. When in a red outfit, she has a 40% chance of changing into a blue outfit, a 30% chance of changing into a green outfit, and a 30% chance of not changing. If she is currently wearing a green outfit, what is the likelihood, in five hours, of her wearing each color outfit?

Perhaps in life more practical examples will find you, but for our purposes here, how can the answer be found?

The answer is: Markov chains. This is a mathematical tool that considers a current state, the likelihood of changing to another state given the current state, and gives the likelihood of being in each state after a certain number of discrete steps. Python is wonderful language to use to create mathematical tools, so let’s consider how to craft a python script that can create and handle Markov chains.

```## This script uses the numpy moduel
import numpy as np

## The names of the states
s_name = []

## Gets how many states in the Markov chain
hold = 1
while hold == 1:
try:
states = int(raw_input('How many states? '))
hold = 0
except:
print 'Dude, seriously.  This needs to be a number.'

## Gets the names of the states
for i in range(1, states+1):
s_name.append(raw_input('What is the name of state '+str(i)+'? '))

## Creates an array of 'states' dimensions
P = np.zeros(shape=(states,states))

## Fills P with probabilities for each state
i = 1
while i &lt; states+1:
sumtoone = 0
print ''
j = 1
while j &lt; states+1:
try:
P[i-1][j-1] = float(raw_input('When in state '+s_name[i-1]+', what is the probability of the next state being '+s_name[j-1]+'? '))
sumtoone += P[i-1][j-1]
hold = 0
j += 1
except:
print 'error'
if sumtoone != 1:
print 'The sum of the probabilities you entered for next state when in state '+s_name[i-1]+' does not sum to 1.'
else:
i += 1
print''
print P
print ''
print 'What is the current state?'
for i in range(0,states):
print 'State 1: '+str(s_name[i])
hold = 1
## Gets the current state
while hold == 1:
try:
c_x = int(raw_input('Current state: '))
if (c_x &gt; 0) &amp; (c_x &lt;= states):
hold = 0
else:
print 'Not a valid state.'
except:
print 'error'

## Creates a 1-dimensional array representing the current state
x = np.zeros(shape=(1,states))
x[c_x-1] = 1

print ''
hold = 1
## Gets how many states into the future the user wants to look
while hold == 1:
try:
n = int(raw_input('Look how many steps into the future?  '))
hold = 0
except:
print 'error'

## Creates a copy of P
## If I simply said P_mult = P, then when I edited P_mult
## P would be edited as well, and we need P to remain what it is
P_mult = np.zeros(shape=(states,states))
for i in range(0,states):
for j in range(0,states):
element = P[i][j]
P_mult[i][j] = element

## Uses the numpy.dot function to matrix multiply P by itself n times
i = 1
while i &lt; n:
P_mult = np.dot(P,P_mult)
i += 1

## Matrix multiplies x by P^n
f_x = np.dot(x,P_mult)
print ''
print 'The probability distribution after '+str(n)+' steps is'
print f_x
```

## Hunter 