Forward algorithm

Why you need this

Two questions come up the moment you have a hidden Markov model. "How plausible is this level under my model?" lets a tool reject a generated level that is technically possible but unlike anything a designer approved. "Given what the player just did, what mode are they probably in?" lets a game adapt without asking. Both are the same computation, and the forward algorithm does it in one pass over the sequence.

The idea

The model is $(\pi, A, B)$: $N$ hidden states, initial distribution $\pi_i = P(X_1 = i)$, transitions $A_{ij} = P(X_{t+1} = j \mid X_t = i)$, emissions $B_{jk} = P(O_t = k \mid X_t = j)$. You are handed observations $O_1, \dots, O_T$ and want

$$P(O_1, \dots, O_T \mid \pi, A, B),$$

the probability that the model produces exactly this sequence by any hidden path at all.

The direct route is to list every hidden path, compute its probability times the probability it emits these observations, and add them up. There are $N^T$ paths. For a modest $N = 6$ section types and a $T = 100$ column level that is $6^{100}$, about $10^{78}$ terms. Not slow: impossible.

The escape is dynamic programming. Define the forward variable

$$\alpha_t(j) = P(O_1, \dots, O_t, X_t = j),$$

the probability of having seen the first $t$ observations and being in state $j$ at time $t$. It bundles together every path that ends in $j$ at $t$, and that bundle is all you need to continue, because of the Markov property.

Initialization. At $t = 1$ there is no previous state. Start in $j$ and emit $O_1$:

$$\alpha_1(j) = \pi_j \, B_{j O_1}.$$

Recursion. To be in $j$ at $t+1$, you were in some $i$ at $t$ (with all of $O_1..O_t$ already accounted for by $\alpha_t(i)$), moved $i \to j$, and emitted $O_{t+1}$:

$$\alpha_{t+1}(j) = \Big[ \sum_{i=1}^{N} \alpha_t(i) \, A_{ij} \Big] B_{j O_{t+1}}.$$

Termination. The sequence ended in some state, so sum over them:

$$P(O_1..O_T) = \sum_{j=1}^{N} \alpha_T(j).$$

The table has $T$ columns and $N$ rows. Each cell costs a sum of $N$ terms. That is $T N^2$ multiplications instead of $N^T$.

Worked example

Two hidden states, 1 Run and 2 Gap. Two tile symbols, 1 flat and 2 pit.

$$\pi = (0.8, 0.2), \qquad A = \begin{pmatrix} 0.7 & 0.3 \\ 0.6 & 0.4 \end{pmatrix}, \qquad B = \begin{pmatrix} 0.9 & 0.1 \\ 0.2 & 0.8 \end{pmatrix}.$$

Rows of $A$: Run stays Run 0.7, Gap returns to Run 0.6. Rows of $B$: Run shows flat 0.9, Gap shows pit 0.8. Observations: $O = (\text{flat}, \text{pit}, \text{flat})$, $T = 3$.

	$t = 1$ (flat)	$t = 2$ (pit)	$t = 3$ (flat)
$\alpha_t(\text{Run})$	$0.8 \cdot 0.9 = 0.72$	$(0.72 \cdot 0.7 + 0.04 \cdot 0.6) \cdot 0.1 = 0.528 \cdot 0.1 = 0.0528$	$(0.0528 \cdot 0.7 + 0.1856 \cdot 0.6) \cdot 0.9 = 0.14832 \cdot 0.9 = 0.133488$
$\alpha_t(\text{Gap})$	$0.2 \cdot 0.2 = 0.04$	$(0.72 \cdot 0.3 + 0.04 \cdot 0.4) \cdot 0.8 = 0.232 \cdot 0.8 = 0.1856$	$(0.0528 \cdot 0.3 + 0.1856 \cdot 0.4) \cdot 0.2 = 0.09008 \cdot 0.2 = 0.018016$
Column sum	0.76	0.2384	0.151504

Termination: $P(O) = 0.133488 + 0.018016 = 0.151504$.

You can check this against the brute-force sum. There are $2^3 = 8$ hidden paths. The largest single term is Run, Gap, Run: $0.8 \cdot 0.9 \cdot 0.3 \cdot 0.8 \cdot 0.6 \cdot 0.9 = 0.093312$. Adding all eight gives $0.151504$, the same number.

Filtering: what state are we in right now?

Divide a column by its sum and you get the probability of each hidden state given everything observed so far:

$$P(X_t = j \mid O_1, \dots, O_t) = \frac{\alpha_t(j)}{\sum_i \alpha_t(i)}.$$

From the table: after "flat" the model is 0.72 / 0.76 = 0.947 sure of Run. After "flat, pit" it is 0.1856 / 0.2384 = 0.779 sure of Gap. After "flat, pit, flat" it is 0.133488 / 0.151504 = 0.881 sure of Run again. The belief tracks the evidence one step behind, which is exactly what you want from a game reading a player's actions.

Log space and scaling

Each $\alpha_t(j)$ is a product of about $2t$ probabilities. If the typical factor is $0.5$, then by $t = 100$ the numbers are around $0.5^{200} \approx 10^{-60}$, and a few hundred columns later they fall below the smallest positive double (about $10^{-308}$), so the whole column becomes zero (see logarithms and underflow). The standard fix is to scale each column. After computing the raw column, sum it to get $c_t$, divide every entry by $c_t$, and remember $\log c_t$. The scaled column is the filtered distribution from the previous section, so it always sums to 1 and never underflows. The total probability is the product of the scale factors, so

$$\log P(O_1..O_T) = \sum_{t=1}^{T} \log c_t.$$

On the example: $c_1 = 0.76$, and the scaled column is $(0.9474, 0.0526)$. Feeding that scaled column into the recursion gives raw values $(0.069474, 0.244211)$, so $c_2 = 0.313684$. One more step gives $c_3 = 0.635503$. Then $\log P = \log 0.76 + \log 0.313684 + \log 0.635503 = -0.2744 - 1.1594 - 0.4533 = -1.8871$, and $e^{-1.8871} = 0.1515$. Same answer, no tiny numbers anywhere.

Report scores as log-probabilities and compare them per column ($\log P / T$) so that levels of different lengths are comparable.

JavaScript

// pi[i], A[i][j], B[j][k] are arrays of numbers; O is an array of symbol indices (0-based).
// Returns logP and filtered[t][j] = P(X_t = j | O_1..O_t).
function forward(pi, A, B, O) {
  const N = pi.length, T = O.length;
  const filtered = [];
  let logP = 0, prev = null;
  for (let t = 0; t < T; t++) {
    const cur = new Array(N).fill(0);
    for (let j = 0; j < N; j++) {
      let s = 0;
      if (t === 0) s = pi[j];
      else for (let i = 0; i < N; i++) s += prev[i] * A[i][j];
      cur[j] = s * B[j][O[t]];
    }
    const c = cur.reduce((a, b) => a + b, 0);       // scale factor for this column
    if (c === 0) return { logP: -Infinity, filtered }; // model says this sequence is impossible
    for (let j = 0; j < N; j++) cur[j] /= c;
    logP += Math.log(c);
    filtered.push(cur);
    prev = cur;
  }
  return { logP, filtered };
}
const pi = [0.8, 0.2], A = [[0.7, 0.3], [0.6, 0.4]], B = [[0.9, 0.1], [0.2, 0.8]];
console.log(forward(pi, A, B, [0, 1, 0]));   // logP -1.8871..., filtered[2] = [0.881, 0.119]

In a map generator

• Filtering the player. Hidden states are player modes (exploring, fighting, stuck); observations are bucketed actions (move, jump, attack, idle). Run one forward step per action and read the filtered column. When "stuck" passes 0.7 for several steps, offer a hint or nudge the regime chain toward Calm. The cost is $N^2$ multiplications per action, which is nothing.
• Scoring a candidate level. After the layout stage emits a column sequence, compute $\log P / T$ under the fitted model and reject or resample anything more than a chosen margin below the training levels' typical score. This is how the tool keeps a fitted generator from wandering into tile combinations the designer never approved.
• Comparing models. Two candidate parameter files, one set of approved levels: the one with the higher total log-probability fits better. This is also the quantity Baum-Welch drives upward.

Common mistakes

• Skipping the scaling. Works on 20-column test levels, returns $\log 0 = -\infty$ for every real level. Scale from day one.
• Comparing raw log-probabilities of different lengths. A 200-column level always scores lower than a 50-column one. Divide by $T$.
• A zero in $B$ for a tile that actually appears. The whole sequence gets probability 0 and every downstream score is $-\infty$. Give every emission a small floor unless it is truly forbidden by the tile catalog.
• Multiplying by $A$ on the wrong side. $\sum_i \alpha_t(i) A_{ij}$ walks down column $j$ of $A$. Using row $j$ instead gives numbers that still sum to something plausible and are wrong.
• Reading the filtered column as the most likely path. The most likely state at each time, read independently, can form a path that has zero probability as a whole. For the best single path use Viterbi.

Cost

Time is $O(T \cdot N^2)$: $T$ observations, and each of the $N$ cells in a column sums over $N$ predecessors. Memory is $O(N)$ if you only keep the current column and the running log sum, or $O(T \cdot N)$ if you keep the filtered table for later inspection or for Baum-Welch. With $N$ under 20 and $T$ in the hundreds this runs in microseconds; it only starts to hurt when Baum-Welch calls it thousands of times over hundreds of sequences.

Going further

• The Viterbi algorithm, the same table with max in place of sum.
• Baum-Welch, which pairs the forward table with a backward table to fit the model.
• The backward variable $\beta_t(i)$ and smoothing, for the probability of a past state given the whole sequence.
• Log-sum-exp, an alternative to scaling that works entirely in log space.