Gaussian Emissions in the Baum-Welch M-Step
In a Gaussian-emission hidden Markov model, each hidden state has its own mean and covariance. The observations do not reveal the states, so the emission parameters are not updated by assigning each observation to its most likely state. Instead, each observation contributes to each state's update in proportion to that state's posterior probability after seeing the full sequence. This produces weighted means and weighted covariances.
Let \((\Omega,\mathcal{F},\mathbb{P})\) support a discrete-time hidden Markov model \((Z_t,Y_t)_{t=1}^T\). The hidden state \(Z_t\) takes values in the finite set \(S=\{1,\ldots,K\}\), and the observation \(Y_t\) takes values in \(\mathbb{R}^d\). The state process has initial distribution \(\mu\) and row-stochastic transition matrix \(P=(p_{ij})\), where \(p_{ij}=\mathbb{P}(Z_{t+1}=j\mid Z_t=i)\). Given the state sequence, observations are conditionally independent, with \(Y_t\mid Z_t=i\) having a Gaussian distribution with mean \(\boldsymbol{\mu}_i\) and covariance matrix \(\boldsymbol{\Sigma}_i\). Each covariance matrix is required to be symmetric positive definite. Write \(\mathcal{F}_t=\sigma(Z_1,\ldots,Z_t,Y_1,\ldots,Y_t)\) for the joint history filtration.
For the observed sequence \(y_{1:T}=(y_1,\ldots,y_T)\), the smoothed state probability is \(\gamma_t(i)=\mathbb{P}(Z_t=i\mid Y_1=y_1,\ldots,Y_T=y_T)\). This conditional probability is well-defined when the observed sequence has positive probability in the discrete case, or when the observed-data density used for conditioning is positive in the continuous case. The forward-backward calculation supplies these probabilities under the current parameter values. At each time \(t\), \(\gamma_t(i)\ge 0\) and \(\sum_{i=1}^K\gamma_t(i)=1\).
The mean and covariance updates below maximize the emission part of the Baum-Welch auxiliary function: the conditional expectation, given all observations under the current parameters, of the complete-data log likelihood. In this maximization the probabilities \(\gamma_t(i)\) are held fixed. The updates do not by themselves update the initial distribution or transition matrix.
Why the Updates Are Weighted
For a \(d\)-dimensional Gaussian, the log density at \(y_t\) in state \(i\) is \(-\tfrac{d}{2}\log(2\pi)-\tfrac{1}{2}\log\det\boldsymbol{\Sigma}_i-\tfrac{1}{2}(y_t-\boldsymbol{\mu}_i)^\mathsf{T}\boldsymbol{\Sigma}_i^{-1}(y_t-\boldsymbol{\mu}_i)\). The expected complete-data emission log likelihood for state \(i\), with terms not involving its parameters omitted, is therefore
The factor \(\gamma_t(i)\) is the expected indicator that the hidden state at time \(t\) equals \(i\). Thus a time point with a larger smoothed probability contributes more to that state's fit, while one with probability zero contributes nothing.
Proof. Fix state \(i\), and abbreviate its weights by \(w_t=\gamma_t(i)\), its total weight by \(W\), its weighted mean by \(\bar y\), and its weighted covariance by \(C\). First, for every vector \(\mu\), expand \(y_t-\mu=(y_t-\bar y)+(\bar y-\mu)\) inside the weighted sum of outer products. The two cross terms vanish because \(\sum_t w_t(y_t-\bar y)=0\). Hence
For any positive-definite \(\Sigma\), the quadratic term in \(Q_i\) is the trace of \(\Sigma^{-1}\) times this matrix. Its value at \(\mu\) exceeds its value at \(\bar y\) by \(W(\mu-\bar y)^\mathsf{T}\Sigma^{-1}(\mu-\bar y)\), which is nonnegative and is zero only when \(\mu=\bar y\). Thus the unique maximizing mean for each fixed covariance is \(\bar y\).
Now suppose \(C\) is positive definite and set \(\mu=\bar y\). The remaining quantity to minimize over positive-definite \(\Sigma\) is \(W(\log\det\Sigma+\operatorname{tr}(\Sigma^{-1}C))/2\). Let \(\lambda_1,\ldots,\lambda_d\) be the positive eigenvalues of \(C^{-1/2}\Sigma C^{-1/2}\). Determinant and trace identities give \(\log\det\Sigma+\operatorname{tr}(\Sigma^{-1}C)=\log\det C+\sum_{r=1}^d(\log\lambda_r+1/\lambda_r)\). For \(x>0\), the function \(\log x+1/x\) has derivative \((x-1)/x^2\), so it decreases up to \(1\), increases after \(1\), and has its unique minimum \(1\) at \(x=1\). The sum is uniquely minimized when every \(\lambda_r=1\), which means \(\Sigma=C\). Together with the mean argument, this proves the claimed unique joint maximizer.
If \(W=0\), all weights are zero and \(Q_i\) is constant in the emission parameters. If \(W>0\) and \(C\) is singular, choose a unit vector \(v\) with \(v^\mathsf{T}Cv=0\), and use covariance matrices that have variance \(\varepsilon\) in direction \(v\) and fixed positive variances on its orthogonal complement. The weighted residual quadratic term has no contribution in direction \(v\), while \(-W\log\det\Sigma/2\) tends to \(+\infty\) as \(\varepsilon\) decreases to zero. Thus no positive-definite covariance maximizes \(Q_i\). \(\square\)
This result explains both the familiar update formulas and an important limitation. They are ordinary maximum-likelihood M-step updates when the state has positive total weight and its weighted observations have nonzero variation in every direction. A singular covariance cannot be repaired by dividing by \(W_i\); it requires an explicit modeling choice, such as a covariance constraint or regularization.
Worked Calculations
Worked Example: Two Scalar State Means and Variances
Consider observations \(y_1=0,y_2=2,y_3=4,y_4=6\) and two hidden states \(S=\{1,2\}\). Let the transition matrix be \(P=\begin{bmatrix}\tfrac12&\tfrac12\\ \tfrac13&\tfrac23\end{bmatrix}\); its first row sums to \(\tfrac12+\tfrac12=1\), and its second row sums to \(\tfrac13+\tfrac23=1\). Take initial distribution \(\mu=(\tfrac12,\tfrac12)\), whose entries sum to \(1\). Suppose the forward-backward E-step gives the following smoothed probabilities, with the second state's probability at each time equal to one minus the first state's probability.
| Observation | \(\gamma_t(1)\) | \(\gamma_t(2)\) |
|---|---|---|
| 0 | \(\tfrac34\) | \(\tfrac14\) |
| 2 | \(\tfrac12\) | \(\tfrac12\) |
| 4 | \(\tfrac14\) | \(\tfrac34\) |
| 6 | \(\tfrac12\) | \(\tfrac12\) |
For state 1, the total weight is \(W_1=\tfrac34+\tfrac12+\tfrac14+\tfrac12=2\). Its weighted sum is \(\tfrac34(0)+\tfrac12(2)+\tfrac14(4)+\tfrac12(6)=0+1+1+3=5\), so \(\bar y_1=5/2\). The weighted squared deviations are \(\tfrac34(0-\tfrac52)^2=\tfrac{75}{16}\), \(\tfrac12(2-\tfrac52)^2=\tfrac18=\tfrac{2}{16}\), \(\tfrac14(4-\tfrac52)^2=\tfrac{9}{16}\), and \(\tfrac12(6-\tfrac52)^2=\tfrac{49}{8}=\tfrac{98}{16}\). Their sum is \(184/16=23/2\), so the updated variance is \((23/2)/2=23/4\).
For state 2, \(W_2=\tfrac14+\tfrac12+\tfrac34+\tfrac12=2\). Its weighted sum is \(\tfrac14(0)+\tfrac12(2)+\tfrac34(4)+\tfrac12(6)=0+1+3+3=7\), giving \(\bar y_2=7/2\). The squared-deviation sum is \(\tfrac14(0-\tfrac72)^2=\tfrac{49}{16}\), \(\tfrac12(2-\tfrac72)^2=\tfrac{18}{16}\), \(\tfrac34(4-\tfrac72)^2=\tfrac{3}{16}\), and \(\tfrac12(6-\tfrac72)^2=\tfrac{50}{16}\). The total is \(120/16=15/2\), so the updated variance is \((15/2)/2=15/4\).
Check: For state 1, the weighted raw second moment divided by its weight is \(\{\tfrac12(2^2)+\tfrac14(4^2)+\tfrac12(6^2)\}/2=(2+4+18)/2=12\). Subtracting the squared mean gives \(12-(5/2)^2=23/4\), as above. For state 2 it is \(\{\tfrac12(2^2)+\tfrac34(4^2)+\tfrac12(6^2)\}/2=(2+12+18)/2=16\); subtracting \((7/2)^2\) gives \(15/4\).
Worked Example: A Two-Dimensional State Covariance
Now let \(Y_t\in\mathbb{R}^2\), with observations \(y_1=(0,0)^\mathsf{T}\), \(y_2=(2,0)^\mathsf{T}\), and \(y_3=(0,2)^\mathsf{T}\). Use two states with transition matrix \(P=\begin{bmatrix}\tfrac12&\tfrac12\\ \tfrac12&\tfrac12\end{bmatrix}\); each row sums to \(\tfrac12+\tfrac12=1\). Let the initial distribution be \((\tfrac12,\tfrac12)\), which sums to \(1\), and suppose \(\gamma_t(1)=\tfrac12\) for all three observations. Each complementary state-2 probability is also \(\tfrac12\).
For state 1, \(W_1=\tfrac12+\tfrac12+\tfrac12=\tfrac32\), and the weighted observation sum is \(\tfrac12(0,0)^\mathsf{T}+\tfrac12(2,0)^\mathsf{T}+\tfrac12(0,2)^\mathsf{T}=(1,1)^\mathsf{T}\). Therefore \(\bar y_1=(2/3,2/3)^\mathsf{T}\). The three weighted outer products of the centered observations are
Their sum is \(\begin{bmatrix}\tfrac43&-\tfrac23\\-\tfrac23&\tfrac43\end{bmatrix}\). Dividing by \(W_1=\tfrac32\) gives \(\boldsymbol{\Sigma}_1=\begin{bmatrix}\tfrac89&-\tfrac49\\-\tfrac49&\tfrac89\end{bmatrix}\). Its determinant is \(\tfrac{64}{81}-\tfrac{16}{81}=\tfrac{48}{81}=\tfrac{16}{27}>0\), and its top-left entry is \(\tfrac89>0\), so this symmetric matrix is positive definite.
Check: The weighted centered vectors sum to \(\tfrac12(-\tfrac23,-\tfrac23)^\mathsf{T}+\tfrac12(\tfrac43,-\tfrac23)^\mathsf{T}+\tfrac12(-\tfrac23,\tfrac43)^\mathsf{T}=(0,0)^\mathsf{T}\), as required for the weighted mean. The covariance is symmetric, and the positive determinant and positive first diagonal entry confirm that it is admissible.
Worked Example: When the Scalar Variance Update Is Degenerate
Take a one-state model with \(S=\{1\}\), \(P=[1]\), and initial distribution \([1]\). The only transition row sums to \(1\), and every smoothed state probability is \(\gamma_t(1)=1\). Suppose the three scalar observations are \(2,2,2\). Then \(W_1=1+1+1=3\), the weighted mean is \((2+2+2)/3=2\), and the weighted variance is \(\{(2-2)^2+(2-2)^2+(2-2)^2\}/3=0\).
A variance of zero is not a valid Gaussian covariance under the positive-variance model. To see why the ordinary maximization has no solution, let the candidate variance be \(v>0\). At the mean \(2\), the weighted log-likelihood contribution is \(-\tfrac32\log(2\pi v)\), since all residuals are zero. As \(v\) decreases toward zero, this quantity increases without bound. Thus the formal weighted-variance formula reaches a boundary value, but there is no maximizing positive variance. A fitted model must impose an additional constraint or regularization to prevent this degeneracy.
Common Mistakes and Practical Interpretation
Common mistake: dividing by the number of observations instead of the state's total weight. In the first example, state 1 has weighted squared-deviation sum \(23/2\), total weight \(2\), and four observations. Dividing by the observation count would give \((23/2)/4=23/8\), which is not the M-step variance. Its mean would also be wrong: dividing its weighted sum \(5\) by \(4\) gives \(5/4\), rather than the correct \(5/2\). The correction is to divide both weighted sums by \(W_1=\sum_t\gamma_t(1)=2\), not by \(T=4\).
A related error is replacing smoothed probabilities by a Viterbi path. The path assigns each time to one state, while \(\gamma_t(i)\) retains uncertainty and uses the full observation sequence. These lead to different parameter updates in general. Also, the smoothed weights are computed using the current parameters and then held fixed during this M-step; they are not recomputed after each individual observation or after each parameter update.
For several independent observation sequences, sum the weights and weighted sums over both sequence and time indices before taking the mean and covariance. The formulas remain the same: divide the aggregate weighted sum by the aggregate state weight, and divide the aggregate centered outer-product sum by that same weight. In practice, inspect each state's total weight and covariance before interpreting the updated Gaussian. A very small weight can make estimates unstable, and a singular weighted covariance means the unconstrained Gaussian update is not an admissible positive-definite covariance.
Check Your Understanding
For each question, keep the state probabilities fixed as E-step weights and check the arithmetic against the defining weighted sums.
- For observations \(1,3\) and state weights \(\tfrac14,\tfrac34\), calculate the weighted mean and variance for that state.
- Why does the Gaussian emission update use \(\gamma_t(i)\), rather than a single most likely hidden state at each time?
- In the two-dimensional example, what feature of the weighted data makes the fitted covariance positive definite?
- What happens to the emission auxiliary function if \(W_i=0\)?
- For three identical observations in one state, explain why the weighted variance formula does not produce a valid positive-variance Gaussian M-step.