HMM经典介绍论文【Rabiner 1989】翻译（八）——学习问题

3.3 问题3的求解（学习问题）

HMM的第三个问题——学习问题是最困难的一个，需要通过最大化观测序列的概率来调整模型参数$(A, B, \pi)$。暂时还没有解析法来解决这个问题。事实上，给定一个有限观测序列作为训练数据，并不存在最优方法得到模型参数。但是，我们可以通过迭代法比如Baum-Welch算法（其实就是EM算法），或者使用梯度法，选择使得$P(O|\lambda)$局部最大的$\lambda = (A, B, \pi)$。这一节，我们讨论选择模型参数的一个迭代方法，这个方法主要基于Baum和他同事的经典工作。

为了描述HMM参数估计过程，首先定义给定模型和观测序列的条件下，时刻$t$的状态为$S_i$且$t+1$的状态为$S_j$的概率为$\xi_t(i, j)$：

$$\xi_t(i, j) = P(q_t = S_i, q_{t+1}= S_j | O, \lambda). \tag{36}$$

(36)的计算结构如图6所示。显然，根据前向变量和后向变量的定义，我们可以得到

$$\begin{align} \xi_t(i, j) &= \frac{\alpha_t(i) a_{ij} b_j(O_{t+1}) \beta_{t+1}(j)}{P(O|\lambda)}\nonumber \\\\ &= \frac{\alpha_t(i) a_{ij} b_j(O_{t+1}) \beta_{t+1}(j)}{\sum_{i=1}^N \sum_{j=1}^N \alpha_t(i) a_{ij} b_j(O_{t+1}) \beta_{t+1}(j) }. \tag{37} \end{align}$$

我们前面已经定义了给定模型和状态序列的条件下，时刻$t$的状态为$S_i$的概率为$\gamma_t(i)$。通过在$j$上求和，有

$$\gamma_t(i) = \sum_{j=1}^N \xi_t(i, j). \tag{38}$$

如果对$\gamma_t(i)$在$t$上求和，得到的和可以理解为状态$S_i$出现次数的期望值，如果去掉$t=T$项，那么得到的和可以理解为从状态$S_i$进行转移的期望次数。类似地，在$t$上（从$t=1$到$t=T-1$）对$\xi_t(i,j)$进行求和得到的结果可以理解为从状态$S_i$转移到$S_j$次数的期望值。即

$$\sum_{t=1}^{T-1}\gamma_t(i) = 从S_i转移的期望次数 \tag{39a}$$

$$\sum{t=1}^{T-1}\xi_t(i, j)=从S_i转移到S_j的期望次数. \tag{39b}$$

利用上面的公式，我们可以给出估计HMM参数的方法。对$\pi, A, B$合理的估计可以是

$$\bar {\pi_i} = 在t=1处于状态S_i的概率 = \gamma_t(i) \tag{40a}$$

$$\bar{a_{ij}} = \frac{从S_i转移到S_j的期望次数}{从S_i转移的期望次数} = \frac{\sum_{t=1}^{T-1}\xi_t(i, j)}{\sum_{t=1}^{T-1}\gamma_t(i)} \tag{40b}$$

$$\bar{b_j}(k) = \frac{处于状态S_j且观测值为v_k的期望次数}{状态S_j的期望次数} = \frac{\sum_{t=1, O_t=v_k}^T \gamma_t(j)} {\sum_{t=1}^{T}\gamma_t(j)}. \tag{40c}$$

如果我们定义当前模型为$\lambda = (A, B, \pi)$，并且用它来计算(40a)-(40c)右边的表达式，然后定义新模型为$\bar {\lambda} = (\bar A, \bar B, \bar{\pi})$，即(40a)-(40c)左边的部分，Baum和他同事证明了1）初始模型$\lambda$定义了似然函数的临界点；2）模型$\bar{\lambda}$比$\lambda$更好，因为$P(O|\bar{\lambda}) > P(O|\lambda)$，即我们找到了新模型$\bar{\lambda}$，更有可能生成观测序列。

基于上面的步骤，如果我们迭代地更新$\lambda$为$\bar {\lambda}$，那么$P(O|\lambda)$会不断增加，知道到达某个极限值。这个过程的结果称为HMM的极大似然估计。注意，前向-后向算法得到的是局部最大值，在大多数问题中，优化表面是非常复杂的并且有很多局部最大值。

更新公式(40a)-(40c)也可以通过最大化Baum的辅助函数推导出来:

$$Q(\lambda, \bar{\lambda}) = \sum_{Q} P(Q|O, \lambda) \log [P(O, Q| \bar{\lambda})] \tag{41}.$$

Baum和他同事证明了最大化$Q(\lambda, \bar{\lambda})$可以增加似然，即

$$\max_{\bar \lambda} [Q[\lambda, \bar{\lambda}]] => P(O|\bar{\lambda}) \geq P(O|\lambda). \tag{42}$$

最终似然函数会收敛到一个临界点。

上面的估计公式可以理解为EM算法的具体实现，其中E步是计算辅助函数$Q(\lambda, \bar{\lambda})$，M步是最大化$Q(\lambda, \bar{\lambda})$。于是，Baum-Welch估计方程本质上等价于EM算法的具体实现。

估计过程重要的一个方面是，HMM参数的随机约束，即

$$\sum_{i=1}^N \bar{\pi}_i = 1 \tag{43a}$$

$$\sum_{j=1}^N {\bar{a}}_{ij} = 1, \quad 1\leq i \leq N \tag{43b}$$

$$\sum_{k=1}^M {\bar{b}}_j(k) = 1, \quad 1 \leq j \leq N \tag{43c}$$

在每次迭代时是自动满足的。通过把参数估计问题看作一个带约束的优化问题，可以用拉格朗日乘子法来得到使$P$（$P=P(O|\lambda)$）最大的$\pi_i, a_{ij}, b_j(k)$值。通过标准拉格朗日优化，可以知道当

$$\pi_i = \frac{\pi_i \frac{\partial P}{\partial \pi_i}}{\sum_{k=1}^N \frac{\partial P}{\partial \pi_k} } \tag{44a}$$

$$a_{ij}= \frac{a_{ij} \frac{\partial P}{\partial a_{ij}} }{ \sum_{k=1}^N a_{ik} \frac{\partial P}{\partial a_{ik}} } \tag{44b}$$

$$b_j(k) = \frac{ b_j(k)\frac{\partial P}{\partial b_j(k)} }{ \sum_{l=1}^M b_j(l) \frac{\partial P}{\partial b_j(l)} } \tag{44c}$$

时，$P$是最大的。通过计算，可以发现公式(44)右边的表达式其实就是(40)右边的表达式。所以，估计公式实际上就是$P$的临界点。

最后，我们注意到由于整个问题可以看作一个优化问题，标准的梯度算法也可以用于求模型参数的最优值。这种方法的结果可以和标准估计过程媲美。