Duda, Hart, Pattern Classification, Section 3.10, Exercise 11, Consider the use of hidden Markov models for classifying sequences of four visible states, A-D....

markov_duda_hart_char_sequence.py

import numpy as np
from hmmlearn import hmm
from matplotlib import pyplot as plt
from matplotlib import rcParams

rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16

# Trying to implement a solution for Duda-Hart Pattern Recognition book, Section 3.10, Exercise 11
# Consider the use of hidden Markov models for classifying sequences of four visible 
# states, A-D. Train two hidden Markov models, each consisting of three hidden states 
# (plus a null initial state and a null final state), fully connected, with the following 
# data. Assume that each sequence starts with a null symbol and ends with an end null 
# symbol (not listed). 


X = ["A", "B", "C", "D"]
Y = ["A", "B", "C", "D"]

m = len(Y)
n = len(X)

def convert_sequence(list):
    x = []
    for u in list:
        if u == 'A':
            x.append([0])
        elif u == 'B':
            x.append([1])
        elif u == 'C':
            x.append([2])
        elif u == 'D':
            x.append([3])
    return x


def plot_transition_matrix(matrix, title):
    plt.figure(figsize=(n + 1, m + 1))
    for krow, row in enumerate(matrix):
        plt.text(5, 10 * krow + 15, Y[krow],
                 horizontalalignment='center',
                 verticalalignment='center')
        for kcol, num in enumerate(row):
            if krow == 0:
                plt.text(10 * kcol + 15, 5, X[kcol],
                         horizontalalignment='center',
                         verticalalignment='center')
            plt.text(10 * kcol + 15, 10 * krow + 15, format(num, '.3f'),
                     horizontalalignment='center',
                     verticalalignment='center')

    plt.axis([0, 10 * (n + 1), 10 * (m + 1), 0])
    plt.title(title)
    plt.xticks(np.linspace(0, 10 * (n + 1), n + 2), [])
    plt.yticks(np.linspace(0, 10 * (m + 1), m + 2), [])
    plt.grid(linestyle="solid")
    # plt.savefig("Table.png", dpi=300)
    plt.show()


original_values = {0: 'A', 1: 'B', 2: 'C', 3: 'D'}

# W1 sequences
val_w1_x1 = ['A', 'A', 'B', 'B', 'C', 'C', 'D', 'D']
val_w1_x2 = ['A', 'B', 'B', 'C', 'B', 'B', 'D', 'D']
val_w1_x3 = ['A', 'C', 'B', 'C', 'B', 'C', 'D']
val_w1_x4 = ['A', 'D']
val_w1_x5 = ['A', 'C', 'B', 'C', 'B', 'A', 'B', 'C', 'D', 'D']
val_w1_x6 = ['B', 'A', 'B', 'A', 'A', 'D', 'D', 'D']
val_w1_x7 = ['B', 'A', 'B', 'C', 'D', 'C', 'C']
val_w1_x8 = ['A', 'B', 'D', 'B', 'B', 'C', 'C', 'D', 'D']
val_w1_x9 = ['A', 'B', 'A', 'A', 'A', 'C', 'D', 'C', 'C', 'D']
val_w1_x10 = ['A', 'B', 'D']

# W2 sequences
val_w2_x1 = ['D', 'D', 'C', 'C', 'B', 'B', 'A', 'A']
val_w2_x2 = ['D', 'D', 'A', 'B', 'C', 'B', 'A']
val_w2_x3 = ['C', 'D', 'C', 'D', 'C', 'B', 'A', 'B', 'A']
val_w2_x4 = ['D', 'D', 'B', 'B', 'A']
val_w2_x5 = ['D', 'A', 'D', 'A', 'C', 'B', 'B', 'A', 'A']
val_w2_x6 = ['C', 'D', 'D', 'C', 'C', 'B', 'A']
val_w2_x7 = ['B', 'D', 'D', 'B', 'C', 'A', 'A', 'A', 'A']
val_w2_x8 = ['B', 'B', 'A', 'B', 'B', 'D', 'D', 'D', 'C', 'D']
val_w2_x9 = ['D', 'D', 'A', 'D', 'D', 'B', 'C', 'A', 'A']
val_w2_x10 = ['D', 'D', 'C', 'A', 'A', 'A']

# Test sequences
T1_sequence = ['A', 'B', 'B', 'B', 'C', 'D', 'D', 'D']
T2_sequence = ['D', 'A', 'D', 'B', 'C', 'B', 'A', 'A']
T3_sequence = ['C', 'D', 'C', 'B', 'A', 'B', 'A']
T4_sequence = ['A', 'D', 'B', 'B', 'B', 'C', 'D']
T5_sequence = ['B', 'A', 'D', 'B', 'D', 'C', 'B', 'A']

# List of test sequences
Test_sequences = [T1_sequence, T2_sequence, T3_sequence, T4_sequence, T5_sequence]

print('--------------------------------------------')
print('--------------------------------------------')
print('--------------------------------------------')
print('-----------------   W1   -------------------')

# Convert W1 sequences
W1_X1 = convert_sequence(val_w1_x1)
W1_X2 = convert_sequence(val_w1_x2)
W1_X3 = convert_sequence(val_w1_x3)
W1_X4 = convert_sequence(val_w1_x4)
W1_X5 = convert_sequence(val_w1_x5)
W1_X6 = convert_sequence(val_w1_x6)
W1_X7 = convert_sequence(val_w1_x7)
W1_X8 = convert_sequence(val_w1_x8)
W1_X9 = convert_sequence(val_w1_x9)
W1_X10 = convert_sequence(val_w1_x10)

X_W1 = np.concatenate([W1_X1, W1_X2, W1_X3, W1_X4, W1_X5, W1_X6, W1_X7, W1_X8, W1_X9, W1_X10])
lengths_W1 = [len(W1_X1), len(W1_X2), len(W1_X3), len(W1_X4), len(W1_X5), len(W1_X6), len(W1_X7), len(W1_X8),
              len(W1_X9),
              len(W1_X10)]

remodelW1 = hmm.GaussianHMM(algorithm='viterbi', n_components=4)
# Train W1
remodelW1.fit(X_W1, lengths_W1)

T1 = convert_sequence(T1_sequence)
T2 = convert_sequence(T2_sequence)
T3 = convert_sequence(T3_sequence)
T4 = convert_sequence(T4_sequence)
T5 = convert_sequence(T5_sequence)

Test_converted_sequences = [T1, T2, T3, T4, T5]

print('Transition matrix:')
print(remodelW1.transmat_)
plot_transition_matrix(remodelW1.transmat_, 'W1 Transition Matrix')

counter_i = 1
for t in Test_converted_sequences:
    print('-------------------------')
    print(f'Evaluating T{counter_i}: ')
    hidden_states = remodelW1.predict(t)
    print('Expected:')
    print(Test_sequences[counter_i - 1])
    print('Hidden state sequence:')
    print([i for i in hidden_states])
    print('Observation sequence:')
    print([original_values[i] for i in hidden_states])
    counter_i = counter_i + 1

print("Means and vars of each hidden state")
for i in range(remodelW1.n_components):
    print("{0}th hidden state".format(i))
    print("mean = ", remodelW1.means_[i])
    print("var = ", np.diag(remodelW1.covars_[i]))
    print()

print('--------------------------------------------')
print('--------------------------------------------')
print('--------------------------------------------')
print('-----------------   W2   -------------------')

# Convert W2 sequences
W2_X1 = convert_sequence(val_w2_x1)
W2_X2 = convert_sequence(val_w2_x2)
W2_X3 = convert_sequence(val_w2_x3)
W2_X4 = convert_sequence(val_w2_x4)
W2_X5 = convert_sequence(val_w2_x5)
W2_X6 = convert_sequence(val_w2_x6)
W2_X7 = convert_sequence(val_w2_x7)
W2_X8 = convert_sequence(val_w2_x8)
W2_X9 = convert_sequence(val_w2_x9)
W2_X10 = convert_sequence(val_w2_x10)

X_W2 = np.concatenate([W2_X1, W2_X2, W2_X3, W2_X4, W2_X5, W2_X6, W2_X7, W2_X8, W2_X9, W2_X10])
lengths_W2 = [len(W2_X1), len(W2_X2), len(W2_X3), len(W2_X4), len(W2_X5), len(W2_X6), len(W2_X7), len(W2_X8),
              len(W2_X9),
              len(W2_X10)]

remodelW2 = hmm.GaussianHMM(algorithm='viterbi', n_components=4)
# Train W2
remodelW2.fit(X_W2, lengths_W2)

print('Transition matrix:')
print(remodelW2.transmat_)
plot_transition_matrix(remodelW2.transmat_, 'W2 Transition Matrix')

i = 1
for t in Test_converted_sequences:
    print('-------------------------')
    print(f'Evaluating T{i}: ')
    predicted = remodelW2.predict(t)
    print('Expected:')
    print(Test_sequences[i - 1])
    print('Hidden state sequence:')
    print([i for i in predicted])
    print('Observation sequence:')
    print([original_values[i] for i in predicted])
    i = i + 1

print("Means and vars of each hidden state")
for i in range(remodelW2.n_components):
    print("{0}th hidden state".format(i))
    print("mean = ", remodelW2.means_[i])
    print("var = ", np.diag(remodelW2.covars_[i]))
    print()