unfair_casino_viterbi.py

"""
========================================
Unfair Casino Problem using Discrete HMM
========================================

This script shows how to use Decoding in a Discrete(Multinomial) HMM, which means
given the model parameters and an observed sequence we wish to find the most likely hidden state sequence
generating the same.
It uses the model given in http://www.rose-hulman.edu/~shibberu/MA490/MA490HMM.html#Example_Dishonest_Casino_
One may also refer to a great lecture series at http://vimeo.com/7175217

"""
print __doc__
import datetime
import numpy as np
import pylab as pl
from sklearn.hmm import MultinomialHMM
import random

###############################################################################
################################################################################
#Specifying HMM parameters

#We can obtain the starting probability by finding the stationary distribution from the transition matrix
#This denotes that the probabilities of starting in Fair state and Biased state are 2/3 and 1/3 respectively
start_prob = np.array([.67, 0.33])

#The Transition Matrix
#This means that given that current state is Fair, probability of remaining in fair state is .95 and transitioning to
#Biased state is .05. Similarly, given that current state is Biased, probability of transitioning to fair state is .1
#and remaining in the same state is .9
trans_mat = np.array([[0.95, 0.05],
                      [.1, 0.9]])

#The Emission Matrix
#This means that if the dice is fair, all 6 outcomes (1 to 6) have same probability,
#whereas, for the biased dice, 6 is observed with probability 0.5 and the remaining 
#with a probability 0.1 each                  
emiss_prob=np.array([[1.0/6,1.0/6,1.0/6,1.0/6,1.0/6,1.0/6],[.1,.1,.1,.1,.1,.5]])

#Initializing the model with 2 state(n_components) and specifying the transition matrix and number of iterations
model = MultinomialHMM(n_components=2,transmat=trans_mat,n_iter=10)

#Number of emitted symbols
model.n_symbols=6
#Assigning emission probability, start probabilty and random seed
model.emissionprob_=emiss_prob
model.startprob_=start_prob
model.prng=np.random.RandomState(10)

#dice_obs is a random sequence of observations. Note that some extra 6's have been added in middle of the sequence
#which might correspond to use of Biased Dice 
dice_obs=[random.randint(1,6) for i in range(0,20)]+[6,6,6,6,6,6]+[random.randint(1,6) for i in range(0,20)]
print "Dice Observations\n",dice_obs

#Note that the emitted symobols must vary from 0 through 5. Thus we subtract 1 from each element of dice_obs
obs=[x-1 for x in dice_obs]

#Finding the likely state sequence and log probabilty of the path using Viterbi algorithm and parameters specified in the model
logprob, state_sequence = model.decode(obs)
print "\nExpected State Sequence- F for fair and B for Biased"
state_sequence_symbols=['F' if x==0 else 'B' for x in state_sequence]
print state_sequence_symbols

print "Log probability of this most likely path is: ",logprob
print "Probability of this most likely path is    : ",np.exp(logprob)