Introduction

A Hidden Markov Model (HMM) is a probabilistic model used to represent systems with observable sequences and latent (hidden) states, such as speech recognition, biological sequence analysis, or financial market prediction. The Baum-Welch algorithm, an Expectation-Maximization (EM) variant, is utilized for training HMMs, enabling the estimation of model parameters when the data is incomplete or has hidden variables.

Explanation of Baum-Welch Algorithm

The Baum-Welch algorithm consists of two primary steps executed iteratively: the Expectation step (E-step) and the Maximization step (M-step). In these steps, we iteratively update the HMM parameters to maximize the likelihood of observed data.

Components of an HMM

An HMM is defined by the following parameters:

• $N$: Number of states in the model.
• $M$: Number of distinct observation symbols per state.
• $\pi$: Initial state distribution vector. $\pi_i$ is the probability of starting in state $i$.
• $A$: State transition probability matrix, where $A_{ij}$ is the probability of moving from state $i$ to state $j$.
• $B$: Observation probability matrix, where $B_{jk}$ is the probability of observing symbol $k$ from state $j$.

Expectation Step (E-step)

In this step, we calculate forward and backward probabilities, which help compute the expected occurrences of transitions and emissions.

Forward Procedure ($\alpha$):

$\alpha_t(i)$ is the probability of observing the sequence $O_1, O_2, \ldots, O_t$ and being in state $i$ at time $t$. The recurrence is:

$$\alpha_t(i) = \left(\sum_{j=1}^{N} \alpha_{t-1}(j) A_{ji}\right) B_i(O_t)$$

with initialization $\alpha_1(i) = \pi_i B_i(O_1)$.

Backward Procedure ($\beta$):

$\beta_t(i)$ is the probability of observing the remaining sequence $O_{t+1}, \ldots, O_T$ given that we are in state $i$ at time $t$. The recurrence is:

$$\beta_t(i) = \sum_{j=1}^{N} A_{ij} B_j(O_{t+1}) \beta_{t+1}(j)$$

with initialization $\beta_T(i) = 1$ for all states.

Derived quantities:

• $\xi_t(i, j)$: Probability of being in state $i$ at time $t$ and state $j$ at time $t+1$:

$$\xi_t(i, j) = \frac{\alpha_t(i) A_{ij} B_j(O_{t+1}) \beta_{t+1}(j)}{P(O|\lambda)}$$

where $P(O|\lambda) = \sum_{i=1}^{N} \alpha_T(i)$.

• $\gamma_t(i)$: Probability of being in state $i$ at time $t$:

$$\gamma_t(i) = \sum_{j=1}^{N} \xi_t(i, j)$$

Maximization Step (M-step)

Using the probabilities from the E-step, we update the model parameters:

Initial state probabilities:

$$\pi_i = \gamma_1(i)$$

Transition probabilities:

$$A_{ij} = \frac{\sum_{t=1}^{T-1} \xi_t(i, j)}{\sum_{t=1}^{T-1} \gamma_t(i)}$$

Observation probabilities:

$$B_j(k) = \frac{\sum_{t=1}^{T} \mathbf{1}(O_t = k) \gamma_t(j)}{\sum_{t=1}^{T} \gamma_t(j)}$$

where $\mathbf{1}(O_t = k)$ is 1 when the observation at time $t$ equals symbol $k$, and 0 otherwise.

Example: Weather and Mood

Consider a simple sequence of observed weather: "Rainy", "Sunny", "Rainy", "Cloudy". We model this with an HMM where the hidden states are "Happy" and "Sad", and the observable symbols are weather conditions.

Initial Parameters

The transition matrix $A$ says: if the person is Happy, there is a 0.7 probability of staying Happy and 0.3 of becoming Sad. The emission matrix $B$ says: if Happy, the probability of Rainy weather is 0.1, Sunny is 0.4, Cloudy is 0.5.

Iteration Process

1. E-step: Compute $\alpha$, $\beta$, $\xi$, and $\gamma$ using the initial $\pi$, $A$, and $B$.
2. M-step: Update $\pi$, $A$, and $B$ based on $\gamma$ and $\xi$ derived from the E-step.
3. Repeat until convergence (when the log-likelihood change between iterations falls below a threshold) or for a specified number of iterations.

Converged Parameters

After several iterations, the parameters might converge to:

The converged model shows that the "Sad" state has a stronger association with rainy weather (0.5) while the "Happy" state is more associated with sunny weather (0.5).

Implementation in Python

1import numpy as np
2
3def baum_welch(observations, N, M, max_iter=100, tol=1e-6):
4    T = len(observations)
5
6    # Random initialization
7    pi = np.random.dirichlet(np.ones(N))
8    A = np.array([np.random.dirichlet(np.ones(N)) for _ in range(N)])
9    B = np.array([np.random.dirichlet(np.ones(M)) for _ in range(N)])
10
11    for iteration in range(max_iter):
12        # E-step: Forward
13        alpha = np.zeros((T, N))
14        alpha[0] = pi * B[:, observations[0]]
15        for t in range(1, T):
16            alpha[t] = (alpha[t-1] @ A) * B[:, observations[t]]
17
18        # E-step: Backward
19        beta = np.zeros((T, N))
20        beta[T-1] = 1
21        for t in range(T-2, -1, -1):
22            beta[t] = A @ (B[:, observations[t+1]] * beta[t+1])
23
24        # Compute xi and gamma
25        xi = np.zeros((T-1, N, N))
26        for t in range(T-1):
27            denom = alpha[t] @ A @ (B[:, observations[t+1]] * beta[t+1])
28            for i in range(N):
29                xi[t, i] = alpha[t, i] * A[i] * B[:, observations[t+1]] * beta[t+1] / denom
30
31        gamma = np.sum(xi, axis=2)  # T-1 x N
32        gamma = np.vstack([gamma, (alpha[T-1] * beta[T-1]) / np.sum(alpha[T-1])])
33
34        # M-step
35        pi_new = gamma[0]
36        A_new = np.sum(xi, axis=0) / np.sum(gamma[:-1], axis=0).reshape(-1, 1)
37        B_new = np.zeros((N, M))
38        for k in range(M):
39            mask = (np.array(observations) == k).astype(float)
40            B_new[:, k] = np.sum(gamma * mask.reshape(-1, 1), axis=0) / np.sum(gamma, axis=0)
41
42        # Check convergence
43        if np.max(np.abs(A_new - A)) < tol:
44            break
45
46        pi, A, B = pi_new, A_new, B_new
47
48    return pi, A, B

Key Points Summary

Summary

The Baum-Welch algorithm is fundamental for training Hidden Markov Models when the state sequence is not directly observed. It works by alternating between computing expected state occupancies and transitions (E-step) and re-estimating model parameters to maximize likelihood (M-step). The algorithm is guaranteed to converge to a local maximum of the likelihood function, though not necessarily the global maximum. In practice, running from multiple random initializations and selecting the best result helps mitigate this limitation.