机器学习(8)-隐马尔可夫模型

这篇文章将会讲解隐马尔可夫模型，我将主要从以下四个方面展开：

模型的定义

隐马尔可夫模型其实是由三个分布确定的模型：即初始概率分布、状态转移概率分布以及观测概率分布，隐马尔可夫模型可以理解为由一个隐藏的马尔科夫链随机生成不可观测的状态随机序列，再由各个状态生成一个观测随机序列的过程，状态序列由状态转移概率分布随机产生，是我们不可知的一个随机序列，观测序列由观测概率分布随机产生，是我们最终得到的观测数据。下面用数学符号来表示这个模型：
我们设Q是所有可能的状态集合，V是所有可能的观测集合：
Q = { q 1 , q 2 , . . . , q N } , V = { v 1 , v 2 , . . . , v M } Q=\{q_1,q_2,...,q_N\},V=\{v_1,v_2,...,v_M\} Q={q1​,q2​,...,qN​},V={v1​,v2​,...,vM​}
N是可能的状态数，M是可能的观测数。假设产生的观测序列长度为T，I是状态序列，O是对应的观测序列：
$$I=(i_1,i_2,...,i_T),O=(o_1,o_2,...,o_T)$$
A是状态转移矩阵：
$$A=[a_{ij}{]}_{N\times N}，其中a_{ij}=P(i_{t+1}=q_j|i_t=q_i)$$
表示从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)$$
表示t时刻处于状态$q_j$的条件下生成观测$v_k$的概率，$\pi$是初始转态概率向量：
$$\pi =({\pi }_i),其中{\pi }_i=P(i_1=q_i)$$
即t=1时刻处于状态$q_i$的概率，因为HMM模型由三个分布确定，因此HMM模型可以表示为：
$$\lambda =(A,B,\pi )$$
除此之外隐马尔可夫模型还有两个基本假设：
（1）齐次马尔可夫性假设，即假设任意时刻t的转态只依赖于前一时刻的状态，与其他任意时刻无关：
$$P(i_t|i_{t-1},o_{t-1},...,i_1,o_1)=P(i_t|i_{t-1})$$
(2)观测独立性假设,即假设任意时刻的观测只依赖于当前时刻的状态，与其他时刻的状态和观测无关：
$$P(o_t|i_T,o_T,i_{T-1},o_{T-1},....,i_{t+1},o_{t+1},i_t,i_{t-1},o_{t-1},...,i_1,o_1)=P(o_t|i_t)$$
下面介绍一下隐马尔可夫模型常见的三个问题：
（1）概率计算问题，给定模型$\lambda =(A,B,\pi )$和观测序列$O=(o_1,o_2,...,o_T)$,计算在模型$\lambda$下观测序列O出现的概率$P(O|\lambda )$
（2）学习问题，已知观测序列$O=(o_1,o_2,...,o_T)$，估计模型参数$\lambda =(A,B,\pi )$，使得在该模型下观测序列概率$P(O|\lambda )$最大
（3）预测问题，已知模型$\lambda =(A,B,\pi )$和观测序列$O=(o_1,o_2,...,o_T)$,求在给定观测序列的条件下，最有可能的状态序列，即求条件概率$P(I|O)$最大的状态序列$I=(i_1,i_2,...,i_T)$
下面就分别介绍如何求解上述三个问题

概率计算问题

给定模型和观测序列，求观测序列的概率主要有三种方法（1）直接计算法（2）前向算法（3）后向算法
直接计算法
思路：按概率公式直接计算，列举出所有可能的长度为T的状态序列，求在不同状态序列下得到观测序列的概率，然后对所有的可能情况求和：
$$P(O|\lambda )=\sum _{I}P(O|I,\lambda )P(I|\lambda )$$
上述方法的时间复杂度为$O(T{N}^T)$，这在计算上显然是不可行的
前向算法
前向算法的本质是利用动态规划的思想，减少了很多重复计算的步骤，时刻t的结果直接从时刻t-1得到，而不需要从头计算，具体步骤如下：
首先定义前向概率${\alpha }_t(i)=P(o_1,o_2,...,o_t,i_t=q_i|\lambda )$表示到时刻t部分的观测序列为$(o_1,o_2,...,o_t)$且时刻t的状态为$q_i$的概率
根据上述定义可以得到，初始时刻t=1的前向概率为：${\alpha }_1(i)={\pi }_i{b}_i(o_1),i=1,2,...,N$,由此可以递推得到任意时刻的前向概率：
$${\alpha }_{t+1}(i)=\left[\sum _{j=1}^N{\alpha }_t(j)a_{ji}\right]b_i(o_{t+1})i=1,2,...,N,t=1,2,...,T-1$$
由此可以得出最终的概率：$P(O|\lambda )={\sum }_{i=1}^N{\alpha }_T(i)$
后向算法
思路和前向算法是一样的，只不过后向算法是从后往前递推，首先定义后向概率${\beta }_t(i)=P(o_{t+1},o_{t+2},...,o_T|i_t=q_i,\lambda )$,表示在时刻t状态为$q_i$的条件下，从t+1到T的部分观测序列为$(o_{t+1},o_{t+2},o_T)$的概率
由定义可以得到初始后向概率为：${\beta }_T(i)=1,i=1,2,...,N$
递推可以得到任意时刻的后向概率为：
$${\beta }_t(i)=\sum _{j=1}^N{a}_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j)i=1,2,...,N,t=T-1,T-2,...,1$$
由此可以得到最终的条件概率$P(O|\lambda )={\sum }_{i=1}^N{\pi }_i{b}_i(o_1){\beta }_1(i)$

学习问题

学习问题根据是否给出状态序列分为监督学习方法和无监督学习方法。
监督学习方法
监督学习方法，是已知T个长度相同的观测序列以及对应的状态序列 { ( O 1 , I 1 ) , ( O 2 , I 2 ) , . . . , ( O T , I T ) } \{(O_1,I_1),(O_2,I_2),...,(O_T,I_T)\} {(O1​,I1​),(O2​,I2​),...,(OT​,IT​)}的情况下，利用极大似然估计法来估计隐马尔可夫模型$\lambda =(A,B,\pi )$的参数
1.转移概率$a_{ij}$的估计
设样本中时刻t状态为i时刻t+1转移到状态j的频数为$A_{ij}$，则状态转移概率$a_{ij}$的估计为：
$${\hat{a}}_{ij}=\frac{A_{ij}}{\sum _{j=1}^N{A}_{ij}}$$
2.观测概率$b_j(k)$的估计
设样本中状态为j并观测为k的频数为$B_{jk}$,则观测概率$b_j(k)$的估计为：
$${\hat{b}}_j(k)=\frac{B_{jk}}{\sum _{k=1}^M{B}_{jk}}$$
3.初始状态概率${\pi }_i$的估计
设样本中初始状态为i的频数为$I_i$，则初始状态概率${\pi }_i$的估计为：
$${\hat{\pi }}_i=\frac{I_i}{\sum _{i=1}^N{I}_i}$$
无监督学习方法
无监督学习是在没有给定状态序列只有观测序列的情况下，估计隐马尔可夫模型的参数，通过Baum-Welch算法对参数进行估计，也就是EM算法。
首先定义两个变量，（1）给定模型$\lambda$和观测O,在时刻t处于状态i的概率
$${\gamma }_t(i)=\frac{{\alpha }_t(i){\beta }_t(i)}{\sum _{j=1}^N{\alpha }_j(t){\beta }_j(t)}$$
(2)给定模型$\lambda$和观测$O$，在时刻t处于状态$q_i$在时刻t+1处于状态$q_j$的概率：
$${\xi }_t(i,j)=\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)}$$
${\alpha }_t(i)，{\beta }_t(i)$是学习问题中的前向概率和后向概率
通过EM算法估计隐马尔可夫模型参数的步骤如下：
（1）初始化，对n=0选取，$a_{ij}^{(0)},b_j(k{)}^{(0)},{\pi }_i^{(0)}$，得到模型${\lambda }^{(0)}=(A^{(0)},B^{(0)},{\pi }^{(0)})$
(2)迭代更新，对$n=0,1,...$
$$\begin{array}{rl}{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)\end{array}$$
等式右边的值都是根据模型参数${\lambda }^{(n)}=(A^{(n)},B^{(n)},{\pi }^{(n)})$计算得到
（3）停止迭代，得到最终参数模型${\lambda }^{(n+1)}=(A^{(n+1)},B^{(n+1)},{\pi }^{(n+1)})$

预测问题

给定模型参数和观测序列预测状态序列两种方法：（1）近似法（2）维特比算法
近似法
近似法的想法是每一个时刻都选择在该时刻最有可能出现的状态$i_t^{\ast }$作为预测状态，最终得到预测状态序列$I^{\ast }=(i_1^{\ast },i_2^{\ast },..,i_T^{\ast })$
在每一个时刻t处于转态$q_i$的概率为${\gamma }_t(i)$,，则
$$i_t^{\ast }=\underset{1\le i\le N}{\mathrm{argmax}}[{\gamma }_t(i)],t=1,2,...,T$$
近似法虽然计算简单，但是通过这样的方法求出来的预测状态序列不能保证在整体上是最有可能的转态序列，因为转态转移概率有为零的部分，会导致预测的状态序列有可能有实际不会发生的部分
维特比算法
维特比算法是通过动态规划的方法来求解最大概率路径，维特比算法基于这样一种思想：假设最优路径在时刻t的转态为$i_t^{\ast }$此时可以将最优路径分为两部分路径，一部分是从$i_1^{\ast }$到$i_t^{\ast }$,另一部分是从$i_t^{\ast }$到$i_T^{\ast }$,那么从$i_1^{\ast }$到$i_t^{\ast }$一定是所有从时刻1到时刻t的状态路径中最优的，如果不是的话我们可以将更优的这条路径和从$i_t^{\ast }$到$i_T^{\ast }$这个部分连接起来形成一条更优的路径，那这和从$i_1^{\ast }$到$i_T^{\ast }$为最优路径这个假设矛盾。
下面我们设${\delta }_t(i)$为时刻t状态为i的所有路径中概率最大值：
$${\delta }_t(i)=\underset{i_1,i_2,...,i_{t-1}}{\max }P(i_t=i,i_{t-1},...,i_1,o_t,...,o_1|\lambda ),i=1,2,...,N$$
由此可以得到$\delta$的转态转移方程为：
$${\delta }_{t+1}(i)=\underset{1\le j\le N}{\max }[{\delta }_t(j)a_{ji}]b_i(o_{t+1})$$
定义${\psi }_t(i)$为时刻t状态为i的所有路径中概率最大对应的路径的第t-1个结点状态：
$${\psi }_t(i)=\underset{1\le j\le N}{\mathrm{argmax}}{\delta }_{t-1}(j)a_{ji},i=1,2,...,N$$
初始条件为：
$${\delta }_1(i)={\pi }_i{b}_i(o_1),i=1,2,...,N$$
$${\psi }_1(i)=0,i=1,2,...,N$$
通过转移方程可以得到最优路径的概率为：
$$P^{\ast }=\underset{1\le j\le N}{\max }{\delta }_T(j)$$
以及最优路径对应的时刻T的状态为：
$$i_T^{\ast }=\underset{1\le j\le N}{\mathrm{argmax}}[{\delta }_T(j)]$$
最后通过最优路径回溯就可以得到最优的状态序列：
$$i_t^{\ast }={\psi }_{t+1}(i_{t+1}^{\ast }),t=T-1,T-2,...,1$$

代码实现

下面的代码实现了如何通过维特比算法对状态序列进行预测

	def predict(self, testfile):
		pi, A, B = self.pi, self.A, self.B

		with open(testfile, 'r', encoding='utf-8') as f:
			lines = f.readlines()

		
		for line in lines:
			line = line.strip()
			print('origin line is {}'.format(line))
			lenl = len(line)
			#维特比算法，进行译码
			delta = np.zeros((lenl, 4))
			psi = np.zeros((lenl, 4)).astype('int')
			for i in range(4):
				delta[0][i] = pi[i] + B[i][ord(line[0])]
			for t in range(1, lenl):
				for i in range(4):
					tmpdelta = list(delta[t-1] + A[:,i] + B[i][ord(line[t])])
					delta[t][i] = max(tmpdelta)
					maxindex = tmpdelta.index(delta[t][i]) #到达当前状态的最优路径
					psi[t][i] = maxindex
			predictlabel = []
			iop_t = np.argmax(delta[lenl-1])
			#状态回溯
			predictlabel.append(iop_t)
			for i in range(lenl-1, 0, -1):
				# print(iop_t)
				iop_t = psi[i][iop_t]
				predictlabel.append(iop_t)
			predictlabel.reverse()
			pre_line = ''
			# print(predictlabel)
			for i in range(lenl):
				pre_line += line[i]
				if (predictlabel[i] == 3 or predictlabel[i] == 2) and i != lenl - 1:
					pre_line += '|'
			print('after participling, line is {}'.format(pre_line))

完整的代码可以在我GitHub项目Machine-Learning-code中进行下载，如有错误或者问题欢迎提出，如果觉得对你有用的话，欢迎star~