# Fit Hidden Markov Model

> Fit hidden Markov models using the Baum-Welch (EM) algorithm with model selection, Viterbi decoding for state sequences, and forward-backward probabilities. Use when observations are generated by unobservable latent states, you need to segment a time series into latent regimes (market regimes, speech phonemes, biological sequences), compute sequence probabilities, decode the most likely hidden state path, or compare models with different numbers of hidden states.

- Skill: `pjt222/fit-hidden-markov-model-10` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pjt222/fit-hidden-markov-model-10`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pjt222/fit-hidden-markov-model-10/raw
- Safety review: pending (external: skill-scanner PASS, skillspector PASS)
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: AI & ML
- License: MIT
- Author: pjt222 (https://skillmd.com/u/pjt222)
- Updated: 2026-09-21
- Page: https://skillmd.com/skills/pjt222/fit-hidden-markov-model-10

---


# Fit Hidden Markov Model

Fit HMM via Baum-Welch EM, decode most likely hidden state sequence via Viterbi, select optimal N hidden states via information criteria.

## Use When

- Observe sequence emissions but underlying generative states not observable
- Data generated by system switching between finite regimes
- Segment time series into latent phases (market regimes, speech phonemes, biological annotation)
- Compute prob of observed sequence under generative model
- Most likely sequence hidden states given observations (decoding)
- Compare models w/ diff N hidden states → complexity-fit trade-off

## In

### Required

| Input | Type | Desc |
|---|---|---|
| `observations` | sequence/matrix | Observed data (univariate/multivariate) |
| `n_hidden_states` | integer | N hidden states (or range for selection) |
| `emission_type` | string | `"gaussian"`, `"discrete"`, `"poisson"`, `"multinomial"` |

### Optional

| Input | Type | Default | Desc |
|---|---|---|---|
| `initial_params` | dict | random/heuristic | Init transition matrix, emission params, start probs |
| `n_restarts` | integer | 10 | Random restarts to mitigate local optima |
| `max_iterations` | integer | 500 | Max EM iterations per restart |
| `convergence_tol` | float | 1e-6 | Log-likelihood convergence threshold |
| `state_range` | list of ints | `[n_hidden_states]` | Range state counts for selection |
| `covariance_type` | string | `"full"` | Gaussian: `"full"`, `"diagonal"`, `"spherical"` |
| `regularization` | float | 1e-6 | Diagonal constant preventing singularity |

## Do

### Step 1: Define Hidden States + Obs Model

1.1. Specify N hidden states `K` (or candidate range Step 5).

1.2. Emission distribution by data type:
   - Continuous: Gaussian (uni/multivariate)
   - Count: Poisson or negative binomial
   - Categorical: discrete/multinomial

1.3. Components:
   - **Transition matrix** `A` size `K x K`: `A[i,j] = P(z_t = j | z_{t-1} = i)`
   - **Emission params** `theta_k` each `k`: distribution-specific (mean + covariance Gaussian)
   - **Initial distribution** `pi`: `pi[k] = P(z_1 = k)`

1.4. Verify data: no missing, consistent dim, length `T >> K^2`.

→ HMM arch w/ `K` states, chosen emission family, clean data `T >> K^2`.

If err: missing → impute or remove. `T` too small → reduce `K` or get more data.

### Step 2: Initialize Params

2.1. Gen initial each of `n_restarts`:
   - **Transition**: Random stochastic (Dirichlet rows) or perturbed uniform
   - **Emission**: K-means clustering → init means; cluster variances Gaussian
   - **Initial distribution**: Uniform or proportional to cluster sizes

2.2. First restart: K-means-informed (strongest). Subsequent: random perturbations.

2.3. Verify valid:
   - Transition rows sum 1, positive
   - Emission in valid domain (PD covariance)
   - Initial sums 1

→ `n_restarts` sets of valid params, ≥1 data-driven.

If err: K-means fails → purely random w/ more restarts. Singular covariance → add regularization to diagonal.

### Step 3: Baum-Welch EM

3.1. **E-step (Forward-Backward)**:
   - Forward `alpha[t,k]` = P(o_1,...,o_t, z_t=k | model):
     - `alpha[1,k] = pi[k] * b_k(o_1)`
     - `alpha[t,k] = sum_j(alpha[t-1,j] * A[j,k]) * b_k(o_t)`
   - Backward `beta[t,k]` = P(o_{t+1},...,o_T | z_t=k, model):
     - `beta[T,k] = 1`
     - `beta[t,k] = sum_j(A[k,j] * b_j(o_{t+1}) * beta[t+1,j])`
   - State posterior `gamma[t,k]` = P(z_t=k | O, model):
     - `gamma[t,k] = alpha[t,k] * beta[t,k] / P(O | model)`
   - Transition posterior `xi[t,i,j]` = P(z_t=i, z_{t+1}=j | O, model).

3.2. **M-step (re-estimate)**:
   - Transition: `A[i,j] = sum_t(xi[t,i,j]) / sum_t(gamma[t,i])`
   - Emission weighted sufficient stats:
     - Gaussian mean: `mu_k = sum_t(gamma[t,k] * o_t) / sum_t(gamma[t,k])`
     - Gaussian covariance: weighted scatter matrix + regularization
     - Discrete: `b_k(v) = sum_t(gamma[t,k] * I(o_t=v)) / sum_t(gamma[t,k])`
   - Initial: `pi[k] = gamma[1,k]`

3.3. Log-likelihood: `log P(O | model) = log sum_k(alpha[T,k])`. Log-sum-exp → prevent underflow.

3.4. **Scaling**: Scaled forward-backward → prevent underflow long sequences. Normalize `alpha` each step + accumulate log scaling factors.

3.5. Repeat E + M until log-likelihood change < `convergence_tol` or `max_iterations`.

3.6. Across restarts → keep params w/ highest final log-likelihood.

→ Monotonically non-decreasing log-likelihood, converge w/in max. Final valid (stochastic matrices, PD covariances).

If err: log-likelihood decreases → bug E/M, verify formulas. Very slow → better init or increase max. Singular covariance → increase regularization.

### Step 4: Viterbi Decoding

4.1. Init:
   - `delta[1,k] = log(pi[k]) + log(b_k(o_1))`
   - `psi[1,k] = 0` (no predecessor)

4.2. Recurse `t = 2,...,T`:
   - `delta[t,k] = max_j(delta[t-1,j] + log(A[j,k])) + log(b_k(o_t))`
   - `psi[t,k] = argmax_j(delta[t-1,j] + log(A[j,k]))`

4.3. Terminate:
   - `z*_T = argmax_k(delta[T,k])`
   - Best path log-prob: `max_k(delta[T,k])`

4.4. Backtrace `t = T-1,...,1`:
   - `z*_t = psi[t+1, z*_{t+1}]`

4.5. Output decoded sequence `z* = (z*_1, ..., z*_T)` + log-prob.

4.6. Compare Viterbi path prob to total sequence prob from forward → dominance.

→ Single most-likely sequence length `T`, each in `{1,...,K}`. Viterbi log-prob ≤ total log-likelihood.

If err: Viterbi -inf log-prob → transition/emission prob zero where shouldn't. Add floor values preventing log(0).

### Step 5: Model Selection (BIC/AIC)

5.1. Each candidate `K` in `state_range` → fit full HMM (Steps 2-4).

5.2. Free params `p`:
   - Transition: `K * (K - 1)` (rows simplex)
   - Emission: family-dependent (Gaussian full covariance `d` dim: `K * (d + d*(d+1)/2)`)
   - Initial: `K - 1`

5.3. Information criteria:
   - `BIC = -2 * log_likelihood + p * log(T)`
   - `AIC = -2 * log_likelihood + 2 * p`
   - `AICc = AIC + 2*p*(p+1) / (T - p - 1)` (small-sample)

5.4. Select lowest BIC (consistency) or AIC (prediction). Report both.

5.5. Tabulate each `K`: log-likelihood, # params, BIC, AIC, convergence.

5.6. Optimal `K` at boundary → extend range + re-fit.

→ Clear min BIC/AIC → optimal N hidden states. Selected converged + interpretable.

If err: no clear min (monotonically decreasing BIC) → misspecified, try diff emission family. All poor log-likelihood → data may not follow HMM structure.

### Step 6: Validate Held-Out + Posterior

6.1. Split training/validation (80/20 or multiple sequences).

6.2. Fit training. Compute held-out log-likelihood via forward (no re-fit).

6.3. **Posterior decoding** (alt to Viterbi):
   - Each step → state w/ highest posterior: `z^_t = argmax_k(gamma[t,k])`
   - Maximizes expected # correctly decoded (vs Viterbi maximizing joint path).

6.4. Compare Viterbi + posterior:
   - Agreement rate between sequences
   - Disagreement regions → ambiguous assignments

6.5. State interpretability:
   - Examine emission params each state (means, variances, discrete)
   - Verify states correspond meaningful regimes in domain
   - Dwell times (diagonal `A`) reasonable

6.6. Held-out log-likelihood per observation + compare across orders → confirm training selection.

→ Held-out reasonably close to training (no severe overfit). Viterbi + posterior agree 90%+. States distinct + interpretable.

If err: held-out much worse than training → overfit, reduce `K` or increase regularization. States not interpretable → diff init or emission family.

## Check

- Log-likelihood monotonically non-decreasing Baum-Welch each restart
- Transition row-stochastic (rows sum 1, non-negative)
- Emission in valid domain (PD covariances, valid prob distributions)
- Viterbi log-prob ≤ total log-prob
- BIC/AIC clear min at selected order
- Held-out confirms generalization
- Forward + backward agree: `P(O) = sum_k(alpha[T,k]) = sum_k(pi[k] * b_k(o_1) * beta[1,k])`

## Traps

- **Local optima EM**: Baum-Welch → local max not global. Always multiple restarts + pick best.
- **Numerical underflow**: Forward-backward probs shrink exponential w/ length. Log-space or scaled variables.
- **Overfit too many states**: Each adds `O(K + d^2)` params. Use BIC not likelihood + validate held-out.
- **Label switching**: States identifiable only up to permutation. Compare across restarts → match by emission params not index.
- **Degenerate states**: State collapses to explain single observation (Gaussian near-zero variance). Regularization prevents.
- **Confuse Viterbi + posterior**: Viterbi = single best joint path; posterior = best marginal state each step. Different questions, can disagree significantly.
- **Ignore dwell times**: Geometric dwell-time in standard HMM may be poor fit for long regime durations. Consider hidden semi-Markov if non-geometric.

## →

- [Model Markov Chain](../model-markov-chain/SKILL.md) — prereq for transition structure underlying hidden layer
- [Simulate Stochastic Process](../simulate-stochastic-process/SKILL.md) — gen synthetic HMM data + simulate fitted for posterior predictive