re-estimation of emission probabilities in HMM

Question

I am confused about the re-erstimation procedure for emissions in HMMs with Baum-Welch (still). I posted two questions concerning this general topic already and I thought I had cleared up my confusion, but not so. At least not in its entirety. The remaining issue I have is this:

Say I have a HMM with some number of states and a set of 4 emissions. Assume for state $s_k$ I have the following distribution of emissions:

e1    e2    e3    e4
.2    .1    .3    .4

Now suppose the next training sequence is e3 e1

Then the re-estimations for the emissions of state $s_k$ are the following:

$$b_{s_k}(e_1) = \frac{\gamma_1(k)}{\gamma_1(k) + \gamma_2(k)}$$

Let's just say this was equal to .3

$$b_{s_k}(e_2) = \frac{0}{\gamma_1(k) + \gamma_2(k)} = 0$$

$$b_{s_k}(e_3) = \frac{\gamma_2(k)}{\gamma_1(k) + \gamma_2(k)}$$

Let's assume this to be .7 then.

$$b_{s_k}(e_4) = \frac{0}{\gamma_1(k) + \gamma_2(k)} = 0$$

Ok so the re-estimated distribution looks like this then:

e1    e2    e3    e4
.3    .0    .7    .0

After doing this re-estimation for all states with the current sequence I assign the newly estimated ditributions to my HMM by swapping them out for the old distributions and go to the next training sequence.

So, for emissions not observed in the current training sequence I have just wiped out all emission probabilities I had estimated so far. If the next sequence would be for example e2 e4 I would end up wipig the estimations for e1 and e3 again. This cannot be correct, so I must (still) be missing something here. What is that though?

Edit 1

The way I now went about this is to just replace the non-zero re-estimates and keep the rest and finally re-normalise all emissions. No idea if that actually does the job or messes with the guanrantee of convergence, but I see no other way.

Edit 2

It does seem to mess with with the rest of the algorithm, because this leads to all transition probabilities per state to converge to 0 with one transition converging to 1.

Why has nobody an answer to this? Is the example confusing? Is my general problem unclear?

Edit 3

Ok, so just to make sure I got this right: I first build expected values like so: I go over the data set and for each sequence I compute $\sum_{t=1}^{|y|} \xi_t(i, j)$ and accumulate these values for instance in a matirx cell $\xi\text{-acc}[i, j]$

So I end up with

$$\forall i, j : \text{expect-trans}[i, j] \leftarrow \sum_{y \in Y} \sum_{t=1}^{|y|} \xi_t(i, j)$$

These are all the expected values for all transitions $i \to j$.

Then I do the same for the gamma values:

$$\forall i : \text{expect-state}[i] \leftarrow \sum_{y \in Y} \sum_{t=1}^{|y|} \gamma_t(i)$$

These are all the expected values for visiting state $i$.

Then I also still need to do this for the emissions

$$\forall i, e : \text{expect-emiss}[i, e] \leftarrow \sum_{y \in Y} \sum_{t=1}^{|y|} \gamma_t(i) \cdot \chi_{[y_t = e]}$$

These are all teh expected counts for emitting $e$ from $i$.

Once I have all of these expected values I can go over to the maximization step: \begin{align*} &\forall i : \text{HMM.init-prob(i)} & \leftarrow &&\frac{\text{expect-state}[i]}{\sum_{j=1}^N \text{expect-state}[j]}\\ \\ &\forall i, j : \text{HMM.trans-prob(i, j)} & \leftarrow &&\frac{\text{expect-trans}[i, j]}{\text{expect-state}[i]}\\ \\ &\forall i, e : \text{HMM.emiss-prob(i, e)} & \leftarrow &&\frac{\text{expect-emiss}[i, e]}{\sum_{j=1}^N \text{expect-state}[j]}\\ \end{align*}

Have I got this right now? :/

---

Or is it rather this:

\begin{align*} &\forall y \in Y, i, j, t : \text{expect-trans}[i, j, t] &\leftarrow && \xi_t(i, j) + \text{expect-trans}[i, j, t]\\ \\ &\forall y \in Y i, t : \text{expect-state}[i, t] &\leftarrow && \gamma_t(i) + \text{expect-state}[i, t]\\ \\ &\forall y \in Y, i, e, t : \text{expect-emiss}[i, e, t] &\leftarrow && \gamma_t(i) \cdot \chi_{[y_t = e]} + \text{expect-emiss}[i, e, t] \end{align*}

with the max $t$ depending on the length of the current $y$.

\begin{align*} &\forall i : \text{HMM.init-prob(i)} & \leftarrow &&\frac{\text{expect-state}[i, 1]}{\sum_{j=1}^N \text{expect-state}[j, 1]}\\ \\ &\forall i, j : \text{HMM.trans-prob(i, j)} & \leftarrow &&\frac{\sum_{t=1}^T\text{expect-trans}[i, j, t]}{\sum_{t=1}^T\text{expect-state}[i, t]}\\ \\ &\forall i, e : \text{HMM.emiss-prob(i, e)} & \leftarrow &&\frac{\sum_{t=1}^T \text{expect-emiss}[i, e, t]}{\sum_{t=1}^T \sum_{j=1}^N \text{expect-state}[j, t]}\\ \end{align*}

with $T$ the max index of the rows $\text{expect-trans}[i, j, t]$ and $\text{expect-emiss}[i, e]$ respectively.

Answer 1

You write:

After doing this re-estimation for all states with the current sequence I assign the newly estimated ditributions to my HMM by swapping them out for the old distributions and go to the next training sequence.

in standard Baum-Welch, or any "vanilla" EM algorithm, you gather expected counts from all the training sequences using one set of parameters $\mathbf{b}^{(i)}$ for the current iteration $i$, and only re-estimate the next iteration's $\mathbf{b}^{(i+1)}$ after you have expected counts from the entire dataset.

There are incremental versions of EM that update parameters after examining only one data item or a "minibatch" of several data items. One incremental EM algorithm due to Neal and Hinton (1999) starts with just an initial guess at the expected counts for data item, and then sweeps through the dataset, potentially many times, to optimize the expected counts for each item.

For HMMs, you would start with an initial guess for each $\gamma^{(0)}_{jl}(s_k)$ for each emission type $e_l$ in each sequence $j$, and compute your initial transition $\mathbf{a}^{(0)}$ and emission $\mathbf{b}^{(0)}$ probabilities using those counts. Next, you would visit each sequence, and estimate your new counts $\gamma^{(i+1)}_{jl}(s_k)$ using your current parameters $\mathbf{a}^{(i)}$ and $\mathbf{b}^{(i)}$, and re-estimate the next value for $\mathbf{b}^{(i+1)}$ by replacing your previous expected counts (for the current sequence only!) with the new expected counts: \begin{align} \gamma^{(i+1)}_{jl}( s_k ) & = \gamma^{(i)}_{jl}( s_k ) & \mbox{for }j\neq\mbox{ the current sequence} \\ \gamma^{(i+1)}_{jl}( s_k ) & = \mbox{E}_{\mathbf{b}^{(i)},\mathbf{a}^{(i)}}\left[ \mbox{count}\left( s_k \rightarrow e_l, j \right) \right] & \mbox{for }j = \mbox{ the current sequence} \\ b^{(i+1)}_{sk}( e_l ) & = \frac{ \sum_{j=1}^n \gamma^{(i+1)}_{jl} }{ \sum_{l'=1}^m\sum_{j=1}^n \gamma^{(i+1)}_{jl'} } \end{align} where $m$ is the number of emission types and $n$ is the number of training sequences. Note that you'll also need to re-estimate transition probabilities incrementally in exactly the same fashion.

Neal and Hinton (1999) point out that you can cheaply update the numerator and denominator by storing the sums in summary variables $\tilde{\gamma}^{(i)}_{l}$ and then subtracting the old counts and adding the new counts. For example, here is the numerator: \begin{align} \tilde{\gamma}^{(0)}_{l} & = \sum_{j=1}^n \gamma^{(0)}_{jl} \\ \tilde{\gamma}^{(i+1)}_{l} & = \tilde{\gamma}^{(i)}_l - \gamma^{(i)}_{jl} + \gamma^{(i+1)}_{jl} \end{align}