This is an old version, view current version.

2.6 Hidden Markov models

A hidden Markov model (HMM) generates a sequence of \(T\) output variables \(y_t\) conditioned on a parallel sequence of latent categorical state variables \(z_t \in \{1,\ldots, K\}\). These “hidden” state variables are assumed to form a Markov chain so that \(z_t\) is conditionally independent of other variables given \(z_{t-1}\). This Markov chain is parameterized by a transition matrix \(\theta\) where \(\theta_k\) is a \(K\)-simplex for \(k \in \{ 1, \dotsc, K \}\). The probability of transitioning to state \(z_t\) from state \(z_{t-1}\) is \[ z_t \sim \textsf{categorical}(\theta_{z[t-1]}). \] The output \(y_t\) at time \(t\) is generated conditionally independently based on the latent state \(z_t\).

This section describes HMMs with a simple categorical model for outputs \(y_t \in \{ 1, \dotsc, V \}\). The categorical distribution for latent state \(k\) is parameterized by a \(V\)-simplex \(\phi_k\). The observed output \(y_t\) at time \(t\) is generated based on the hidden state indicator \(z_t\) at time \(t\), \[ y_t \sim \textsf{categorical}(\phi_{z[t]}). \] In short, HMMs form a discrete mixture model where the mixture component indicators form a latent Markov chain.

Supervised parameter estimation

In the situation where the hidden states are known, the following naive model can be used to fit the parameters \(\theta\) and \(\phi\).

data {
  int<lower=1> K;             // num categories
  int<lower=1> V;             // num words
  int<lower=0> T;             // num instances
  array[T] int<lower=1, upper=V> w;  // words
  array[T] int<lower=1, upper=K> z;  // categories
  vector<lower=0>[K] alpha;   // transit prior
  vector<lower=0>[V] beta;    // emit prior
}
parameters {
  array[K] simplex[K] theta;        // transit probs
  array[K] simplex[V] phi;          // emit probs
}
model {
  for (k in 1:K) {
    theta[k] ~ dirichlet(alpha);
  }
  for (k in 1:K) {
    phi[k] ~ dirichlet(beta);
  }
  for (t in 1:T) {
    w[t] ~ categorical(phi[z[t]]);
  }
  for (t in 2:T) {
    z[t] ~ categorical(theta[z[t - 1]]);
  }
}

Explicit Dirichlet priors have been provided for \(\theta_k\) and \(\phi_k\); dropping these two statements would implicitly take the prior to be uniform over all valid simplexes.

Start-state and end-state probabilities

Although workable, the above description of HMMs is incomplete because the start state \(z_1\) is not modeled (the index runs from 2 to \(T\)). If the data are conceived as a subsequence of a long-running process, the probability of \(z_1\) should be set to the stationary state probabilities in the Markov chain. In this case, there is no distinct end to the data, so there is no need to model the probability that the sequence ends at \(z_T\).

An alternative conception of HMMs is as models of finite-length sequences. For example, human language sentences have distinct starting distributions (usually a capital letter) and ending distributions (usually some kind of punctuation). The simplest way to model the sequence boundaries is to add a new latent state \(K+1\), generate the first state from a categorical distribution with parameter vector \(\theta_{K+1}\), and restrict the transitions so that a transition to state \(K+1\) is forced to occur at the end of the sentence and is prohibited elsewhere.

Calculating sufficient statistics

The naive HMM estimation model presented above can be sped up dramatically by replacing the loops over categorical distributions with a single multinomial distribution.

The data are declared as before. The transformed data block computes the sufficient statistics for estimating the transition and emission matrices.

transformed data {
  array[K, K] int<lower=0> trans;
  array[K, V] int<lower=0> emit;
  for (k1 in 1:K) {
    for (k2 in 1:K) {
      trans[k1, k2] = 0;
    }
  }
  for (t in 2:T) {
    trans[z[t - 1], z[t]] += 1;
  }
  for (k in 1:K) {
    for (v in 1:V) {
      emit[k, v] = 0;
    }
  }
  for (t in 1:T) {
    emit[z[t], w[t]] += 1;
  }
}

The likelihood component of the model based on looping over the input is replaced with multinomials as follows.

model {
  // ...
  for (k in 1:K) {
    trans[k] ~ multinomial(theta[k]);
  }
  for (k in 1:K) {
    emit[k] ~ multinomial(phi[k]);
  }
}

In a continuous HMM with normal emission probabilities could be sped up in the same way by computing sufficient statistics.

Analytic posterior

With the Dirichlet-multinomial HMM, the posterior can be computed analytically because the Dirichlet is the conjugate prior to the multinomial. The following example illustrates how a Stan model can define the posterior analytically. This is possible in the Stan language because the model only needs to define the conditional probability of the parameters given the data up to a proportion, which can be done by defining the (unnormalized) joint probability or the (unnormalized) conditional posterior, or anything in between.

The model has the same data and parameters as the previous models, but now computes the posterior Dirichlet parameters in the transformed data block.

transformed data {
  vector<lower=0>[K] alpha_post[K];
  vector<lower=0>[V] beta_post[K];
  for (k in 1:K) {
    alpha_post[k] = alpha;
  }
  for (t in 2:T) {
    alpha_post[z[t - 1], z[t]] += 1;
  }
  for (k in 1:K) {
    beta_post[k] = beta;
  }
  for (t in 1:T) {
    beta_post[z[t], w[t]] += 1;
  }
}

The posterior can now be written analytically as follows.

model {
  for (k in 1:K) {
    theta[k] ~ dirichlet(alpha_post[k]);
  }
  for (k in 1:K) {
    phi[k] ~ dirichlet(beta_post[k]);
  }
}

Semisupervised estimation

HMMs can be estimated in a fully unsupervised fashion without any data for which latent states are known. The resulting posteriors are typically extremely multimodal. An intermediate solution is to use semisupervised estimation, which is based on a combination of supervised and unsupervised data. Implementing this estimation strategy in Stan requires calculating the probability of an output sequence with an unknown state sequence. This is a marginalization problem, and for HMMs, it is computed with the so-called forward algorithm.

In Stan, the forward algorithm is coded as follows. First, two additional data variable are declared for the unsupervised data.

data {
  // ...
  int<lower=1> T_unsup;             // num unsupervised items
  array[T_unsup] int<lower=1, upper=V> u;  // unsup words
  // ...
}

The model for the supervised data does not change; the unsupervised data are handled with the following Stan implementation of the forward algorithm.

model {
  // ...
  array[K] real acc;
  array[T_unsup, K] real gamma;
  for (k in 1:K) {
    gamma[1, k] = log(phi[k, u[1]]);
  }
  for (t in 2:T_unsup) {
    for (k in 1:K) {
      for (j in 1:K) {
        acc[j] = gamma[t - 1, j] + log(theta[j, k]) 
                 + log(phi[k, u[t]]);
      }
      gamma[t, k] = log_sum_exp(acc);
    }
  }
  target += log_sum_exp(gamma[T_unsup]);
}

The forward values gamma[t, k] are defined to be the log marginal probability of the inputs u[1],...,u[t] up to time t and the latent state being equal to k at time t; the previous latent states are marginalized out. The first row of gamma is initialized by setting gamma[1, k] equal to the log probability of latent state k generating the first output u[1]; as before, the probability of the first latent state is not itself modeled. For each subsequent time t and output j, the value acc[j] is set to the probability of the latent state at time t-1 being j, plus the log transition probability from state j at time t-1 to state k at time t, plus the log probability of the output u[t] being generated by state k. The log_sum_exp operation just multiplies the probabilities for each prior state j on the log scale in an arithmetically stable way.

The brackets provide the scope for the local variables acc and gamma; these could have been declared earlier, but it is clearer to keep their declaration near their use.

Predictive inference

Given the transition and emission parameters, \(\theta_{k, k'}\) and \(\phi_{k,v}\) and an observation sequence \(u_1, \dotsc, u_T \in \{ 1, \dotsc, V \}\), the Viterbi (dynamic programming) algorithm computes the state sequence which is most likely to have generated the observed output \(u\).

The Viterbi algorithm can be coded in Stan in the generated quantities block as follows. The predictions here is the most likely state sequence y_star[1], ..., y_star[T_unsup] underlying the array of observations u[1], ..., u[T_unsup]. Because this sequence is determined from the transition probabilities theta and emission probabilities phi, it may be different from sample to sample in the posterior.

generated quantities {
  array[T_unsup] int<lower=1, upper=K> y_star;
  real log_p_y_star;
  {
    array[T_unsup, K] int back_ptr;
    array[T_unsup, K] real best_logp;
    real best_total_logp;
    for (k in 1:K) {
      best_logp[1, k] = log(phi[k, u[1]]);
    }
    for (t in 2:T_unsup) {
      for (k in 1:K) {
        best_logp[t, k] = negative_infinity();
        for (j in 1:K) {
          real logp;
          logp = best_logp[t - 1, j]
                  + log(theta[j, k]) + log(phi[k, u[t]]);
          if (logp > best_logp[t, k]) {
            back_ptr[t, k] = j;
            best_logp[t, k] = logp;
          }
        }
      }
    }
    log_p_y_star = max(best_logp[T_unsup]);
    for (k in 1:K) {
      if (best_logp[T_unsup, k] == log_p_y_star) {
        y_star[T_unsup] = k;
      }
    }
    for (t in 1:(T_unsup - 1)) {
      y_star[T_unsup - t] = back_ptr[T_unsup - t + 1,
                                      y_star[T_unsup - t + 1]];
    }
  }
}

The bracketed block is used to make the three variables back_ptr, best_logp, and best_total_logp local so they will not be output. The variable y_star will hold the label sequence with the highest probability given the input sequence u. Unlike the forward algorithm, where the intermediate quantities were total probability, here they consist of the maximum probability best_logp[t, k] for the sequence up to time t with final output category k for time t, along with a backpointer to the source of the link. Following the backpointers from the best final log probability for the final time t yields the optimal state sequence.

This inference can be run for the same unsupervised outputs u as are used to fit the semisupervised model. The above code can be found in the same model file as the unsupervised fit. This is the Bayesian approach to inference, where the data being reasoned about is used in a semisupervised way to train the model. It is not “cheating” because the underlying states for u are never observed — they are just estimated along with all of the other parameters.

If the outputs u are not used for semisupervised estimation but simply as the basis for prediction, the result is equivalent to what is represented in the BUGS modeling language via the cut operation. That is, the model is fit independently of u, then those parameters used to find the most likely state to have generated u.