4
$\begingroup$

Is there a closed-form solution for the M-step of the expectation maximization (EM) algorithm for a multivariate Gaussian hidden Markov model (HMM) where observations are missing at random? $\newcommand\given[1][]{\:#1\vert\:}$


My HMM notation, largely borrowed from https://www.princeton.edu/~rvan/orf557/hmm080728.pdf, is:

  • $X_t \in \mathrm{E}$ is a discrete hidden state with transition matrix $P$ and initial distribution $\mu$
  • The state space $\mathrm{E} \equiv \left\{1, 2, \ldots, d\right\}$ is finite
  • The observations are $\left(\,Y_t \given X_t = i\,\right)\sim \mathcal{N}\left(\mathbf{m}_i, \mathbf{v}_i\right)$
  • If the observations are scalars, $\mathbf{m}$ and $\mathbf{v}$ are simply vectors of means and variances (one for each hidden state)
  • If the observations are multivariate normal, say $Y_t \in \mathbb{R}^k$, each $\mathbf{m}_i$ is a $k$-vector of means and each $\mathbf{v}_i$ is a $k$-by-$k$ covariance matrix

I'm interested in estimating $P$, $\mu$, $\mathbf{m}$ and $\mathbf{v}$ using the EM algorithm. Let $\mathbf{y}$ denote an $n$-by-$k$ matrix of observed data. Let $\pi_{it}$ denote the posterior probability over hidden state $i$ for observation $t$, computed under the previous iteration's parameters, i.e. $\pi_{it} = \Pr[X_t = i \given \mathbf{y}]$.

The notes I linked to earlier give the M step when $k$ = 1, i.e. scalar Gaussian observations (section 6.2, page 87):

$$\begin{equation*} \mathbf{m}_i = \frac{\sum^n_{t=1} \mathbf{y}_t \, \pi_{it}}{\sum^n_{t=1} \pi_{it}} \end{equation*}$$ $$\begin{equation*} \mathbf{v}_i = \frac{\sum^n_{t=1} \left(\mathbf{y}_t - \mathbf{m}_i\right)^2 \, \pi_{it}}{\sum^n_{t=1} \pi_{it}} \end{equation*}$$

I'm interested mainly in the updates to $\mathbf{m}$ and $\mathbf{v}$ so I'll omit the updates to $P$ and $\mu$ for clarity.

The M-step updates for multivariate Gaussian observations are quite similar. The update to $\mathbf{m}_i$ is identical (the only difference being that $y_t$ is now a $k$-vector), and the update for the covariance matrices is:

$$\begin{equation*} \mathbf{v}_i = \frac{\sum^n_{t=1} \left(\mathbf{y}_t - \mathbf{m}_i\right) \, \left(\mathbf{y}_t - \mathbf{m}_i\right)^{\intercal} \, \pi_{it}}{\sum^n_{t=1} \pi_{it}} \end{equation*}$$

This equation appears on slide 9 of http://personal.ee.surrey.ac.uk/Personal/P.Jackson/tutorial/hmm_tut4.pdf, for example, with slightly different notation.


What is the M-step if elements of $\mathbf{y}$ (the $n$-by-$k$ matrix of observations) are missing at random? My guess is to first define a matrix $\tilde{\mathbf{y}}$ containing 1s where $\mathbf{y}$ is observed and 0s where it is missing (i.e. a matrix of indicators for non-missing data). I then recode the missing observations in $\mathbf{y}$ as 0.

My updates (after the recoding described above) are:

$$\begin{equation*} \mathbf{m}_i = \frac{\sum^n_{t=1} \mathbf{y}_t \, \pi_{it}}{\sum^n_{t=1} \tilde{\mathbf{y}}_t \, \pi_{it}} \end{equation*}$$ $$\begin{equation*} \mathbf{v}_i = \frac{\sum^n_{t=1} \left(\mathbf{y}_t - \mathbf{m}_i\right) \, \left(\mathbf{y}_t - \mathbf{m}_i\right)^{\intercal} \, \pi_{it}}{\sum^n_{t=1} \tilde{\mathbf{y}}_t \, \tilde{\mathbf{y}}_t^{\intercal} \pi_{it}} \end{equation*}$$

In the equation above, both the numerator and denominator are $k$-by-$k$ matrices. The fraction should be read as an element-wise division.

Is this correct -- is it even guaranteed to produce a positive definite covariance matrix? Does a correct analytical update exist? Any pointers to papers or lecture notes that discuss multivariate Gaussian HMMs with observations missing at random would be much appreciated.


A relevant paragraph from https://en.wikipedia.org/wiki/Estimation_of_covariance_matrices:

In the general case, the unbiased estimate of the covariance matrix provides an acceptable estimate when the data vectors in the observed data set are all complete: that is they contain no missing elements. One approach to estimating the covariance matrix is to treat the estimation of each variance or pairwise covariance separately, and to use all the observations for which both variables have valid values. Assuming the missing data are missing at random this results in an estimate for the covariance matrix which is unbiased. However, for many applications this estimate may not be acceptable because the estimated covariance matrix is not guaranteed to be positive semi-definite. This could lead to estimated correlations having absolute values which are greater than one, and/or a non-invertible covariance matrix.

Do any references discuss how to deal with this problem in the context of multivariate Gaussian HMMs?


Here's a little R example showing covariance matrix estimates with negative eigenvalues:

library(mixtools)  # For rmvnorm

get_sigma_hat <- function(n_obs, sigma, pr_missing) {
    ## Generate random multivariate data y with observations missing at random
    ## Return estimated variance-covariance matrix
    stopifnot(0 <= pr_missing && pr_missing < 1)
    stopifnot(all(eigen(sigma)$values > 0.0))
    stopifnot(identical(sigma, t(sigma)))
    y <- rmvnorm(n_obs, sigma=sigma)
    stopifnot(all(dim(y) == c(n_obs, nrow(sigma))))
    y_non_missing <- matrix(1*(runif(n_obs * nrow(sigma)) > pr_missing), n_obs, nrow(sigma))
    mean(y_non_missing)  # Should be close to (1 - pr_missing)
    observed_y <- ifelse(y_non_missing, y, 0)  # Careful, NAs recoded to zero
    sigma_hat <- t(observed_y) %*% observed_y / (t(y_non_missing) %*% y_non_missing)
    return(list(sigma_hat=sigma_hat,
                sigma_hat_positive_definite=all(eigen(sigma_hat)$values > 0.0)))
}

set.seed(543987)

sigma <- matrix(c(1.0, 0.5, 0.6,
                  0.5, 2.0, 0.7,
                  0.6, 0.7, 3.0), 3, 3)  # True variance-covariance matrix for observations in R^3
n_reps <- 1000
reps <- replicate(n_reps, get_sigma_hat(n_obs=500, sigma=sigma, pr_missing=0.10), simplify=FALSE)

mean(sapply(reps, function(x) x$sigma_hat_positive_definite))  # 1.0 i.e. all positive definite
mean(sapply(reps, function(x) x$sigma_hat[1, 1]))  # Close to sigma[1, 1]
mean(sapply(reps, function(x) x$sigma_hat[3, 2]))  # Close to sigma[3, 2]

reps <- replicate(n_reps, get_sigma_hat(n_obs=400, sigma=sigma, pr_missing=0.80), simplify=FALSE)
mean(sapply(reps, function(x) x$sigma_hat_positive_definite))  #  Some negative eigenvalues!
sigma_hat <- Filter(function(x) !x$sigma_hat_positive_definite, reps)[[1]]$sigma_hat
eigen(sigma_hat)$values  # Third eigenvalue is negative
cov2cor(sigma_hat)  # ...but all correlations are in [-1, 1]

another_sigma_hat <- Filter(function(x) !x$sigma_hat_positive_definite, reps)[[2]]$sigma_hat
eigen(another_sigma_hat)$values  # Again, third eigenvalue is negative
cov2cor(another_sigma_hat)  # This one has some correlations above 1

The code simulates a simplified version of the problem I describe above, where the $Y$s have mean zero and there are no weights $\pi$.

$\endgroup$
1
$\begingroup$

Typically, when using the EM algorithm you would consider missing observations as hidden variables that would be estimated. Instead of setting missing observations to 0, you would have to estimate them in the M-step.

However, for the case of HMM, given the current state, each observation is independent. Because of this, it turns out that we do not need to estimate missing observations, as long as we estimate all the hidden states.

The update equations for the M-step turn out to be what you've guessed! This is explained in detail in the below paper:

Roberta Paroli and Luigi Spezia, "Parameter estimation of Gaussian hidden Markov models when missing observations occur", Metron-International Journal of Statistics, 2002 (https://www.researchgate.net/publication/5182349_Parameter_estimation_of_Gaussian_hidden_Markov_models_when_missing_observations_occur)

$\endgroup$
  • $\begingroup$ Thank you for your answer. The paper you posted is almost what I'm looking for, but I'd like something that explicitly discusses multivariate observations. If I am reading them correctly, Paroli and Spezia only consider scalar observations, right? $\endgroup$ – Adrian Aug 13 '16 at 8:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.