A first hidden Markov model - P8149 Computational Vignettes

We ask one question. How can a sequence of noisy observations tell us about a sequence of hidden states?

A hidden Markov model (HMM) has three parts.

Hidden state. $Z_t\in\{0,1\}$ is the state at position $t$ .
Transition. $P(Z_t\neq Z_{t-1})=s$ controls how often the state changes.
Emission. $P(X_t\neq Z_t)=e$ controls how often the observation differs from the state.

The forward--backward algorithm uses dynamic programming to sum over all $2^T$ possible hidden paths without listing them. It computes $P(Z_t=1\mid X_1,\ldots,X_T)$ at every position.

One fixed sequence¶

We use one simulated hidden path and four deliberately flipped observations. The hidden path is shown because this is a teaching example. In real data, it is unknown.

hidden_state <- c(rep(0, 8), rep(1, 9), rep(0, 7), rep(1, 8))
observed <- hidden_state
flip_positions <- c(4, 12, 20, 27)
observed[flip_positions] <- 1 - observed[flip_positions]
position <- seq_along(observed)

forward_backward <- function(observed, switch_probability, emission_error) {
  n_position <- length(observed)
  transition <- matrix(
    c(1 - switch_probability, switch_probability,
      switch_probability, 1 - switch_probability),
    nrow = 2, byrow = TRUE
  )
  emission_probability <- function(x) {
    ifelse(c(0, 1) == x, 1 - emission_error, emission_error)
  }

  forward <- matrix(0, n_position, 2)
  forward[1, ] <- 0.5 * emission_probability(observed[1])
  forward[1, ] <- forward[1, ] / sum(forward[1, ])
  for (t in 2:n_position) {
    forward[t, ] <- as.vector(forward[t - 1, ] %*% transition) *
      emission_probability(observed[t])
    forward[t, ] <- forward[t, ] / sum(forward[t, ])
  }

  backward <- matrix(1, n_position, 2)
  for (t in (n_position - 1):1) {
    backward[t, ] <- transition %*%
      (emission_probability(observed[t + 1]) * backward[t + 1, ])
    backward[t, ] <- backward[t, ] / sum(backward[t, ])
  }

  posterior <- forward * backward
  posterior <- posterior / rowSums(posterior)
  posterior[, 2]
}

sequence_data <- data.frame(
  position, hidden_state, observed,
  flipped = position %in% flip_positions
)
sequence_data[sequence_data$flipped, ]

Transition probability controls persistence¶

Move the slider for $s$ . A small value says neighboring positions probably share a state, so the posterior smooths over isolated observations. A larger value permits more frequent state changes, so the posterior follows local observations more closely.

switch_grid <- c(0.01, 0.03, 0.05, 0.10, 0.20, 0.35, 0.50)
switch_frames <- do.call(rbind, lapply(switch_grid, function(s) {
  data.frame(
    position = position,
    hidden_state = hidden_state,
    observed = observed,
    posterior = forward_backward(observed, s, emission_error = 0.15),
    setting = sprintf("%.2f", s)
  )
}))
switch_frames$setting <- factor(
  switch_frames$setting, levels = sprintf("%.2f", switch_grid)
)

invisible(switch_frames)

A posterior probability is not a decoded path. It measures uncertainty at each position. Values near $1/2$ mean that the two hidden states remain difficult to distinguish.

Emission error controls how much the observations are trusted¶

Now hold $s=0.08$ fixed and move the slider for $e$ . When $e$ is small, an observation strongly favors the matching state. As $e$ approaches $1/2$ , the observations carry less information and the posterior is driven mainly by state persistence.

error_grid <- c(0.01, 0.05, 0.10, 0.20, 0.35, 0.49)
error_frames <- do.call(rbind, lapply(error_grid, function(e) {
  data.frame(
    position = position,
    hidden_state = hidden_state,
    observed = observed,
    posterior = forward_backward(
      observed, switch_probability = 0.08, emission_error = e
    ),
    setting = sprintf("%.2f", e)
  )
}))
error_frames$setting <- factor(
  error_frames$setting, levels = sprintf("%.2f", error_grid)
)

invisible(error_frames)

What to remember¶

The transition probability describes how states persist or change along a sequence.
The emission probability describes how a hidden state produces an observation.
Forward--backward uses dynamic programming to sum over all possible hidden paths and return a posterior probability at each position.

The Li–Stephens model uses the same structure. The hidden state is the copied reference haplotype, a transition changes the copying source, and an emission mismatch represents imperfect copying.