Hidden Markov Models Viterbi Baum-Welch Sequential Inference
# Hidden Markov Models: Viterbi, Baum-Welch & Sequential Inference
## Introduction & Motivation
Hidden Markov Models (HMMs) are fundamental probabilistic models for sequences where an underlying sequence of hidden (latent) states evolves according to a Markov process, and observed data is generated stochastically from those hidden states. This simple yet powerful framework captures the essence of many real-world phenomena: speech evolves through phonetic states we don't directly observe; protein sequences fold through intermediate configurations; disease progresses through clinical stages inferred from symptoms.
Historical Context: Developed in the 1960s and 1970s, HMMs were the dominant approach in speech recognition and natural language processing before deep learning. They remain widely deployed in production systems and offer interpretable alternatives to neural methods.
Core Motivation: Many practical problems involve inferring hidden structure from noisy observations. An HMM formalizes this by answering three canonical questions:
1. Decoding: Given observations, what is the most likely hidden state sequence? (Viterbi algorithm)
2. Likelihood: What is the probability of observed data? (Forward algorithm)
3. Learning: Given observations, what model parameters best explain the data? (Baum-Welch EM algorithm)
HMMs excel when the underlying process is truly sequential, dependencies are primarily local (Markovian), and interpretability of hidden states matters. Applications span speech recognition, named entity recognition (NER), part-of-speech (POS) tagging, bioinformatics sequence alignment, and anomaly detection in time series.
---
## Core Concepts & Theory
### Markov Assumption & Conditional Independence
A Markov chain models a sequence of states s_1, s_2, \ldots, s_T where the next state depends only on the current state, not the full history:
$$P(s_t | s_{t-1}, s_{t-2}, \ldots, s_1) = P(s_t | s_{t-1})$$
This first-order Markov assumption enables tractable inference via dynamic programming.
### HMM Components
An HMM is defined by:
- Hidden States: S = \{s_1, \ldots, s_N\} (e.g., phonemes in speech)
- Observations: O = \{o_1, \ldots, o_M\} (e.g., acoustic features)
- Initial distribution: \boldsymbol{\pi} = [\pi_i] where \pi_i = P(s_1 = i)
- Transition probabilities: \mathbf{A} = [a_{ij}] where a_{ij} = P(s_t = j | s_{t-1} = i)
- Emission probabilities: \mathbf{B} = [b_i(o_t)] where b_i(o_t) = P(o_t | s_t = i)
Together, (\boldsymbol{\pi}, \mathbf{A}, \mathbf{B}) defines the HMM.
### Generative Process
An HMM generates a sequence (o_1, \ldots, o_T) by:
1. Start: sample s_1 \sim \boldsymbol{\pi}
2. For t = 1, \ldots, T:
- Emit o_t \sim b_{s_t}(o_t)
- Transition: s_{t+1} \sim a_{s_t, \cdot} (if t < T)
---
## Mathematical Formulation
### Forward Algorithm (Computing Likelihood)
The forward probability \alpha_t(i) is the probability of observing o_1, \ldots, o_t and being in state i at time t:
$$\alpha_t(i) = P(o_1, \ldots, o_t, s_t = i)$$
Recursion (dynamic programming):
$$\alpha_t(j) = b_j(o_t) \sum_{i=1}^{N} \alpha_{t-1}(i) a_{ij}$$
Base case: \alpha_1(i) = \pi_i b_i(o_1)
Total likelihood: P(o_1, \ldots, o_T) = \sum_{i=1}^{N} \alpha_T(i)
Complexity: O(N^2 T) instead of O(N^T) (brute force).
### Backward Algorithm
The backward probability \beta_t(i) is the probability of observing o_{t+1}, \ldots, o_T given state i at time t:
$$\beta_t(i) = P(o_{t+1}, \ldots, o_T | s_t = i)$$
Recursion (backward in time):
$$\beta_t(i) = \sum_{j=1}^{N} a_{ij} b_j(o_{t+1}) \beta_{t+1}(j)$$
Base case: \beta_T(i) = 1 (no more observations)
### Viterbi Algorithm (Finding Most Likely State Path)
The Viterbi algorithm finds the hidden state sequence s^* = \arg\max_{s_1, \ldots, s_T} P(s_1, \ldots, s_T | o_1, \ldots, o_T) via dynamic programming:
$$\delta_t(j) = \max_{s_1, \ldots, s_{t-1}} P(s_1, \ldots, s_{t-1}, s_t = j, o_1, \ldots, o_t)$$
Recursion:
$$\delta_t(j) = b_j(o_t) \max_i (\delta_{t-1}(i) a_{ij})$$
Backpointer: \psi_t(j) = \arg\max_i (\delta_{t-1}(i) a_{ij})
Termination: s_T^* = \arg\max_j \delta_T(j)
Backtrack: s_{t}^* = \psi_{t+1}(s_{t+1}^*) for t = T-1, \ldots, 1
Complexity: O(N^2 T), same as forward but single best path.
### Baum-Welch Algorithm (Expectation-Maximization)
Given observations, learn HMM parameters via EM:
E-Step: Compute posterior probabilities:
$$\gamma_t(i) = P(s_t = i | o_1, \ldots, o_T) = \frac{\alpha_t(i) \beta_t(i)}{P(\mathbf{o})}$$
Transition posterior:
$$\xi_t(i, j) = P(s_t = i, s_{t+1} = j | o_1, \ldots, o_T) = \frac{\alpha_t(i) a_{ij} b_j(o_{t+1}) \beta_{t+1}(j)}{P(\mathbf{o})}$$
M-Step: Update parameters:
$$\pi_i = \gamma_1(i)$$
$$a_{ij} = \frac{\sum_{t=1}^{T-1} \xi_t(i, j)}{\sum_{t=1}^{T-1} \gamma_t(i)}$$
$$b_i(o_k) = \frac{\sum_{t: o_t = o_k} \gamma_t(i)}{\sum_{t=1}^{T} \gamma_t(i)}$$
Repeat E and M steps until convergence.
---
## Advanced Theory & Extensions
### Higher-Order Markov Models
Replace first-order assumption with P(s_t | s_{t-1}, s_{t-2}) for longer-range dependencies. Inference complexity increases to O(N^3 T).
### Factorial HMMs
Multiple independent Markov chains generate observations jointly. Decoding requires belief propagation or sampling (intractable exact inference).
### Input-Output HMMs (Discriminative HMMs)
Condition on input: P(s_t | s_{t-1}, \mathbf{x}_t). Learn via CRF-style discriminative objective instead of generative likelihood.
### Hierarchical HMMs
Hierarchical structure with high-level and low-level states. Models multi-scale temporal dependencies.
### Duration Modeling
Standard HMM implicitly models exponential duration distributions (due to Markov assumption). Explicit duration models capture more realistic lengths via semi-Markov HMMs.
---
## Computational Considerations
### Time Complexity
- Forward/Backward: O(N^2 T) (quadratic in states, linear in sequence length)
- Viterbi: O(N^2 T) (same complexity, single path)
- Baum-Welch (single iteration): O(N^2 T) (forward + backward + parameter update)
### Space Complexity
- Forward tables: O(NT) (store \alpha_t(i) for all i, t)
- Backpointers (Viterbi): O(NT) (store \psi_t(i) for traceback)
### Numerical Stability
Direct computation of \alpha_t(i) = b_i(o_t) \sum_{i} \alpha_{t-1}(i) a_{ij} causes underflow for long sequences. Solution: use log-space computation (log-forward, log-Viterbi).
$$\log \alpha_t(j) = \log b_j(o_t) + \log \sum_i \exp(\log \alpha_{t-1}(i) + \log a_{ij})$$
Use log-sum-exp trick: \log(e^a + e^b) = \max(a, b) + \log(1 + e^{-|a-b|})
---
## Practical Implementation Strategies
### Initializing Parameters
Random: Initialize \mathbf{A}, \mathbf{B}, \boldsymbol{\pi} uniformly or with small random noise.
Data-Driven: Initialize b_i(o_k) from empirical observation frequencies per state (from rough manual labeling or clustering).
Warm Start: Use parameters from related task or pre-trained model.
### Choosing Number of Hidden States
- Domain Knowledge: Use linguistic or domain intuition (e.g., phoneme inventory size).
- Cross-Validation: Try different N, evaluate on held-out likelihood.
- Bayesian Model Selection: Compute marginal likelihood via nested sampling or Laplace approximation.
### Handling Sparse Observations
In practice, many observation symbols never occur in training. Smooth emission probabilities via Laplace smoothing (add pseudo-count):
$$b_i(o_k) = \frac{count(s_t = i, o_t = o_k) + \alpha}{count(s_t = i) + \alpha M}$$
where M is vocabulary size, \alpha is smoothing strength.
### Multiple Sequences
Given multiple observation sequences \mathbf{o}^{(1)}, \ldots, \mathbf{o}^{(D)}, Baum-Welch generalizes naturally: accumulate statistics across all sequences before M-step update.
### Decoding with Confidence
Output both Viterbi path and posterior \gamma_t(i) to estimate confidence per state. Low \gamma_t(s_t^*) indicates uncertainty.
---
## Benchmark Datasets & Evaluation
### Synthetic Datasets
- Toy Examples: Small state spaces (N = 2–5), short sequences (T = 10–50). Useful for verifying algorithm correctness.
- Synthetic Data Generator: Sample from known HMM; evaluate parameter recovery.
### Real Benchmark Datasets
- Penn Treebank POS Tagging: ~1M words, 45 POS tags (hidden states). Standard NLP benchmark.
- BioSeq: Protein sequences, secondary structure labels (hidden states).
- Speech Datasets (e.g., TIMIT): Speech waveforms; phoneme labels (hidden states); low-level acoustic features (observations).
- Stock Market: Daily stock prices; market regimes (bull/bear/crash) as hidden states.
### Evaluation Metrics
Supervised (labels available):
- Accuracy: Fraction of correctly predicted hidden states.
- F1-Score (per hidden state): Balance precision and recall.
Unsupervised:
- Likelihood: \log P(\mathbf{o}) on test set. Higher is better.
- Perplexity: 2^{-\frac{1}{T} \log P(\mathbf{o})}. Lower is better.
---
## Key Challenges & Limitations
### Markov Assumption
The first-order assumption P(s_t | s_{t-1}) is often violated. Long-range dependencies require higher-order models (computationally expensive).
### Discrete State Space
HMMs assume discrete hidden states. Continuous state spaces require switching to Kalman filters or particle filters.
### Training Data Requirements
Learning N^2 + NM parameters requires sufficient data. Sparse states or rare observation-state pairs lead to poor estimation. Smoothing is essential.
### Computational Scaling
O(N^2 T) complexity becomes prohibitive for very large N or very long sequences. Solutions: state pruning, approximations, or sampling-based methods.
### Limited Expressive Power
HMMs cannot capture complex hierarchical or deeply sequential patterns (unlike RNNs/LSTMs). For intricate dynamics, neural models often outperform.
---
## Hyperparameter Tuning
### Number of Hidden States N
Cross-validate over N \in \{2, 5, 10, 20, 50\}. Use likelihood on validation set or downstream task performance (e.g., POS tagging accuracy).
### Observation Model
- Discrete: If observations are categorical (e.g., words, phonemes), use categorical emission model b_i(o).
- Continuous: If observations are real-valued (e.g., acoustic features), use Gaussian emission \mathcal{N}(\boldsymbol{\mu}_i, \boldsymbol{\Sigma}_i) or mixture of Gaussians.
### Smoothing Strength
Laplace smoothing parameter \alpha: try \{0.1, 1.0, 10.0\}. Balance between covering unseen data and overfitting.
### Convergence Criteria (Baum-Welch)
Stop when |L_{ ext{new}} - L_{ ext{old}}| < \epsilon (e.g., \epsilon = 10^{-6}) or max iterations reached. Typically 10–100 iterations suffice.
Tuning Strategy: Grid search or random search over N, smoothing, observation model; evaluate on validation likelihood and downstream task.
---
## Real-World Applications & Case Studies
### Speech Recognition
Scenario: Convert audio waveform to text.
HMM Approach: Hidden states represent phonemes. Observations are MFCC features (acoustic features). Train HMM for each phoneme; combine into word/sentence HMMs. Decode using Viterbi.
Outcome: Dominated speech recognition for decades; still baseline in hybrid systems. Modern end-to-end neural models (e.g., attention-based) often outperform but lack interpretability.
### Part-of-Speech Tagging
Scenario: Assign grammatical tags (noun, verb, etc.) to words in a sentence.
HMM Approach: Hidden states = POS tags. Observations = words. Emission model P( ext{word} | ext{tag}); transition model P( ext{tag}_t | ext{tag}_{t-1}) captures grammatical constraints (e.g., adjectives rarely follow verbs).
Outcome: Bilingual tagger achieves ~95% accuracy. Standard baseline for NLP.
### Bioinformatics: Gene Finding
Scenario: Identify genes in DNA sequence.
HMM Approach: Hidden states represent exons, introns, intergenic regions. Observations are nucleotides. Train on known genes; decode unknown sequences.
Outcome: Enables automated genome annotation. Combined with more complex models (phylogenetic HMMs) for multi-sequence alignment.
### Anomaly Detection in Time Series
Scenario: Detect equipment failures or fraud in sensor data.
HMM Approach: Train HMM on normal-behavior data. Low likelihood or unusual state transitions on new data flag anomalies.
Outcome: Interpretable alternative to deep learning for domain-specific domains with limited training data.
---
## Integration with Other Methods
### HMM + CRF (Conditional Random Field)
CRF is discriminative version: condition on observations, learn weights of state and transition features. Often outperforms HMM for tagging tasks.
### HMM + Neural Emission Model
Replace hand-crafted emission model with neural network: b_i(o_t) = ext{Neural}(o_t; heta_i). More expressive; still interpretable via Viterbi.
### HMM + Language Model
In speech recognition, combine acoustic HMM with n-gram language model: P( ext{word}_t | ext{word}_{t-1}) for decoding constraints.
### Hierarchical HMM + Segmentation
Segment sequences first (via change-point detection), then fit separate HMMs. Improves accuracy on multi-regime data.
---
## Future Research Directions
### Neural HMMs & Structured Prediction
Hybrid models combining HMM structure (transparency, interpretability) with neural components (expressive power). Growing interest in "structured" deep learning.
### Continuous State Spaces
Extend to continuous latent dynamics (via Gaussian processes or VAEs). Bridges HMMs and modern generative models.
### Online & Streaming Learning
Incremental Baum-Welch or online EM for data that arrives sequentially without fixed endpoint.
### Approximate Inference
Variational inference, expectation propagation, or particle filters for models where exact Baum-Welch intractable.
### Causality & HMMs
Understanding causal structure in hidden states; intervening on states for counterfactual reasoning.
---
## Summary & Key Takeaways
Hidden Markov Models are foundational probabilistic models for sequential data. Three canonical algorithms solve key tasks: Forward (likelihood), Viterbi (decoding), Baum-Welch (learning). HMMs remain invaluable for interpretable sequence modeling, especially when domain structure is partially known.
Key Principles:
1. HMMs formalize the notion of hidden structure in sequential data via Markovian assumptions.
2. Viterbi algorithm efficiently finds the most likely hidden path using dynamic programming.
3. Baum-Welch EM learns parameters from unlabeled or partially labeled sequences.
4. Numerical stability (log-space computation) essential for long sequences.
5. HMMs are interpretable but limited by Markov assumption; neural methods often outperform for complex dependencies.
6. Hybrid approaches (neural emission models, structured neural networks) combine HMM interpretability with deep learning expressiveness.
HMMs remain essential in production NLP, speech, bioinformatics, and anomaly detection systems.
---
---
## Appendix: Practical Labs
### Lab 1: Forward Algorithm and Likelihood Computation
Implement forward algorithm; verify likelihood computation on toy HMM.
import numpy as np
class SimpleHMM:
def __init__(self, pi, A, B):
"""
Initialize HMM.
Args:
pi: initial state distribution (N,)
A: transition matrix (N, N)
B: emission matrix (N, M) where B[i, k] = P(o=k | s=i)
"""
self.pi = pi
self.A = A
self.B = B
self.N = len(pi) # number of states
def forward(self, observations):
"""
Forward algorithm: compute likelihood P(observations).
Args:
observations: sequence of observation indices (T,)
Returns:
alpha: forward probabilities (T, N)
likelihood: P(observations)
"""
T = len(observations)
N = self.N
alpha = np.zeros((T, N))
# Base case: t=0
alpha[0, :] = self.pi * self.B[:, observations[0]]
# Recursion: t=1..T-1
for t in range(1, T):
for j in range(N):
alpha[t, j] = self.B[j, observations[t]] * np.sum(alpha[t-1, :] * self.A[:, j])
# Total likelihood
likelihood = np.sum(alpha[T-1, :])
return alpha, likelihood
# Toy HMM: 2 states, 2 observations
pi = np.array([0.6, 0.4]) # Initial: more likely to start in state 0
A = np.array([[0.7, 0.3], # From state 0: 70% stay, 30% go to state 1
[0.4, 0.6]]) # From state 1: 40% go to state 0, 60% stay
B = np.array([[0.9, 0.1], # State 0: emit obs 0 with prob 0.9
[0.2, 0.8]]) # State 1: emit obs 1 with prob 0.8
hmm = SimpleHMM(pi, A, B)
# Observation sequence: [0, 0, 1, 1, 0]
observations = np.array([0, 0, 1, 1, 0])
alpha, likelihood = hmm.forward(observations)
print("Forward Algorithm Test:")
print(f"Observations: {observations}")
print(f"Likelihood: {likelihood:.6f}")
print(f"Alpha shape: {alpha.shape}")
# Tests
assert likelihood > 0, "Likelihood must be positive"
assert likelihood <= 1, "Likelihood must be at most 1"
print("✓ Likelihood is valid probability")
# Verify alpha values are non-negative and decreasing likelihood
assert np.all(alpha >= 0), "Alpha values must be non-negative"
print("✓ Alpha values are non-negative")
# Test: alpha final row should sum to likelihood
final_sum = np.sum(alpha[-1, :])
assert np.isclose(final_sum, likelihood), "Final alpha sum should equal likelihood"
print(f"✓ Final alpha sum matches likelihood: {final_sum:.6f}")
if __name__ == "__main__":
print("Lab 1: Forward Algorithm - PASSED")### Lab 2: Viterbi Algorithm and Decoding
Implement Viterbi algorithm; find most likely hidden state sequence.
import numpy as np
class SimpleHMM:
def __init__(self, pi, A, B):
self.pi = pi
self.A = A
self.B = B
self.N = len(pi)
def forward(self, observations):
T = len(observations)
N = self.N
alpha = np.zeros((T, N))
alpha[0, :] = self.pi * self.B[:, observations[0]]
for t in range(1, T):
for j in range(N):
alpha[t, j] = self.B[j, observations[t]] * np.sum(alpha[t-1, :] * self.A[:, j])
return alpha, np.sum(alpha[T-1, :])
def viterbi(self, observations):
"""
Viterbi algorithm: find most likely hidden state sequence.
Args:
observations: sequence of observation indices (T,)
Returns:
best_path: most likely state sequence (T,)
best_prob: probability of best path
"""
T = len(observations)
N = self.N
# DP table: delta[t, i] = max prob of path ending in state i at time t
delta = np.zeros((T, N))
psi = np.zeros((T, N), dtype=int) # Backpointers
# Base case: t=0
delta[0, :] = self.pi * self.B[:, observations[0]]
# Recursion: t=1..T-1
for t in range(1, T):
for j in range(N):
# Find best predecessor state
temp = delta[t-1, :] * self.A[:, j]
psi[t, j] = np.argmax(temp)
delta[t, j] = self.B[j, observations[t]] * np.max(temp)
# Termination: best final state
best_final_state = np.argmax(delta[T-1, :])
best_prob = np.max(delta[T-1, :])
# Backtrack
best_path = np.zeros(T, dtype=int)
best_path[T-1] = best_final_state
for t in range(T-2, -1, -1):
best_path[t] = psi[t+1, best_path[t+1]]
return best_path, best_prob
# Toy HMM
pi = np.array([0.6, 0.4])
A = np.array([[0.7, 0.3], [0.4, 0.6]])
B = np.array([[0.9, 0.1], [0.2, 0.8]])
hmm = SimpleHMM(pi, A, B)
observations = np.array([0, 0, 1, 1, 0])
best_path, best_prob = hmm.viterbi(observations)
print("Viterbi Algorithm Test:")
print(f"Observations: {observations}")
print(f"Best path: {best_path}")
print(f"Best prob: {best_prob:.6f}")
# Tests
assert len(best_path) == len(observations), "Path length mismatch"
print("✓ Path length matches observations")
assert best_prob > 0, "Probability must be positive"
assert best_prob <= 1, "Probability must be at most 1"
print("✓ Probability is valid")
# Compare to forward: Viterbi path prob should be <= forward likelihood
alpha, forward_likelihood = hmm.forward(observations)
assert best_prob <= forward_likelihood, "Viterbi prob should be <= forward likelihood"
print(f"✓ Viterbi prob {best_prob:.6f} <= Forward likelihood {forward_likelihood:.6f}")
if __name__ == "__main__":
print("Lab 2: Viterbi Algorithm - PASSED")### Lab 3: Backward Algorithm and Posterior Computation
Implement backward algorithm; compute posterior state probabilities.
import numpy as np
class SimpleHMM:
def __init__(self, pi, A, B):
self.pi = pi
self.A = A
self.B = B
self.N = len(pi)
def forward(self, observations):
T = len(observations)
N = self.N
alpha = np.zeros((T, N))
alpha[0, :] = self.pi * self.B[:, observations[0]]
for t in range(1, T):
for j in range(N):
alpha[t, j] = self.B[j, observations[t]] * np.sum(alpha[t-1, :] * self.A[:, j])
return alpha, np.sum(alpha[T-1, :])
def backward(self, observations):
"""
Backward algorithm: compute backward probabilities.
Args:
observations: sequence of observation indices (T,)
Returns:
beta: backward probabilities (T, N)
"""
T = len(observations)
N = self.N
beta = np.zeros((T, N))
# Base case: t=T-1
beta[T-1, :] = 1.0
# Recursion: t=T-2..0 (backward in time)
for t in range(T-2, -1, -1):
for i in range(N):
beta[t, i] = np.sum(self.A[i, :] * self.B[:, observations[t+1]] * beta[t+1, :])
return beta
def posterior(self, observations):
"""
Compute posterior state probabilities: P(s_t | observations).
"""
alpha, likelihood = self.forward(observations)
beta = self.backward(observations)
T = len(observations)
gamma = np.zeros((T, self.N))
for t in range(T):
gamma[t, :] = (alpha[t, :] * beta[t, :]) / likelihood
return gamma
# Toy HMM
pi = np.array([0.6, 0.4])
A = np.array([[0.7, 0.3], [0.4, 0.6]])
B = np.array([[0.9, 0.1], [0.2, 0.8]])
hmm = SimpleHMM(pi, A, B)
observations = np.array([0, 0, 1, 1, 0])
gamma = hmm.posterior(observations)
print("Posterior State Probabilities:")
print(f"Gamma shape: {gamma.shape}")
print("Posterior probabilities:")
for t in range(len(observations)):
print(f" t={t}, obs={observations[t]}: P(state=0|obs)={gamma[t, 0]:.3f}, P(state=1|obs)={gamma[t, 1]:.3f}")
# Tests
assert gamma.shape == (len(observations), hmm.N), "Gamma shape mismatch"
print("✓ Gamma shape correct")
# Each row should sum to 1 (valid probability)
row_sums = np.sum(gamma, axis=1)
assert np.allclose(row_sums, 1.0), "Posterior rows must sum to 1"
print("✓ Posterior rows sum to 1 (valid probabilities)")
# All values should be in [0, 1]
assert np.all(gamma >= 0) and np.all(gamma <= 1), "Posteriors must be in [0,1]"
print("✓ All posterior values in [0, 1]")
if __name__ == "__main__":
print("Lab 3: Backward Algorithm and Posterior - PASSED")### Lab 4: Baum-Welch EM Algorithm
Implement Baum-Welch; learn HMM parameters from observations.
import numpy as np
class SimpleHMM:
def __init__(self, pi, A, B):
self.pi = pi.copy()
self.A = A.copy()
self.B = B.copy()
self.N = len(pi)
def forward(self, observations):
T = len(observations)
N = self.N
alpha = np.zeros((T, N))
alpha[0, :] = self.pi * self.B[:, observations[0]]
for t in range(1, T):
for j in range(N):
alpha[t, j] = self.B[j, observations[t]] * np.sum(alpha[t-1, :] * self.A[:, j])
return alpha, np.sum(alpha[T-1, :])
def backward(self, observations):
T = len(observations)
N = self.N
beta = np.zeros((T, N))
beta[T-1, :] = 1.0
for t in range(T-2, -1, -1):
for i in range(N):
beta[t, i] = np.sum(self.A[i, :] * self.B[:, observations[t+1]] * beta[t+1, :])
return beta
def baum_welch(self, observations, max_iter=50, tol=1e-6):
"""
Baum-Welch EM algorithm: learn HMM parameters.
Args:
observations: sequence of observation indices (T,)
max_iter: max iterations
tol: convergence tolerance
Returns:
likelihoods: log-likelihood at each iteration
"""
likelihoods = []
for iteration in range(max_iter):
alpha, likelihood = self.forward(observations)
beta = self.backward(observations)
likelihoods.append(likelihood)
T = len(observations)
N = self.N
M = self.B.shape[1] # Number of observation symbols
# E-step: compute posterior probabilities
gamma = np.zeros((T, N))
for t in range(T):
gamma[t, :] = (alpha[t, :] * beta[t, :]) / likelihood
# M-step: update parameters
# Update initial distribution
self.pi = gamma[0, :]
# Update transition matrix
xi_sum = np.zeros((N, N))
gamma_sum = np.zeros(N)
for t in range(T-1):
for i in range(N):
for j in range(N):
xi = (alpha[t, i] * self.A[i, j] * self.B[j, observations[t+1]] * beta[t+1, j]) / likelihood
xi_sum[i, j] += xi
gamma_sum[i] += gamma[t, i]
for i in range(N):
for j in range(N):
self.A[i, j] = xi_sum[i, j] / gamma_sum[i] if gamma_sum[i] > 0 else 1.0 / N
# Update emission matrix
for i in range(N):
for k in range(M):
numerator = np.sum(gamma[observations == k, i])
denominator = np.sum(gamma[:, i])
self.B[i, k] = numerator / denominator if denominator > 0 else 1.0 / M
# Check convergence
if iteration > 0 and abs(likelihoods[-1] - likelihoods[-2]) < tol:
print(f"Converged at iteration {iteration}")
break
return np.array(likelihoods)
# Initialize random HMM
np.random.seed(42)
pi_init = np.random.dirichlet([1, 1])
A_init = np.random.dirichlet([1, 1], size=2)
B_init = np.random.dirichlet([1, 1], size=2)
hmm = SimpleHMM(pi_init, A_init, B_init)
# Generate observations from true HMM
true_pi = np.array([0.6, 0.4])
true_A = np.array([[0.7, 0.3], [0.4, 0.6]])
true_B = np.array([[0.9, 0.1], [0.2, 0.8]])
true_hmm = SimpleHMM(true_pi, true_A, true_B)
observations = np.array([0, 0, 1, 1, 0, 0, 1, 1, 1, 0])
# Train
likelihoods = hmm.baum_welch(observations, max_iter=50)
print("Baum-Welch EM Algorithm:")
print(f"Initial likelihood: {likelihoods[0]:.6f}")
print(f"Final likelihood: {likelihoods[-1]:.6f}")
print(f"Likelihood improvement: {(likelihoods[-1] - likelihoods[0]):.6f}")
# Tests
assert len(likelihoods) > 0, "No likelihoods computed"
print("✓ Likelihoods computed")
assert np.all(likelihoods > 0), "Likelihoods must be positive"
print("✓ All likelihoods positive")
# Likelihood should be non-decreasing (or very close due to numerical issues)
diffs = np.diff(likelihoods)
assert np.all(diffs >= -1e-6), "Likelihood not non-decreasing"
print("✓ Likelihood is non-decreasing")
# Verify parameters are valid probabilities
assert np.allclose(hmm.pi.sum(), 1.0), "Initial distribution doesn't sum to 1"
assert np.allclose(hmm.A.sum(axis=1), 1.0), "Transition rows don't sum to 1"
assert np.allclose(hmm.B.sum(axis=1), 1.0), "Emission rows don't sum to 1"
print("✓ Learned parameters are valid probabilities")
if __name__ == "__main__":
print("Lab 4: Baum-Welch EM Algorithm - PASSED")