隐马尔科夫模型

基本概念​

隐马尔可夫模型是关于时序的概率模型，描述由一个隐藏的马尔可夫链随机生成不可观察的状态随机序列，再由各个状态生成一个观测而产生随机序列的过程。隐藏的马尔可夫链随机生成的状态的序列，称为状态序列；每个状态生成一个观测，而由此产生的观测的随机序列，称为观测序列。序列的每一个未知又可以看做是一个时刻。​

用途​

隐马尔科夫模型是可用于标注问题的统计学习模型，描述由隐藏的马尔可夫链随机生成观测序列的过程，属于生成模型。

前提

随机变量

随机变量（Random Variable），通常用大写字母来表示一个随机事件。比如看下面的例子：

$X$：河水是咸的；$Y$：井水是甜的

很显然，$X$与$Y$两个随机事件是没有关系的。也就是说和之间是相互独立的。记作：

$$X \perp Y$$

随机过程

对于一类随机变量来说，它们之间存在着某种关系。

比如： ：表示在 时刻某支股票的价格，那么$S_{t+1}$和$S_t$之间一定是有关系的，至于具体什么样的关系，这里原先不做深究，但有一点可以确定，两者之间一定存在的一种关系。随着时间$t$的变化，可以写出下面的形式：

$$\cdots S_t,S_{t+1},S_{t+2}\cdots$$

这样就生成了一组随机变量，它们之间存在着一种相当复杂的关系，也就是说，各个随机变量之间存在着关系，即不相互独立。由此，我们会把按照某个时间或者次序上的一组不相互独立的随机变量的这样一个整体作为研究对象。这样的话，也就引出了另外的一个概念：随机过程（Stochastic Process）。也就是说随机过程的研究对象不在是单个的随机变量，而是一组随机变量，并且这一组随机变量之间存在着一种非常紧密的关系（不相互独立）。记作：

$$\{S_t\}^{\infty}_{t=1}$$

马尔可夫链/过程

**马尔科夫链（Markov Chain）**即马尔可夫过程，是一种特殊的随机过程——具备马尔可夫性的随机过程。

• 马尔可夫性：（Markov Property）: 如果满足$P(S_{t+1}|S_t,S_{t-1}\cdots S_1)=P(S_{t+1}|S_t)$，即具备了马尔可夫性。简单来说， 和之间存在关系，和以前的时刻的没有关系，即只和“最近的状态” 有关系。

• 现实例子：**下一个时刻仅依赖于当前时刻，跟过去无关。**比如：一个老师讲课，明天的讲课状态一定和今天的状态最有关系，和过去十年的状态基本就没关系了。

• 最主要考量：为了简化计算。$P(S_{t+1}|S_t,S_{t-1}\cdots S_1)=P(S_{t+1}|S_t)$ ，如果$S_{t+1}$和 $S_t,S_{t-1}\cdots S_1$都有关系的话，计算就会太大了，可能造成无法计算。

状态空间模型

状态空间模型（State Space Model），常应用于 HMM,Kalman Filterm Particle Filter。在这里就是指马尔可夫链 + 观测变量，即Markov Chain + Obervation。

image-20210125210024542

如上图所示，s1-s2-s3为马尔可夫链，a1, a2, a3为观测变量，以a2为例，a2只和s2有关和s1, s3无关。状态空间模型可以说是由马尔可夫链演化而来的模型。记作：

$$P(S_{t+1}|S_t) = P(S_{t+1}|S1,\cdots,S_t)$$

马尔可夫奖励过程

马尔可夫奖励过程（Markov Reward Process），即马尔可夫链+奖励，即：Markov Chain + Reward。如下图：

image-20210125210341550

举个例子，比如说你买了一支股票，然后你每天就会有“收益”，当然了这里的收益是泛化的概念，收益有可能是正的，也有可能是负的，有可能多，有可能少，总之从今天的状态$S_t$到明天的状态 $s_{t+1}$，会有一个reward。记作：

$$R_s = E[R_{t+1}|S_t]$$

模型定义​

image-20210104222033435

隐马尔可夫模型由初始概率分布、状态转移概率分布以及观测概率分布确定。隐马尔可夫模型的形式定义如下：

​ 设$Q$是所有可能的状态集合，$V$是所有可能的观测集合。

$$Q = \{q_1,q_2,\cdots,q_N\}，V = \{v_1,v_2,\cdots,v_M\}$$

其中，$N$是可能的状态数，$M$是可能的观测数。

​ $I$是长度为$T$的状态序列，$O$是对应的观测序列。

$$I = (i_1,i_2,\cdots,i_T)，O = (o_1,o_2,\cdots,o_T)$$

​ $A$是状态转移概率矩阵：

$$A = [a_{ij}]_{N\times N}$$

其中，

$$a_{ij} = 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_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$的概率。

$(i)$ 隐马尔可夫模型$\lambda = (\pi,A,B)$

隐马尔可夫模型由初始状态概率向量$\pi$、状态转移概率矩阵$A$以及观测概率矩阵$B$决定。$\pi$和$A$决定状态序列，$B$决定观测序列。因此，隐马尔可夫模型$\lambda$由如下表示：

$$\lambda = (\pi,A,B)$$

$(ii)$ 从定义可知，隐马尔可夫模型作了两个基本假设：

1. 齐次马尔可夫性假设：假设隐藏的马尔科夫链在任意时刻$t$的状态只依赖于其前一时刻的状态，与其他的时刻的状态及观测无关，也与时刻$t$无关。

$$P(i_{t+1}|i_{1\cdots t},o_{1\cdots t}) = P(i_{t+1}|i_t)$$

image-20210104222912358

1. 观测独立性假设：假设任意时刻的观测只依赖于该时刻的马尔可夫链的状态，与其他观测及状态无关。

$$P(O_t|i_{1\cdots t},o_{1\cdots t}) = P(O_t|i_t)$$

image-20210104223203109

$(iii)$ 三个问题

$(1).$ 概率计算问题。给定$\lambda$，求$P(O|\lambda)$

$(2).$ 学习问题

$(3).$ 预测问题（解码问题）

概率计算问题

直接计算法

状态序列$I=(i_1,i_2,\cdots,i_T)$的概率是：

$$P(I|\lambda)= P(i_{1\cdots T}|\lambda) = P(i_{1\cdots T-1},i_T|\lambda)=P(i_T|i_{T-1})\cdot P(i_{1\cdots T}|\lambda) = \pi_{i_1}a_{i_1i_2}a_{i_2i_3}\cdots a_{i_{T-1}i_T} = \pi_{i_1}\prod_{t=2}^Ta_{i_{t-1}i_t}$$

状态$i_t$的状态转移概率分布{$a_{i_t,i_{t+1}}$}产生状态$i_{t+1}$

对固定的状态序列$I=(i_1,i_2,\cdots,i_T)$，观测序列$O=(o_1,o_2,\cdots,o_T)$的概率是$P(O|I,\lambda)$：

$$P(O|I,\lambda) = b_{i_1}(o_1)b_{i_2}(o_2)\cdots b_{i_T}(o_T)$$

状态$i_t$的观测概率分布$b_{i_t}(k)$生成$o_t$

然后，对所有可能的状态序列$I$求和，得到观测序列$O$的概率$P(O|\lambda)$，即

$$P(O|\lambda)=\sum_{I} P(I,O|\lambda)=\sum_{I} P(I|\lambda)P(O|I,\lambda)=\sum_{i_1,i_2,\cdots,i_T}\pi_{i_1}\prod_{t=2}^T{a_{i_{t-1},i_t}}\prod_{t=1}^Tb_{i_t}(o_t)$$

由于空间复杂度为$O(TN^T)$，故该算法理论上可行，实际不可行。

前向算法

image-20210105100709601

前向概率定义：给定隐马尔可夫模型$\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_i|\lambda)$$

有，

$$\alpha_T(i) = P(O,i_t=q_i|\lambda)$$

算法过程：

输入：隐马尔可夫模型$\lambda$，观测序列$O$；

输出：观测序列概率$P(O|\lambda)$。

1. 初值

$$\alpha_1(i) = \pi_ib_i(o_1)，i=1,2,\cdots,N$$

1. 递推，对$t=1,2,\cdots,T-1$

$$\begin{align} \alpha_{t+1}(j) &= P(o_1,\cdots,o_t,o_{t+1},i_{t+1}=q_j|\lambda) \\ &= \sum_{i=1}^N P(o_1,\cdots,o_t,o_{t+1},i_{t+1}=q_j,i_t=q_i|\lambda) \\ &= \sum_{i=1}^N P(o_{t+1}|o_1,\cdots,o_t,i_t=q_i,i_{t+1}=q_j,\lambda)\cdot P(o_1,\cdots,o_t,i_t=q_i,i_{t+1}=q_j|\lambda) \\ &= \sum_{i=1}^N P(o_{t+1}|i_{t+1}=q_j,\lambda)\cdot P(i_{t+1}=q_j|o_1,\cdots,o_t,i_t=q_i,\lambda)\cdot P(o_1,\cdots,o_t,i_t=q_i|\lambda) \\ &= \sum_{i=1}^N P(o_{t+1}|i_{t+1}=q_j,\lambda)\cdot P(i_{t+1}=q_j|i_t=q_i,\lambda)\cdot P(o_1,\cdots,o_t,i_t=q_i|\lambda) \\ &= \sum_{i=1}^N b_j(o_{t+1})\cdot a_{ij}\cdot \alpha_t(i) \\ &= \sum_{i=1}^N [\alpha_t(i)a_{ij}]b_j(o_{t+1}) \end{align}$$

1. 终止

$$P(O|\lambda) = \sum_{i=1}^N \alpha_T(i)$$

后向算法

image-20210105134637569

后向概率定义：给定隐马尔可夫模型$\lambda$，定义在时刻$t$状态为$q_i$的条件下，从$t+1$到$T$的部分观测序列为$o_{t+1},o_{t+2},\cdots,o_T$的概率为后向概率，记作：

$$\beta_t(i) = P(o_{t+1},o_{t+2},\cdots,o_T|i_t=q_i,\lambda)$$

算法过程：

输入：隐马尔可夫模型$\lambda$，观测序列$O$；

输出：观测序列概率$P(O|\lambda)$。

1. 初值

$$\beta_T(i) = 1，i=1,2,\cdots,N$$

1. 递推，对$t=T-1,T-2,\cdots,1$

$$\begin{align} \beta_t(i) &= P(o_{t+1},o_{t+2},\cdots,o_T|i_t=q_i,\lambda) \\ &= \sum_{j=1}^N P(o_{t+1},o_{t+2},\cdots,o_T,i_{t+1}=q_j|i_t=q_i,\lambda) \\ &= \sum_{j=1}^N P(o_{t+1},o_{t+2},\cdots,o_T|i_t=q_i,i_{t+1}=q_j,\lambda)\cdot P(i_{t+1}=q_j|i_t=q_i,\lambda) \\ &= \sum_{j=1}^N P(o_{t+2},\cdots,o_T|i_{t+1}=q_j,\lambda)\cdot P(o_{t+1}|o_{t+2},\cdots,o_T,i_{t+1}=q_j,i_t=q_i,\lambda)\cdot P(i_{t+1}=q_j|i_t=q_i,\lambda) \\ &= \sum_{j=1}^N P(o_{t+1}|i_{t+1}=q_j,\lambda)\cdot P(i_{t+1}=q_j|i_t=q_i,\lambda)\cdot P(o_{t+2},\cdots,o_T|i_{t+1}=q_j,\lambda) \\ &= \sum_{j=1}^N b_j(o_{t+1})a_{ij}\beta_{t+1}(j) \end{align}$$

1. 终止

$$\begin{align} P(O|\lambda) &= P(o_1,\cdots,o_T|\lambda) \\ &= \sum_{i=1}^N P(o_1,\cdots,o_T,i_1=q_i|\lambda) \\ &= \sum_{i=1}^N P(o_1,\cdots,o_T|i_1=q_i,\lambda)\cdot P(i_1=q_i|\lambda) \\ &= \sum_{i=1}^N P(o_1|o_2,\cdots,o_T,i_1=q_i,\lambda)\cdot P(o_2,\cdots,o_T|i_1=q_i,\lambda) \cdot \pi_i \\ &= \sum_{i=1}^N P(o_1|i_1=q_i,\lambda)\cdot P(o_2,\cdots,o_T|i_1=q_i,\lambda) \cdot \pi_i \\ &= \sum_{i=1}^N b_i(o_1)\beta_1(i)\pi_i \end{align}$$

其他概率与期望值的计算

利用前向概率和后向概率，可以得到关于单个状态和两个状态概率的计算公式。

1. 由前向概率$\alpha_t(i)$和后向概率$\beta_t(i)$定义可知：

$$\begin{align} \alpha_t(i)\beta_t(i) &= P(o_1,\cdots,o_t,i_t=q_i|\lambda)\cdot P(o_{t+1},\cdots,o_T|i_t=q_i,\lambda) \\ &= P(o_1,\cdots,o_t,i_t=q_i|\lambda)\cdot P(o_{t+1},\cdots,o_T|i_t=q_i,o_1,\cdots,o_t,\lambda) \\ &= P(o_1,\cdots,o_t,o_{t+1},\cdots,o_T,i_t=q_i|\lambda) \\ &= P(O,i_t=q_i|\lambda) \end{align}$$

1. 给定模型$\lambda$和观测$O$，在时刻$t$处于状态$q_i$的概率，记为：

$$\begin{align} \gamma_t(i) &= P(i_t=q_i|O,\lambda) \\ &=\frac{P(i_t=q_i,O|\lambda)}{P(O|\lambda)} \\ &=\frac{P(i_t=q_i,O|\lambda)}{\sum_{j=1}^N P(i_t=q_j,O|\lambda)} \\ &=\frac{\alpha_t(i)\beta_t(i)}{\sum_{j=1}^N \alpha_t(j)\beta_t(j)} \end{align}$$

1. 给定模型$\lambda$和观测$O$，在时刻$t$处于状态$q_i$且在时刻$t+1$处于状态$q_j$的概率。记为：

$$\begin{align} \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{P(i_t=q_i,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)} \end{align}$$

其中，

$$\begin{align} P(i_t=q_i,i_{t+1}=q_j,O|\lambda) &= P(i_t=q_i,o_1,\cdots,o_t,i_{t+1}=q_j,o_{t+1},\cdots,o_T | \lambda) \\ &= P(i_t=q_i,o_1,\cdots,o_t|\lambda)\cdot P(i_{t+1}=q_j,o_{t+1},\cdots,o_T|i_t=q_i,o_1,\cdots,o_t,\lambda) \\ &= \alpha_t(i)\cdot P(i_{t+1}=q_j|i_t=q_i,O,\lambda)\cdot P(o_{t+1},\cdots,o_T|i_{t+1}=q_j,i_t=q_i,o_1,\cdots,o_t,\lambda) \\ &= \alpha_t(i)\cdot a_{ij}\cdot P(o_{t+1},o_{t+2},\cdots,o_T|i_{t+1}=q_j,i_t=q_i,o_1,\cdots,o_t,\lambda) \\ &= \alpha_t(i)\cdot a_{ij}\cdot P(o_{t+1}|i_{t+1}=q_j,i_t=q_i,o_1,\cdots,o_t)\cdot P(o_{t+2},\cdots,o_T|i_{t+1}=q_j,i_t=q_i,o_1,\cdots,o_t,o_{t+1}) \\ &= \alpha_t(i)\cdot a_{ij}\cdot P(o_{t+1}|i_{t+1})\cdot P(o_{t+2},\cdots,o_T|i_{t+1}=q_j) \\ &= \alpha_t(i) a_{ij} b_j(o_{t+1})\beta_{t+1}(j) \end{align}$$

1. 将$\gamma_t(i)$和$\xi_t(i,j)$对各个时刻的$t$求和，可以得到一些有用的期望值。

$(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 \xi_t(i,j)$$

学习问题

隐马尔可夫模型的学习，根据训练数据是包括观测序列和对应的状态序列还是只有观测序列，可以分别由监督学习与非监督学习实现。

监督学习算法

假设已给训练数据包含$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_{i=1}^N A_{ij}}，i=1,2,\cdots,N;j=1,2,\cdots,N$$

1. 观测概率$b_j(k)$的估计

设样本中状态为$j$并观测为$k$的频数是$B_{jk}$，那么状态为$j$观测为$k$的概率$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$$

1. 初始化状态概率$\pi_i$的估计$\hat{\pi_i}$为$S$个样本中初始状态为$q_i$的概率

由于监督学习需要使用训练数据，而人工标注训练数据往往代价很高，有时就会利用非监督学习的方法。

非监督学习算法： $Baum-Welch$算法

假设给定训练数据只包含$S$个长度为$T$的观测序列${O_1,O_2,\cdots,O_S}$而没有对应的状态序列，目标是学习隐马尔可夫模型$\lambda=(A,B,\pi)$的参数。将观测序列数据看做观测数据$O$，状态序列数据看做不可预测的隐数据$I$，那么因马尔可夫模型事实上是一个含有隐变量的概率模型：

$$P(O|\lambda) = \sum_I P(O|I,\lambda)P(I|\lambda)$$

该模型的参数可以由$EM$算法实现。

算法过程：

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

所有观测数据写成$O=(o_1,o_2,\cdots,o_T)$，所有隐数据写成$I=(i_1,i_2,\cdots,i_T)$，完全数据是$(O,I)=(o_1,o_2,\cdots,o_T,i_1,i_2,\cdots,i_T)$。完全数据的对数似然函数是$log\ P(O,I|\lambda)$。

1. $EM$算法的$E$步：求$Q$函数$Q(\lambda,\lambda^{(t)})$

$$Q(\lambda,\lambda^{(t)}) = \sum_I lop\ P(O,I|\lambda)P(O,I|\lambda^{(t)})$$

其中，$\lambda^{(t)}$是隐马尔可夫模型参数的当前估计值，$\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)=\pi_{i_1}\prod_{t=2}^Ta_{i_{t-1},i_t}\prod_{t=1}^T b_{i_t}(o_t)$$

于是：

$$\begin{align} Q(\lambda,\lambda^{(t)}) &= \sum_I[log \ \pi_{i_1}\prod_{t=2}^Ta_{i_{t-1},i_t}\prod_{t=1}^T b_{i_t}(o_t)]P(O,I|\lambda^{(t)}) \\ &= \sum_I[log \ \pi_{i_1}+ \sum_{t=2}^T log\ a_{i_{t-1},i_t} + \sum_{t=1}^T log\ b_{i_t}(o_t)]P(O,I|\lambda^{(t)}) \\ \end{align}$$

1. $EM$算法的$M$步：极大化$Q$函数$Q(\lambda,\lambda^{(t)})$求模型参数$A$，$B$，$\pi$

$$\lambda^{(t+1)} = {argmax}_{\lambda}Q(\lambda,\lambda^{(t)})= argmax_{\lambda} \sum_I[\textcolor{blue}{log \ \pi_{i_1}}+ \textcolor{green}{\sum_{t=2}^T log\ a_{i_{t-1},i_t}} + \textcolor{orange}{\sum_{t=1}^T log\ b_{i_t}(o_t)}]P(O,I|\lambda^{(t)})$$

由上可知，三个参数分别单独在三个部分，在求解其中一个时，可以忽略掉其他两个，单独讨论。

对于$\pi$：

$$\sum_I \pi_{i_1}P(O,I|\lambda) = \sum_{i=1}^N log\ \pi_iP(O,i_1=i|\lambda^{(t)})$$

由于$\pi$是一个概率分布，因此有一个约束条件$\sum_{i=1}^N \pi_i=1$，于是利用拉格朗日乘子法，构造拉格朗日函数，有：

$$\sum_{i=1}^N log\ \pi_iP(O,i_1=i|\lambda^{(t)})+\gamma(\sum_{i=1}^N \pi_i-1)$$

对其求偏导数并令结果为0，有：

$$\frac{\partial}{\partial\pi_i}[\sum_{i=1}^N log\ \pi_iP(O,i_1=i|\lambda^{(t)})+\gamma(\sum_{i=1}^N \pi_i-1)] = 0$$

得：

$$\begin{align} \frac{1}{\pi_i}P(O,i_1=i|\lambda^{(t)}) + \gamma &= 0 \\ P(O,i_1=i|\lambda^{(t)}) + \gamma\pi_i&=0 \end{align}$$

同时进行连加，得到$\gamma$：

$$\sum_{i=1}^N P(O,i_1=i|\lambda^{(t)}+\gamma\pi_i) = P(O|\lambda^{(t)}) + \gamma = 0$$

有：

$$\gamma = -P(O|\lambda^{(t)})$$

带入$P(O,i_1=i|\lambda^{(t)}) + \gamma\pi_i=0$得：

$$\pi_i^{(t+1)} = \frac{P(O,i_1=i|lambda^{(t)})}{P(O|\lambda^{(t)})}$$

对于$a_{ij}$：

$$\sum_I \sum_{t=2}^T log\ a_{i_{t-1},i_t}P(O,I|\lambda) = \sum_{i=1}^N\sum_{j=1}^N \sum_{t=1}^{T-1} log\ a_{ij} log\ P(O,i_t=i,i_{t+1}=j|\lambda^{(t)})$$

由于$A$是一个概率分布的矩阵，其约束条件为$\sum_{j=1}^N a_{ij}=1$，与上述相同，利用拉格朗日乘子法，构造拉格朗日函数，求出得：

$$a_{ij}^{(t+1)} = \frac{P(O,i_t=i,i_{t+1}=j|\lambda^{(t)}}{\sum_{t=1}^{T-1}P(O,i_t=i|\lambda^{(t)})}$$

对于$b_j(k)$：

$$\sum_I[\sum_{t=1}^T log\ b_{i_t}(o_t)]P(O,I|\lambda^{(t)}) = \sum_{j=1}^N\sum_{t=1}^T log\ b_{j}(o_t)P(O,i_t=j|\lambda^{(t)})$$

由于$B$是一个概率分布的矩阵，其约束条件为$\sum_{k=1}^M b_j(o_k)=1$，与上述相同，利用拉格朗日乘子法，构造拉格朗日函数，求出得：

$$\begin{align} \frac{\partial}{\partial b_j(k)}[\sum_{j=1}^N\sum_{t=1}^T log\ b_{j}(o_t)P(O,i_t=j|\lambda^{(t)})+\sum_{i=1}^N\eta_i(\sum_{k=1}^M (b_j(k)-1)] &= 0 \\ \sum_{t=1}^T \frac{1}{b_j(k)}P(O,i_t=j|\lambda^{(t)})\textcolor{blue}{I(o_t=v_k)}+\sum_{i=1}\eta_i &= 0 \end{align}$$

注意，只有在$o_t=v_k$时，$b_j(o_t)$对$b_j(k)$的偏导数才不为0，以$I(o_t=v_k)$表示，该函数的作用是满足条件值为1，否则为0。

于是：

$$b_j(k)^{(t+1)} = -\frac{\sum_{t=1}^T P(O,i_t=j|\lambda^{(t)}){I(o_t=v_k)}}{\sum_{i=1}\eta_i}$$

$Baum-Welch$模型的参数估计公式

利用其它概率与期望值的计算的公式，进行改写为如下式子：

$$\begin{align} a_{ij} &= \frac{\sum_{t=1}^T\xi_t(i,j)}{\sum_{t=1}^T\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) \end{align}$$

预测算法

求概率最大的状态序列。

维特比算法

利用动态规划求解隐马尔可夫模型预测问题，即用动态规划求概率最大路径（最优路径）。这时一条路径对应着一个状态序列。

定义在时刻$t$状态为$i$的所有单个路径$(i_1,i_2,\cdots,i_t)$中概率最大值为：

$$\delta_t(i) = max_{i_1,i_2,\cdots,i_{t-1}} P(i_t=i,i_{t-1},\cdots,i_1,o_t,\cdots,o_1|\lambda)，i=1,2,\cdots,N$$

由定义可得变量$\delta$的递推公式为：

$$\begin{align} \delta_{t+1}(i) &= max_{i_1,i_2,\cdots,i_{t}} P(i_{t+1}=i,i_{t},\cdots,i_1,o_{t+1},\cdots,o_1|\lambda) \\ &= max_{1\leq j \leq N}[\delta_t(j)a_{ji}]b_i(o_{t+1})，i=1,2,\cdots,N;t=1,2,\cdots,T-1 \end{align}$$

定义在时刻$t$状态为$i$的所有单个路径$(i_1,i_2,\cdots,i_{t-1},i)$中概率最大的路径的第$t-1$个结点为：

$$\psi_t(i) = argmax_{1 \leq j \leq N}[\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)$。

1. 初始化

$$\delta_1(i) = \pi_ib_i(o_1)，i=1,2,\cdots,N \\ \psi_1(i) = 0，i=1,2,\cdots,N$$

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

$$\delta_{t}(i) = max_{1\leq j \leq N}[\delta_{t-1}(j)a_{ji}]b_i(o_{t})，i=1,2,\cdots,N \\ \psi_t(i) = argmax_{1 \leq j \leq N}[\delta_{t-1}(j)a_{ji}]，i=1,2,\cdots,N$$

1. 终止

$$P^* = max_{1\leq i \leq N} \delta_T(i) \\ i^*_T = argmax_{1\leq i \leq N} [\delta_T(i)]$$

1. 最有路径回溯，对$t=T-1,T-2,\cdots,1$

$$i^*_t = \psi_{t+1}(i^*_{t+1})$$

求得最优路径$I^=(i^_1,i^_2,\cdots,i^_T)$。

代码实现

import logging
import numpy as np
np.set_printoptions(suppress=True, threshold=np.inf,precision=10)

class HMM:
    def __init__(self,M,N,epoches=10):
        self.A  = np.abs(np.random.randn(M,M))
        self.B  = np.abs(np.random.randn(M,N))
        self.pi = np.ones(M)
        self.T  = None
        self.O  = None
        self.M  = M # 可能的状态数
        self.N  = N # 可能的观察数
        self.alpha = None
        self.beta = None
        self.epoches = epoches
        self.logger = logging.getLogger(__name__)
        self.logger.setLevel(logging.DEBUG)
    
    def forward(self,):
        alpha = np.zeros((self.T,self.M))
        for t in range(self.T):
            if t == 0:
                alpha[t,:] = self.B[:,O[t]]*self.pi # 初始化
                continue
            alpha[t,:] = np.array([self.A[:,i].dot(alpha[t-1,:])*self.B[i,O[t]] for i in range(self.M)])
        return alpha

    def backward(self,):
        beta = np.zeros((self.T,self.M))
        beta[-1,:] = np.ones(self.M)
        for t in list(range(self.T-2,-1,-1)):
            beta[t,:] = np.array([np.sum(np.array([self.A[i,j]*self.B[j,O[t+1]]*beta[t+1,j] for j in range(self.M)])) for i in range(self.M)])
        return beta

    def calc_gamma(self,t,i):
        gamma = (self.alpha[t,i]*self.beta[t,i]) / np.sum(np.array([self.alpha[t,j]*self.beta[t,j] for j in range(self.M)]))
        return gamma

    def calc_kesai(self,t,i,j):
        kesai = (self.alpha[t,i]*self.A[i,j]*self.B[j,self.O[t+1]]*self.beta[t+1,j]) / np.sum(np.array([[self.alpha[t,i]*self.A[i,j]*self.B[j,self.O[t+1]]*self.beta[t+1,j] for j in range(M)] for i in range(M)]))
        return kesai
    
    def init_params(self,O,A,B,pi):
        self.T = len(O)
        self.O = O
        if A is not None:
            self.A = A
        if B is not None:
            self.B = B
        if pi is not None:
            self.pi = pi
        
        self.M,self.N = self.B.shape
        
        if A is None:
            self.A = self.A / np.sum(self.A)
            self.B = self.B / np.sum(self.B)
            self.pi = self.pi / np.sum(self.pi)
            

    def cal_probality(self,select=None):
        p = 0
        if select != "backward":
            self.logger.info("前向算法")
        else:
            self.logger.info("后向算法")
        self.alpha = self.forward()
        self.beta = self.backward()
        p = np.sum(self.alpha[-1,:])
        return p

    def baum_welch(self,):
        for e in range(self.epoches):
            
            self.alpha = self.forward()
            self.beta = self.backward()
            
            self.logger.info("第{}次迭代".format(e))
            
            A_ = np.zeros((self.M,self.M))
            B_ = np.zeros((self.M,self.N))
            pi_ = np.zeros(self.M)
                        
            # a_{ij}
            for i in range(self.M):
                for j in range(self.M):
                    molecular = 0.0
                    denominator = 0.0
                    for t in range(self.T-1):
                        molecular += self.calc_kesai(t,i,j)
                        denominator += self.calc_gamma(t,i)
                    A_[i,j] = molecular / denominator
                    
            # b_{jk}
            for j in range(self.M):
                for k in range(self.N):
                    molecular = 0.0
                    denominator = 0.0
                    for t in range(self.T):
                        if k == self.O[t]:
                            molecular += self.calc_gamma(t,j)
                        denominator += self.calc_gamma(t,j)
                    B_[j,k] = molecular / denominator
                    
            # pi_{i}
            for i in range(self.M):
                pi_[i] = self.calc_gamma(0,i)
            
            # 更新
            self.A  = A_/np.sum(A_)
            self.B  = np.array([self.B[i,:]/np.sum(self.B[i,:]) for i in range(self.M)])
            self.pi = pi_/np.sum(pi_)
        self.logger.info("更新完毕:\n A:\n{}\nB:\n{}\npi:\n{}\n".format(self.A.round(5),
                                                                                self.B.round(5),self.pi.round(5)))
        
    def viterbi(self,):
        delta = np.zeros((self.T,self.M))
        delta[0,:] = np.array([self.pi[i]*self.B[i,self.O[0]] for i in range(self.M)])
        psi = np.zeros((self.T,self.M))
        
        for t in range(1,self.T):
            for i in range(self.M):
                max_delta = np.array([delta[t-1,j]*self.A[j,i] for j in range(self.M)])
                delta[t,i] = np.max(max_delta)*self.B[i,self.O[t]]
                psi[t,i] = np.argmax(max_delta)
        print("delta:\n{},\npsi:\n{}".format(delta,psi))
        # 终止
        max_delta = delta[self.T-1,:]
        I = np.zeros(self.T,dtype=int)
        P = np.max(max_delta)
        I[-1] = np.argmax(max_delta)
        print("P(*)-最优路径的概率：{}；最优路径的终点是：{}".format(P.round(5),I[-1]))
        
        
        # 回溯
        for t in range(self.T-2,-1,-1):
            I[t] = psi[t+1,I[t+1]]
        print("最优路径I(*):\n{}".format(I))
        
        
    def fit(self,O,A=None,B=None,pi=None,select=None):
        if A is None:
            select = "bw"
        self.init_params(O,A,B,pi)
        self.logger.info("初始化:\n A:\n{}\nB:\n{}\npi:\n{}\n状态数:{},观测数:{}".format(self.A.round(5),
                                                                                self.B.round(5),self.pi.round(5),self.M,self.N))
        if select == "bw" or select == "baum_welch":
            self.baum_welch()
        self.p = self.cal_probality()
        self.logger.info("当前p(o|lambda)={}".format(self.p.round(5)))
        self.p = self.cal_probality("backward")
        self.logger.info("当前p(o|lambda)={}".format(self.p.round(5)))
        print("p(o|lambda)的概率为:{}".format(self.p))
        
A = np.array([ # 状态转移概率分布
    [.5,.2,.3],
    [.3,.5,.2],
    [.2,.3,.5]
])
B = np.array([ # 观测概率分布
    [.5,.5],
    [.4,.6],
    [.7,.3]
])
pi = np.array([ # 初始概率分布 => 状态的分布
    .2,.4,.4
])

O = [0,1,0] # 1表示红，0表示白 # 观测序列


M,N = A.shape[1],B.shape[0]    
hmm = HMM(10,4)
hmm.fit(O,A,B,pi)
hmm.viterbi()

image-20210105165436434

O = [0,1,0] # 1表示红，0表示白 # 观测序列

M,N = A.shape[1],B.shape[0]    
hmm = HMM(3,2)
hmm.fit(O)
hmm.viterbi()

image-20210105165453785