MP5: HMM

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.