Hidden Markov Model in Python

I am going to implement a hidden markov model (HMM) in this tutorial, this model can be used to predict something based on evidence in the current state and the previous state. Hidden markov models is probabilty networks of observable states, hidden states and transitions between hidden states. A transition can only be from one time period to another (markov process).

HMM:s can be used in partially observable environments, in which an agent only have a limited knowledge about the world. When the environment is partially observable, an agent can at best predict how the world will evolve in the next time step. An agent uses a sensor model to perceive the world and a transition model to predict the next state in the world.

Hidden markov models assumes that the future only is dependent on observations in current time and transitions from a previous time, this means that the model is fast and that it does not need to store a lot of historical information. A process that only depends on observations in current time and transitions from a previous time is called a markov process.

I using the viterbi algorithm to make predictions in this tutorial. The viterbi algorithm finds the most likely sequence of predictions given an HMM and a sequence of observations.

Problem and libraries

A security guard is stationed in an underground installation and tries to guess if it is raining based on observing if the director comes to work with an umbrella or not. I am using the following libraries in this tutorial: numpy, pandas and pprint.

Hidden Markov Model

The following code is used to model the problem with probability matrixes. A probability matrix is created for umbrella observations and the weather, another probability matrix is created for the weather on day 0 and the weather on day 1 (transitions between hidden states). The output from a run is shown below the code.

# Import libraries
import numpy as np
import pandas as pd
import pprint

# Get markov edges
def get_markov_edges(df):

    # Create a dictionary
    edges = {}

    # Loop columns
    for column in df.columns:
        # Loop rows
        for row in df.index:
            edges[(row,column)] = df.loc[row,column]

    # Return edges
    return edges

# Viterbi algorithm for shortest path
# https://github.com/alexsosn/MarslandMLAlgo/blob/master/Ch16/HMM.py
def viterbi(pi, a, b, obs):
    nStates = np.shape(b)[0]
    T = np.shape(obs)[0]
    path = np.zeros(T)
    delta = np.zeros((nStates, T))
    phi = np.zeros((nStates, T))
    delta[:, 0] = pi * b[:, obs[0]]
    phi[:, 0] = 0

    for t in range(1, T):
        for s in range(nStates):
            delta[s, t] = np.max(delta[:, t-1] * a[:, s]) * b[s, obs[t]] 
            phi[s, t] = np.argmax(delta[:, t-1] * a[:, s])
    path[T-1] = np.argmax(delta[:, T-1])
    for t in range(T-2, -1, -1):
        path[t] = phi[int(path[t+1]),int(t+1)]
    return path, delta, phi

# The main entry point for this module
def main():

    # Observation states
    # The director can have an umbrella or not have an umbrella (equally likely)
    observation_states = ['Umbrella', 'No umbrella']

    # Create hidden states with probabilities (250 rainy days per year)
    p = [0.32, 0.68]
    #p = [0.5, 0.5]
    #p = [0.7, 0.3]
    hidden_states = ['Sunny', 'Rainy']
    state_space = pd.Series(p, index=hidden_states, name='states')

    # Print hidden states
    print('--- Hidden states ---')

    # Create a hidden states transition matrix with probabilities
    hidden_df = pd.DataFrame(columns=hidden_states, index=hidden_states)
    hidden_df.loc[hidden_states[0]] = [0.75, 0.25]
    hidden_df.loc[hidden_states[1]] = [0.25, 0.75]

    # Print transition matrix
    print('--- Transition matrix for hidden states ---')

    # Create matrix of observations with sensor probabilities
    observations_df = pd.DataFrame(columns=observation_states, index=hidden_states)
    observations_df.loc[hidden_states[0]] = [0.1, 0.9]
    observations_df.loc[hidden_states[1]] = [0.8, 0.2]

    # Print observation matrix
    print('--- Sensor matrix ---')

    # Create graph edges and weights
    hidden_edges = get_markov_edges(hidden_df)
    observation_edges = get_markov_edges(observations_df)

    # Print edges
    print('--- Hidden edges ---')
    print('--- Sensor edges ---')

    # Observations
    observations_map = {0:'Umbrella', 1:'No umbrella'}
    observations = np.array([1,1,1,0,1,1,1,0,0,0])
    observerations_path = [observations_map[v] for v in list(observations)]

    # Get predictions with the viterbi algorithm
    path, delta, phi = viterbi(p, hidden_df.values, observations_df.values, observations)
    state_map = {0:'Sunny', 1:'Rainy'}
    state_path = [state_map[v] for v in path]
    state_delta = np.amax(delta, axis=0)

    # Print predictions
    print('--- Predictions ---')

# Tell python to run main method
if __name__ == "__main__": main()
--- Hidden states ---
Sunny    0.32
Rainy    0.68
Name: states, dtype: float64

--- Transition matrix for hidden states ---
      Sunny Rainy
Sunny  0.75  0.25
Rainy  0.25  0.75

Sunny    1.0
Rainy    1.0
dtype: float64

--- Sensor matrix ---
      Umbrella No umbrella
Sunny      0.1         0.9
Rainy      0.8         0.2

Sunny    1.0
Rainy    1.0
dtype: float64

--- Hidden edges ---
{('Rainy', 'Rainy'): 0.75,
 ('Rainy', 'Sunny'): 0.25,
 ('Sunny', 'Rainy'): 0.25,
 ('Sunny', 'Sunny'): 0.75}

--- Sensor edges ---
{('Rainy', 'No umbrella'): 0.2,
 ('Rainy', 'Umbrella'): 0.8,
 ('Sunny', 'No umbrella'): 0.9,
 ('Sunny', 'Umbrella'): 0.1}

--- Predictions ---
   Observation Prediction     Delta
0  No umbrella      Sunny  0.288000
1  No umbrella      Sunny  0.194400
2  No umbrella      Sunny  0.131220
3     Umbrella      Sunny  0.026244
4  No umbrella      Sunny  0.006643
5  No umbrella      Sunny  0.004484
6  No umbrella      Sunny  0.003027
7     Umbrella      Rainy  0.000605
8     Umbrella      Rainy  0.000363
9     Umbrella      Rainy  0.000218

Leave a Reply

Your email address will not be published. Required fields are marked *