xbe.at

Statistics and maths / Time series and stochastic

Viterbi Algorithm for Hidden Markov Models in Python

10 min read, 15 sections

Contents
  1. Introduction to Viterbi Algorithm for Hidden Markov Models
  2. Hidden Markov Models: The Foundation
  3. Viterbi Algorithm: Core Concept
  4. Initialization Step
  5. Recursion Step
  6. Backtracking Step
  7. Handling Numerical Underflow
  8. Real-Life Example: Part-of-Speech Tagging
  9. Real-Life Example: Speech Recognition
  10. Implementing the Forward Algorithm
  11. Comparing Viterbi and Forward Algorithms
  12. Implementing the Backward Algorithm
  13. Baum-Welch Algorithm: HMM Parameter Estimation
  14. Scaling in Forward-Backward Algorithms
  15. Additional Resources

Introduction to Viterbi Algorithm for Hidden Markov Models

The Viterbi algorithm is a dynamic programming approach used to find the most likely sequence of hidden states in a Hidden Markov Model (HMM). It's widely used in various fields, including speech recognition, natural language processing, and bioinformatics. This slideshow will explore the algorithm's concepts, implementation, and applications using Python.

import numpy as np

class HiddenMarkovModel:
    def __init__(self, states, observations, start_prob, trans_prob, emit_prob):
        self.states = states
        self.observations = observations
        self.start_prob = start_prob
        self.trans_prob = trans_prob
        self.emit_prob = emit_prob

# Example initialization
states = ['Rainy', 'Sunny']
observations = ['Walk', 'Shop', 'Clean']
start_prob = {'Rainy': 0.6, 'Sunny': 0.4}
trans_prob = {
    'Rainy': {'Rainy': 0.7, 'Sunny': 0.3},
    'Sunny': {'Rainy': 0.4, 'Sunny': 0.6}
}
emit_prob = {
    'Rainy': {'Walk': 0.1, 'Shop': 0.4, 'Clean': 0.5},
    'Sunny': {'Walk': 0.6, 'Shop': 0.3, 'Clean': 0.1}
}

hmm = HiddenMarkovModel(states, observations, start_prob, trans_prob, emit_prob)

Hidden Markov Models: The Foundation

Hidden Markov Models (HMMs) are statistical models where the system being modeled is assumed to be a Markov process with unobserved (hidden) states. The Viterbi algorithm operates on these models to decode the most probable sequence of hidden states given a sequence of observations.

def generate_sequence(hmm, length):
    current_state = np.random.choice(hmm.states, p=list(hmm.start_prob.values()))
    sequence = [current_state]
    observations = []
    
    for _ in range(length - 1):
        observation = np.random.choice(hmm.observations, p=list(hmm.emit_prob[current_state].values()))
        observations.append(observation)
        current_state = np.random.choice(hmm.states, p=list(hmm.trans_prob[current_state].values()))
        sequence.append(current_state)
    
    observation = np.random.choice(hmm.observations, p=list(hmm.emit_prob[current_state].values()))
    observations.append(observation)
    
    return sequence, observations

# Generate a sequence
hidden_sequence, observed_sequence = generate_sequence(hmm, 5)
print(f"Hidden sequence: {hidden_sequence}")
print(f"Observed sequence: {observed_sequence}")

Viterbi Algorithm: Core Concept

The Viterbi algorithm finds the most likely sequence of hidden states by using dynamic programming. It calculates the probability of being in each state at each time step, considering both the transition probabilities between states and the emission probabilities of observations.

def viterbi(hmm, observations):
    V = [{}]
    path = {}

    # Initialize base cases (t == 0)
    for state in hmm.states:
        V[0][state] = hmm.start_prob[state] * hmm.emit_prob[state][observations[0]]
        path[state] = [state]

    # Run Viterbi for t > 0
    for t in range(1, len(observations)):
        V.append({})
        newpath = {}

        for state in hmm.states:
            (prob, prev_state) = max((V[t-1][prev_state] * hmm.trans_prob[prev_state][state] * hmm.emit_prob[state][observations[t]], prev_state) for prev_state in hmm.states)
            V[t][state] = prob
            newpath[state] = path[prev_state] + [state]

        path = newpath

    # Find the most likely final state
    (prob, state) = max((V[len(observations) - 1][state], state) for state in hmm.states)
    return prob, path[state]

# Example usage
observations = ['Walk', 'Shop', 'Clean']
prob, state_sequence = viterbi(hmm, observations)
print(f"Most likely state sequence: {state_sequence}")
print(f"Probability: {prob}")

Initialization Step

The Viterbi algorithm begins by initializing the probability of being in each state at the first time step. This is done by multiplying the start probability of each state with the emission probability of the first observation.

def initialize_viterbi(hmm, first_observation):
    V = {}
    for state in hmm.states:
        V[state] = hmm.start_prob[state] * hmm.emit_prob[state][first_observation]
    return V

# Example usage
first_observation = 'Walk'
initial_probabilities = initialize_viterbi(hmm, first_observation)
print(f"Initial probabilities: {initial_probabilities}")

Recursion Step

The core of the Viterbi algorithm lies in its recursive step. For each subsequent observation, it calculates the probability of being in each state by considering all possible previous states, the transition probabilities, and the emission probabilities.

def viterbi_step(hmm, V_prev, observation):
    V_current = {}
    backpointer = {}
    for current_state in hmm.states:
        max_prob, best_prev_state = max(
            (V_prev[prev_state] * hmm.trans_prob[prev_state][current_state] * hmm.emit_prob[current_state][observation], prev_state)
            for prev_state in hmm.states
        )
        V_current[current_state] = max_prob
        backpointer[current_state] = best_prev_state
    return V_current, backpointer

# Example usage
observations = ['Walk', 'Shop']
V_prev = initialize_viterbi(hmm, observations[0])
V_current, backpointer = viterbi_step(hmm, V_prev, observations[1])
print(f"Current probabilities: {V_current}")
print(f"Backpointers: {backpointer}")

Backtracking Step

After computing probabilities for all observations, the Viterbi algorithm performs a backtracking step to determine the most likely sequence of hidden states. It starts from the most probable final state and follows the backpointers to construct the optimal path.

def viterbi_backtrack(V, backpointers):
    last_step = max(V[-1], key=V[-1].get)
    path = [last_step]
    
    for t in range(len(V) - 2, -1, -1):
        path.insert(0, backpointers[t+1][path[0]])
    
    return path

# Example usage
V = [
    {'Rainy': 0.06, 'Sunny': 0.24},
    {'Rainy': 0.0384, 'Sunny': 0.0432},
    {'Rainy': 0.01344, 'Sunny': 0.002592}
]
backpointers = [
    None,
    {'Rainy': 'Sunny', 'Sunny': 'Sunny'},
    {'Rainy': 'Rainy', 'Sunny': 'Rainy'}
]

path = viterbi_backtrack(V, backpointers)
print(f"Most likely state sequence: {path}")

Handling Numerical Underflow

When dealing with long observation sequences, the probability values can become extremely small, leading to numerical underflow. To mitigate this, we can work in the log space, converting multiplications to additions.

import math

def viterbi_log(hmm, observations):
    V = [{}]
    path = {}

    # Initialize base cases (t == 0)
    for state in hmm.states:
        V[0][state] = math.log(hmm.start_prob[state]) + math.log(hmm.emit_prob[state][observations[0]])
        path[state] = [state]

    # Run Viterbi for t > 0
    for t in range(1, len(observations)):
        V.append({})
        newpath = {}

        for state in hmm.states:
            (log_prob, prev_state) = max(
                (V[t-1][prev_state] + math.log(hmm.trans_prob[prev_state][state]) + math.log(hmm.emit_prob[state][observations[t]]), prev_state)
                for prev_state in hmm.states
            )
            V[t][state] = log_prob
            newpath[state] = path[prev_state] + [state]

        path = newpath

    # Find the most likely final state
    (log_prob, state) = max((V[len(observations) - 1][state], state) for state in hmm.states)
    return math.exp(log_prob), path[state]

# Example usage
observations = ['Walk', 'Shop', 'Clean', 'Walk', 'Shop']
prob, state_sequence = viterbi_log(hmm, observations)
print(f"Most likely state sequence: {state_sequence}")
print(f"Probability: {prob}")

Real-Life Example: Part-of-Speech Tagging

One common application of the Viterbi algorithm is in natural language processing for part-of-speech tagging. Given a sequence of words, we can use the Viterbi algorithm to determine the most likely sequence of parts of speech.

# Simplified POS tagging HMM
pos_states = ['Noun', 'Verb', 'Adjective']
pos_observations = ['The', 'cat', 'catches', 'the', 'mouse']
pos_start_prob = {'Noun': 0.3, 'Verb': 0.3, 'Adjective': 0.4}
pos_trans_prob = {
    'Noun': {'Noun': 0.3, 'Verb': 0.6, 'Adjective': 0.1},
    'Verb': {'Noun': 0.5, 'Verb': 0.1, 'Adjective': 0.4},
    'Adjective': {'Noun': 0.8, 'Verb': 0.1, 'Adjective': 0.1}
}
pos_emit_prob = {
    'Noun': {'The': 0.4, 'cat': 0.4, 'catches': 0.1, 'the': 0.1, 'mouse': 0.4},
    'Verb': {'The': 0.01, 'cat': 0.01, 'catches': 0.97, 'the': 0.01, 'mouse': 0.01},
    'Adjective': {'The': 0.5, 'cat': 0.1, 'catches': 0.1, 'the': 0.5, 'mouse': 0.1}
}

pos_hmm = HiddenMarkovModel(pos_states, pos_observations, pos_start_prob, pos_trans_prob, pos_emit_prob)

sentence = ['The', 'cat', 'catches', 'the', 'mouse']
prob, pos_sequence = viterbi(pos_hmm, sentence)
print(f"Sentence: {' '.join(sentence)}")
print(f"Most likely POS sequence: {pos_sequence}")

Real-Life Example: Speech Recognition

Another application of the Viterbi algorithm is in speech recognition. In this simplified example, we'll use the algorithm to determine the most likely sequence of phonemes given a series of acoustic observations.

# Simplified speech recognition HMM
phoneme_states = ['S', 'IH', 'L', 'EH', 'N', 'T']
acoustic_observations = ['high_freq', 'low_freq', 'medium_freq', 'low_freq', 'high_freq', 'medium_freq']
phoneme_start_prob = {state: 1/len(phoneme_states) for state in phoneme_states}
phoneme_trans_prob = {
    state: {next_state: 0.1 if state != next_state else 0.5 for next_state in phoneme_states}
    for state in phoneme_states
}
phoneme_emit_prob = {
    'S': {'high_freq': 0.6, 'medium_freq': 0.3, 'low_freq': 0.1},
    'IH': {'high_freq': 0.1, 'medium_freq': 0.7, 'low_freq': 0.2},
    'L': {'high_freq': 0.1, 'medium_freq': 0.2, 'low_freq': 0.7},
    'EH': {'high_freq': 0.3, 'medium_freq': 0.6, 'low_freq': 0.1},
    'N': {'high_freq': 0.2, 'medium_freq': 0.3, 'low_freq': 0.5},
    'T': {'high_freq': 0.7, 'medium_freq': 0.2, 'low_freq': 0.1}
}

speech_hmm = HiddenMarkovModel(phoneme_states, acoustic_observations, phoneme_start_prob, phoneme_trans_prob, phoneme_emit_prob)

acoustic_sequence = ['high_freq', 'low_freq', 'medium_freq', 'low_freq', 'high_freq', 'medium_freq']
prob, phoneme_sequence = viterbi(speech_hmm, acoustic_sequence)
print(f"Acoustic observations: {acoustic_sequence}")
print(f"Most likely phoneme sequence: {phoneme_sequence}")

Implementing the Forward Algorithm

The Forward algorithm is closely related to the Viterbi algorithm and calculates the probability of an observation sequence given the HMM. It can be used to evaluate how well a model fits the observed data.

def forward(hmm, observations):
    F = [{} for _ in range(len(observations))]
    
    # Initialize base cases (t == 0)
    for state in hmm.states:
        F[0][state] = hmm.start_prob[state] * hmm.emit_prob[state][observations[0]]
    
    # Run Forward algorithm for t > 0
    for t in range(1, len(observations)):
        for state in hmm.states:
            F[t][state] = sum(F[t-1][prev_state] * hmm.trans_prob[prev_state][state] * hmm.emit_prob[state][observations[t]]
                              for prev_state in hmm.states)
    
    # Final probability is the sum of the terminal probabilities
    return sum(F[-1].values())

# Example usage
observations = ['Walk', 'Shop', 'Clean']
probability = forward(hmm, observations)
print(f"Probability of the observation sequence: {probability}")

Comparing Viterbi and Forward Algorithms

While both the Viterbi and Forward algorithms work with HMMs, they serve different purposes. Viterbi finds the most likely state sequence, while Forward calculates the overall probability of the observation sequence.

def compare_viterbi_forward(hmm, observations):
    viterbi_prob, viterbi_path = viterbi(hmm, observations)
    forward_prob = forward(hmm, observations)
    
    print(f"Observations: {observations}")
    print(f"Viterbi - Most likely state sequence: {viterbi_path}")
    print(f"Viterbi - Probability of the most likely path: {viterbi_prob}")
    print(f"Forward - Probability of the observation sequence: {forward_prob}")
    
    if viterbi_prob > forward_prob:
        print("Viterbi probability is higher (as expected)")
    else:
        print("Unexpected result: Forward probability is higher")

# Example usage
observations = ['Walk', 'Shop', 'Clean']
compare_viterbi_forward(hmm, observations)

Implementing the Backward Algorithm

The Backward algorithm is another important algorithm for HMMs. It calculates the probability of the ending partial sequence given a starting state. When combined with the Forward algorithm, it can be used for parameter estimation and smoothing.

def backward(hmm, observations):
    B = [{} for _ in range(len(observations))]
    
    # Initialize base cases (t == len(observations)-1)
    for state in hmm.states:
        B[-1][state] = 1.0
    
    # Run Backward algorithm for t < len(observations)-1
    for t in range(len(observations) - 2, -1, -1):
        for state in hmm.states:
            B[t][state] = sum(hmm.trans_prob[state][next_state] * 
                              hmm.emit_prob[next_state][observations[t+1]] * 
                              B[t+1][next_state]
                              for next_state in hmm.states)
    
    # Initial probability
    initial_prob = sum(hmm.start_prob[state] * 
                       hmm.emit_prob[state][observations[0]] * 
                       B[0][state] for state in hmm.states)
    
    return initial_prob, B

# Example usage
observations = ['Walk', 'Shop', 'Clean']
initial_prob, backward_probs = backward(hmm, observations)
print(f"Backward probability: {initial_prob}")

Baum-Welch Algorithm: HMM Parameter Estimation

The Baum-Welch algorithm is an iterative method used to estimate the parameters of an HMM given a set of observations. It uses the Forward-Backward algorithm as a subroutine and is a special case of the Expectation-Maximization (EM) algorithm.

def baum_welch(observations, states, n_iter=100):
    # Initialize HMM parameters randomly
    n_states = len(states)
    n_observations = len(set(observations))
    
    start_prob = np.random.dirichlet(np.ones(n_states))
    trans_prob = np.random.dirichlet(np.ones(n_states), size=n_states)
    emit_prob = np.random.dirichlet(np.ones(n_observations), size=n_states)
    
    for _ in range(n_iter):
        # E-step: Compute forward and backward probabilities
        alpha = forward(observations, start_prob, trans_prob, emit_prob)
        beta = backward(observations, start_prob, trans_prob, emit_prob)
        
        # M-step: Update parameters
        start_prob, trans_prob, emit_prob = update_parameters(observations, alpha, beta, trans_prob, emit_prob)
    
    return start_prob, trans_prob, emit_prob

# Note: The forward, backward, and update_parameters functions are not implemented here for brevity
# In a real implementation, these functions would need to be defined

Scaling in Forward-Backward Algorithms

When dealing with long observation sequences, both the Forward and Backward algorithms can suffer from numerical underflow. Scaling is a technique used to prevent this issue by normalizing the probabilities at each step.

def scaled_forward(hmm, observations):
    F = [{} for _ in range(len(observations))]
    scale = [0] * len(observations)
    
    # Initialize base cases (t == 0)
    for state in hmm.states:
        F[0][state] = hmm.start_prob[state] * hmm.emit_prob[state][observations[0]]
    scale[0] = sum(F[0].values())
    for state in hmm.states:
        F[0][state] /= scale[0]
    
    # Run Forward algorithm for t > 0
    for t in range(1, len(observations)):
        for state in hmm.states:
            F[t][state] = sum(F[t-1][prev_state] * hmm.trans_prob[prev_state][state] * hmm.emit_prob[state][observations[t]]
                              for prev_state in hmm.states)
        scale[t] = sum(F[t].values())
        for state in hmm.states:
            F[t][state] /= scale[t]
    
    return F, scale

# Example usage
observations = ['Walk', 'Shop', 'Clean', 'Walk', 'Shop']
scaled_forward_probs, scaling_factors = scaled_forward(hmm, observations)
print(f"Scaling factors: {scaling_factors}")

Additional Resources

For those interested in diving deeper into the Viterbi algorithm and Hidden Markov Models, here are some valuable resources:

  1. "A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition" by Lawrence R. Rabiner (1989) ArXiv: https://arxiv.org/abs/cs/0702084
  2. "An Introduction to Hidden Markov Models and Bayesian Networks" by Zoubin Ghahramani (2001) ArXiv: https://arxiv.org/abs/1301.6725
  3. "The Viterbi Algorithm" by G. David Forney Jr. (1973) While not available on ArXiv, this seminal paper can be found in many digital libraries.

These resources provide in-depth explanations of the theoretical foundations and practical applications of HMMs and the Viterbi algorithm in various fields of study.

More in Time series and stochastic

  1. Understanding Data Distribution in Python5 min

    Data distribution refers to the way data is spread or organized within a dataset. Understanding data distribution is crucial for effective data analysis...

  2. Tree Attention Topology-aware Decoding for Long-Context on GPUs8 min

    Tree Attention is a novel approach to handling long-context attention in large language models. It addresses the quadratic complexity problem of traditional...

  3. Mathematical Definitions in Data Science22 min

    Gradient Descent is an optimization algorithm used to minimize a function by iteratively moving in the direction of steepest descent. In data science, it's...

  4. Forward Algorithm for Hidden Markov Models in Python7 min

    The Forward Algorithm is a crucial component in Hidden Markov Models (HMMs), used to efficiently calculate the probability of observing a sequence of events....

  5. Estimating Portfolio Losses with Python's Monte Carlo Simulation7 min

    Monte Carlo simulation is a powerful statistical technique used to model complex systems and estimate probabilities. It's named after the famous casino in...

  6. Eliminating Skewness in Data Using Log Transform in Python8 min

    Log transformation is a powerful technique used to handle skewed data distributions. It can help normalize data, reduce the impact of outliers, and make...