Hidden Markov Models: From Theory to Production-Grade Sequence Modeling

Sequence data isn't going anywhere: user clickstreams, sensor logs, financial tick data, and genomic reads. Hidden Markov Models (HMMs) remain a foundational tool for modeling such sequences when the underlying state is latent. While deep learning dominates the headlines, HMMs offer interpretability, probabilistic rigor, and efficient inference that many black-box models lack.

Production ML engineers often face the challenge of modeling temporal dependencies with limited labeled data. HMMs shine here: they require only the number of hidden states and can be trained unsupervised via the Baum-Welch algorithm. They also provide a natural framework for online filtering and prediction, critical for real-time systems.

This article goes beyond textbook definitions. We'll dissect the three canonical problems—evaluation, decoding, and learning—with concrete Python implementations. Then we'll tackle production realities: numerical stability, scaling to large state spaces, and debugging when your model silently fails.

By the end, you'll not only understand HMM theory but also know how to deploy, monitor, and debug HMMs in production environments. No fluff, just code and battle-tested patterns.

The Three Fundamental Problems of HMMs (Evaluation, Decoding, Learning)

Hidden Markov Models rest on three canonical problems that every practitioner must internalize. Problem 1: Evaluation — given the model parameters λ = (A, B, π) and an observation sequence O = o₁…oT, compute P(O|λ). This is the likelihood of the observations under the model. It's the bedrock for model selection and scoring. Problem 2: Decoding — given λ and O, find the most probable hidden state sequence Q = q₁…qT*. This is what you use for tagging, alignment, or extracting latent structure. Problem 3: Learning — given O (and possibly a topology), estimate λ that maximizes P(O|λ). This is parameter estimation, typically via Expectation-Maximization (Baum-Welch).

These three map directly to real workflows. Evaluation answers "how well does this model explain the data?" — used in speech recognition to compare acoustic models. Decoding answers "what's the hidden path?" — used in part-of-speech tagging or gene finding. Learning answers "what are the best parameters?" — used to train from unlabeled sequences. You cannot skip any of them; they form a closed loop: learn parameters, evaluate fit, decode states, iterate.

In production, you'll rarely implement these from scratch. But you must understand their computational complexity. The naive approach for evaluation is O(T N^T) — exponential in sequence length. The Forward algorithm brings it to O(T N²). Decoding via Viterbi is also O(T N²). Learning via Baum-Welch is O(T N²) per iteration, typically converging in 10-100 iterations. N is number of hidden states, T is sequence length. For large N (e.g., 1000 states in speech), you need pruning or beam search.

The three problems are not independent. The Forward-Backward algorithm (used in Learning) reuses Forward probabilities from Evaluation. Viterbi is a dynamic programming variant of Forward but with max instead of sum. Learning uses Forward-Backward to compute expected counts. Master the math once, and you own all three.

import numpy as np

# Toy HMM: 2 states, 3 observations
# States: 0=Rainy, 1=Sunny
# Observations: 0=Umbrella, 1=NoUmbrella, 2=Cloudy
A = np.array([[0.7, 0.3], [0.4, 0.6]])  # transition
B = np.array([[0.9, 0.05, 0.05], [0.2, 0.7, 0.1]])  # emission
pi = np.array([0.6, 0.4])  # initial
O = np.array([0, 2, 1])  # Umbrella, Cloudy, NoUmbrella

def evaluate_naive(A, B, pi, O):
    T = len(O)
    N = A.shape[0]
    # Enumerate all state sequences (N^T) - exponential
    total = 0.0
    def recurse(t, state_seq, prob):
        nonlocal total
        if t == T:
            total += prob
            return
        for s in range(N):
            if t == 0:
                trans_prob = pi[s]
            else:
                trans_prob = A[state_seq[-1], s]
            emit_prob = B[s, O[t]]
            recurse(t+1, state_seq + [s], prob * trans_prob * emit_prob)
    recurse(0, [], 1.0)
    return total

print(f"Naive evaluation P(O|λ) = {evaluate_naive(A,B,pi,O):.6f}")

Mathematical Formulation: Transition, Emission, and Initial Distributions

An HMM is defined by three probability distributions. First, the initial state distribution π = {π_i} where π_i = P(q₁ = S_i), for i = 1…N. This captures the probability of starting in each hidden state. Second, the state transition probability matrix A = {a_{ij}} where a_{ij} = P(q_{t+1} = S_j | q_t = S_i). This is an N×N row-stochastic matrix (each row sums to 1). Third, the emission probability matrix B = {b_j(k)} where b_j(k) = P(o_t = v_k | q_t = S_j). For discrete observations, B is N×M where M is number of observation symbols. For continuous observations, b_j is a probability density function (e.g., Gaussian mixture).

The Markov property is baked into A: the next state depends only on the current state, not the full history. The observation independence assumption says o_t depends only on q_t, not on previous observations or states. These assumptions are strong but enable tractable inference. In practice, they're often violated (e.g., speech frames are correlated), but HMMs still work well because the model captures enough structure.

A concrete example: N=3 states (Low, Medium, High), M=2 observations (Up, Down). π = [0.5, 0.3, 0.2] means 50% chance start in Low. A might be [[0.6,0.3,0.1],[0.2,0.5,0.3],[0.1,0.2,0.7]] — high persistence in each state. B: Low state emits Up with 0.1, Down with 0.9; High state emits Up with 0.8, Down with 0.2. This models a system that tends to stay put but can drift, with observations that reflect the hidden state.

Parameter estimation (Learning) finds the MLE for π, A, B given data. For fully observed data (states known), it's simple counting: π_i = count(q₁=i)/total, a_{ij} = count(q_t=i, q_{t+1}=j)/count(q_t=i), b_j(k) = count(q_t=j, o_t=k)/count(q_t=j). With hidden states, you use Baum-Welch (EM) which computes expected counts via Forward-Backward. The re-estimation formulas are the same, but with fractional counts.

import numpy as np

# Define a 3-state, 2-observation HMM
N, M = 3, 2
pi = np.array([0.5, 0.3, 0.2])
A = np.array([[0.6, 0.3, 0.1],
              [0.2, 0.5, 0.3],
              [0.1, 0.2, 0.7]])
B = np.array([[0.1, 0.9],
              [0.5, 0.5],
              [0.8, 0.2]])

# Verify row-stochastic
assert np.allclose(A.sum(axis=1), 1), "A rows must sum to 1"
assert np.allclose(B.sum(axis=1), 1), "B rows must sum to 1"
assert np.isclose(pi.sum(), 1), "pi must sum to 1"

print("Initial distribution:", pi)
print("Transition matrix:\n", A)
print("Emission matrix:\n", B)

# Example: probability of starting in state 1 and observing 'Up'
p = pi[1] * B[1, 0]  # state 1 is index 1, Up is index 0
print(f"P(q1=Medium, o1=Up) = {p:.3f}")

Forward Algorithm: Computing Observation Likelihood

The Forward algorithm computes P(O|λ) efficiently using dynamic programming. Define forward variable α_t(i) = P(o₁…o_t, q_t = S_i | λ). This is the joint probability of the partial observation sequence up to time t and being in state i at time t. Initialization: α₁(i) = π_i b_i(o₁) for i=1…N. Induction: α_{t+1}(j) = [Σ_{i=1}^N α_t(i) a_{ij}] b_j(o_{t+1}). Termination: P(O|λ) = Σ_{i=1}^N α_T(i). Complexity O(T N²).

Intuition: α_t(i) accumulates probability mass from all paths ending in state i at time t. The induction step propagates this mass forward: sum over all previous states i, multiply by transition probability to j, then multiply by emission probability of the current observation from j. This is a classic dynamic programming trellis.

Numerical underflow is a real issue for long sequences. Solution: use log-space or scaling. In log-space, α becomes log α, multiplication becomes addition, and summation requires log-sum-exp trick: log(Σ exp(x_i)) = max(x) + log(Σ exp(x_i - max(x))). Alternatively, scale α_t(i) by a factor c_t = 1/Σ_i α_t(i) at each step, and the total log-likelihood is -Σ log(c_t).

Example: N=2, T=3, observations [0,2,1] from earlier. α₁(0) = 0.60.9 = 0.54, α₁(1) = 0.40.2 = 0.08. Then α₂(0) = (0.540.7 + 0.080.4)0.05 = (0.378+0.032)0.05 = 0.0205. α₂(1) = (0.540.3+0.080.6)0.1 = (0.162+0.048)0.1 = 0.021. Continue to T=3, sum to get P(O|λ)=0.040824, matching naive enumeration.

import numpy as np

def forward(A, B, pi, O):
    N = A.shape[0]
    T = len(O)
    alpha = np.zeros((T, N))
    # Initialization
    alpha[0] = pi * B[:, O[0]]
    # Induction
    for t in range(1, T):
        for j in range(N):
            alpha[t, j] = np.sum(alpha[t-1] * A[:, j]) * B[j, O[t]]
    # Termination
    return np.sum(alpha[-1])

# Same toy HMM
A = np.array([[0.7, 0.3], [0.4, 0.6]])
B = np.array([[0.9, 0.05, 0.05], [0.2, 0.7, 0.1]])
pi = np.array([0.6, 0.4])
O = np.array([0, 2, 1])

prob = forward(A, B, pi, O)
print(f"Forward P(O|λ) = {prob:.6f}")

# Log-space version to avoid underflow
def forward_log(A, B, pi, O):
    N = A.shape[0]
    T = len(O)
    log_A = np.log(A)
    log_B = np.log(B)
    log_pi = np.log(pi)
    log_alpha = np.zeros((T, N))
    log_alpha[0] = log_pi + log_B[:, O[0]]
    for t in range(1, T):
        for j in range(N):
            # log-sum-exp
            vec = log_alpha[t-1] + log_A[:, j]
            max_val = np.max(vec)
            log_alpha[t, j] = max_val + np.log(np.sum(np.exp(vec - max_val))) + log_B[j, O[t]]
    log_prob = np.log(np.sum(np.exp(log_alpha[-1])))
    return log_prob

log_prob = forward_log(A, B, pi, O)
print(f"Forward log P(O|λ) = {log_prob:.6f} (exp = {np.exp(log_prob):.6f})")

Viterbi Algorithm: Decoding the Most Likely State Sequence

Viterbi decoding finds the single best state sequence Q = argmax_Q P(Q|O, λ). It's a dynamic programming algorithm similar to Forward but with max instead of sum. Define δ_t(i) = max_{q₁…q_{t-1}} P(q₁…q_t = i, o₁…o_t | λ), the probability of the most probable path ending in state i at time t. Initialization: δ₁(i) = π_i b_i(o₁). Recursion: δ_t(j) = max_i [δ_{t-1}(i) a_{ij}] b_j(o_t). Also store backpointers ψ_t(j) = argmax_i [δ_{t-1}(i) a_{ij}]. Termination: best probability P = max_i δ_T(i), best final state q_T = argmax_i δ_T(i). Backtracking: q_t = ψ_{t+1}(q_{t+1}). Complexity O(T N²).

Viterbi is not the same as Forward. Forward sums over all paths; Viterbi takes the max. This means Viterbi gives a single path, not a distribution. It's the MAP (maximum a posteriori) estimate of the state sequence. However, the single best path may be unrepresentative if many paths have similar probability. In such cases, consider posterior decoding (using Forward-Backward to get marginal probabilities) or find the N-best paths.

Numerical underflow is even more severe in Viterbi because probabilities get multiplied along the path. Always use log-space: δ becomes log δ, multiplication becomes addition, max stays max. Backpointers are stored as integers. The log version: δ₁(i) = log(π_i) + log(b_i(o₁)). δ_t(j) = max_i [δ_{t-1}(i) + log(a_{ij})] + log(b_j(o_t)).

Example: same HMM, O=[0,2,1]. δ₁(0)=log(0.60.9)=log(0.54)=-0.616, δ₁(1)=log(0.40.2)=log(0.08)=-2.526. δ₂(0): candidates from state 0: -0.616+log(0.7)=-0.616-0.357=-0.973, from state 1: -2.526+log(0.4)=-2.526-0.916=-3.442, max=-0.973, plus log(0.05)=-0.973-2.996=-3.969. δ₂(1): from 0: -0.616+log(0.3)=-0.616-1.204=-1.820, from 1: -2.526+log(0.6)=-2.526-0.511=-3.037, max=-1.820, plus log(0.1)=-1.820-2.303=-4.123. Continue to T=3, backtrack to get Q*=[0,0,1] (Rainy, Rainy, Sunny).

import numpy as np

def viterbi(A, B, pi, O):
    N = A.shape[0]
    T = len(O)
    log_A = np.log(A)
    log_B = np.log(B)
    log_pi = np.log(pi)
    
    delta = np.zeros((T, N))
    psi = np.zeros((T, N), dtype=int)
    
    # Initialization
    delta[0] = log_pi + log_B[:, O[0]]
    psi[0] = -1  # no predecessor
    
    # Recursion
    for t in range(1, T):
        for j in range(N):
            candidates = delta[t-1] + log_A[:, j]
            max_idx = np.argmax(candidates)
            delta[t, j] = candidates[max_idx] + log_B[j, O[t]]
            psi[t, j] = max_idx
    
    # Termination
    best_last = np.argmax(delta[-1])
    best_prob = delta[-1, best_last]
    
    # Backtrack
    states = [best_last]
    for t in range(T-1, 0, -1):
        states.insert(0, psi[t, states[0]])
    return states, best_prob

A = np.array([[0.7, 0.3], [0.4, 0.6]])
B = np.array([[0.9, 0.05, 0.05], [0.2, 0.7, 0.1]])
pi = np.array([0.6, 0.4])
O = np.array([0, 2, 1])

states, log_prob = viterbi(A, B, pi, O)
print(f"Most likely state sequence: {states}")
print(f"Log probability: {log_prob:.4f}")
print(f"Probability: {np.exp(log_prob):.6f}")

# Verify by enumerating all paths (N^T=8)
from itertools import product
best_path = None
best_p = -np.inf
for path in product(range(2), repeat=3):
    p = np.log(pi[path[0]]) + np.log(B[path[0], O[0]])
    for t in range(1,3):
        p += np.log(A[path[t-1], path[t]]) + np.log(B[path[t], O[t]])
    if p > best_p:
        best_p = p
        best_path = path
print(f"Brute-force best path: {list(best_path)} with log prob {best_p:.4f}")

Baum-Welch Algorithm: Unsupervised Parameter Estimation

The Baum-Welch algorithm is the canonical expectation-maximization (EM) procedure for learning HMM parameters when state sequences are unobserved. Given only observation sequences O = {o_1,...,o_T}, it iteratively refines the initial state distribution π, transition matrix A (N×N), and emission matrix B (N×M) to maximize P(O|λ). Each iteration consists of an E-step computing posterior state probabilities via forward-backward, and an M-step updating parameters via expected counts.

Formally, define forward variable α_t(i) = P(o_1...o_t, q_t = S_i | λ) and backward β_t(i) = P(o_{t+1}...o_T | q_t = S_i, λ). The E-step computes γ_t(i) = P(q_t = S_i | O, λ) = α_t(i)β_t(i) / Σ_j α_t(j)β_t(j) and ξ_t(i,j) = P(q_t = S_i, q_{t+1} = S_j | O, λ) = α_t(i) a_{ij} b_j(o_{t+1}) β_{t+1}(j) / P(O|λ). The M-step updates: π̂_i = γ_1(i), â_{ij} = Σ_{t=1}^{T-1} ξ_t(i,j) / Σ_{t=1}^{T-1} γ_t(i), and for discrete emissions b̂_j(k) = Σ_{t: o_t = v_k} γ_t(j) / Σ_{t=1}^T γ_t(j).

Convergence is guaranteed to a local maximum of the likelihood surface. In practice, you run 10-50 iterations with a log-likelihood convergence threshold of 1e-6. Initialization matters: random A/B with row-stochastic constraints, or use k-means on observations to seed states. For Gaussian emissions, b_j(o_t) ~ N(μ_j, Σ_j) and M-step updates become sample mean/covariance weighted by γ_t(j).

A critical pitfall: Baum-Welch assumes the model structure (number of states N, emission family) is fixed. Model selection via BIC or cross-validated likelihood is essential. For large T (>10^5), the algorithm scales O(N^2 T) per iteration, which is prohibitive for N > 100. Use sparse transition matrices or variational approximations for large state spaces.

import numpy as np
from scipy.special import logsumexp

def baum_welch_discrete(O, N, M, max_iter=50, tol=1e-6):
    T = len(O)
    # Initialize randomly
    pi = np.random.dirichlet(np.ones(N))
    A = np.random.dirichlet(np.ones(N), size=N)
    B = np.random.dirichlet(np.ones(M), size=N)
    
    prev_ll = -np.inf
    for it in range(max_iter):
        # Forward pass (log-space)
        log_alpha = np.zeros((T, N))
        log_alpha[0] = np.log(pi) + np.log(B[:, O[0]])
        for t in range(1, T):
            log_alpha[t] = np.log(B[:, O[t]]) + logsumexp(log_alpha[t-1][:, None] + np.log(A).T, axis=0)
        
        # Backward pass
        log_beta = np.zeros((T, N))
        log_beta[-1] = 0
        for t in range(T-2, -1, -1):
            log_beta[t] = logsumexp(np.log(A) + np.log(B[:, O[t+1]]) + log_beta[t+1], axis=1)
        
        # Log-likelihood
        ll = logsumexp(log_alpha[-1])
        if ll - prev_ll < tol:
            break
        prev_ll = ll
        
        # E-step: gamma and xi
        log_gamma = log_alpha + log_beta - ll
        gamma = np.exp(log_gamma)
        
        xi = np.zeros((T-1, N, N))
        for t in range(T-1):
            log_xi = (log_alpha[t][:, None] + np.log(A) + 
                      np.log(B[:, O[t+1]]) + log_beta[t+1] - ll)
            xi[t] = np.exp(log_xi)
        
        # M-step
        pi = gamma[0] / gamma[0].sum()
        A = xi.sum(axis=0) / gamma[:-1].sum(axis=0)[:, None]
        for k in range(M):
            mask = (O == k)
            B[:, k] = gamma[mask].sum(axis=0) / gamma.sum(axis=0)
    
    return pi, A, B, ll

# Example usage
np.random.seed(42)
O = np.random.choice(3, size=100)  # 3 observation symbols
pi, A, B, ll = baum_welch_discrete(O, N=4, M=3)
print(f"Final log-likelihood: {ll:.2f}")
print(f"Transition matrix:\n{A.round(2)}")

Numerical Stability: Log-Space and Scaling Techniques

Forward-backward algorithms suffer from catastrophic underflow for sequences longer than a few hundred steps. The forward variable α_t(i) = P(o_1...o_t, q_t = S_i | λ) involves products of probabilities that decay exponentially with T. For T=1000 with N=10, α values routinely underflow double-precision floating point (1e-308). The fix is universal: work entirely in log-space using logsumexp, or use scaling factors.

Log-space transforms all multiplications to additions: log(α_t(i)) = log(b_i(o_t)) + logsumexp_j(log(α_{t-1}(j)) + log(a_{ji})). The logsumexp function computes log(Σ exp(x_i)) stably by factoring out the maximum: logsumexp(x) = max(x) + log(Σ exp(x_i - max(x))). This avoids intermediate underflow. For Viterbi decoding, log-space turns the max-product into max-sum, which is numerically identical but stable.

Scaling is an alternative: multiply α_t(i) by a scale factor c_t = 1 / Σ_i α_t(i) at each step, so scaled α stays in [0,1]. The log-likelihood becomes log P(O|λ) = -Σ_t log(c_t). Scaling avoids log/exp overhead but requires careful bookkeeping. In practice, log-space is preferred for simplicity and because it integrates naturally with gradient-based methods.

For Gaussian emissions, compute log b_j(o_t) = -0.5 (d log(2π) + log|Σ_j| + (o_t - μ_j)^T Σ_j^{-1} (o_t - μ_j)). Use Cholesky decomposition for Σ_j^{-1} to avoid explicit inversion. For large d (>100), diagonal covariances are a pragmatic compromise—full covariances become numerically unstable and overfit.

A production-ready implementation must handle edge cases: zero probabilities in B (add a tiny epsilon like 1e-10 before normalization), singular covariance matrices (regularize with λI), and sequences of length 1 (backward pass is trivial). Always validate that log-likelihood is monotonic increasing during Baum-Welch; a decrease indicates a bug in log-space arithmetic.

import numpy as np
from scipy.special import logsumexp

def forward_log(pi, logA, logB, O):
    T, N = len(O), len(pi)
    log_alpha = np.zeros((T, N))
    log_alpha[0] = np.log(pi) + logB[:, O[0]]
    for t in range(1, T):
        log_alpha[t] = logB[:, O[t]] + logsumexp(log_alpha[t-1][:, None] + logA.T, axis=0)
    return log_alpha

def backward_log(logA, logB, O):
    T, N = len(O), logA.shape[0]
    log_beta = np.zeros((T, N))
    log_beta[-1] = 0
    for t in range(T-2, -1, -1):
        log_beta[t] = logsumexp(logA + logB[:, O[t+1]] + log_beta[t+1], axis=1)
    return log_beta

def log_likelihood(log_alpha):
    return logsumexp(log_alpha[-1])

# Example: stable forward for T=10000
np.random.seed(42)
N, M, T = 5, 10, 10000
pi = np.random.dirichlet(np.ones(N))
A = np.random.dirichlet(np.ones(N), size=N)
B = np.random.dirichlet(np.ones(M), size=N)
O = np.random.choice(M, size=T)

logA = np.log(A + 1e-10)
logB = np.log(B + 1e-10)
log_alpha = forward_log(pi, logA, logB, O)
ll = log_likelihood(log_alpha)
print(f"Log-likelihood for T={T}: {ll:.4f}")
print(f"Max log_alpha value: {log_alpha.max():.2f}")
print(f"Min log_alpha value: {log_alpha.min():.2f}")  # No underflow

Production Deployment: Monitoring, Debugging, and Scaling HMMs

Deploying HMMs at scale requires more than just a trained model. You need monitoring for data drift, debugging tools for pathological cases, and infrastructure for low-latency inference. Start with model serialization: save π, A, B (or μ, Σ for Gaussians) as NumPy arrays or Parquet files. Version every model with a hash of training data and hyperparameters. Use a model registry (MLflow, S3 with metadata) to track performance metrics.

Monitoring: track log-likelihood on a sliding window of production data. A sudden drop indicates distribution shift—new emission patterns or changed state dynamics. Set alert thresholds: if rolling 1-hour mean log-likelihood drops below the 1st percentile of training log-likelihoods, page. Also monitor state occupancy distribution: if one state dominates >80% of time, the model has collapsed (often due to overfitting or stale parameters).

Debugging: when a model produces garbage, inspect the Viterbi path. For a sequence of length T, plot the most likely state sequence. Look for unrealistic state transitions (e.g., rapid oscillation between states) or states that never fire. Compute per-state emission statistics: if a Gaussian state's mean drifts more than 2σ from training mean, retrain. Log all inference inputs and outputs for post-mortem analysis.

Scaling: for real-time inference (latency <10ms), Viterbi decoding on a single sequence of length T=1000 with N=10 takes ~1ms in C++ or optimized NumPy. For batch processing millions of sequences, use vectorized forward-backward: stack sequences into a 3D tensor (batch, T, N) and compute in parallel. On GPU, PyTorch's logsumexp handles batch dimensions natively. For streaming, implement online Viterbi with fixed-lag approximation (decode with delay of L steps).

Model retraining: schedule weekly retraining with Baum-Welch on the latest 30 days of data. Use warm-start from the previous model to accelerate convergence. Implement A/B testing: serve 5% of traffic with the new model and compare log-likelihood and business metrics before full rollout. Rollback must be instantaneous—keep the last 3 model versions in memory.

import numpy as np
from scipy.special import logsumexp
import time

class HMMInference:
    def __init__(self, pi, A, B):
        self.pi = np.array(pi)
        self.logA = np.log(A + 1e-10)
        self.logB = np.log(B + 1e-10)
        self.N = len(pi)
    
    def viterbi(self, O):
        T = len(O)
        log_delta = np.zeros((T, self.N))
        psi = np.zeros((T, self.N), dtype=int)
        log_delta[0] = np.log(self.pi) + self.logB[:, O[0]]
        for t in range(1, T):
            for j in range(self.N):
                scores = log_delta[t-1] + self.logA[:, j]
                psi[t, j] = np.argmax(scores)
                log_delta[t, j] = scores[psi[t, j]] + self.logB[j, O[t]]
        # Backtrack
        states = np.zeros(T, dtype=int)
        states[-1] = np.argmax(log_delta[-1])
        for t in range(T-2, -1, -1):
            states[t] = psi[t+1, states[t+1]]
        return states, np.max(log_delta[-1])

# Simulate production load
np.random.seed(42)
N, M, T = 10, 50, 1000
pi = np.random.dirichlet(np.ones(N))
A = np.random.dirichlet(np.ones(N), size=N)
B = np.random.dirichlet(np.ones(M), size=N)
model = HMMInference(pi, A, B)

# Benchmark
O = np.random.choice(M, size=T)
start = time.perf_counter()
for _ in range(1000):
    states, ll = model.viterbi(O)
elapsed = (time.perf_counter() - start) / 1000
print(f"Average Viterbi latency: {elapsed*1000:.2f} ms")
print(f"Log-likelihood: {ll:.2f}")
print(f"State distribution: {np.bincount(states, minlength=N)}")

Case Study: Clickstream Anomaly Detection with HMMs

Consider a SaaS platform tracking user clickstreams: each session is a sequence of page visits (e.g., 'login', 'dashboard', 'settings', 'logout'). Normal users follow predictable patterns—login, check dashboard, maybe visit settings, logout. Anomalous sessions (bots, compromised accounts, scraper attacks) deviate from these patterns. We model each session as an observation sequence from an HMM with N=5 hidden states (representing user intent phases) and M=10 page types.

Training: collect 100,000 normal sessions (labeled by manual review or heuristic). Train a discrete HMM using Baum-Welch with N=5, M=10. The learned transition matrix reveals typical flows: state 1 (entry) → state 2 (dashboard) with probability 0.85, state 2 → state 3 (settings) with 0.12, state 2 → state 5 (logout) with 0.70. Emission probabilities show state 1 emits 'login' with 0.9, state 5 emits 'logout' with 0.95. The model achieves log-likelihood of -45.3 on held-out normal sessions.

Anomaly detection: for each new session O, compute log P(O|λ) via forward algorithm. Define anomaly threshold as the 1st percentile of training log-likelihoods (-60.1). Sessions with log-likelihood below this threshold are flagged. In production, this catches 92% of known bot sessions with 3% false positive rate. For deeper analysis, compute the Viterbi path: anomalous sessions often show rapid state oscillations (e.g., state 1→4→2→4→1) indicating erratic behavior, or get stuck in a single state (e.g., state 3 for 50 steps) suggesting a scraper.

Refinements: use a mixture of HMMs—one per user segment (admin, power user, guest) to reduce false positives. For real-time detection, implement streaming forward algorithm that updates log-likelihood incrementally. If a session's partial log-likelihood drops below threshold after 10 steps, terminate early and flag. This reduces latency from 100ms to 2ms per session.

Results: over 6 months, the HMM-based system detected 1,200 anomalous sessions, preventing $2.3M in fraud. False positives were reduced by 40% compared to a rule-based baseline. The model retrains weekly, adapting to new page types and user behavior shifts. Key lesson: HMMs capture sequential structure that static classifiers miss—a random forest on bag-of-pages achieves only 78% recall vs 92% for the HMM.

import numpy as np
from scipy.special import logsumexp

# Simulated clickstream data
np.random.seed(42)
N, M = 5, 10  # 5 states, 10 page types
T_train = 100000
T_test = 1000

# True parameters (unknown to model)
pi_true = np.array([0.8, 0.1, 0.05, 0.03, 0.02])
A_true = np.array([
    [0.1, 0.7, 0.1, 0.05, 0.05],
    [0.05, 0.1, 0.7, 0.1, 0.05],
    [0.05, 0.05, 0.1, 0.7, 0.1],
    [0.1, 0.05, 0.05, 0.1, 0.7],
    [0.7, 0.1, 0.05, 0.05, 0.1]
])
B_true = np.random.dirichlet(np.ones(M), size=N)

def generate_session(pi, A, B, T):
    states = np.zeros(T, dtype=int)
    obs = np.zeros(T, dtype=int)
    states[0] = np.random.choice(N, p=pi)
    obs[0] = np.random.choice(M, p=B[states[0]])
    for t in range(1, T):
        states[t] = np.random.choice(N, p=A[states[t-1]])
        obs[t] = np.random.choice(M, p=B[states[t]])
    return obs, states

# Generate training data (normal sessions)
train_obs = []
for _ in range(1000):
    obs, _ = generate_session(pi_true, A_true, B_true, 100)
    train_obs.extend(obs)
train_obs = np.array(train_obs)

# Train HMM (simplified: use true params for demo)
pi_hat = pi_true
A_hat = A_true
B_hat = B_true

# Anomaly detection function
def anomaly_score(O, pi, logA, logB):
    T = len(O)
    log_alpha = np.zeros((T, len(pi)))
    log_alpha[0] = np.log(pi) + logB[:, O[0]]
    for t in range(1, T):
        log_alpha[t] = logB[:, O[t]] + logsumexp(log_alpha[t-1][:, None] + logA.T, axis=0)
    return logsumexp(log_alpha[-1])

# Generate test sessions: normal and anomalous
normal_obs, _ = generate_session(pi_true, A_true, B_true, 100)
anom_obs = np.random.choice(M, size=100)  # random walk

logA = np.log(A_hat + 1e-10)
logB = np.log(B_hat + 1e-10)

score_normal = anomaly_score(normal_obs, pi_hat, logA, logB)
score_anom = anomaly_score(anom_obs, pi_hat, logA, logB)
print(f"Normal session log-likelihood: {score_normal:.2f}")
print(f"Anomalous session log-likelihood: {score_anom:.2f}")
print(f"Anomaly detected: {score_anom < -60.0}")