So far, we’ve considered state space models with a continuous latent space \boldsymbol{z}. However, we may often want to run our HMM models in situations were we are moving between discrete classes in the latent space. For example, we may be observing the fMRI scans of a patient that is waking up and falling back asleep. In that case our latent space consists of two discrete states - awake and asleep - that determine the fMRI patterns we will observe. In this class, we will call these Discrete HMM models¹. In the discrete case, we will still treat \boldsymbol{z}_{t} as a vector of the one-hot encoding. That is:

z_{t, k}= \begin{cases}1 & \text { if in class } k \text { at time } t \tag{7.22}\\ 0 & \text { otherwise }\end{cases}

If we assume that we have K classes, then our index k can go from 1, \ldots, K. If we only have K discrete classes, then our transition probability is fully described by specifying², for all k, k^{\prime}, the probability of transitioning from state k to k^{\prime} :

Since A_{k, k^{\prime}} is a probability we can state that:

From Equation 7.25 we can say that the matrix \boldsymbol{A} has K(K-1) degrees of freedom. We can write our full transition probability function as:

The initial latent node \boldsymbol{z}_{0} has a prior that is fully determined by a vector of probabilities \pi :

Again, we have the restriction that:

\begin{align} & 0 \leq \pi_{k} \leq 1 \tag{7.28}\\ & \sum_{k=1}^{K} \pi_{k}=1 \tag{7.29} \end{align}

meaning that there are only K-1 degrees of freedom to our vector \boldsymbol{\pi}.

The only remaining distribution we need to specify in our model is p\left(\boldsymbol{x}_{t} \mid \boldsymbol{z}_{t}\right). Until we want to optimize the parameters of p\left(\boldsymbol{x}_{t} \mid \boldsymbol{z}_{t}\right), we will not need to assume any form for the distribution. We will only require that we are able to evaluate it (the same assumption we made for particle filtering).

Our first goal is to infer the latent states of our discrete HMM model. The full posterior, p\left(\boldsymbol{z}_{1: T} \mid \boldsymbol{x}_{1: T}\right), requires us evaluating K^{T} possible configurations. For large latent states or long sequences, this quickly becomes computationally intractable. As with the LDS system, we’ll start instead by focusing on the exact marginals, p\left(\boldsymbol{z}_{t} \mid \boldsymbol{x}_{1: T}\right) and joint p\left(\boldsymbol{z}_{t}, \boldsymbol{z}_{t+1} \mid \boldsymbol{x}_{1: T}\right). These will be the quantities we need for future prediction and for optimizing the parameters of our model.

When we introduced our HMM models, we started by dividing our marginal distribution into two terms:

We utilize this refactoring of the equation because both \alpha and \beta can be calculated recursively through a forward and backward pass. For \alpha we have:

For \beta we have:

Note that we can write our normalizing constant from Equation 7.30 as:

p\left(\boldsymbol{x}_{1: T}\right)=\sum_{\boldsymbol{z}_{t}} \alpha\left(\boldsymbol{z}_{t}\right) \beta\left(\boldsymbol{z}_{t}\right) \tag{7.35}

and we can write our joint distribution as:

If we calculate the \alpha and \beta terms we have everything we need to evaluate the marginal and joint distribution we’re interested in. Equation 7.33 and Equation 7.34 are both recursive, so we need a starting point for our equation. For \beta, this is easy since:

because there is no data beyond T. For \alpha we can take advantage of our prior:

Before we work through an example, there’s an additional implementation caveat for the forward and backward pass of the \alpha and \beta calculation.

Specifically, as we go deeper into our forward pass, the values of \alpha are being multiplied by more and more factors of the transition and observation probability. Since these values are less than 1, as our sequence gets longer we expect \alpha to tend towards zero. This will eventually run into the floating point limit of our computer, and so it will be important to work with rescaled variables. The standard choice is:

\hat{\alpha}\left(\boldsymbol{z}_{t}\right)=p\left(\boldsymbol{z}_{t} \mid \boldsymbol{x}_{1: t}\right)=\frac{\alpha\left(\boldsymbol{z}_{t}\right)}{p\left(\boldsymbol{x}_{1: t}\right)} \tag{7.43}

It’s useful to define a scaling factor c_{t} as:

which is useful since we can write:

These scaling factors can be used to write our rescaled \hat{\alpha} in terms of the original \alpha :

This rescaled variable is normalized to one by construction, so it is guaranteed to be numerically stable. However, we now need a recursion relation for \hat{\alpha} :

In practice, c_{t} is just the normalization of the right-hand side of Equation 7.34.

\begin{align} c_{t} \sum_{\boldsymbol{z}_{t}} \hat{\alpha}\left(\boldsymbol{z}_{t}\right) & =\sum_{\boldsymbol{z}_{t}} p\left(\boldsymbol{x}_{t} \mid \boldsymbol{z}_{t}\right) \sum_{\boldsymbol{z}_{t-1}} \hat{\alpha}\left(\boldsymbol{z}_{t-1}\right) p\left(\boldsymbol{z}_{t} \mid \boldsymbol{z}_{t-1}\right) \tag{7.51}\\ c_{t} & =\sum_{\boldsymbol{z}_{t}} p\left(\boldsymbol{x}_{t} \mid \boldsymbol{z}_{t}\right) \sum_{\boldsymbol{z}_{t-1}} \hat{\alpha}\left(\boldsymbol{z}_{t-1}\right) p\left(\boldsymbol{z}_{t} \mid \boldsymbol{z}_{t-1}\right) \tag{7.52} \end{align}

Notice that \hat{\alpha} is the Kalman filtering distribution we calculated for the LDS model.

Given that we have a rescaled \alpha we will also want a rescaled \beta. As we will see, a good choice is:

This gives us the recursion relation:

Finally, we’ll leave it as an exercise for the reader to use Equation 7.27 and Equation 7.38 to verify that:

Let’s work through an example that uses these rescaled values. Imagine we have a latent space with K=2 and we have two observations such that T=2. Let’s start by specifying our prior and transition matrix:

\begin{align} \boldsymbol{\pi} & =\left[\begin{array}{l} 0.2 \\ 0.8 \end{array}\right] \tag{7.58}\\ \boldsymbol{A} & =\left[\begin{array}{ll} 0.5 & 0.5 \\ 0.5 & 0.5 \end{array}\right] \tag{7.59} \end{align}

Our prior preferences the k=2 state, but our transition matrix says that we are equally likely to move to either state no matter what state we are in. The last piece we need is p\left(\boldsymbol{x}_{t} \mid \boldsymbol{z}_{t}\right). Let’s keep this example simple by asserting that:

\begin{align} & p\left(\boldsymbol{x}_{1} \left\lvert\, \boldsymbol{z}_{1}=\left[\begin{array}{l} 1.0 \\ 0.0 \end{array}\right]\right.\right)=0.2 \tag{7.60}\\ & p\left(\boldsymbol{x}_{1} \left\lvert\, \boldsymbol{z}_{1}=\left[\begin{array}{l} 0.0 \\ 1.0 \end{array}\right]\right.\right)=0.8 \tag{7.61}\\ & p\left(\boldsymbol{x}_{2} \left\lvert\, \boldsymbol{z}_{2}=\left[\begin{array}{l} 1.0 \\ 0.0 \end{array}\right]\right.\right)=0.9 \tag{7.62}\\ & p\left(\boldsymbol{x}_{2} \left\lvert\, \boldsymbol{z}_{2}=\left[\begin{array}{l} 0.0 \\ 1.0 \end{array}\right]\right.\right)=0.1 \tag{7.63} \end{align}

That’s all we need to conduct our calculations. Let’s adopt the notation:

p\left(\boldsymbol{z}_{t}=\left[\begin{array}{l} 0.0 \tag{7.64}\\ 1.0 \end{array}\right] \left\lvert\, \boldsymbol{z}_{t-1}=\left[\begin{array}{l} 1.0 \\ 0.0 \end{array}\right]\right.\right)=p\left(z_{t, 2} \mid z_{t-1,1}\right)

and calculate \hat{\alpha} :

\begin{align} \hat{\alpha}\left(z_{0,1}\right) & =\pi_{1}^{z_{0,1}} \tag{7.65}\\ & =0.2 \tag{7.66}\\ \hat{\alpha}\left(z_{0,2}\right) & =0.8 \tag{7.67}\\ \hat{\alpha}\left(z_{1,1}\right) & =\frac{1}{c_{1}} p\left(\boldsymbol{x}_{1} \mid z_{1,1}\right) \sum_{\boldsymbol{z}_{0}} \hat{\alpha}\left(\boldsymbol{z}_{0}\right) p\left(z_{1,1} \mid \boldsymbol{z}_{0}\right) \tag{7.68}\\ & =\frac{1}{c_{1}} p\left(\boldsymbol{x}_{1} \mid z_{1,1}\right)\left[\hat{\alpha}\left(z_{0,1}\right) p\left(z_{1,1} \mid z_{0,1}\right)+\hat{\alpha}\left(z_{0,2}\right) p\left(z_{1,1} \mid z_{0,2}\right)\right] \tag{7.69}\\ & =\frac{1}{c_{1}} \times 0.2 \times[0.2 \times 0.5+0.8 \times 0.5] \tag{7.70}\\ & =\frac{1}{c_{1}} \times 0.1 \tag{7.71}\\ \hat{\alpha}\left(z_{1,2}\right) & =\frac{1}{c_{1}} p\left(\boldsymbol{x}_{1} \mid z_{1,2}\right) \sum_{\boldsymbol{z}_{0}} \hat{\alpha}\left(\boldsymbol{z}_{0}\right) p\left(z_{1,2} \mid \boldsymbol{z}_{0}\right) \tag{7.72}\\ & =\frac{1}{c_{1}} p\left(\boldsymbol{x}_{1} \mid z_{1,2}\right)\left[\hat{\alpha}\left(z_{0,1}\right) p\left(z_{1,2} \mid z_{0,1}\right)+\hat{\alpha}\left(z_{0,2}\right) p\left(z_{1,2} \mid z_{0,2}\right)\right] \tag{7.73}\\ & =\frac{1}{c_{1}} \times 0.8 \times[0.2 \times 0.5+0.8 \times 0.5] \tag{7.74}\\ & =\frac{1}{c_{1}} \times 0.4 \tag{7.75}\\ c_{1} & =0.5 \tag{7.76} \end{align}

Notice how we only get our value for c_{1} after we’ve calculated both of our \hat{\alpha} terms. The final two terms are:

\begin{align} \hat{\alpha}\left(z_{2,1}\right) & =\frac{1}{c_{2}} p\left(\boldsymbol{x}_{2} \mid z_{2,1}\right) \sum_{z_{1}} \hat{\alpha}\left(\boldsymbol{z}_{1}\right) p\left(z_{2,1} \mid \boldsymbol{z}_{1}\right) \tag{7.77}\\ & =\frac{1}{c_{2}} \times 0.45 \tag{7.78}\\ \hat{\alpha}\left(z_{2,2}\right) & =\frac{1}{c_{2}} p\left(\boldsymbol{x}_{2} \mid z_{2,2}\right) \sum_{\boldsymbol{z}_{1}} \hat{\alpha}\left(\boldsymbol{z}_{1}\right) p\left(z_{2,2} \mid \boldsymbol{z}_{1}\right) \tag{7.79}\\ & =\frac{1}{c_{2}} \times 0.05 \tag{7.80}\\ c_{2} & =0.5 \tag{7.81} \end{align}

We leave it as an exercise to the reader to calculate the \hat{\beta} values, but notice that you already have all the constants c_{t} you will need for that calculation.

We did not need to know the explicit form of p\left(\boldsymbol{x}_{t} \mid \boldsymbol{z}_{t}\right) to calculate the exact values of \hat{\alpha}, we only needed to be able to evaluate it. As we will see in next week’s lecture, we will only be interested in the parameterization of p\left(\boldsymbol{x}_{t} \mid \boldsymbol{z}_{t}\right) when we want to update those parameters.

As with out continuous latent space, we want to build the machinery to infer the parameters of our model:

\begin{equation*} p(\boldsymbol{\theta} \mid \boldsymbol{X}) \propto p(\boldsymbol{X} \mid \boldsymbol{\theta}) p(\boldsymbol{\theta}), \tag{7.82} \end{equation*}

where \theta=\{\boldsymbol{\pi}, \boldsymbol{A}, \ldots\}. We have not yet specified the parameters of our observation probability function. The EM algorithm we derived for a continuous latent space (see Algorithm 1) translates perfectly to our discrete case. We already know the choice we want to make for q_{n} :

\begin{equation*} q_{n}\left(\boldsymbol{z}_{1: T}\right)=p\left(\boldsymbol{z}_{1: T} \mid \boldsymbol{x}_{1: T}, \boldsymbol{\theta}_{n-1}\right) . \tag{7.83} \end{equation*}

As with the continuous case, we do not have access to this full joint posterior. We only have the marginal and the joint over two neighboring latent states. However, as with our continuous case, the structure of our model will come to our rescue. Once we have specified q_{n} we need to maximize:

\begin{align*} \mathbb{E}_{q_{n}}\left[\log p\left(\boldsymbol{x}_{1: T}, \boldsymbol{z}_{1: T} \mid \boldsymbol{\theta}\right)\right] &= \mathbb{E}_{q_{n}}\left[\log p\left(\boldsymbol{z}_{0}\right)+\sum_{t=1}^{T} \log p\left(\boldsymbol{z}_{t} \mid \boldsymbol{z}_{t-1}\right)\right. \\ & \quad \left.+\sum_{t=1}^{T} \log p\left(\boldsymbol{x}_{t} \mid \boldsymbol{z}_{t}\right)\right] \tag{7.84} \end{align*}

with respect to \boldsymbol{\theta}. All of the terms in our equation above either depend on a single latent state or a pair of latent states, so taking an expectation of the full joint boils down to summing over the probabilities for the marginals and joints of neighboring variables:

\begin{align*} & \mathbb{E}_{q_{n}}\left[\log p\left(\boldsymbol{x}_{1: T}, \boldsymbol{z}_{1: T} \mid \boldsymbol{\theta}\right)\right]=\sum_{k=1}^{K} p\left(z_{0, k} \mid \boldsymbol{x}_{1: T}\right) \log \pi_{k} \\ & +\sum_{t=1}^{T} \sum_{k^{\prime}=1}^{K} \sum_{k=1}^{K} p\left(z_{t, k^{\prime}}, z_{t-1, k} \mid \boldsymbol{x}_{1: T}\right) \log A_{k k^{\prime}} \\ & \quad+\sum_{k=1}^{K} \sum_{t=1}^{T} p\left(z_{t, k} \mid \boldsymbol{x}_{1: T}\right) \log p\left(\boldsymbol{x}_{t} \mid z_{t, k}\right) \tag{7.85} \end{align*}

We have some constraints on \boldsymbol{\pi} and \boldsymbol{A} that must be enforced with Lagrange multipliers, but otherwise we can maximize these equations with respect to \boldsymbol{\pi} and \boldsymbol{A} to find:

where the n subscript encodes that these are the parameters corresponding to EM step n where we used distribution q_{n}.

We can also solve for the parameters of our observation probability function so long as we specify its form. One choice that will often come up (including in lab) is the form:

\begin{equation*} p\left(\boldsymbol{x}_{t} \mid z_{t, k}\right)=\mathcal{N}\left(\boldsymbol{x}_{t} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right) \tag{7.88} \end{equation*}

This is a model where each discrete latent state k represents a multivariate normal observational probability function with its own mean, \boldsymbol{\mu}_{k}, and covariance, \boldsymbol{\Sigma}_{k}. If we plug this choice into Equation 7.85 we get:

---

Viterbi algorithm with discrete HMMs

The discrete nature of our latent space makes it possible for us to ask what the most likely sequence of latent states is. This is not the same as the MAP estimate for \boldsymbol{z}_{t} at each t because we are now considering the full joint p\left(\boldsymbol{z}_{1: T} \mid \boldsymbol{x}_{1: T}\right). We do not want to calculate the full joint posterior since this would require K^{T} terms. However, if we are only interested in the MAP sequence and not the full joint, we can do a much more efficient dynamic programming calculation known as the Viterbi algorithm.

To start, imagine that at time step t you have the probability of the most likely path where z_{t, k}=1. We will call this w_{t, k} :

Imagine that we also have the index m_{t-1, k} of the class for z_{t-1} that maximized the expression in Equation 7.91. First, notice the dimensionality of these two objects. The weights, w_{t, k}, must be calculated for our K classes for each latent state, and therefore has dimensionality (T+1) \times K. Since each weight except \boldsymbol{w}_{0} has a corresponding m_{t-1, k}, the indices will have dimensionality T \times K. This is much more manageable than K^{T}.

If we already had all of these weights, then:

\begin{align*} k_{T}^{\max } & =\underset{k}{\operatorname{argmax}} w_{T, k} \tag{7.92}\\ & =\underset{k}{\operatorname{argmax}} \max _{\boldsymbol{z}_{1: T-1}} \log p\left(\boldsymbol{z}_{1: T-1}, z_{T, k} \mid \boldsymbol{x}_{1: T}\right), \tag{7.93} \end{align*}

which is the MAP class for the latent state \boldsymbol{z}_{T}. The MAP class for the latent state \boldsymbol{z}_{T-1}, k_{T-1}^{\max }, would then be:

\begin{equation*} k_{T-1}^{\max }=m_{T-1, k_{T}^{\max }} \tag{7.94} \end{equation*}

and more generally:

\begin{equation*} k_{t-1}^{\max }=m_{t-1, k_{t}^{\max }} \tag{7.95} \end{equation*}

So now we just need to show that we can calculate w_{t, k} and m_{t-1, k} recursively. If we have the maximum probability of the most likely path for each class k at timestep t-1 then:

And so the Viterbi algorithm allows you to make one pass through the data to calculate the MAP sequence for our model and observations. In our framing, the initial weight distribution of the Viterbi algorithm is defined by:

\begin{equation*} w_{0, k}=\log p\left(z_{0, k}\right) \tag{7.98} \end{equation*}