Baum-Welch

Why you need this

Hand-authoring a 6 by 6 transition matrix and a 6 by 12 emission matrix is possible. Hand-authoring them so the generator's output resembles the twenty levels a designer already approved is guesswork. Baum-Welch turns the approved levels into the tables. The same procedure fits a regime chain from play logs, once you have some.

The idea

You have one or more observation sequences $O_1, \dots, O_T$ (tile symbols per column) and a chosen number of hidden states $N$. You want $\pi$, $A$, and $B$: $\pi_i = P(X_1 = i)$, $A_{ij} = P(X_{t+1} = j \mid X_t = i)$, $B_{jk} = P(O_t = k \mid X_t = j)$.

If the hidden states were visible, this would be counting. $A_{ij}$ is the number of times $i$ was followed by $j$, divided by the number of times $i$ occurred. $B_{jk}$ is the number of times $j$ emitted $k$, divided by the number of times $j$ occurred. But the states are hidden, so you cannot count them.

Expectation-maximization gets around this in plain words: guess the model, use it to compute how often each transition and emission probably happened (expected counts rather than actual counts), re-estimate the tables from those expected counts as if they were real, and repeat. Each round the model explains the data at least as well as the round before.

The backward variable

The forward algorithm gives $\alpha_t(i) = P(O_1..O_t, X_t = i)$, the probability of the observations up to $t$ ending in state $i$. You also need the probability of the observations after $t$ given state $i$:

$$\beta_t(i) = P(O_{t+1}, \dots, O_T \mid X_t = i).$$

It is computed right to left. At the end there is nothing left to observe, so $\beta_T(i) = 1$ for every $i$. Then

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

From $i$ at time $t$, move to $j$, emit $O_{t+1}$, and account for the rest with $\beta_{t+1}(j)$. It needs the same scaling as the forward pass to avoid underflow.

Two posterior quantities

Multiply forward and backward at the same time step and divide by $P(O) = \sum_j \alpha_T(j)$, and you have the probability of being in state $i$ at time $t$ given the whole sequence:

$$\gamma_t(i) = P(X_t = i \mid O) = \frac{\alpha_t(i) \, \beta_t(i)}{P(O)}.$$

Insert one transition and one emission between them, and you have the probability of the pair $(i, j)$ at times $(t, t+1)$:

$$\xi_t(i, j) = P(X_t = i, X_{t+1} = j \mid O) = \frac{\alpha_t(i) \, A_{ij} \, B_{j O_{t+1}} \, \beta_{t+1}(j)}{P(O)}.$$

Summing $\xi_t(i, j)$ over $j$ gives $\gamma_t(i)$ back, which is a useful check.

Re-estimation

Now treat $\gamma$ and $\xi$ as soft counts.

$$\hat{\pi}_i = \gamma_1(i)$$

is the expected number of times the sequence started in $i$ (with one sequence, just the probability).

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

is the expected number of transitions from $i$ to $j$ divided by the expected number of times in $i$ with a step still to come.

$$\hat{B}_{jk} = \frac{\sum_{t : O_t = k} \gamma_t(j)}{\sum_{t=1}^{T} \gamma_t(j)}$$

is the expected number of times state $j$ emitted symbol $k$ divided by the expected number of times in $j$.

With several training sequences, add the numerators across sequences and the denominators across sequences before dividing. Then replace $(\pi, A, B)$ with the hats and go again.

Worked example

One iteration on the smallest useful case. Same model as the forward and Viterbi pages: states 1 Run and 2 Gap, 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},$$

and the single sequence $O = (\text{flat}, \text{pit}, \text{flat})$, $T = 3$. The forward page computed $\alpha$ and $P(O) = 0.151504$.

	$t = 1$	$t = 2$	$t = 3$
$\alpha_t(\text{Run})$	0.72	0.0528	0.133488
$\alpha_t(\text{Gap})$	0.04	0.1856	0.018016
$\beta_t(\text{Run})$	0.1971	0.69	1
$\beta_t(\text{Gap})$	0.2398	0.62	1
$\gamma_t(\text{Run})$	0.93669	0.24047	0.88109
$\gamma_t(\text{Gap})$	0.06331	0.75953	0.11891

Backward, right to left. $\beta_2(\text{Run}) = 0.7 \cdot 0.9 \cdot 1 + 0.3 \cdot 0.2 \cdot 1 = 0.69$ (next observation is flat). $\beta_2(\text{Gap}) = 0.6 \cdot 0.9 + 0.4 \cdot 0.2 = 0.62$. $\beta_1(\text{Run}) = 0.7 \cdot 0.1 \cdot 0.69 + 0.3 \cdot 0.8 \cdot 0.62 = 0.0483 + 0.1488 = 0.1971$ (next observation is pit). $\beta_1(\text{Gap}) = 0.6 \cdot 0.1 \cdot 0.69 + 0.4 \cdot 0.8 \cdot 0.62 = 0.0414 + 0.1984 = 0.2398$. Check: $0.72 \cdot 0.1971 + 0.04 \cdot 0.2398 = 0.151504 = P(O)$.

Gamma, for instance $\gamma_2(\text{Run}) = 0.0528 \cdot 0.69 / 0.151504 = 0.036432 / 0.151504 = 0.24047$.

Xi, numerators then divided by $P(O)$:

$\xi_t(i, j)$	Run to Run	Run to Gap	Gap to Run	Gap to Gap
$t = 1$ (next is pit)	$0.72 \cdot 0.7 \cdot 0.1 \cdot 0.69 = 0.034776 \to 0.22954$	$0.72 \cdot 0.3 \cdot 0.8 \cdot 0.62 = 0.107136 \to 0.70715$	$0.04 \cdot 0.6 \cdot 0.1 \cdot 0.69 = 0.001656 \to 0.01093$	$0.04 \cdot 0.4 \cdot 0.8 \cdot 0.62 = 0.007936 \to 0.05238$
$t = 2$ (next is flat)	$0.0528 \cdot 0.7 \cdot 0.9 \cdot 1 = 0.033264 \to 0.21956$	$0.0528 \cdot 0.3 \cdot 0.2 \cdot 1 = 0.003168 \to 0.02091$	$0.1856 \cdot 0.6 \cdot 0.9 \cdot 1 = 0.100224 \to 0.66153$	$0.1856 \cdot 0.4 \cdot 0.2 \cdot 1 = 0.014848 \to 0.09800$

Check one row: $0.22954 + 0.70715 = 0.93669 = \gamma_1(\text{Run})$.

Re-estimate:

• $\hat{\pi} = (0.93669, 0.06331)$.
• $\hat{A}_{\text{Run},\text{Run}} = (0.22954 + 0.21956) / (0.93669 + 0.24047) = 0.44910 / 1.17716 = 0.3815$, and $\hat{A}_{\text{Run},\text{Gap}} = 0.72806 / 1.17716 = 0.6185$.
• $\hat{A}_{\text{Gap},\text{Run}} = (0.01093 + 0.66153) / (0.06331 + 0.75953) = 0.67246 / 0.82284 = 0.8172$, and $\hat{A}_{\text{Gap},\text{Gap}} = 0.1828$.
• $\hat{B}_{\text{Run},\text{flat}} = (0.93669 + 0.88109) / (0.93669 + 0.24047 + 0.88109) = 1.81778 / 2.05825 = 0.8832$, so $\hat{B}_{\text{Run},\text{pit}} = 0.1168$.
• $\hat{B}_{\text{Gap},\text{flat}} = (0.06331 + 0.11891) / 0.94175 = 0.1935$, so $\hat{B}_{\text{Gap},\text{pit}} = 0.8065$.

Every row still sums to 1. Under the new tables the same sequence has $P(O) = 0.331$ (log $-1.105$), up from $0.1515$ (log $-1.887$). It went up, as promised. It also shows the danger: from three observations the model has decided that Run almost always turns into Gap. That is the tiny dataset talking, not the world.

Convergence, local optima, initialization, and data

The log-likelihood $\log P(O \mid \pi, A, B)$ never decreases from one iteration to the next. Stop when it improves by less than some small amount, or after a fixed number of iterations (a few dozen is typical).

It converges to a local optimum, not necessarily the best one. Different starting tables land in different places, and a poor start can converge to states that mean nothing. The cheap cure is a good start: decode your training levels with Viterbi under a hand-guessed model, count the labeled transitions and emissions, and use those normalized counts as the initial $A$ and $B$. Run Baum-Welch from several starts and keep the one with the highest final log-likelihood.

Be honest about data size. Twenty approved levels of 100 columns is 2,000 observations, which sounds like a lot until you divide it among 36 transition entries and 72 emission entries. Some will be estimated from three or four soft occurrences. Regularize with pseudo-counts: add a small constant (0.5 or 1) to every expected count in the numerators, and the corresponding total to each denominator, before dividing. No entry goes to zero, and rare transitions stay possible instead of vanishing because they happened not to appear in the sample. In the worked example, adding 1 to each of the four expected transition counts pulls $\hat{A}$ back toward uniform: $\hat{A}_{\text{Run},\text{Run}} = (0.4491 + 1) / (1.1772 + 2) = 0.456$ instead of $0.3815$.

In a map generator

• Fitting the side-scroller HMM. The design tool holds levels the designer approved. Symbolize each column, initialize from Viterbi labels under a hand-guessed model, run Baum-Welch with pseudo-counts, and write $(\pi, A, B)$ into the parameter file. Generation then produces levels statistically like the approved ones. Re-run when the approved set grows.
• Fitting the regime chain. Play logs record which events a region produced each tick (spawns, resource pulls, floods). Treat the events as observations and the regime as hidden; the fitted $A$ replaces the hand-tuned one on the regime chains page.
• Validation. Hold out a few approved levels, and score them with the forward algorithm under the fitted model. If the held-out score is far below the training score, the model has memorized rather than learned. Add pseudo-counts or reduce $N$.

Common mistakes

• Not scaling $\beta$. The forward pass is scaled and the backward pass is not, so $\gamma$ is garbage after 300 columns. Use the same $c_t$ factors for both.
• Zero entries that never recover. A zero in $A$ or $B$ stays zero forever, because its expected count is always zero. Initialize with no exact zeros unless the zero is a rule you mean to enforce.
• Too many hidden states. With $N = 12$ and twenty levels, several states collapse into near-copies of each other and the log-likelihood barely moves. Start at the number of section types the designer can name.
• Fitting to one long sequence. Concatenating levels puts a fake transition from the end of one to the start of the next. Fit as separate sequences and sum the counts.
• Trusting one run. The first local optimum is not the best. Several restarts, compare log-likelihoods.

Cost

One iteration is a forward pass, a backward pass, and the $\xi$ sums, each $O(T \cdot N^2)$ per sequence, with $T$ the sequence length and $N$ the number of hidden states. Over $S$ sequences and $I$ iterations the total is $O(I \cdot S \cdot T \cdot N^2)$; with one sequence, $O(I \cdot T \cdot N^2)$. Memory is $O(T \cdot N)$ for the $\alpha$ and $\beta$ tables of the sequence being processed. For $N = 6$, twenty levels of 100 columns, and 50 iterations, that is about 3.6 million multiplications, well under a second. It hurts when $N$ climbs past 50 or when the training set reaches tens of thousands of sequences, neither of which a level tool reaches.

Going further

• Expectation-maximization in general, of which Baum-Welch is the HMM special case.
• Viterbi training (hard EM), which replaces $\gamma$ and $\xi$ with the counts from the decoded path: faster, cruder.
• Model selection, for choosing $N$ from held-out likelihood.
• Dirichlet priors, the principled version of pseudo-counts.