统计学习方法 隐马尔可夫模型

读李航的《统计学习方法》时，关于隐马尔可夫模型的笔记

隐马尔可夫模型（hidden Markov model, HMM）是可用于标注问题的统计学习模型，属于生成模型。

基本概念

隐马尔可夫模型：在 Markov 链的基础上，随机生成不可观测的状态序列，再由每个状态生成一个观测，产生观测序列，序列的每个位置可看成一个时刻。

隐马尔可夫模型由初始概率分布 $\pi$ 、状态转移概率分布 $A$ 以及观测概率分布 $B$ 确定，符号定义如下：

设 $Q$ 是所有可能状态的集合，$V$ 是所有可能观测的集合：
Q = {   q 1 ,   q 2 ,   ⋯   ,   q N   } , V = {   v 1 ,   v 2 ,   ⋯   ,   v M   } Q=\set{q_1,\,q_2,\,\cdots,\,q_N},\quad V=\set{v_1,\,v_2,\,\cdots,\,v_M} Q={q1​,q2​,⋯,qN​},V={v1​,v2​,⋯,vM​}
$N$ 是可能的状态数，$M$ 是可能的观测数。$I$ 是长度为 $T$ 的状态序列，$O$ 是对应的观测序列：
I = {   i 1 ,   i 2 ,   ⋯   ,   i T   } , O = ( o 1 ,   o 2 ,   ⋯   ,   o T ) I=\set{i_1,\,i_2,\,\cdots,\,i_T},\quad O=(o_1,\,o_2,\,\cdots,\,o_T) I={i1​,i2​,⋯,iT​},O=(o1​,o2​,⋯,oT​)
$A$ 是状态转移概率矩阵，跟 Markov 链的定义是一样的：
$$A=[a_{ij}{]}_{N\times N}$$
其中 $a_{ij}$ 代表从状态 $q_i$ 转移到状态 $q_j$ 的概率：
$$a_{ij}=P(i_{t+1}=q_j|i_t=q_i),i=1,2,\cdots ,N;j=1,2,\cdots ,N$$
$B$ 是观测概率矩阵：
$$B=[b_j(k){]}_{N\times M}$$
其中 $b_j(k)$ 代表状态 $q_j$ 产生观测结果 $v_k$ 的概率：
$$b_j(k)=P(o_t=v_k|i_t=q_j),k=1,2,\cdots ,M;j=1,2,\cdots ,N;$$
$\pi$ 是初始概率向量：
$$\pi =({\pi }_i)$$
其中 ${\pi }_i$ 代表初始时刻（$t=1$）处于状态 $q_i$ 的概率：
$${\pi }_i=P(i_1=q_i),i=1,2,\cdots ,N$$
即隐马尔可夫模型 $\lambda$ 可以表示为三元组：
$$\lambda =(\pi ,A,B)$$
隐马尔可夫模型的基本假设：有两个基本假设

1. 齐次 Markov 性：任意时刻 $t$ 的状态只依赖于前一时刻状态，而与其他时刻的状态及观测无关，与时刻 $t$ 也无关：

$$P(i_t|i_{t-1},o_{t-1},\cdots ,i_1,o_1)=P(i_t|i_{t-1}),i=1,2,\cdots ,T$$

2. 观测独立性：任意时刻的观测只依赖于该时刻所处的 Markov 链上的状态，与其他任何观测及状态无关：

$$P(o_t|i_T,o_T,\cdots ,i_{t+1},o_{t+1},i_{t-1},o_{t-1},\cdots ,i_1,o_1)=P(o_t|i_t)$$

标注问题：对于标注问题，我们可以认为标注问题的数据是由 HMM 生成的，并且状态对应着标记，只要通过学习和预测算法就可以进行标注。

例（盒子和球模型）：有四个盒子，每个盒子里装有若干红、白两种颜色的球：
$$\begin{array}{ccccc}& \text{box}1& \text{box}2& \text{box}3& \text{box}4\\ \text{red ball}& 5& 3& 6& 8\\ \text{white ball}& 5& 7& 4& 2\end{array}$$
按照以下抽球方式，产生一个颜色的观测序列：

• 开始时，从四个盒子中等概率地选择一个，有放回地取出一个球并记录颜色；
• 然后转移到下一个盒子：如果是盒子 1，则必转移到盒子 2；如果是盒子 2 或盒子 3，则分别以 0.4 和 0.6 的概率转移到左边和右边的盒子；如果是盒子 4，则各以 0.5 的概率留在盒子 4 或转移到盒子 3；
• 重复，比如进行 5 次，得到一个观测序列：

O = {   r ,   r ,   w ,   w ,   r   } O=\set{r,\,r,\,w,\,w,\,r} O={r,r,w,w,r}

这个过程中，观察者只能观测到球的颜色，而并不知道球是从哪个盒子里取出的。因此盒子对应状态，即：
Q = {   1 ,   2 ,   3 ,   4   } , N = 4 Q=\set{1,\,2,\,3,\,4},\quad N=4 Q={1,2,3,4},N=4
球的颜色为观测集合：
V = {   r ,   w   } , M = 2 V=\set{r,\,w},\quad M=2 V={r,w},M=2
初始概率分布为：
$$\pi =(0.25,0.25,0.25,0.25{)}^T$$
状态转移概率分布为：
$$A=\left[\begin{array}{cccc}0& 1& 0& 0\\ 0.4& 0& 0.6& 0\\ 0& 0.4& 0& 0.6\\ 0& 0& 0.5& 0.5\end{array}\right]$$
观测概率分布为：
$$B=\left[\begin{array}{cc}0.5& 0.5\\ 0.3& 0.7\\ 0.6& 0.4\\ 0.8& 0.2\end{array}\right]$$
三个基本问题：给定 HMM，随机生成一个观测序列是很容易的。但还有以下三个不那么明显的问题：

1. 概率计算问题：给定模型和观测序列，计算该模型产生该观测序列的概率 $P(O|\lambda )$ ；
2. 学习问题：给定观测序列，估计模型参数，使得$P(O|\lambda )$ 最大，即用极大似然估计的方法估计参数；
3. 预测问题：也称解码问题，已知模型参数和观测序列，求使得条件概率 $P(I|O)$ 最大的状态序列 $I=(i_1,i_2,\cdots ,i_T)$ ，即给定观测序列，求出最有可能的对应的状态序列。

概率计算问题

概率计算问题是指，给定模型和观测序列，计算该模型产生该观测序列的概率。

直接计算法

算法：直接枚举所有可能的长度为 $T$ 的状态序列，然后求出各状态序列产生该观测序列的概率，最后求和。

出现状态序列 $I=(i_1,i_2,\cdots ,i_T)$ 的概率为：
$$P(I|\lambda )={\pi }_{i_1}{a}_{i_1{i}_2}{a}_{i_2{i}_3}\cdots a_{i_{T-1}{i}_T}$$
对于该状态序列，产生观测 $O=(o_1,o_2,\cdots o_T)$ 的概率为：
$$P(O|I,\lambda )=b_{i_1}(o_1)b_{i_2}(o_2)\cdots b_{i_T}(o_T)$$
则该观测序列出现的总的概率为：
$$\begin{array}{rl}P(O|\lambda )=& \sum _{I}P(O|I,\lambda )P(I|\lambda )\\ =& \sum _{i_1,i_2,\cdots ,i_T}{\pi }_{i_1}{b}_{i_1}(o_1)a_{i_1{i}_2}{b}_{i_2}(o_2)a_{i_2{i}_3}\cdots a_{i_{T-1}{i}_T}{b}_{i_T}(o_T)\end{array}$$
但是这个算法的时间复杂度为 $O(T{N}^T)$ ，是指数级别的。

前向算法

前向概率：给定模型 $\lambda$ ，定义前向概率为，到时刻 $t$ 时，前面出现的观测序列为 $O[0,t]=(0_1,o_2,\cdots ,o_t)$ ，而此时状态正好为 $q_i$ 的概率，记为：
$${\alpha }_t(i)=P(o_1,o_2,\cdots ,o_t,i_t=q_t|\lambda )$$
可以想到，如果知道了前一时刻 $t-1$ 所有状态的前向概率，则可以得到下一时刻 $t$ 任一状态 $q_i$ 的前向概率，只需要前向概率乘以状态转移概率，再乘以产生 $o_t$ 的概率，最后求和即可：

并且，知道了 $T$ 时刻所有状态的前向概率以后，将所有状态的前向概率求和即可得到

算法：观测序列概率的前向算法

• 输入：HMM $\lambda$ ，观测序列 $O$ ；
• 输出：产生该观测序列的概率 $P(O|\lambda )$ ；

1. 初值：

$${\alpha }_1(i)={\pi }_i{b}_i(o_1),i=1,2,\cdots ,N$$

2. 递推：对于 $t=1,2,\cdots ,T-1$ ，递推地计算后一时刻的前向概率：

$${\alpha }_{t+1}(i)=b_i(o_{t+1})\sum _{j=1}^N{\alpha }_t(j)a_{ji}$$

3. 终止：求 $T$ 时刻所有状态的前向概率的和：

$$P(O|\lambda )=\sum _1^N{\alpha }_T(i)$$

整个计算过程即如下的 DAG（有点像 Bellman-Ford 的计算路径），想象一下，其实直接算法的计算路径就如同一颗以状态数 $N$ 为子节点数的 $N$ 叉树，而前向算法则是将每一层重复节点进行了合并（就如同 DP 优化暴力算法一样），因此得到的时间复杂度为 $O(N^{2}T)$ ：

例（盒子和球模型）：状态集合为 Q = {   1 ,   2 ,   3   } Q=\set{1,\,2,\,3} Q={1,2,3}，观测集合为 V = {   r ,   w   } V=\set{r,\,w} V={r,w} ，
$$A=\left[\begin{array}{ccc}0.5& 0.2& 0.3\\ 0.3& 0.5& 0.2\\ 0.2& 0.3& 0.5\end{array}\right]B=\left[\begin{array}{cc}0.5& 0.5\\ 0.4& 0.6\\ 0.7& 0.3\end{array}\right]\pi =\left[\begin{array}{c}0.2\\ 0.4\\ 0.4\end{array}\right]$$
假设得到的观测序列为 $O=(r,w,r)$ ，则前向概率为：
$$\begin{array}{cccc}T& 1& 2& 3\\ \text{box}1& 0.10& 0.077& 0.04187\\ \text{box}2& 0.16& 0.1104& 0.03551\\ \text{box}3& 0.28& 0.0606& 0.05284\end{array}$$
这个表格的计算顺序是一列一列计算的，最后一列的和即为所求结果：
$$P(O|\lambda )=0.13022$$

后向算法

后向概率：给定模型 $\lambda$ ，后向概率定义为，时刻 $t$ 为状态 $q_i$ 的条件下，后续的观测序列正好是 $O[t+1,T]=(o_{t+1},o_{t+2},\cdots ,o_T)$ 的概率：
$${\beta }_t(i)=P(o_{t+1},o_{t+1},\cdots ,o_T|i_t=q_i,\lambda )$$
可以想到，如果知道 $t+1$ 时刻所有状态的后向概率，那么就可以求得 $t$ 时刻任一状态的后向概率，只需要将后向概率乘以状态转移概率，再乘上产生观测为 $o_{t+1}$ 的概率，最后求和即可：

注意：后向概率 ${\beta }_t(i)$ 中并没有包含 $t$ 时刻的观测 $o_t$ 出现的概率。

算法：观测序列概率的后向算法

• 输入：HMM $\lambda$ ，观测序列 $O$ ；
• 输出：产生该观测序列的概率 $P(O|\lambda )$ ；

1. 初值：时刻 $T$ 以后没有观测了，所以可以认为处于每个状态的后向概率是 1：

$${\beta }_T(i)=1,i=1,2,\cdots ,N$$

2. 递归：对于 $t=T-1,T-2,\cdots ,1$ 递归地计算前一时刻的后向概率：

$${\beta }_t(i)=\sum _{j=1}^N{a}_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j),i=1,2,\cdots ,N$$

3. 求和：最后时刻 $t=1$ 时，后向概率要乘以初始分布位于该状态的概率，以及产生第一个观测的概率，并求和：

$$P(O|\lambda )=\sum _{i=1}^N{\pi }_i{b}_i{\beta }_1(i)$$

前向概率和后向概率

状态转移的过程可以看成是如下图的某条路径：

经过某个点的概率：利用前向和后向概率的定义，可以将观测序列概率 $P(O|\lambda )$ 统一写成：
$$P(O|\lambda )=\sum _{j=1}^N{\alpha }_t(j){\beta }_t(j),t=1,2,\cdots ,T-1$$
可以理解成，我们任意选择其中一层，其中共有 $N$ 个节点。对于每个节点，$a_t(i)$ 代表产生前面的观测序列后到达该节点的概率，${\beta }_t(j)$ 代表从该节点出发得到剩余的观测序列的概率。所有这些乘在一起则代表：给定模型 $\lambda$ ，在 $t$ 时刻正好经过该节点并且最终也正好得到 $O$ 的概率

即：
$$P(i_t=q_i,O|\lambda )={\alpha }_t(j){\beta }_t(j)$$
而给定模型 $\lambda$ 和观测序列 $O$ ，在时刻 $t$ 处于状态 $q_i$ 的概率，记为：
$$\begin{array}{rl}{\gamma }_t(i)=& P(i_t=q_i|O,\lambda )\\ =& \frac{P(i_t=q_i,O|\lambda )}{P(O|\lambda )}\\ =& \frac{{\alpha }_t(j){\beta }_t(j)}{\sum _{j=1}^N{\alpha }_t(j){\beta }_t(j)}\end{array}$$
经过某条边的概率：利用前向和后向概率的定义，可以将观测序列概率 $P(O|\lambda )$ 统一写成：
$$P(O|\lambda )=\sum _{i=1}^N\sum _{j=1}^N{a}_t(i)a_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j),t=1,2,\cdots ,T-1$$
可以理解成，我们任意选择相邻的某两层节点，其中共有 $N^2$ 条边。对于每条边，$a_t(i)$ 代表产生前面的观测序列后到达左节点的概率，$a_{ij}$ 代表从左节点转移到右节点的概率，$b_j(o_{t+1}){\beta }_{t+1}(j)$ 代表从右节点出发得到剩余的观测序列的概率。所有这些乘在一起则代表：给定模型 $\lambda$ ，在 $t\to t+1$ 时刻正好经过这条边并且最终也正好得到 $O$ 的概率，即：
$$P(i_t=q_i,i_{t+1}=q_j,O|\lambda )=a_t(i)a_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j)$$
对相邻的某两层节点之间所有边，将经过每条边的概率求和，即得到 $P(O|\lambda )$ 。

而给定模型 $\lambda$ 和观测序列 $O$ ， 在 $t\to t+1$ 时刻正好经过某条边的概率为：
$$\begin{array}{rl}{\xi }_t(i,j)=& P(i_t=q_i,i_{t+1}=q_j|O,\lambda )\\ =& \frac{P(i_t=q_i,i_{t+1}=q_j,O|\lambda )}{P(O|\lambda )}\\ =& \frac{a_t(i)a_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j)}{\sum _{i=1}^N\sum _{j=1}^N{a}_t(i)a_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j)}\end{array}$$

学习问题

学习问题是指，给定观测序列，估计模型参数，使得$P(O|\lambda )$ 最大，即用极大似然估计的方法估计参数。

监督学习算法

就是使用完全数据，假设训练数据包含 $S$ 个长度相同的观测序列和对应的状态序列，那么可以使用极大似然估计法来估计 HMM 的参数：
train_data = {   ( O 1 ,   I 1 ) ,   ( O 1 ,   I 2 ) ,   ⋯   ,   ( O S ,   I S )   } \text{train\_data}=\set{(O_1,\,I_1),\,(O_1,\,I_2),\,\cdots,\,(O_S,\,I_S)} train_data={(O1​,I1​),(O1​,I2​),⋯,(OS​,IS​)}

1. 转移概率 $a_{ij}$ 的估计：就是从状态 $q_i$ 出发的转移中，转移到 $q_j$ 的频率；设样本中从状态 $q_i$ 转移到状态 $q_j$ 的频数为$A_{ij}$ ，则状态转移概率 $a_{ij}$ 的估计是：

$${\hat{a}}_{ij}=\frac{A_{ij}}{\sum _{j=1}^N{A}_{ij}},i=1,2,\cdots ,N;j=1,2,\cdots ,N$$

2. 观测概率 $b_j(k)$ 的估计：就是处于状态 $q_j$ 时，产生观测 $v_k$ 的频率；设样本中处于状态 $q_j$ 时，产生观测 $v_k$ 的频数为 $B_{jk}$ ，则观测概率 $b_j(k)$ 的估计为：

$${\hat{b}}_j(k)=\frac{B_{jk}}{\sum _{k=1}^M{B}_{jk}},j=1,2,\cdots ,N;k=1,2,\cdots ,M$$

3. 初始状态概率 ${\pi }_i$ 的估计：就是 $S$ 个样本中初始状态为 $q_i$ 的频率；

Baum-Welch 算法

Baum-Welch 算法是 EM 算法在 HMM 的参数估计中的具体实现。

E 步

这里 $O=(o_1,o_2,\cdots ,o_T)$ 相当于 EM 算法中的 $Y$ ，$I=(i_1,i_2,\cdots ,i_T)$ 相当于 EM 算法中的 $Z$ ，而 $\lambda =(\pi ,A,B)$ 相当于 EM 算法中的 $\theta$ 。因此，完全数据的对数似然函数为：
$$\log P(O,I|\lambda )$$
$Q$ 函数为：
$$\begin{array}{rl}Q(\lambda |\hat{\lambda })=& E_I[\log P(O,I|\lambda )|O,\hat{\lambda }]\\ =& \sum _I\log P(O,I|\lambda )P(I|\hat{\lambda })\\ =& \sum _I\log P(O,I|\lambda )\frac{P(O,I|\hat{\lambda })}{P(O|\hat{\lambda })}\end{array}$$
因为 $Q$ 函数的变量为 $\lambda$ ， $\frac{1}{P(O|\hat{\lambda })}$ 对于 $Q$ 函数来说是个常数，不影响极值的求解，因此我们将其省去，得到：
$$Q(\lambda |\hat{\lambda })=\sum _I\log P(O,I|\lambda )P(O,I|\hat{\lambda })$$
而：
$$P(O,I|\lambda )={\pi }_{i_1}{b}_{i_1}(o_1)a_{i_1{i}_2}\cdots a_{i_{T-1}{i}_T}{b}_{i_T}(o_T)$$
因此 $Q$ 函数可以写成：
$$\begin{array}{rl}Q(\lambda |\hat{\lambda })=& \sum _I\log {\pi }_{i_1}P(O,I|\hat{\lambda })\\ +& \sum _I\left(\sum _{t=1}^{T-1}\log a_{i_t{i}_{t+1}}\right)P(O,I|\hat{\lambda })\\ +& \sum _I\left(\sum _{t=1}^T\log b_{i_t}(o_t)\right)P(O,I|\hat{\lambda })\end{array}$$

M 步

由于 ${\pi }_i$ 、$a_{ij}$ 和 $b_j(k)$ 分别出现在 $Q$ 函数的三项中，因此只需要对各项分别极大化：

1. 求 ${\pi }_i$ 极值点：可以进行一些变换，使用了类似于导出三硬币的 EM 算法时的技巧：

$$\begin{array}{rl}\sum _I\log {\pi }_{i_1}P(O,I|\hat{\lambda })=& \sum _{i_1}\sum _{i_t,t\ne 1}\log {\pi }_{i_1}P(O,I|\hat{\lambda })\end{array}$$

注意到在固定 $i_1$ 的取值时，有：
$$\sum _{i_t,t\ne 1}\log {\pi }_{i_1}P(O,I|\hat{\lambda })=\log {\pi }_{i_1}\sum _{i_t,t\ne 1}P(O,I|\hat{\lambda })=\log {\pi }_{i_1}P(O,i_1=i_1|\hat{\lambda })$$
$P(O,i_1=i_1|\hat{\lambda })$ 这个记号写得不好。。。其实就是代表给定参数 $\hat{\lambda }$ 的情况下，得到观测 $O$ 并且第一个状态正好是刚刚固定的 $i_1$ 的概率。因此：
$$\sum _I\log {\pi }_{i_1}P(O,I|\hat{\lambda })=\sum _{i=1}^N\log {\pi }_{i}P(O,i_1=q_i|\hat{\lambda })$$
这里 ${\pi }_i$ 满足约束条件 $\sum _{i=1}^N{\pi }_i=1$ ，使用 Lagrange 乘子法，Lagrangian 为：
$$L(\pi )=\sum _{i=1}^N\log {\pi }_{i}P(O,i_1=q_i|\hat{\lambda })-\gamma (1-\sum _{i=1}^N{\pi }_i)$$
FOC 为：
$$\begin{array}{rl}\frac{\partial L(\pi )}{\partial {\pi }_i}=& \frac{P(O,i_1=q_i|\hat{\lambda })}{{\pi }_i}+\gamma =0\\ \Rightarrow & P(O,i_1=q_i|\hat{\lambda })+{\pi }_i\gamma =0\end{array}$$
对所有 FOC 求和，解得 $\gamma$ 为：
$$\gamma =-\sum _{i=1}^{N}P(O,i_1=q_i|\hat{\lambda })=-P(O|\hat{\lambda })$$
代回解得：
$${\pi }_i=\frac{P(O,i_1=q_i|\hat{\lambda })}{P(O|\hat{\lambda })}$$
这实际上就是给定参数 $\hat{\lambda }$ 和观测 $O$ ，1 时刻经过状态 $q_i$ 的概率。

2. 求 $a_{ij}$ 极值点：同样地，对目标函数进行一些变换：

$$\begin{array}{rl}& \sum _I\left(\sum _{t=1}^{T-1}\log a_{i_t{i}_{t+1}}\right)P(O,I|\hat{\lambda })\\ =& \sum _I\sum _{t=1}^{T-1}\log a_{i_t{i}_{t+1}}P(O,I|\hat{\lambda })\text{(因为P(O,I∣λ^)与t无关)}\\ =& \sum _{t=1}^{T-1}\sum _I\log a_{i_t{i}_{t+1}}P(O,I|\hat{\lambda })\text{(可以交换求和顺序)}\\ =& \sum _{t=1}^{T-1}\sum _{i_t}\sum _{i_{t+1}}\sum _{i_{\text{else}}}\log a_{i_t{i}_{t+1}}P(O,I|\hat{\lambda })\\ =& \sum _{t=1}^{T-1}\sum _{i_t}\sum _{i_{t+1}}\log a_{i_t{i}_{t+1}}\sum _{i_{\text{else}}}P(O,I|\hat{\lambda })\text{(因为logaitit+1与ielse无关)}\\ =& \sum _{t=1}^{T-1}\sum _{i=1}^N\sum _{j=1}^N\log a_{ij}P(O,i_t=q_i,i_{t+1}=q_j|\hat{\lambda })\\ =& \sum _{i=1}^N\sum _{j=1}^N\sum _{t=1}^{T-1}\log a_{ij}P(O,i_t=q_i,i_{t+1}=q_j|\hat{\lambda })\text{(可以交换求和顺序)}\end{array}$$

约束条件同样是 $\sum _{j=1}^N{a}_{ij}=1$ ，Lagrangian 为：
$$L(a)=\sum _{i=1}^N\sum _{j=1}^N\sum _{t=1}^{T-1}\log a_{ij}P(O,i_t=q_i,i_{t+1}=q_j|\hat{\lambda })-\sum _{i=1}^N{\gamma }_i(1-\sum _{j=1}^N{a}_{ij})$$
FOC 为：
$$\begin{array}{rl}\frac{\partial L(a)}{\partial a_{ij}}=& \frac{\sum _{t=1}^{T-1}P(O,i_t=q_i,i_{t+1}=q_j|\hat{\lambda })}{a_{ij}}+{\gamma }_i=0\\ \Rightarrow & \sum _{t=1}^{T-1}P(O,i_t=q_i,i_{t+1}=q_j|\hat{\lambda })+a_{ij}{\gamma }_i=0\end{array}$$
将所有 $i$ 相同的 FOC 相加，得到：
$${\gamma }_i=-\sum _{t=1}^{T-1}P(O,i_t=q_i|\hat{\lambda })$$
代回解得：
$$a_{ij}=\frac{\sum _{t=1}^{T-1}P(O,i_t=q_i,i_{t+1}=q_j|\hat{\lambda })}{\sum _{t=1}^{T-1}P(O,i_t=q_i|\hat{\lambda })}$$

3. 求 $b_j(k)$ 极值点：同样地，对目标函数进行一些变换：

$$\begin{array}{rl}& \sum _I\left(\sum _{t=1}^T\log b_{i_t}(o_t)\right)P(O,I|\hat{\lambda })\\ =& \sum _I\sum _{t=1}^T\log b_{i_t}(o_t)P(O,I|\hat{\lambda })\text{(因为P(O,I∣λ^)与t无关)}\\ =& \sum _{t=1}^T\sum _I\log b_{i_t}(o_t)P(O,I|\hat{\lambda })\text{(可以交换求和顺序)}\\ =& \sum _{t=1}^T\sum _{i_t}\sum _{i_{\text{else}}}\log b_{i_t}(o_t)P(O,I|\hat{\lambda })\\ =& \sum _{t=1}^T\sum _{i_t}\log b_{i_t}(o_t)\sum _{i_{\text{else}}}P(O,I|\hat{\lambda })\text{(因为logbit(ot)与ielse无关)}\\ =& \sum _{t=1}^T\sum _{j=1}^N\log b_j(o_t)P(O,i_t=q_j|\hat{\lambda })\\ =& \sum _{j=1}^N\sum _{t=1}^T\log b_j(o_t)P(O,i_t=q_j|\hat{\lambda })\end{array}$$

约束条件为 $\sum _{k=1}^M{b}_i(k)=1$ ，Lagrangian 为：
$$L(b)=\sum _{j=1}^N\sum _{t=1}^T\log b_j(o_t)P(O,i_t=q_j|\hat{\lambda })-\sum _{j=1}^N{\gamma }_j(1-\sum _{k=1}^M{b}_i(k))$$
FOC 为（注意，只有 $o_t=v_k$ 时 $b_j(o_t)$ 对 $b_j(k)$ 的偏导数才不为 0，这里用指示函数来表示）：
$$\begin{array}{rl}\frac{\partial L(b)}{\partial b_j(k)}=& \frac{\sum _{t=1}^{T}P(O,i_t=q_j|\hat{\lambda })I(o_t=v_k)}{b_j(k)}+{\gamma }_j=0\\ \Rightarrow & \sum _{t=1}^{T}P(O,i_t=q_j|\hat{\lambda })I(o_t=v_k)+b_j(k){\gamma }_j=0\end{array}$$
将所有 $j$ 相同的 FOC 相加，得到（因为有且仅有一个 $k$ 使得 $I(o_t=v_k)=1$）：
$${\gamma }_j=-\sum _{t=1}^{T}P(O,i_t=q_j|\hat{\lambda })$$
代回解得：
$$b_j(k)=\frac{\sum _{t=1}^{T}P(O,i_t=q_j|\hat{\lambda })I(o_t=v_k)}{\sum _{t=1}^{T}P(O,i_t=q_j|\hat{\lambda })}$$

参数估计公式

以上参数估计的公式可以使用前向概率和后向概率表示，即：
$$a_{ij}=\frac{\sum _{t=1}^{T-1}{\xi }_t(i,j)}{\sum _{t=1}^{T-1}{\gamma }_t(i)},b_j(k)=\frac{\sum _{t=1,o_t=v_k}^T{\gamma }_t(i)}{\sum _{t=1}^T{\gamma }_t(i)},{\pi }_i={\gamma }_1(i)$$

算法描述

输入：观测数据 $O=(o_1,o_2,\cdots ,o_T)$ ；

输出：HMM 参数

1. 初始化，对 $n=0$ ，选取模型参数初始值 $a_{ij}^{(0)}$ ，$b_j(k{)}^{(0)}$ ，${\pi }_i^{(0)}$ ；
2. 递推，对于 $n=1,2,\cdots$ ，使用前一次递推得到的参数计算前向概率和后向概率，并更新参数：

$$a_{ij}=\frac{\sum _{t=1}^{T-1}{\xi }_t(i,j)}{\sum _{t=1}^{T-1}{\gamma }_t(i)},b_j(k)=\frac{\sum _{t=1,o_t=v_k}^T{\gamma }_t(i)}{\sum _{t=1}^T{\gamma }_t(i)},{\pi }_i={\gamma }_1(i)$$

3. 中止，得到模型参数 ${\lambda }^{(n+1)}=(A^{(n+1)},B^{(n+1)},{\pi }^{(n+1)})$ ；

解码问题

解码问题是指，已知模型参数，给定观测序列，求出最有可能的对应的状态序列。

近似算法

回忆一下前面 ${\gamma }_t(i)$ 的定义是：给定模型 $\lambda$ 和观测序列 $O$ ，在时刻 $t$ 处于状态 $q_i$ 的概率。

近似算法是直接使用每一个时刻出现概率最大的状态作为输出，即：
$$i_t^{\ast }=\arg \underset{1\le i\le N}{\max }[{\gamma }_t(i)],t=1,2,\cdots ,T$$
从而得到状态序列 $I^{\ast }=(i_1^{\ast },i_2^{\ast },\cdots ,i_T^{\ast })$ 。

近似算法顾名思义，就是近似的，没有更多地考虑到状态转移概率 $a_{ij}$ ，但是它计算简单。

Viterbi 算法

维特比算法实际上是动态规划，属于精确算法。原理为：最优路径的子路径一定是最优的（这里的最优指的是概率最大）。例如，若最优路径在 $t$ 时刻经过节点 $i_t^{\ast }$ ，则从 $i_t^{\ast }$ 到 $i_T^{\ast }$ 的子路径一定是所有从 $i_t^{\ast }$ 到 $i_T^{\ast }$ 的路径中最优的。否则我们可以选择更优的从 $i_t^{\ast }$ 到 $i_T^{\ast }$ 的路径去替换原来的子路径，得到更优的总的路径。因此，我们从 $t=1$ 开始，递推地计算时刻 $t$ 状态为 $i$ 的各部分路径的最大概率，直至时刻 $t=T$ 。

我们定义变量 $\delta$ ，代表时刻 $t$ 状态为 $i$ 的所有单个路径 $(i_1,i_2,\cdots ,i_t)$ 中概率的最大值：
$${\delta }_t(i)=\underset{(i_1,i_2,\cdots ,i_t)}{\max }P(i_t=i,o_t,\cdots ,o_1|\lambda ),i=1,2,\cdots ,N$$
假设一个共同的起始节点为 $q_0$ ，其到每个状态的转移概率为每个状态对应的初始概率，即 $a_{0i}={\pi }_i$ 。

假设我们前面已经计算了时刻 $[0,t]$ ，得到了 $0$ 时刻从 $i_0=q_0$ 节点出发，$t$ 时刻到达各个节点的最优路径和对应的概率 ${\delta }_t(j)$ 。则计算 $t+1$ 时刻时，对于某个状态 $q_i$ ，我们只需求得 $j=\arg \underset{j}{\max }{\delta }_t(j)a_{ji}{b}_j(o_{t+1})$ ，则说明 ${\delta }_{t+1}(i)={\delta }_t(j)a_{ji}{b}_j(o_{t+1})$ ，并且 $0$ 时刻从 $i_0$ 节点出发，$t+1$ 时刻到达状态 $q_i$ 的最优路径中，前一个节点的状态就是 $q_j$ 。

相应的递推公式为：
$${\delta }_{t+1}(i)=\underset{1\le j\le N}{\max }{\delta }_t(j)a_{ji}{b}_j(o_{t+1}),i=1,2,\cdots ,N;t=1,2,\cdots ,T-1$$
算法中需要记录这个 $q_j$ ，我们定义变量 ${\Psi }_t(i)$ ，代表时刻 $t$ 状态为 $q_i$ 的所有单个路径 $(i_1,i_2,\cdots ,i_{t-1},i_t)$ 中概率最大的路径的第 $t-1$ 个节点为：
$${\Psi }_t(i)=\arg \underset{1\le j\le N}{\max }{\delta }_t(j)a_{ji}{b}_j(o_{t+1}),i=1,2,\cdots ,N$$
（我觉得算法中计算的东西可以叫做最优前缀子路径）

最后，当计算完所有最优前缀子路径后，我们取 $i=\arg \underset{i}{\max }{\delta }_t(i)$ ，则说明我们所要求的总的最优路径的最后一个节点的状态是 $q_i$ ，此时只需按照所记录的前一个节点的状态递推回去，就能得到整个最优路径。

算法：维比特算法

• 输入：模型 $\lambda =(A,B,\pi )$ 和观测 $O=(o_1,o_2,\cdots ,o_T)$ ；
• 输出：最优路径 $I^{\ast }=(i_1^{\ast },i_2^{\ast },\cdots ,i_T^{\ast })$ ；

1. 初始化：

$$\begin{array}{c}{\delta }_1(i)={\pi }_i{b}_i(o_1),i=1,2,\cdots ,N\\ {\Psi }_1(i)=0,i=1,2,\cdots ,N\end{array}$$

2. 递推，对于 $t=2,3,\cdots ,T$ ：

$$\begin{array}{c}{\delta }_t(i)=\underset{1\le j\le N}{\max }{\delta }_{t-1}(j)a_{ji}{b}_i(o_t),i=1,2,\cdots ,N\\ {\Psi }_t(i)=\arg {\max }_{1\le j\le N}{\delta }_{t-1}(j)a_{ji}{b}_i(o_t),i=1,2,\cdots ,N\end{array}$$

注意，这里 ${\Psi }_t(i)$ 的表达式也可以写成以下形式，因为 $b_i(o_t)$ 对于第 $t-1$ 层而言是个常数，乘不乘进去都可以，反正对大家来说都是一样的效果：
$${\Psi }_t(i)=\arg \underset{1\le j\le N}{\max }{\delta }_{t-1}(j)a_{ji},i=1,2,\cdots ,N$$

3. 中止：

$$\begin{array}{c}{P}^{\ast }=\underset{1\le i\le N}{\max }{\delta }_T(i)\\ i_T^{\ast }=\arg \underset{1\le i\le N}{\max }{\delta }_T(i)\end{array}$$

4. 最优路径回溯：对于 $t=T-1,T-2,\cdots ,1$ ：

$$i_t^{\ast }={\Psi }_{t+1}(i_{t+1}^{\ast })$$

即得到最优路径 $I^{\ast }=(i_1^{\ast },i_2^{\ast },\cdots ,i_T^{\ast })$ ；

例（盒子和球模型）：状态集合为 Q = {   1 ,   2 ,   3   } Q=\set{1,\,2,\,3} Q={1,2,3}，观测集合为 V = {   r ,   w   } V=\set{r,\,w} V={r,w} ，
$$A=\left[\begin{array}{ccc}0.5& 0.2& 0.3\\ 0.3& 0.5& 0.2\\ 0.2& 0.3& 0.5\end{array}\right]B=\left[\begin{array}{cc}0.5& 0.5\\ 0.4& 0.6\\ 0.7& 0.3\end{array}\right]\pi =\left[\begin{array}{c}0.2\\ 0.4\\ 0.4\end{array}\right]$$
假设得到的观测序列为 $O=(r,w,r)$ ，则计算过程为：

最优路径为 $I^{\ast }=(3,3,3)$ ；