机器学习之隐马尔科夫模型（HMM）

1、隐马尔科夫模型介绍
2、隐马尔科夫数学原理
3、Python代码实现隐马尔科夫模型
4、总结

隐马尔可夫模型介绍

马尔科夫模型（hidden Markov model，HMM）是关于时序的概率模型，描述由一个隐藏的马尔科夫随机生成不可观测的状态随机序列，再由各个状态生成一个观测从而产生观测随机序列的过程，属于一个生成模型。

下面我们来从概率学角度定义马尔科夫模型，从一个典型例子开始：

假设有4个盒子，每个盒子里面有不同数量的红、白两种颜色的球，具体如下表：

盒子编号	1	2	3	4
红球数	5	3	6	8
白球数	5	7	4	2

现在从这些盒子中取出T个球，取样规则为每次选择一个盒子取出一个球，记录其颜色，放回。在这个过程中，我们只能观测到球的颜色的序列，观测不到球是从哪个盒子中取出来的，即观测不到盒子的序列，这里有两个随机序列，一个是盒子的序列（状态序列），一个是球的颜色的观测序列（观测序列），前者是隐藏的，只有后者是可观测的。这里就构成了一个马尔科夫的例子。
定义 $Q$ 是所有的可能的状态集合，V是所有的可能的观测的集合：
$Q = \{q_1,q_2,\cdots,q_N\},　　　　　V = \{v_1,v_2,\cdots,v_M\}$
其中，Ｎ是可能的状态数，Ｍ是可能的观测数，例如上例中Ｎ＝４，Ｍ＝２。

$I$ 是长度为T的状态序列， $O$ 是对应的观测序列：
$I = (i_1,i_2,\cdots,i_T),　　　　　　O = (o_1,o_2,\cdots,o_T)$
A是状态转移概率矩阵：
$A = [a_{i,j}]_{N \times N}$
其中， $a_{i,j} = P(i_{t+1} = q_j|i_t = q_i),　i=1,2,\cdots,N;j=1,2,\cdots,N$ 是指在时刻 $t$ 处于状态 $q_i$ 的条件下在时刻 $t+1$ 转移到状态 $q_j$ 的概率。

B是观测概率矩阵：
$B = [b_j(k)]_{N \times M}$
其中， $b_j(k) = P(o_t = v_k|i_t = q_j),　k=1,2,\cdots,M;j=1,2,\cdots,N$ 是指在时刻 $t$ 处于状态 $q_j$ 的条件下生成观测 $v_k$ 的概率。

$\pi$ 是初始状态概率向量：
$\pi = (\pi_i)$
其中， $\pi_i = P(i_1 = q_i),　i=1,2,\cdots,N$ 是指在时刻 $t$ =1处于状态 $q_i$ 的概率。

由此可得到，隐马尔可夫模型 $\lambda$ 的三元符号表示，即
$\lambda = (A,B,\pi)$
$A,B,\pi$ 称为隐马尔可夫模型的三要素。

由定义可知隐马尔可夫模型做了两个基本假设：

(1)齐次马尔科夫性假设，即假设隐藏的马尔科夫链在任意时刻 $t$ 的状态只和 $t$ -1状态有关；
$P(i_t|i_{t-1},o_{t-1},\cdots,i_1,o_1) = P(i_t|i_{t-1}),　　　　t=1,2,\cdots,T$
(2)观测独立性假设，观测只和当前时刻状态有关；
$P(o_t|i_T,i_{T-1},o_{T-1},\cdots,i_{t+1},o_{t+1},i_t,i_{t-1},o_{t-1},\cdots,i_1,o_1) = P(O_t|i_t)$
仍以上面的盒子取球为例，假设我们定义盒子和球模型：

状态集合： $Q$ = {盒子1，盒子2，盒子3，盒子4}， N=4
观测集合： $V$ = {红球，白球} M=2
初始化概率分布：
$\pi = (0.25,0.25,0.25,0.25)^T$
状态转移矩阵：
$A = \left[ \begin{matrix} 0 & 1 & 0 & 0 \\ 0.4 & 0 & 0.6 & 0 \\ 0 & 0.4 & 0 & 0.6 \\ 0 & 0 & 0.5 & 0.5 \end{matrix} \right]$
观测矩阵:
$B = \left[ \begin{matrix} 0.5 & 0.5 \\ 0.3 & 0.7 \\ 0.6 & 0.4 \\ 0.8 & 0.2 \end{matrix} \right]$

隐马尔可夫模型的三个基本问题

1、概率计算问题

给定： $\lambda = (A,B,\pi)　O=(o_1,o_2,\cdots,o_T)$

计算： $P(O|\lambda)$
2、学习问题

已知： $O=(o_1,o_2,\cdots,o_T)$

估计： $\lambda = (A,B,\pi)$ ,使 $P(O|\lambda)$ 最大
3、预测问题（解码）

已知： $\lambda = (A,B,\pi)　O=(o_1,o_2,\cdots,o_T)$

求：使 $P(I|O)$ 最大的状态序列 $I=(i_1,i_2,\cdots,i_T)$

下面我们使用python代码写一个HMM模型生成序列 $O$ 的示例代码

import numpy as np

class HMM(object):
    def __init__(self, N, M, pi=None, A=None, B=None):
        self.N = N
        self.M = M
        self.pi = pi
        self.A = A
        self.B = B

    def get_data_with_distribute(self, dist): # 根据给定的概率分布随机返回数据（索引）
        r = np.random.rand()
        for i, p in enumerate(dist):
            if r < p: return i
            r -= p

    def generate(self, T: int):
        '''
        根据给定的参数生成观测序列
        T: 指定要生成观测序列的长度
        '''
        result = []
        for ind in range(T):        # 依次生成状态和观测数据
            if ind==0:
                i = self.get_data_with_distribute(self.pi)
            else:
                i = self.get_data_with_distribute(self.A[i])
            o = self.get_data_with_distribute(self.B[i])
            result.append(o)
        return result

if __name__ == "__main__":
    pi = np.array([0.25, 0.25, 0.25, 0.25])
    A = np.array([
        [0,  1,  0, 0],
        [0.4, 0, 0.6, 0],
        [0, 0.4, 0, 0.6],
        [0, 0, 0.5, 0.5]])
    B = np.array([
        [0.5, 0.5],
        [0.3, 0.7],
        [0.6, 0.4],
        [0.8, 0.2]])
    hmm = HMM(4, 2, pi, A, B)
    print(hmm.generate(10))  # 生成10个数据
 
# 生成结果如下
[0, 0, 0, 0, 1, 1, 1, 1, 0, 0]   # 0代表红球，1代表白球

隐马尔可夫模型数学原理

下面我们从隐马尔可夫模型的三个基本问题出发，逐个解释问题解法：

1、概率计算问题
- 直接计算法
- 前向算法
- 后向算法
直接计算法：
- 状态序列 $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}$
- 对固定的状态序列 $I$ ，观测序列 $O$ 的概率： $P(O|I,\lambda)$
  $P(O|I,\lambda) = b_{i_1}(o_1) b_{i_2}(o_2) \cdots b_{i_T}(o_T)$
- $O$ 和 $I$ 同时出现的联合概率为：
  $P(O,I|\lambda) = P(O|I,\lambda)P(I|\lambda) = \pi_{i_1}b_{i_1}(o_1)a_{i_1 i_2} b_{i_2}(o_2) \cdots a_{i_{T-1}i_T}b_{i_T}(o_T)$
- 对所有的可能的状态序列 $I$ 求和，得到观测序列 $O$ 的概率：
  $P(O|\pi) = \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) \cdots a_{i_{T-1}i_T}b_{i_T}(o_T)$
  此算法的复杂度为 $O(TN^T)$ ，计算量太大，不可行。
前向算法：
- 前向概率定义：给定隐马尔科夫模型 $\lambda$ ，定义到时刻t部分观测序列为： $o_1,o_2,\cdots,o_t$ ，且状态为 $q_i$ 的概率为前向概率，记作：
  $\alpha_t(i) = P(o_1,o_2,\cdots,o_t,i_t=q_t|\lambda)$
  输入：隐马尔可夫模型 $\lambda$ ，观测序列 $O$ ；
  
  输出：观测序列概率 $P(O|\lambda)$
  
  初值： $\alpha_1(i) = \pi_ib_i(o_1)，i=1,2,\cdots,N$
  
  递推： $\alpha_{t+1}(i) = \left[ \sum_{j=1}^N \alpha_t(j) a_{_ji} \right]b_i(o_{t+1}),i=1,2,\cdots,N$
  
  终止： $P(O|\lambda) = \sum_{i=1}^N \alpha_T(i)$
  
  算法复杂度 $O(N^2T)$ ，算法复杂度大大降低
  
  因为： $\alpha_T(i) = P(o_1,o_2,\cdots,o_T,i_T=q_i|\lambda)$
  
  所以： $P(O|\lambda) = \sum_{i=1}^N \alpha_T(i)$
- 下面以一个具体的例子来演示前向算法计算过程：
  $A = \left[ \begin{matrix} 0.5 & 0.2 & 0.3 \\ 0.3 & 0.5 & 0.2 \\ 0.2 & 0.3 & 0.5 \end{matrix} \right], B = \left[ \begin{matrix} 0.5 & 0.5 \\ 0.4 & 0.6 \\ 0.7 & 0.3 \end{matrix} \right], \pi = (0.2,0.4,0.4)^T$
  假设T = 3，O = (红，白，红)，使用前向算法计算 $P(O|\lambda)$ ，模型 $\lambda = (A,B,\pi)$ ，状态集合 $Q=\{1,2,3\}$ ，观测集合 $V = \{红，白\}$
  
  （1）计算初值
  $\alpha_1(1) = \pi_1 b_1(o_1) = 0.2\times0.5 = 0.1 \\ \alpha_1(2) = \pi_2 b_2(o_1) = 0.4\times0.4 = 0.16 \\ \alpha_1(3) = \pi_3 b_3(o_1) = 0.4\times0.7 = 0.28$
  （2）递推计算
  $\begin{align} \alpha_2(1) & = \left[ \begin{matrix} \sum_{i=1}^3\alpha_1(i)a_{i1} \end{matrix} \right]b_1(o_2) = (0.1\times0.5+0.16\times0.3+0.28\times0.2)\times0.5 = 0.154\times0.5 = 0.077 \\ \alpha_1(2) & = \left[ \begin{matrix} \sum_{i=1}^3\alpha_1(i)a_{i2} \end{matrix} \right]b_2(o_2) = 0.184\times0.6 = 0.1104 \\ \alpha_1(3) & = \left[ \begin{matrix} \sum_{i=1}^3\alpha_1(i)a_{i3} \end{matrix} \right]b_3(o_2) = 0.202\times0.3 = 0.0606 \\ \alpha_3(1) & = \left[ \begin{matrix} \sum_{i=1}^3\alpha_2(i)a_{i1} \end{matrix} \right]b_1(o_3) = 0.04187 \\ \alpha_3(2) & = \left[ \begin{matrix} \sum_{i=1}^3\alpha_2(i)a_{i2} \end{matrix} \right]b_2(o_3) = 0.03551 \\ \alpha_3(3) & = \left[ \begin{matrix} \sum_{i=1}^3\alpha_2(i)a_{i3} \end{matrix} \right]b_3(o_3) = 0.05284 \end{align}$
  （3）终止
  $P(O|\lambda) = \sum_{i=1}^3 \alpha_3(i) = 0.13022$
后向算法：
- 后向概率定义：给定隐马尔可夫模型 $\lambda$ ，定义在时刻t状态为 $q_i$ 的条件下，从t+1到T的部分观测序列为： $o_{i+1},o_{i+2},\cdots,o_T$ 的概率为后向概率，记作：
  $\beta_t(i) = P(o_{t+1},o_{t+2},\cdots,o_T|i_t=q_i,\lambda)$
  仍然可以使用递推的方法求得后向概率 $\beta_t(i)$ 及观测序列概率 $P(O|\lambda)$
  
  输入：隐马尔科夫模型 $\lambda$ ，观测序列O；
  
  输出：观测序列概率 $P(O|\lambda)$ .
  
  （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 $P(O|\lambda)=\sum_{i=1}^N \pi_ib_i(o_1)\beta_1(i)$
  
  前后向统一写为：（t=1和t=T-1分别对应）
  $P(O|\lambda) = \sum_{i=1}^N \sum_{j=1}^N \alpha_t(i)a_{ij}b_j(o_{t+1}) \beta_{t+1}(j), 　t=1,2,\cdots,T-1$
  一些概率和期望值计算方式：
  
  1、给定模型 $\lambda$ 和观测 $O$ ，在时刻 $t$ 处于状态 $q_i$ 的概率，记 $\gamma_t(i)=P(i_t=q_t|O,\lambda)$
  $\gamma_t(i)=P(i_t=q_t|O,\lambda)=\frac{P(i_t=q_t,O|\lambda)}{P(O|\lambda)}=\frac{\alpha_t(i)\beta_t(i)}{\sum_{j=1}^N \alpha_t(j)\beta_t(j)}$
  2、给定模型 $\lambda$ 和观测 $O$ ，在时刻 $t$ 处于状态 $q_i$ 且在时刻 $t+1$ 处于 $q_j$ 概率，记 $\xi_t(i,j)=P(i_t=q_i,i_{t+1}=q_j|O,\lambda)$
  $\xi_t(i,j)=P(i_t=q_i,i_{t+1}=q_j,O|\lambda)=\frac{P(i_t=q_t,i_{t+1}=q_j,O|\lambda)}{P(O|\lambda)}=\frac{P(i_t=q_t,i_{t+1}=q_j,O|\lambda)}{\sum_{i=1}^N \sum_{j=1}^N P(i_t=q_i,i_{t+1}=q_j,O|\lambda)} = \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)}$
  3、由1和2可得，(1)在观测O下状态i出现的期望值： $\sum_{t=1}^T \gamma_t(i)$ ，（2）在观测O下由状态i转移的期望值： $\sum_{t=1}^{T-1} \gamma_t(i)$ ，（3）在观测O下由状态i转移到状态j的期望值： $\sum_{t=1}^{T-1} \xi_t(i,j)$ 。
2、学习问题

监督学习方法：
- 假设训练数据是包括观测序列O和对应的状态序列I， ${(O_1,I_1),(O_2,I_2),\cdots,(O_S,I_S)}$
- 可以利用极大似然估计发来估计隐马尔可夫模型参数。
已知： ${(O_1,I_1)(O_2,I_2),\cdots,(O_S,I_S)}$

（1）转移概率 $a_{ij}$ 的估计：假设样本中时刻t处于状态i，时刻t+1转移到状态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)$ 的估计：设样本中状态为j并观测为k的频数是 $B_j(k)$ 那么状态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$ 的估计 $\hat{\pi}_i$ 为S个样本中初始状态为 $q_i$ 的频率。

非监督学习方法：

Baum-Welch算法

假定训练数据只包括{O1,O2,...Os}，求模型参数 $\lambda=(A,B,\pi)$ ，实质上是有隐变量的概率模型：EM算法 $P(O|\lambda)=\sum_I P(O|I,\lambda)P(I|\lambda)$ 。

（1）确定完全数据的对数似然函数

完全数据 $(O,I)=(o_1,o_2,\cdots,o_T,i_1,i_2,\cdots,i_r)$

完全数据的对数似然函数 $logP(O,I|\lambda)$

（2）EM的E步，求Q函数 $Q(\lambda,\bar{\lambda})=\sum_I logP(O,I|\lambda)P(O,I|\bar{\lambda})$ ， $P(O,I|\lambda)=\pi_{i_1}b_{i_1}(O_1)a_{i_1 i_2}b_{i_2}(o_2)\cdots a_{i_{T-1}i_T}b_{i_T}(O_T)$ 则：
$Q(\lambda,\bar{\lambda})=\sum_I log\pi_{i_1}P(O,I|\hat{\lambda})+\sum_I(\sum_{t=1}^{T-1}log a_{i_t i_{t+1}})P(O,I|\bar{\lambda})+\sum_I(\sum_{t=1}^T logb_{i_t}(o_t))P(O,I|\bar{\lambda})$
对序列总长度T进行

（3）EM算法的M步，极大化 $Q(\lambda,\bar{\lambda})$ 求模型参数 $A,B,\pi$

第一项 $\sum_I log\pi_{i_0}P(O,I|\hat{\lambda})=\sum_{i=1}^N log\pi_i P(O,i_1=i|\bar{\lambda})$ ，由约束条件： $\sum_{i=1}^N \pi_i = 1$ 利用拉格朗日乘子 $\sum_{i=1}^N log\pi_iP(O,i_1=i|\bar{\lambda})+\gamma(\sum_{i=1}^N \pi_i -1)$ 求偏导数，并令结果等于0
$\frac{\partial}{\partial \pi_i}[\sum_{i=1}^N log\pi_i P(O,i_1=i|\bar{\lambda})+\gamma(\sum_{i=1}^N \pi_i -1)]=0$
得： $P(O,i_1 = i|\bar{\lambda})+\gamma\pi_i = 0$ $\gamma=-P(O|\bar{\lambda})$ $\pi_i = \frac{P(O,i_1=i|\bar{\lambda})}{P(O|\bar{\lambda})}$

第二项可写成 $\sum_I(\sum_{t=1}^{T-1}log a_{i_t i_{t+1}})P(O,I|\bar{\lambda})=\sum_{i=1}^N \sum_{j=1}^N \sum_{t=1}^{T-1} log a_{ij}P(O,i_t=i,i_{t+1}=j|\bar{\lambda})$ 由约束条件 $\sum_{j=1}^N a_{ij}=1$ ，拉格朗日乘子法得：
$a_{ij}=\frac{\sum_{t=1}^{T-1}P(O,i_t=i,i_{t+1}=j|\bar{\lambda})}{\sum_{t=1}^{T-1}P(O,i_t=i|\bar{\lambda})}$
第三项 $\sum_I(\sum_{t=1}^T logb_{i_t}(o_t))P(O,I|\bar{\lambda})=\sum_{j=1}^N \sum_{t=1}^T logb_j(o_t)P(O,i_t=j|\bar{\lambda})$ 由约束条件 $\sum_{k=1}^M b_j(k)=1$ ，注意，只有在 $o_t=v_k$ 时 $b_j(o_t)$ 对 $b_j(k)$ 的偏导数才不为0，以 $I(o_t=v_k)$ 表示，求得
$b_j(k)=\frac{\sum_{t=1}^T P(O,i_t=j|\bar{\lambda})I(o_t=v_k)}{\sum_{t=1}^N P(O,i_t=j|\bar{\lambda})}$
将以上得到的概率分别用 $\gamma_t(i)，\xi_t(i,j)$ 表示
$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(j)}{\sum_{t=1}^T \gamma_t(j)} \\ \pi_i = \gamma_1(i)$
Baum-Welch算法流程：

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

输出：隐马尔可夫模型参数。

（1）初始化

对n=0，选取 $a_{ij}^{(0)}，b_j(k)^{(0)}，\pi_i^{(0)}$ 得到模型 $\lambda^{(0)}=(A^{(0)},B^{(0)},\pi^{(0)})$

（2）递推，对n=1,2,...,
$a_{ij}^{(n+1)}=\frac{\sum_{t=1}^{T-1}\xi_t(i,j)}{\sum_{t=1}^{T-1}\gamma_t(i)} \\ b_j(k)^{(n+1)}=\frac{\sum_{t=1,o_t=v_k}^T \gamma_t(j)}{\sum_{t=1}^T \gamma_t(j)} \\ \pi_i^{(n+1)} = \gamma_1(i)$
右侧各值按照观测 $O=(o_1,o_2,\cdots,o_T)$ 和模型 $\lambda^{(n)}=(A^{(n)},B^{(n)},\pi^{(n)})$ 计算

（3）终止，得到模型参数 $\lambda^{(n+1)}=(A^{(n+1)},B^{(n+1)},\pi^{(n+1)})$
3、预测问题

预测算法

（1）近似算法

想法：在每个时刻t选择在该时刻最后可能的状态 $i_t^*$ ，从而得到一个序列状态序列 $I^*=(i_1^*,i_2^*,\cdots,i_T^*)$ ，将它作为预测的结果，在时刻t处于状态 $q_i$ 的概率：
$\gamma_t(i)=\frac{\alpha_t(i) \beta_t(i)}{P(O|\lambda)}=\frac{\alpha_t(i) \beta_t(i)}{\sum_{j=1}^N \alpha_t(j) \beta_t(j)}$
在每一时刻t最有可能的状态是： $i_t^*=\underset{1\leq i \leq N}{\operatorname{argmax}}[\gamma_t(i)],t=1,2,\cdots,T$

从而得到状态序列： $I_*=(i_1^*,i_2^*,\cdots,i_T^*)$

得到的状态有可能实际不发生

（2）维比特算法

用动态规划解概率最大路径，一个路径对应一个状态序列。
最优路径具有这样的特性：如果最优路径在时刻t通过结点 $i_t^*$ ，那么这一路径从结点 $i_t^*$ 到终点 $i_T^*$ 的部分路径，对于从 $i_t^*$ 到 $i_T^*$ 的所有可能的部分路径来说，必须是最优的。
只需从时刻t=1开始，递推地计算在时刻t状态为i的各条部分路径的最大概率，直至得到时刻t=T状态为i的各条路径的最大概率，时刻t=T的最大概率即为最优路径的概率 $P^*$ 最有路径的终结点 $i_T^*$ 也同时得到。
之后，为了找出最优路径的各个结点，从终结点开始，由后向前逐步求得结点 $i_{T-1}^*,\cdots,i_1^*$ ，得到最优路径 $I^*=(i_1^*,i_2^*,\cdots,i_T^*)$

导入两个变量 $\delta$ 和 $\psi$ ，定义在时刻t状态为i的所有单个路径 $(i_1,i_2,\cdots,i_t)$ 中概率最大值为：
$\delta_t(i)=\underset{i_1,i_2,\cdots,i_{t-1}}{\operatorname{max}}P(i_t=i,i_{t-1},\cdots,i_1,o_t,\cdots,o_1|\lambda),i=1,2,\cdots,N$
由定义可得变量 $\delta$ 的递推公式：
$\delta_{t+1}(i)=\underset{i_1,i_2,\cdots,i_t}{\operatorname{max}}P(i_{t+1}=i,i_t,\cdots,i_1,o_{t+1},\cdots,o_1|\lambda)=\underset{1\leq j \leq N}{\operatorname{max}}[\delta_t(j)a_{ji}]b_i(o_{t+1}),　　i=1,2,\cdots,N;t=1,2,\cdots,T-1$
定义在时刻t状态为i的所有单个路径 $(i_1,i_2,\cdots,i_{t-1},i)$ 中概率最大的路径的第t-1个结点为
$\psi_t(i)=\underset{1\leq j \leq N}{\operatorname{argmax}}[\delta_{t-1}(j)a_{ji}],　　i=1,2,\cdots,N$
输入：模型 $\lambda=(A,B,\pi)$ 和观测 $O=(o_1,o_2,\cdots,o_T)$ ；

输出：最优路径 $I^*=(i_1^*,i_2^*,\cdots,i_T^*)$ .

①初始化
$\delta_1(i)=\pi_ib_i(o_1),　　i=1,2,\cdots,N \\ \psi_1(i)=0,　　i=1,2,\cdots,N$
②递推，对 $t=2,3,\cdots,T$
$\delta_t(i)=\underset{1\leq j \leq N}{\operatorname{max}}[\delta_{t-1}(j)a_{ji}]b_i(o_t),　　i=1,2,\cdots,N \\ \psi_t(i)=\underset{1\leq j \leq N}{\operatorname{argmax}}[\delta_{t-1}(j)a_{ji}],　　i=1,2,\cdots,N$
③终止
$P*=\underset{1\leq i \leq N}{\operatorname{max}}\delta_T(i) \\ i_T^*=\underset{1\leq i \leq N}{\operatorname{argmax}}[\delta_T(i)]$
④最优路径回溯，对 $t=T-1,T-2,\cdots,1$
$i_t^*=\psi_{t+1}(i_{t+1}^*)$
求得最优路径 $I^*=(i_1^*,i_2^*,\cdots,i_T^*)$
示例：
$A= \left[\begin{matrix} 0.5 & 0.2 & 0.3 \\ 0.3 & 0.5 & 0.2 \\ 0.2 & 0.3 & 0.5 \end{matrix}\right]， B = \left[\begin{matrix} 0.5 & 0.5 \\ 0.4 & 0.6 \\ 0.7 & 0.3 \end{matrix}\right]， \pi = (0.2,0.4,0.4)^T$
O=(红，白，红)，试求最优状态序列，即最优路径 $I^*=(i_1^*,i_2^*,i_3^*)$

①初始化：在t=1时，对每一个状态i，i=1,2,3，求状态i观测O1为红的概率，记为： $\delta_1(i)=\pi_ib_i(o_1)=\pi_ib_i(红),i=1,2,3$

代入实际数据： $\delta_1(1)=0.10，\delta_1(2)=0.16，\delta_1(3)=0.28$ ，记 $\psi_1(i)=0,i=1,2,3$

②在t=2时，对每一个状态i，i=1,2,3，求在t=1时状态为j观测O1为红并在t=2时状态为i观测O2为白的路径的最大概率，记为 $\delta_2(i)$
$\delta_2(i)=\underset{1\leq j \leq 3}{\operatorname{max}}[\delta_1(j)a_{ji}]b_i(o_2)$
同时，对每个状态i，记录概率最大路径的前一个状态j， $\psi_2(i)=\underset{1 \leq j \leq 3}{\operatorname{argmax}}[\delta_1(j)a_{ji}]，i=1,2,3$

计算：
$\delta_2(1)=\underset{1 \leq j \leq 3}{\operatorname{max}}[\delta_1(j)a_{j1}]b_1(o_2)=\underset{j}{\operatorname{max}}\{0.10\times0.5,0.16\times0.3,0.28\times0.2\}\times0.5=0.028 \\ \psi_2(1)=3 \\ \delta_2(2)=0.0504，\psi_2(2)=3 \\ \delta_2(3)=0.042，\psi_2(3)=3$
同样，在t=3时
$\delta_3(i)=\underset{1 \leq j \leq 3}{\operatorname{max}}[\delta_2(j)a_{ji}]b_i(o_3) \\ \psi_3(i)=\underset{1\leq j \leq 3}{\operatorname{argmax}}[\delta_2(j)a_{ji}] \\ \delta_3(1) = 0.00756，\psi_3(1)=2 \\ \delta_3(2) = 0.01008，\psi_3(2)=2 \\ \delta_3(3) = 0.0147，\psi_3(3)=3$
③以P*表示最优路径的概率：
$P^*=\underset{1 \leq i \leq 3}{\operatorname{max}}\delta_3(i)=0.0147$
最优路径的终点是：
$i_3^* = \underset{i}{\operatorname{argmax}}[\delta_3(i)]=3$
④由最优路径的终点 $i_3^*$ ，逆向找到 $i_2^*,i_1^*$
$在t=2时，i_2^*=\psi_3(i_3^*)=\psi_3(3)=3 \\ 在t=1时，i_1^*=\psi_2(i_2^*)=\psi_2(3)=3 \\$
于是求得最优路径，即最优状态序列：
$I^*=(i_1^*,i_2^*,i_3^*)=(3,3,3)$

隐马尔科夫模型Python实现

import numpy as np

class HMM(object):
    def __init__(self, N, M, pi=None, A=None, B=None):
        self.N = N
        self.M = M
        self.pi = pi
        self.A = A
        self.B = B

    def get_data_with_distribute(self, dist): # 根据给定的概率分布随机返回数据（索引）
        r = np.random.random()
        for i, p in enumerate(dist):
            if r < p: return i
            r -= p

    def generate(self, T: int):
        '''
        根据给定的参数生成观测序列
        T: 指定要生成观测序列的长度
        '''
        result = []
        i = 0
        for ind in range(T):        # 依次生成状态和观测数据
            if ind==0:
                i = self.get_data_with_distribute(self.pi)
            else:
                i = self.get_data_with_distribute(self.A[i,:])
            o = self.get_data_with_distribute(self.B[i,:])
            result.append(o)
        return result
    
    def forward(self,X):
        '''
        根据给定的参数计算条件概率(前向算法)
        X：观测数据
        '''
        alpha = self.pi * self.B[:,X[0]]
        for x in X[1:]:
            #将alpha构造成(N,1)维度，使用广播机制
            alpha = np.sum(self.A * alpha[...,np.newaxis],axis=0) * self.B[:,x]
        
        return alpha.sum()
    
    def backward(self,X):
        '''
        根据给定的参数计算条件概率(后向算法)
        X：观测数据
        '''
        beta = np.ones(self.N)
        for x in X[:0:-1]:
            #beta默认为(1,N)维度，使用广播机制
            beta = np.sum(self.A * self.B[:,x] * beta,axis=1)

        return np.sum(beta * self.pi * self.B[:,X[0]],axis=0)
    
    def get_alpha_beta(self,X):
        '''
        根据给定数据与参数，计算所有时刻的前向概率和后向概率
        '''
        T = len(X)
        alpha = np.zeros((T,self.N))
        alpha[0,:] = self.pi * self.B[:,X[0]]
        for t in range(T-1):
            alpha[t+1,:] = np.sum(self.A * alpha[t,:][...,np.newaxis],axis=0) * self.B[:,X[t+1]]
            
        beta = np.ones((T,self.N))
        for t in range(T-1,0,-1):
            beta[t-1,:] = np.sum(self.A * beta[t,:] * self.B[:,X[t]],axis=1)

        # print(alpha[T-1,:].sum(),np.sum(beta[0,:] * self.pi * self.B[:,X[0]],axis=0))        
        return alpha,beta
    
    def fit(self,X):
        '''
        根据观测序列估计参数
        '''
        #初始化参数pi，A，B
        self.pi = np.random.random(self.N)
        self.pi /= np.sum(self.pi)
        self.A = np.random.random((self.N,self.N))
        self.A /= np.sum(self.A,axis=1)[:,np.newaxis]
        self.B = np.random.random((self.N,self.M))
        self.B /= np.sum(self.B,axis=1)[:,np.newaxis]
        
        T = len(X)
        for epoch in range(10):
            #根据公式计算下一时刻参数值
            alpha,beta = self.get_alpha_beta(X)
            gamma = alpha * beta
            for i in range(self.N):
                for j in range(self.N):
                    self.A[i,j] = np.sum(alpha[:-1,i] * beta[1:,j] * self.A[i,j] * self.B[j,X[1:]],axis=0) / np.sum(gamma[:-1,i],axis=0)
            for i in range(self.N):
                for j in range(self.M):
                    self.B[i,j] = np.sum(gamma[:,i] * (X == j),axis=0) / np.sum(gamma[:,i],axis=0)
                    
            self.pi = gamma[0] / gamma[0].sum()
            
    def decode(self,X):
        '''
        解码问题，参数已知，根据给定的观测序列，返回最大概率的隐藏序列
        '''
        T = len(X)
        x = X[0]
        delta = self.pi[:,np.newaxis] * self.B[:,x]
        varphi = np.zeros((T,self.N),dtype=np.int)
        path = [0] * T
        for t in range(1,T):
            # delta = delta.reshape(-1,1)
            tmp = delta * self.A
            varphi[t,:] = np.argmax(tmp,axis=0)
            delta = np.max(tmp,axis=0) * self.B[:,X[t]]

        path[-1] = np.argmax(delta)
        #回溯最优路径
        for t in range(T-1,0,-1):
            path[t-1] = varphi[t,path[t]]
            
        return path
 
if __name__ == "__main__":
    pi = np.array([0.25, 0.25, 0.25, 0.25])
    A = np.array([
        [0,  1,  0, 0],
        [0.4, 0, 0.6, 0],
        [0, 0.4, 0, 0.6],
        [0, 0, 0.5, 0.5]])
    B = np.array([
        [0.5, 0.5],
        [0.3, 0.7],
        [0.6, 0.4],
        [0.8, 0.2]])
    hmm = HMM(4, 2, pi, A, B)
    print(hmm.generate(10))  # 生成10个数据
    X = [0,0,1,1,0] #观测序列为红，红，白，白，红
    print(hmm.forward(X))
    print(hmm.backward(X))

    import matplotlib.pyplot as plt
    def triangle_data(T):   # 生成三角波形状的序列
        data = []
        for x in range(T):
            x = x % 2
            data.append(x if x <= 1 else 2-x)
        return data
    
    data = np.array(triangle_data(15))
    hmm.fit(data)
    gen_obs = hmm.generate(15)  # 再根据学习的参数生成数据
    x = np.arange(15)
    plt.scatter(x, gen_obs, marker='*', color='r')
    plt.plot(x, data, color='g')
    plt.show()

    pi = np.array([0.2, 0.4, 0.4])
    A = np.array([
        [0.5,  0.2,  0.3],
        [0.3, 0.5, 0.2],
        [0.2, 0.3, 0.5]])
    B = np.array([
        [0.5, 0.5],
        [0.4, 0.6],
        [0.7, 0.3]])
    O = [0,1,0]
    hmm = HMM(3, 2, pi, A, B)
    path = hmm.decode(O)
    print(path) #输出(2,2,2)为索引

总结

隐马尔可夫模型是关于时序的概率模型，由初始状态概率向量 $\pi$ 、状态转移概率矩阵A和观测概率矩阵B决定，可以写成 $\lambda=(A,B,\pi)$ 的形式，隐马尔可夫模型是一个生成模型，它的应用很多，在分词、词性标注、NLP、语音识别等领域都有着广泛的应用，是一个非常经典的机器学习模型。

最后编辑于：2021.06.30 14:06:23

禁止转载，如需转载请通过简信或评论联系作者。

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 204,189评论 6赞 478
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 85,577评论 2赞 381
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 150,857评论 0赞 337
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 54,703评论 1赞 276
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 63,705评论 5赞 366
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,620评论 1赞 281
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 37,995评论 3赞 396
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,656评论 0赞 258
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 40,898评论 1赞 298
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,639评论 2赞 321
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,720评论 1赞 330
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,395评论 4赞 319
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 38,982评论 3赞 307
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,953评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,195评论 1赞 260
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 44,907评论 2赞 349
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,472评论 2赞 342

机器学习之隐马尔科夫模型（HMM）