---

name: fit-hidden-markov-model locale: caveman-ultra source_locale: en source_commit: 82c77053 fence_basis_commit: 82c77053 translator: "Julius Brussee homage — caveman" translation_date: "2026-04-24" description: > 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. license: MIT allowed-tools: Read Write Edit Bash Grep Glob metadata: author: Philipp Thoss version: "1.0" domain: stochastic-processes complexity: advanced language: multi tags: stochastic, hmm, baum-welch, viterbi, em-algorithm

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 — prereq for transition structure underlying hidden layer
• Simulate Stochastic Process — gen synthetic HMM data + simulate fitted for posterior predictive