ECE/CS 434 | MP5: HMM

Objective

In this MP, you will:

Explore two fundamental problems that are of interest to Hidden Markov Models (HMM):

What is the probability of a specific sequence of observations occuring given an HMM?

What is the optimal hidden state sequence given a sequence of observations and an HMM?

Implement either the forward or backward algorithm to answer the first problem.

Implement the Viterbi algorithm to answer the second problem.

References

Material in this MP is adapted from the iconic paper A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition by Lawrence R. Rabiner. You are encouraged to follow the pseudocode in this paper for your implementation.

You are highly encouraged to watch this video by Prof.Patterson from Westmont College before beginning the MP.

Imports & Setup

Installing requirements correctly

First. we will make sure that the correct versions of required modules are installed. This ensures that your local Python environment is consistent with the one running on the Gradescope autograder. Just convert the following cell to code and run:

Note: It’s preferred that your local environment matches the autograder to prevent possible inconsistencies. However, if you’re running into annoying Python version issues but haven’t had any issues getting consistent results on the autograder, there is no need to stress over it. Just skip for now and come back when you do encounter inconsistencies:) Ditto below.

WARNING: ENSURE THE FOLLOWING CELL IS MARKDOWN OR DELETED BEFORE SUBMITTING. THE AUTOGRADER WILL FAIL OTHERWISE.

if name == ‘main’: import sys !{sys.executable} -m pip install -r requirements.txt

Your imports

Write your import statements below. If Gradescope reports an error and you believe it is due to an unsupported import, check with the TA to see if it could be added.

import numpy as np

import random

import scipy

import pandas as pd

# This function is used to format test results. You don’t need to touch it.

def display_table(data):

from IPython.display import HTML, display

html = “<table>”

for row in data:

html += “<tr>”

for field in row:

html += “<td><h4>{}</h4><td>”.format(field)

html += “</tr>”

html += “</table>”

display(HTML(html))

Constructing an HMM

Recall that an HMM consists of three parameters:

A , the state transition probability matrix of size N×N , where N represents the number of hidden states(or true states).

B , the observation probability matrix of size N×M , where M represents the number of observed states.

π , an initial distribution vector of size 1×N

In this MP, we will be randomly generating these parameters. Run the cell below to see what a model looks like and how to access it. Notice you’ll get slightly different numbers each time you run due to the random functions used.

if __name__ == ‘__main__’:

def generate_matrix(R, C, mode):

model = np.zeros((R, C))

for i in range(R):

row = np.zeros(C)

if(mode == 0):

while(row.sum() != 1):

row = np.random.normal(0, 2, C)

row = row – row.min()

row = row / row.sum()

elif(mode == 1):

while(row.sum() != 1):

row = np.random.dirichlet(np.ones(C),size=1)

row = row / row.sum()

elif(mode == 2):

while(row.sum() != 1):

row = np.random.gamma(5, 5, C)

row = row / row.sum()

model[i] = row

return model

example_model = {‘A’: generate_matrix(5, 5, 0),

‘B’: generate_matrix(5, 3, 0),

‘pi’: pi}

print(‘A – State Transition Probability Matrix’)

print(pd.DataFrame(example_model[‘A’]))

print(‘The probability of transitioning from state 2 to state 3 is ‘ + str(example_model[‘A’][2][3]))

print(‘\n\nB – Observation Probability Matrix’)

print(pd.DataFrame(example_model[‘B’]))

print(‘The probability of being in hidden state 1 while observing observable state 2 is ‘ + str(example_model[‘B’][1][2]))

A – State Transition Probability Matrix

0 1 2 3 4

0 0.000000 0.248270 0.127124 0.085380 0.539226

1 0.082389 0.201989 0.231160 0.000000 0.484461

2 0.254655 0.103178 0.248717 0.000000 0.393449

3 0.273374 0.139821 0.024535 0.562270 0.000000

4 0.186978 0.000000 0.312920 0.249575 0.250527

The probability of transitioning from state 2 to state 3 is 0.0

B – Observation Probability Matrix

0 1 2

0 0.000000 0.926460 0.073540

1 0.430079 0.569921 0.000000

2 0.597683 0.402317 0.000000

3 0.481032 0.000000 0.518968

4 0.230072 0.000000 0.769928

The probability of being in hidden state 1 while observing observable state 2 is 0.0

Problem 1. Determining the Best Model

You are given a sequence of observations, O , and a few candidate HMMs, {λ1=(A1,B1,π1),λ2=(A2,B2,π2)…} . Your task is to find the model that best describes the sequence of observations.

You could approach this problem by calculating P(O∣λ) for every λ , and returning the index of the model yielding the highest likelihood.

Hint: The Forward Algorithm generates a table of size T×N where T is the length of the observation sequence. Then, P(O∣λ) is essentially the sum of the last row of this table. The Backward Algorithm works similarly. You may choose one to implement.

Note: You must use dynamic programming in your implementation. Otherwise, the tests will timeout.

Implement your algorithm in the find_model(O, models) function below. Do NOT change the function signature. You are, however, free to define and use helper functions.

def find_model(O, models):

“””Finds the HMM that is mostly likely to have generated a given sequence of observations

Args:

O: A list of observed states.

models: a list of HMM models. E.x.models[0][‘A’] accesses the transition matrix of the first HMM.

Returns: The index corresponding to the best-fit HMM.

“””

return 0

Run & Test

Use the cell below to run and test find_model(O, models).

if __name__ == ‘__main__’:

def test_metadata(N, M, T):

“””Gathers constants that describes a test

Args:

N: Number of hidden states.

M: Number of observable states.

T: Length of observable sequence.

“””

return {‘N’: N, ‘M’: M, ‘T’: T}

tests = [test_metadata(N=5, M=3, T=100),

test_metadata(N=10, M=5, T=200),

test_metadata(N=10, M=3, T=300)]

output = [[‘Test’, ‘% of Correct Trials’, ‘Grade’]]

for idx, test in enumerate(tests):

print(“Running test ” + str(idx) + “….”)

count = 0

for trial in range(50):

# 1. Create Candidate HMMs

models = []

pi = np.ones((test[‘N’]))

pi = pi / pi.sum()

for i in range(3):

models.append({‘A’: generate_matrix(test[‘N’], test[‘N’], i),

‘B’: generate_matrix(test[‘N’], test[‘M’], i),

‘pi’: pi})

# 2. Pick 1 HMM

TRUE_MODEL = random.randint(0,2)

# 3. Create observation based on model

OBS = []

for i in range(test[‘T’]):

cur_h_state = 0

if(i == 0):

cur_h_state = np.random.choice(np.arange(0, test[‘N’]), p=models[TRUE_MODEL][‘pi’])

else:

cur_h_state = np.random.choice(np.arange(0, test[‘N’]), p=models[TRUE_MODEL][‘A’][int(cur_h_state)])

OBS.append(int(np.random.choice(np.arange(0, test[‘M’]), p=models[TRUE_MODEL][‘B’][int(cur_h_state)])))

# 3. Test find_model

est_model = find_model(OBS, models)

if(est_model == TRUE_MODEL):

count = count + 1

score = 10 if(count > (50 / len(tests))) else 0

output.append([idx, ‘{:.1%}’.format(count / 50), “{:0.0f} / 10”.format(score)])

output.append([‘<i>👻 Hidden test 0 👻</i>’,'<i>???</i>’, ‘<i>???</i> / 10’])

output.append([‘<i>👻 Hidden test 1 👻</i>’,'<i>???</i>’, ‘<i>???</i> / 10’])

display_table(output)

Running test 0….

Running test 1….

Running test 2….

Test

% of Correct Trials

Grade

0

26.0%

0 / 10

1

40.0%

10 / 10

2

20.0%

0 / 10

👻 Hidden test 0 👻

???

??? / 10

👻 Hidden test 1 👻

???

??? / 10

Problem 2. Recovering the Hidden States

You are given a sequence of observed states and an HMM. Recover the most likely sequence of hidden states. You are encouraged to implement the Viterbi algorithm here.

def recover_states(O, model):

“”” Recovers the most likely sequence of hidden states given observations and a model

Args:

O: A list of observed states.

model: The HMM.

Returns: A list of hidden states

“””

return []

Run & Test

Use the cell below to run and test recover_states(O, model).

if __name__ == ‘__main__’:

def test_metadata(N, M, T):

“””

Args:

N: Number of hidden states.

M: Number of observable states.

T: Length of observable sequence.

“””

return {‘N’: N, ‘M’: M, ‘T’: T}

def random_guessing(T, M):

return [random.randint(0,M-1) for i in range(T)]

tests = [test_metadata(N=5, M=3, T=100),

test_metadata(N=10, M=3, T=200),

test_metadata(N=12, M=3, T=300)]

output = [[‘Test’, ‘% of States Uncovered’, ‘Random Guessing’, ‘Grade’]]

for idx, test in enumerate(tests):

score = 0

rand_score = 0

print(“Running test ” + str(idx) + “…..”)

for trial in range(100):

# 1. Generate model

pi = np.ones((test[‘N’]))

pi = pi / pi.sum()

model = {‘A’: generate_matrix(test[‘N’], test[‘N’], idx),

‘B’: generate_matrix(test[‘N’], test[‘M’], idx),

‘pi’: pi}

# 2. Generate observations and hidden states

OBS = []

TRUE_STATES = []

for i in range(test[‘T’]):

cur_h_state = 0

if(i == 0):

cur_h_state = np.random.choice(np.arange(0, test[‘N’]), p=model[‘pi’])

else:

cur_h_state = np.random.choice(np.arange(0, test[‘N’]), p=model[‘A’][int(cur_h_state)])

OBS.append(int(np.random.choice(np.arange(0, test[‘M’]), p=model[‘B’][int(cur_h_state)])))

TRUE_STATES.append(int(cur_h_state))

# 3. Test viterbi

est_states = recover_states(OBS, model)

random_states = random_guessing(test[‘T’], test[‘M’])

# 4. Calculate score

count = sum(x == y for x, y in zip(np.array(TRUE_STATES), np.array(est_states)))

score = score + count / test[‘T’]

rand_count = sum(x == y for x, y in zip(np.array(TRUE_STATES), np.array(random_states)))

rand_score = rand_score + rand_count / test[‘T’]

grade = 10 if(score > rand_score) else 0

output.append([idx, ‘{:.1%}’.format(score / 100), ‘{:.1%}’.format(rand_score / 100), “{:0.0f} / 10”.format(grade)])

output.append([‘<i>👻 Hidden test 0 👻</i>’,'<i>???</i>’, ‘<i>???</i>’, ‘<i>???</i> / 10’])

output.append([‘<i>👻 Hidden test 1 👻</i>’,'<i>???</i>’, ‘<i>???</i>’, ‘<i>???</i> / 10’])

display_table(output)

Running test 0…..

Running test 1…..

Running test 2…..

Test

% of States Uncovered

Random Guessing

Grade

0

0.0%

20.7%

0 / 10

1

0.0%

10.3%

0 / 10

2

0.0%

8.3%

0 / 10

👻 Hidden test 0 👻

???

???

??? / 10

👻 Hidden test 1 👻

???

???

??? / 10

Rubric

Since HMM are probabilistic models, we will run your code many times and your grade should stabilize after many trials by the Law of Large Numbers. The script provided runs 50 trials to make testing more pleasant. The actual grading script will run more trials. On Gradescope, we will display your total score, including hidden test cases. You can expect your final grade for this MP to match what you see after submission.

Determing the Best Model (50 points)

You will be graded on the 3 sets of provided cases as well as 2 sets of hidden cases. 10 points will be rewarded if your algorithm performs better than random guessing. For example, if the test provides 5 models, your algorithm will receive full points if it outputs the correct model over 20% of the times.

Recovering the Hidden States (50 points)

You will be graded on the 3 sets of provided cases as well as 2 sets of hidden cases. 10 points will be rewarded if your algorithm recovered more states than random guessing.

Submission Guideline

This Jupyter notebook is the only file you need to submit on Gradescope. If you are working in a pair, make sure your partner is correctly added on Gradescope and that both of your names are filled in at the top of this file.

Make sure any code you added to this notebook, except for import statements, is either in a function or guarded by __main__(which won’t be run by the autograder). Gradescope will give you immediate feedback using the provided test cases. It is your responsibility to check the output before the deadline to ensure your submission runs with the autograder.