文章目录
1.隐马尔可夫模型1.1 简介与定义1.2 观测序列的生成2. 隐马尔可夫模型的3个基本问题2.1 概率计算方法2.1.1 直接计算法2.1.2 前向算法2.1.3 后向算法2.2 学习算法2.2.1 监督学习方法2.2.2 Baum-Welch算法2.2.3 Baum-Welch模型参数估计公式2.3 预测算法2.3.1 近似算法2.3.2 维特比算法3. 隐马尔可夫模型总结参考1.隐马尔可夫模型
1.1 简介与定义
隐马尔可夫模型隐马尔可夫模型是关于时序的概率模型,描述由一个隐藏的马尔可夫链随机生成不同观测的状态随机序列,再由各个状态生成一个观测从而产生观测随机序列的过程。隐藏的马尔可夫链随机生成的状态的序列,称为状态序列;每个状态生成一个观测,而由此产生的观测随机序列,称为观测序列。序列的每一个位置又可以看作是一个时刻。
隐马尔可夫模型由初始概率分布、状态转移概率分布以及观测概率分布确定。即:设QQQ是所有可能的状态的集合,VVV是所有可能的观测的集合:Q={q1,q2,...,qN},V={v1,v2,...,vM}Q=\{q_1,q_2,...,q_N\},\ \ \ V=\{v_1,v_2,...,v_M\}Q={q1,q2,...,qN},V={v1,v2,...,vM}
其中,NNN是可能的状态数,MMM是可能的观测数。III是长度为TTT的状态序列,OOO是对应的观测序列:I=(i1,i2,...,iT),O=(o1,o2,...,oT)I=(i_1,i_2,...,i_T),\ \ \ O=(o_1,o_2,...,o_T)I=(i1,i2,...,iT),O=(o1,o2,...,oT)
AAA是状态转移矩阵:A=[aij]N×NA=\left[a_{ij}\right]_{N×N}A=[aij]N×N
其中:aij=P(it+1=qj∣it=qi),i=1,2,...,N;j=1,2,...,Na_{ij}=P(i_{t+1}=q_j|i_t=q_i),\ \ \ i=1,2,...,N;\ \ \ j=1,2,...,Naij=P(it+1=qj∣it=qi),i=1,2,...,N;j=1,2,...,N
是在时刻ttt处于状态qiq_iqi的条件下再时刻t+1t+1t+1转移到状态qjq_jqj的概率。BBB是观测概率矩阵:B=[bj(k)]N×MB=\left[b_j(k)\right]_{N×M}B=[bj(k)]N×M
其中:bj(k)=P(ot=vk∣it=qj),k=1,2,...,M;j=1,2,...,Nb_j(k)=P(o_t=v_k|i_t=q_j),\ \ \ k=1,2,...,M;\ \ \ j=1,2,...,Nbj(k)=P(ot=vk∣it=qj),k=1,2,...,M;j=1,2,...,N
是在时刻ttt处于状态qjq_jqj的条件下生成观测vkv_kvk的概率。π\piπ是初始状态概率向量:π=(πi)\pi=(\pi_i)π=(πi)
其中:πi=P(i1=qi),i=1,2,...,N\pi_i=P(i_1=q_i),\ \ \ i=1,2,...,Nπi=P(i1=qi),i=1,2,...,N
是时刻t=1t=1t=1处于状态qiq_iqi的概率。隐马尔可夫模型由初始状态概率向量π\piπ、状态转移概率矩阵AAA和观测概率矩阵BBB决定。π\piπ和AAA决定状态序列,BBB决定观测序列。因此,隐马尔可夫模型λ\lambdaλ可由以下表示:λ=(A,B,π)\lambda=(A,B,\pi)λ=(A,B,π)
状态转移矩阵AAA和初始状态概率向量π\piπ确定了隐藏的马尔可夫链,生成不可观测的状态序列。观测概率矩阵BBB确定了如何从状态生成观测,与状态序列综合确定了如何产生观测序列。从定义可知,隐马尔可夫模型作了两个基本的假设:
(1)齐次马尔可夫性假设,即假设隐藏的马尔可夫链再任意时刻ttt的状态只依赖于其前一时刻的状态,与其他时刻的状态及观测无关,也与时刻ttt无关:P(it∣it−1,ot−1,...,i1,o1)=P(it∣it−1),t=1,2,...,TP(i_t|i_{t-1},o_{t-1},...,i_1,o_1)=P(i_t|i_{t-1}),\ \ \ t=1,2,...,TP(it∣it−1,ot−1,...,i1,o1)=P(it∣it−1),t=1,2,...,T
(2)观测独立性假设,即假设任意时刻的观测只依赖于该时刻的马尔可夫链的状态,与其他观测及状态无关:P(ot∣iT,oT,iT−1,oT−1,...,it+1,ot+1,it−1,ot−1,...,i1,o1)=P(ot∣it)P(o_t|i_T,o_T,i_{T-1},o_{T-1},...,i_{t+1},o_{t+1},i_{t-1},o_{t-1},...,i_1,o_1)=P(o_t|i_t)P(ot∣iT,oT,iT−1,oT−1,...,it+1,ot+1,it−1,ot−1,...,i1,o1)=P(ot∣it)
例题假设有444个盒子,每个盒子里都装有红、白两种颜色的球,盒子里的红、白球数如下:
现按照下面的方法抽球,产生一个球的颜色的观测序列:
开始,从444个盒子里以等概率随机选取111个盒子,从这个盒子里随机抽出111个球,记录其颜色后,放回;然后,从当前盒子随机转移到下一个盒子,规则是:如果当前盒子是111,那么下一个盒子一定是盒子222;如果当前盒子是222或333,那么分别以概率0.40.40.4和0.60.60.6转移到左边或右边的盒子;如果当前盒子是444,那么各以0.50.50.5的概率停留在盒子444或转移到盒子333;确定转移的盒子后,再从这个盒子里随机抽出111个球,记录其颜色,放回;如下下去,重复进行555次,得到一个球的颜色的观测序列:O=(红,红,白,白,红)O=(红, 红,白,白,红)O=(红,红,白,白,红)
解
在这个过程中,观察者只能观测到球的颜色的序列,观测不到球是从哪个盒子里取出的,即观测不到盒子的序列。在这个例子中有两个随机序列,一个是盒子的序列,一个是球的颜色的序列。前者是隐藏的,只有后者是可观测的。这是一个隐马尔可夫模型,根据所给条件可以明确状态集合、观测集合、序列长度。盒子对应状态,状态的集合是:Q={盒子1,盒子2,盒子3,盒子4},N=4Q=\{盒子1,盒子2,盒子3,盒子4\},\ \ \ N=4Q={盒子1,盒子2,盒子3,盒子4},N=4
球的颜色对应观测。观测的集合是:V={红,白},M=2V=\{红,白\},\ \ \ M=2V={红,白},M=2
状态序列和观测序列长度T=5T=5T=5。初始概率分布为:π=(0.25,0.25,0.25,0.25)T\pi=(0.25,0.25,0.25,0.25)^{\rm T}π=(0.25,0.25,0.25,0.25)T
状态转移概率分布为:A=[01000.400.6000.400.6000.50.5]A=\begin{gathered} \begin{bmatrix} 0 & 1 & 0 & 0\\ 0.4 & 0 & 0.6 & 0\\ 0 & 0.4 & 0 & 0.6\\ 0 & 0 & 0.5 & 0.5 \end{bmatrix} \end{gathered}A=⎣⎢⎢⎡00.400100.4000.600.5000.60.5⎦⎥⎥⎤
观测概率分布为:B=[0.50.50.30.70.60.40.80.2]B=\begin{gathered} \begin{bmatrix} 0.5 & 0.5\\ 0.3 & 0.7\\ 0.6 & 0.4\\ 0.8 & 0.2 \\ \end{bmatrix} \end{gathered}B=⎣⎢⎢⎡0.50.30.60.80.50.70.40.2⎦⎥⎥⎤
1.2 观测序列的生成
根据隐马尔可夫模型定义,可以将一个长度为TTT的观测序列O=(o1,o2,...,oT)O=(o_1,o_2,...,o_T)O=(o1,o2,...,oT)的生成过程描述如下。
观测序列的生成
输入隐马尔可夫模型λ=(A,B,π)\lambda=(A,B,\pi)λ=(A,B,π),观测序列长度TTT;
输出观测序列O=(o1,o2,...,oT)O=(o_1,o_2,...,o_T)O=(o1,o2,...,oT)。
(1)按照初始状态分布π\piπ产生状态i1i_1i1;
(2)令t=1t=1t=1;
(3)按照状态iti_tit的观测概率分布bit(k)b_{i_t}(k)bit(k)生成oto_tot;
(4)按照状态iti_tit的状态转移概率分布状态{aitit+1}\{a_{i_ti_{t+1}}\}{aitit+1}产生状态it+1i_{t+1}it+1,it+1=1,2,...,Ni_{t+1}=1,2,...,Nit+1=1,2,...,N;
(5)令t=t+1t=t+1t=t+1;如果t<Tt<Tt<T,转向(3);否则,终止。
2. 隐马尔可夫模型的3个基本问题
隐马尔可夫模型有333个基本问题:
(1)概率计算问题给定模型λ=(A,B,π)\lambda=(A,B,\pi)λ=(A,B,π)和观测序列O=(o1,o2,...,oT)O=(o_1,o_2,...,o_T)O=(o1,o2,...,oT),计算在模型λ\lambdaλ下观测序列OOO出现的概率P(O∣λ)P(O|\lambda)P(O∣λ);
(2)学习问题已知观测序列O=(o1,o2,...,oT)O=(o_1,o_2,...,o_T)O=(o1,o2,...,oT),估计模型λ=(A,B,π)\lambda=(A,B,\pi)λ=(A,B,π)参数,使得在该模型下观测序列概率P(O∣λ)P(O|\lambda)P(O∣λ)最大;
(3)预测问题已知模型λ=(A,B,π)\lambda=(A,B,\pi)λ=(A,B,π)和观测序列O=(o1,o2,...,oT)O=(o_1,o_2,...,o_T)O=(o1,o2,...,oT),求对给定观测序列条件概率P(I∣O)P(I|O)P(I∣O)最大的状态序列I=(i1,i2,...,iT)I=(i_1,i_2,...,i_T)I=(i1,i2,...,iT),。即给定观测序列,求最有可能的对应的状态序列。
2.1 概率计算方法
2.1.1 直接计算法
给定模型λ=(A,B,π)\lambda=(A,B,\pi)λ=(A,B,π)和观测序列O=(o1,o2,...,oT)O=(o_1,o_2,...,o_T)O=(o1,o2,...,oT),计算观测序列OOO出现的概率P(O∣λ)P(O|\lambda)P(O∣λ)。最直接的方法是按概率公式直接计算。通过列举所有可能的长度为TTT的状态序列I=(i1,i2,...,iT)I=(i_1,i_2,...,i_T)I=(i1,i2,...,iT),求各个状态序列III与观测序列O=(o1,o2,...,oT)O=(o_1,o_2,...,o_T)O=(o1,o2,...,oT)的联合概率P(O,I∣λ)P(O,I|\lambda)P(O,I∣λ),然后对所有可能的状态序列求和,得到P(O∣λ)P(O|\lambda)P(O∣λ)。状态序列I=(i1,i2,...,iT)I=(i_1,i_2,...,i_T)I=(i1,i2,...,iT)的概率是:P(I∣λ)=πi1ai1i2ai2i3...aiT−1iTP(I|\lambda)=\pi_{i_1}a_{i_1i_2}a_{i_2i_3}...a_{i_{T-1}i_T}P(I∣λ)=πi1ai1i2ai2i3...aiT−1iT
即初始概率乘以每次的转移概率。对固定的状态序列I=(i1,i2,...,iT)I=(i_1,i_2,...,i_T)I=(i1,i2,...,iT),观测序列O=(o1,o2,...,oT)O=(o_1,o_2,...,o_T)O=(o1,o2,...,oT)的概率是:P(O∣I,λ)=bi1(o1)bi2(o2)...biT(oT)P(O|I,\lambda)=b_{i_1}(o_1)b_{i_2}(o_2)...b_{i_T}(o_T)P(O∣I,λ)=bi1(o1)bi2(o2)...biT(oT)
即观测序列中每一项为状态序列中相应元素的概率。则OOO和III同时出现的联合概率为:P(O,I∣λ)=P(O∣I,λ)P(I∣λ)=πi1bi1(o1)ai1i2bi2(o2)ai2i3...aiT−1iTbiT(oT)\begin{aligned} P(O,I|\lambda)&=P(O|I,\lambda)P(I|\lambda)\\&=\pi_{i_1}b_{i_1}(o_1)a_{i_1i_2}b_{i_2}(o_2)a_{i_2i_3}...a_{i_{T-1}i_T}b_{i_T}(o_T) \end{aligned}P(O,I∣λ)=P(O∣I,λ)P(I∣λ)=πi1bi1(o1)ai1i2bi2(o2)ai2i3...aiT−1iTbiT(oT)
然后,对所有可能的状态序列III求和,得到观测序列OOO的概率P(O∣λ)P(O|\lambda)P(O∣λ),即:P(O∣λ)=∑IP(O∣I,λ)P(Iλ)=∑i1,i2,...,iTπi1bi1(o1)ai1i2bi2(o2)ai2i3...aiT−1iTbiT(oT)\begin{aligned} P(O|\lambda)&=\sum_IP(O|I,\lambda)P(I\lambda)\\&=\sum_{i_1,i_2,...,i_T}\pi_{i_1}b_{i_1}(o_1)a_{i_1i_2}b_{i_2}(o_2)a_{i_2i_3}...a_{i_{T-1}i_T}b_{i_T}(o_T) \end{aligned}P(O∣λ)=I∑P(O∣I,λ)P(Iλ)=i1,i2,...,iT∑πi1bi1(o1)ai1i2bi2(o2)ai2i3...aiT−1iTbiT(oT)
但是,上式的计算量很大,是O(TNT)O(TN^{T})O(TNT)阶的,这种算法是不可行的。
2.1.2 前向算法
前向概率给定隐马尔可夫模型λ\lambdaλ,定义到时刻ttt部分观测序列为o1,o2,...,oTo_1,o_2,...,o_To1,o2,...,oT且状态为qiq_iqi的概率为前向概率,记作:αt(i)=P(o1,o2,...ot,it=qi∣λ)\alpha_t(i)=P(o_1,o_2,...o_t,i_t=q_i|\lambda)αt(i)=P(o1,o2,...ot,it=qi∣λ)
可以递推地求得前向概率αt(i)\alpha_t(i)αt(i)及观测序列概率P(O∣λ)P(O|\lambda)P(O∣λ)。
观测序列概率的前向算法
输入隐马尔可夫模型λ\lambdaλ,观测序列OOO;
输出观测序列模型P(O∣λ)P(O|\lambda)P(O∣λ)。
(1)初值:α1(i)=πibi(o1),i=1,2,...,N(1)\alpha_1(i)=\pi_ib_i(o_1),\ \ \ i=1,2,...,N\tag{1}α1(i)=πibi(o1),i=1,2,...,N(1)
(2)递推,对t=1,2,...,T−1t=1,2,...,T-1t=1,2,...,T−1:αt+1(i)=[∑j=1Nαt(i)aji]bi(ot+1),i=1,2,...,N(2)\alpha_{t+1}(i)=\left[\sum_{j=1}^N\alpha_t(i)a_{ji}\right]b_i(o_{t+1}),\ \ \ i=1,2,...,N\tag{2}αt+1(i)=[j=1∑Nαt(i)aji]bi(ot+1),i=1,2,...,N(2)
(3)终止:P(O∣λ)=∑i=1NαT(i)(3)P(O|\lambda)=\sum_{i=1}^N\alpha_T(i)\tag{3}P(O∣λ)=i=1∑NαT(i)(3)
例考虑盒子和球模型λ=(A,B,π)\lambda=(A,B,\pi)λ=(A,B,π),状态集合Q={1,2,3}Q=\{1,2,3\}Q={1,2,3},观测集合V={红,白}V=\{红,白\}V={红,白},A=[0.50.20.30.30.50.20.20.30.5],B=[0.50.50.40.60.70.3],π=[0.20.40.4]\begin{gathered} A=\begin{bmatrix} 0.5 & 0.2 & 0.3\\ 0.3 & 0.5 & 0.2\\ 0.2 & 0.3 & 0.5 \end{bmatrix},\ \ \ B=\begin{bmatrix} 0.5 & 0.5\\ 0.4 & 0.6\\ 0.7 & 0.3 \end{bmatrix},\ \ \ \pi=\begin{bmatrix} 0.2\\ 0.4\\ 0.4 \end{bmatrix} \end{gathered}A=⎣⎡0.50.30.20.20.50.30.30.20.5⎦⎤,B=⎣⎡0.50.40.70.50.60.3⎦⎤,π=⎣⎡0.20.40.4⎦⎤
设T=3T=3T=3,O=(红,白,红)O=(红,白,红)O=(红,白,红),用前向算法计算P(O∣λ)P(O|\lambda)P(O∣λ)。
解
(1)计算初值:α1(1)=π1b1(o1)=0.10\alpha_1(1)=\pi_1b_1(o_1)=0.10α1(1)=π1b1(o1)=0.10
α1(2)=π2b2(o1)=0.16\alpha_1(2)=\pi_2b_2(o_1)=0.16α1(2)=π2b2(o1)=0.16
α1(3)=π3b3(o1)=0.28\alpha_1(3)=\pi_3b_3(o_1)=0.28α1(3)=π3b3(o1)=0.28
(2)递推计算:
α2(1)=[∑i=13α1(i)ai1]b1(o2)=0.07700\alpha_2(1)=\left[\sum_{i=1}^3\alpha_1(i)a_{i1}\right]b_1(o_2)=0.07700α2(1)=[i=1∑3α1(i)ai1]b1(o2)=0.07700
α2(2)=[∑i=13α1(i)ai2]b2(o2)=0.11040\alpha_2(2)=\left[\sum_{i=1}^3\alpha_1(i)a_{i2}\right]b_2(o_2)=0.11040α2(2)=[i=1∑3α1(i)ai2]b2(o2)=0.11040
α2(3)=[∑i=13α1(i)ai3]b3(o2)=0.06060\alpha_2(3)=\left[\sum_{i=1}^3\alpha_1(i)a_{i3}\right]b_3(o_2)=0.06060α2(3)=[i=1∑3α1(i)ai3]b3(o2)=0.06060
α3(1)=[∑i=13α2(i)ai1]b1(o3)=0.04187\alpha_3(1)=\left[\sum_{i=1}^3\alpha_2(i)a_{i1}\right]b_1(o_3)=0.04187α3(1)=[i=1∑3α2(i)ai1]b1(o3)=0.04187
α3(2)=[∑i=13α2(i)ai2]b2(o3)=0.03551\alpha_3(2)=\left[\sum_{i=1}^3\alpha_2(i)a_{i2}\right]b_2(o_3)=0.03551α3(2)=[i=1∑3α2(i)ai2]b2(o3)=0.03551
α3(3)=[∑i=13α2(i)ai3]b3(o3)=0.05284\alpha_3(3)=\left[\sum_{i=1}^3\alpha_2(i)a_{i3}\right]b_3(o_3)=0.05284α3(3)=[i=1∑3α2(i)ai3]b3(o3)=0.05284
(3)终止:P(O∣λ)=∑i=13α3(i)=0.13022P(O|\lambda)=\sum_{i=1}^3\alpha_3(i)=0.13022P(O∣λ)=i=1∑3α3(i)=0.13022
2.1.3 后向算法
后向概率给定隐马尔可夫模型λ\lambdaλ,定义在时刻ttt状态为qiq_iqi的条件下,从t+1t+1t+1到TTT的部分观测序列(ot+1,ot+2,...,oT)(o_{t+1},o_{t+2},...,o_T)(ot+1,ot+2,...,oT)的概率为后向概率,记作:βt(i)=P(ot+1,ot+2,...,oT∣it=qi,λ)\beta_t(i)=P(o_{t+1},o_{t+2},...,o_T|i_t=q_i,\lambda)βt(i)=P(ot+1,ot+2,...,oT∣it=qi,λ)
可以用递推的方法求得后向概率βt(i)\beta_t(i)βt(i)及观测序列概率P(O∣λ)P(O|\lambda)P(O∣λ)。
观测序列概率的后向算法
输入隐马尔可夫模型λ\lambdaλ,观测序列OOO;
输出观测序列模型P(O∣λ)P(O|\lambda)P(O∣λ)。
(1)βT(i)=1,i=1,2,...,N(4)\beta_T(i)=1,\ \ \ i=1,2,...,N\tag{4}βT(i)=1,i=1,2,...,N(4)
(2)对t=T−1,T−2,...,1t=T-1,T-2,...,1t=T−1,T−2,...,1:βt(i)=∑j=1Naijbj(ot+1)βt+1(j),i=1,2,...,N(5)\beta_t(i)=\sum_{j=1}^Na_{ij}b_j(o_{t+1})\beta_{t+1}(j),\ \ \ i=1,2,...,N\tag{5}βt(i)=j=1∑Naijbj(ot+1)βt+1(j),i=1,2,...,N(5)
(3)P(O∣λ)=∑i=1Nπibi(o1)β1(i)(6)P(O|\lambda)=\sum_{i=1}^N\pi_ib_i(o_1)\beta_1(i)\tag{6}P(O∣λ)=i=1∑Nπibi(o1)β1(i)(6)
2.2 学习算法
2.2.1 监督学习方法
假设已给训练数据包含SSS个长度相同的观测序列和对应的状态序列{(O1,I1),(O2,I2),...,(OS,IS)}\{(O_1,I_1),(O_2,I_2),...,(O_S,I_S)\}{(O1,I1),(O2,I2),...,(OS,IS)},那么可以利用极大似然估计法来估计隐马尔可夫模型的参数。
转移概率aija_{ij}aij的估计设样本中时刻ttt处于状态iii时刻t+1t+1t+1转移到状态jjj的频数为AijA_{ij}Aij,那么状态转移概率aija_{ij}aij的估计是:a^ij=Aij∑i=1NAij,i=1,2,...,N;j=1,2,...,N(7)\hat a_{ij}=\frac{A_{ij}}{\sum \limits_{i=1}^N}A_{ij},\ \ \ i=1,2,...,N;\ \ \ j=1,2,...,N\tag{7}a^ij=i=1∑NAijAij,i=1,2,...,N;j=1,2,...,N(7)
观测概率bj(k)b_j(k)bj(k)的估计设样本中状态为jjj并观测为kkk的频数是BjkB_{jk}Bjk,那么状态为jjj观测为kkk的概率bjkb_{jk}bjk的估计是:b^j(k)=Bjk∑k=1MBjk,j=1,2,...,N;K=1,2,...,M(8)\hat b_j(k)=\frac{B_{jk}}{\sum\limits_{k=1}^MB_{jk}},\ \ \ j=1,2,...,N;\ \ \ K=1,2,...,M\tag{8}b^j(k)=k=1∑MBjkBjk,j=1,2,...,N;K=1,2,...,M(8)
初始状态概率πi\pi_iπi的估计π^i\hat\pi_iπ^i为SSS个样本中初始状态为qiq_iqi的概率
2.2.2 Baum-Welch算法
假设给定训练数据只包含SSS个长度为TTT的观测序列{O1,O2,...,OT}\{O_1,O_2,...,O_T\}{O1,O2,...,OT}而没有对应的状态序列,目标是学习马尔可夫模型λ=(A,B,π)\lambda=(A,B,\pi)λ=(A,B,π)的参数。我们将观测序列数据看作观测数据OOO,状态序列数据看作不可观测的隐数据III,那么隐马尔可夫模型事实上是一个含有隐变量的概率模型:P(O∣λ)=∑IP(O∣I,λ)P(I∣λ)(9)P(O|\lambda)=\sum_IP(O|I,\lambda)P(I|\lambda)\tag{9}P(O∣λ)=I∑P(O∣I,λ)P(I∣λ)(9)
它的参数估计可以由EMEMEM算法实现。
确定完全数据的对数似然函数
所有观测数据写成O=(o1,o2,...,oT)O=(o_1,o_2,...,o_T)O=(o1,o2,...,oT),所有隐数据写成I=(i1,i2,...,iT)I=(i_1,i_2,...,i_T)I=(i1,i2,...,iT),完全数据是(O,I)=(o1,o2,...,oT,i1,i2,...,iT)(O,I)=(o_1,o_2,...,o_T,i_1,i_2,...,i_T)(O,I)=(o1,o2,...,oT,i1,i2,...,iT)。完全数据的对数似然函数是logP(O,I∣λ)\log P(O,I|\lambda)logP(O,I∣λ)。
EM算法的E步:求Q函数
Q(λ,λ‾)=∑IlogP(O,I∣λ)P(O,I∣λ‾)Q(\lambda,\overline\lambda)=\sum_I\log P(O,I|\lambda)P(O,I|\overline\lambda)Q(λ,λ)=I∑logP(O,I∣λ)P(O,I∣λ)
其中λ‾\overline\lambdaλ是隐马尔可夫模型参数的当前估计值,λ\lambdaλ是要极大化的隐马尔可夫模型参数。P(O,I∣λ)=πi1bi1(o1)ai1i2bi2(o2)ai2i3...aiT−1iTbiT(oT)P(O,I|\lambda)=\pi_{i_1}b_{i_1}(o_1)a_{i_1i_2}b_{i_2}(o_2)a_{i_2i_3}...a_{i_{T-1}i_T}b_{i_T}(o_T)P(O,I∣λ)=πi1bi1(o1)ai1i2bi2(o2)ai2i3...aiT−1iTbiT(oT)
于是函数Q(λ,λ‾)Q(\lambda,\overline\lambda)Q(λ,λ)可以写成:Q(λ,λ‾)=∑Ilogπi1P(O,I∣λ‾)+∑I(∑t=1T−1logaitit+1)P(O,I∣λ‾)+∑I(∑t=1Tlogbit(ot))P(O,I∣λ‾)Q(\lambda,\overline\lambda)=\sum_I\log\pi_{i_1}P(O,I|\overline\lambda)+\sum_I\left(\sum_{t=1}^{T-1}\log a_{i_ti_{t+1}}\right)P(O,I|\overline\lambda)+\sum_I\left(\sum_{t=1}^T\log b_{i_t}(o_t)\right)P(O,I|\overline\lambda)Q(λ,λ)=I∑logπi1P(O,I∣λ)+I∑(t=1∑T−1logaitit+1)P(O,I∣λ)+I∑(t=1∑Tlogbit(ot))P(O,I∣λ)
式中求和都是对所有数据的序列总长度TTT进行的。
EM算法的M步:极大化Q函数
对上式的三项分别极大化,得到:πi=P(O,i1=i∣λ‾)P(O,λ)\pi_i=\frac{P(O,i_1=i|\overline\lambda)}{P(O,\lambda)}πi=P(O,λ)P(O,i1=i∣λ)
aij=∑t=1T−1P(O,it=i,it+1=j∣λ‾)∑t=1T−1P(O,it=i∣λ‾)a_{ij}=\frac{\sum \limits_{t=1}^{T-1}P(O,i_t=i,i_{t+1}=j|\overline\lambda)}{\sum \limits_{t=1}^{T-1}P(O,i_t=i|\overline\lambda)}aij=t=1∑T−1P(O,it=i∣λ)t=1∑T−1P(O,it=i,it+1=j∣λ)
bj(k)=∑t=1TP(O,it=j∣λ‾)I(ot=vk)∑t=1TP(O,it=j∣‾λ)b_j(k)=\frac{\sum \limits_{t=1}^TP(O,i_t=j|\overline\lambda)I(o_t=v_k)}{\sum \limits_{t=1}^TP(O,i_t=j\overline|\lambda)}bj(k)=t=1∑TP(O,it=j∣λ)t=1∑TP(O,it=j∣λ)I(ot=vk)
2.2.3 Baum-Welch模型参数估计公式
Baum-Welch算法
输入观测数据O=(o1,o2,...,oT)O=(o_1,o_2,...,o_T)O=(o1,o2,...,oT);
输出隐马尔可夫模型参数。
(1)初始化。对n=0n=0n=0,选取aij(0)a_{ij}^{(0)}aij(0),bj(k)(0)b_j(k)^{(0)}bj(k)(0),πi(0)\pi_i^{(0)}πi(0),得到模型λ(0)=(A(0),B(0),π(0))\lambda^{(0)}=(A^{(0)},B^{(0)},\pi^{(0)})λ(0)=(A(0),B(0),π(0));
(2)递推。对n=1,2,...n=1,2,...n=1,2,...:aij=∑t=1T−1P(O,it=i,it+1=j∣λ‾)∑t=1T−1P(O,it=i∣λ‾)a_{ij}=\frac{\sum \limits_{t=1}^{T-1}P(O,i_t=i,i_{t+1}=j|\overline\lambda)}{\sum \limits_{t=1}^{T-1}P(O,i_t=i|\overline\lambda)}aij=t=1∑T−1P(O,it=i∣λ)t=1∑T−1P(O,it=i,it+1=j∣λ)
bj(k)=∑t=1TP(O,it=j∣λ‾)I(ot=vk)∑t=1TP(O,it=j∣‾λ)b_j(k)=\frac{\sum \limits_{t=1}^TP(O,i_t=j|\overline\lambda)I(o_t=v_k)}{\sum \limits_{t=1}^TP(O,i_t=j\overline|\lambda)}bj(k)=t=1∑TP(O,it=j∣λ)t=1∑TP(O,it=j∣λ)I(ot=vk)
πi=P(O,i1=i∣λ‾)P(O,λ)\pi_i=\frac{P(O,i_1=i|\overline\lambda)}{P(O,\lambda)}πi=P(O,λ)P(O,i1=i∣λ)
右端各值按观测O=(o1,o2,...,oT)O=(o_1,o_2,...,o_T)O=(o1,o2,...,oT)和模型λ(n)=(A(n),B(n),π(n))\lambda^{(n)}=(A^{(n)},B^{(n)},\pi^{(n)})λ(n)=(A(n),B(n),π(n))计算。
(3)终止。得到模型参数λ(n+1)=(A(n+1),B(n+1),π(n+1))\lambda^{(n+1)}=(A^{(n+1)},B^{(n+1)},\pi^{(n+1)})λ(n+1)=(A(n+1),B(n+1),π(n+1))。
2.3 预测算法
2.3.1 近似算法
近似算法的想法是,在每个时刻ttt选择在该时刻最有可能出现的状态it∗i_t^*it∗,从而得到一个状态序列I∗=(i1∗,i2∗,...,iT∗)I^*=(i_1^*,i_2^*,...,i_T^*)I∗=(i1∗,i2∗,...,iT∗),将它作为预测的结果。给定隐马尔可夫模型λ\lambdaλ和观测序列OOO,在时刻ttt处于状态qiq_iqi的概率γt(i)\gamma_t(i)γt(i)是:γt(i)=αt(i)βt(i)∑j=1Nαt(i)βt(i)\gamma_t(i)=\frac{\alpha_t(i)\beta_t(i)}{\sum \limits_{j=1}^N\alpha_t(i)\beta_t(i)}γt(i)=j=1∑Nαt(i)βt(i)αt(i)βt(i)
在每一时刻ttt最后可能的状态it∗i_t^*it∗是:it∗=argmax1≤i≤N[γt(i)],t=1,2,...,Ti_t^*=\arg\max_{1\leq i\leq N}[\gamma_t(i)],\ \ \ t=1,2,...,Tit∗=arg1≤i≤Nmax[γt(i)],t=1,2,...,T
从而得到I∗=(i1∗,i2∗,...,iT∗)I^*=(i_1^*,i_2^*,...,i_T^*)I∗=(i1∗,i2∗,...,iT∗)。近似算法的优点是计算简单,其缺点是不能保证预测的状态序列整体是最优可能的状态序列,因为预测的状态序列可能有实际不发生的部分。
2.3.2 维特比算法
维特比算法实际使用动态规划解隐马尔可夫模型预测问题,即用动态规划求概率最大路径。这时的路径对应着一个状态。首先导入两个变量δ\deltaδ和Ψ\PsiΨ。定义在时刻ttt状态为iii的所有单个路径(i1,i2,...,iT)(i_1,i_2,...,i_T)(i1,i2,...,iT)中概率最大值为:δt(i)=maxi1,i2,...,it−1P(it=i,it−1,...,i1,ot,...,o1∣λ),i=1,2,...,N\delta_t(i)=\max_{i_1,i_2,...,i_{t-1}}P(i_t=i,i_{t-1},...,i_1,o_t,...,o_1|\lambda),\ \ \ i=1,2,...,Nδt(i)=i1,i2,...,it−1maxP(it=i,it−1,...,i1,ot,...,o1∣λ),i=1,2,...,N
由定义得到变量δ\deltaδ的递推公式:δt+1(i)=maxi1,i2,...,it−1P(it=i,it−1,...,i1,ot,...,o1∣λ)=max1≤j≤N[δt(j)aji]bi(ot+1),i=1,2,...,N;t=1,2,...,T−1\begin{aligned} \delta_{t+1}(i)&=\max_{i_1,i_2,...,i_{t-1}}P(i_t=i,i_{t-1},...,i_1,o_t,...,o_1|\lambda)\\&=\max_{1\leq j\leq N}[\delta_t(j)a_{ji}]b_i(o_{t+1}),\ \ \ i=1,2,...,N;\ \ \ t=1,2,...,T-1 \end{aligned}δt+1(i)=i1,i2,...,it−1maxP(it=i,it−1,...,i1,ot,...,o1∣λ)=1≤j≤Nmax[δt(j)aji]bi(ot+1),i=1,2,...,N;t=1,2,...,T−1
定义在时刻ttt状态为iii的所有单个路径(i1,i2,...,it−1,i)(i_1,i_2,...,i_{t-1},i)(i1,i2,...,it−1,i)中概率最大的路径的第t−1t-1t−1个结点为:Ψ(i)=argmax1≤j≤N[δt−1(j)aji],i=1,2,...,N\Psi(i)=\arg\max_{1\leq j\leq N}[\delta_{t-1}(j)a_{ji}],\ \ \ i=1,2,...,NΨ(i)=arg1≤j≤Nmax[δt−1(j)aji],i=1,2,...,N
维特比算法
输入模型λ=(A,B,π)\lambda=(A,B,\pi)λ=(A,B,π)和观测O=(o1,o2,...,oT)O=(o_1,o_2,...,o_T)O=(o1,o2,...,oT);
输出最优路径I∗=(i1∗,i2∗,...,iT∗)I^*=(i_1^*,i_2^*,...,i_T^*)I∗=(i1∗,i2∗,...,iT∗)。
(1)初始化:δ1(i)=πibi(o1),i=1,2,...,N\delta_1(i)=\pi_ib_i(o_1),\ \ \ i=1,2,...,Nδ1(i)=πibi(o1),i=1,2,...,N
Ψ1(i)=0,i=1,2,...,N\Psi_1(i)=0,\ \ \ i=1,2,...,NΨ1(i)=0,i=1,2,...,N
(2)递推。对t=2,3,...,Tt=2,3,...,Tt=2,3,...,T:δt(i)=max1≤j≤N[δt−1(j)aji]bi(ot+1),i=1,2,...,N\delta_t(i)=\max_{1\leq j\leq N}[\delta_{t-1}(j)a_{ji}]b_i(o_{t+1}),\ \ \ i=1,2,...,Nδt(i)=1≤j≤Nmax[δt−1(j)aji]bi(ot+1),i=1,2,...,N
Ψ(i)=argmax1≤j≤N[δt−1(j)aji],i=1,2,...,N\Psi(i)=\arg\max_{1\leq j\leq N}[\delta_{t-1}(j)a_{ji}],\ \ \ i=1,2,...,NΨ(i)=arg1≤j≤Nmax[δt−1(j)aji],i=1,2,...,N
(3)终止:P∗=max1≤i≤NδT(i)P^*=\max_{1\leq i\leq N}\delta_T(i)P∗=1≤i≤NmaxδT(i)
iT∗=argmax1≤i≤N[δT(i)]i^*_T=\arg\max_{1\leq i\leq N}[\delta_T(i)]iT∗=arg1≤i≤Nmax[δT(i)]
(3)最优路径回溯。对t=T−1,T−2,...,1t=T-1,T-2,...,1t=T−1,T−2,...,1:it∗=Ψt+1(it+1∗)i_t^*=\Psi_{t+1}(i_{t+1}^*)it∗=Ψt+1(it+1∗)
求得最优路径I∗=(i1∗,i2∗,...,iT∗)I^*=(i_1^*,i_2^*,...,i_T^*)I∗=(i1∗,i2∗,...,iT∗)。
解例题
(1)初始化。在t=1t=1t=1时,对每个状态iii,i=1,2,3i=1,2,3i=1,2,3,求状态为iii观测o1o_1o1为红的概率,记此概率为δ1(i)\delta_1(i)δ1(i),则:δ1(i)=πibi(o1)=πibi(红),i=1,2,3\delta_1(i)=\pi_ib_i(o_1)=\pi_ib_i(红),\ \ \ i=1,2,3δ1(i)=πibi(o1)=πibi(红),i=1,2,3
带入实际数据:δ1(i)=0.10,δ2(i)=0.16,δ3(i)=0.28\delta_1(i)=0.10,\ \ \ \delta_2(i)=0.16,\ \ \ \delta_3(i)=0.28δ1(i)=0.10,δ2(i)=0.16,δ3(i)=0.28
记Ψ1(i)=0\Psi_1(i)=0Ψ1(i)=0,i=1,2,3i=1,2,3i=1,2,3。
(2)在t=2t=2t=2时,对每个状态iii,i=1,2,3i=1,2,3i=1,2,3,求在t=1t=1t=1时状态为jjj观测为红并在t=2t=2t=2时状态为iii观测o2o_2o2为白的路径的最大概率,记次最大概率为δ2(i)\delta_2(i)δ2(i),则:δ2(i)=max1≤j≤3[δ1(j)aji]bi(o2)\delta_2(i)=\max_{1\leq j\leq 3}[\delta_1(j)a_{ji}]b_i(o_2)δ2(i)=1≤j≤3max[δ1(j)aji]bi(o2)
同时,对每个状态iii,i=1,2,3i=1,2,3i=1,2,3,记录概率最大路径的前一个状态jjj:Ψ2(i)=argmax1≤j≤3[δ1(j)aji],i=1,2,3\Psi_2(i)=\arg\max_{1\leq j\leq3}[\delta_1(j)a_{ji}],\ \ \ i=1,2,3Ψ2(i)=arg1≤j≤3max[δ1(j)aji],i=1,2,3
计算:δ2(1)=max1≤j≤3[δ1(j)aj1]b1(o2)=maxj{0.10×0.5,0.16×0.3,0.28×0.2}×0.5=0.028\begin{aligned} \delta_2(1)&=\max_{1\leq j\leq3}[\delta_1(j)a_{j1}]b_1(o_2)\\&=\max_j\{0.10×0.5,0.16×0.3,0.28×0.2\}×0.5\\&=0.028 \end{aligned}δ2(1)=1≤j≤3max[δ1(j)aj1]b1(o2)=jmax{0.10×0.5,0.16×0.3,0.28×0.2}×0.5=0.028
Ψ2(1)=argmax1≤j≤3[δ1(j)aj1]=argmaxj{0.10×0.5,0.16×0.3,0.28×0.2}=3\begin{aligned} \Psi_2(1)&=\arg\max_{1\leq j\leq3}[\delta_1(j)a_{j1}]\\&=\arg\max_j\{0.10×0.5,0.16×0.3,0.28×0.2\}\\&=3 \end{aligned}Ψ2(1)=arg1≤j≤3max[δ1(j)aj1]=argjmax{0.10×0.5,0.16×0.3,0.28×0.2}=3
δ2(2)=max1≤j≤3[δ1(j)aj2]b2(o2)=maxj{0.10×0.2,0.16×0.5,0.28×0.3}×0.6=0.0504\begin{aligned} \delta_2(2)&=\max_{1\leq j\leq3}[\delta_1(j)a_{j2}]b_2(o_2)\\&=\max_j\{0.10×0.2,0.16×0.5,0.28×0.3\}×0.6\\&=0.0504 \end{aligned}δ2(2)=1≤j≤3max[δ1(j)aj2]b2(o2)=jmax{0.10×0.2,0.16×0.5,0.28×0.3}×0.6=0.0504
Ψ2(2)=argmax1≤j≤3[δ1(j)aj2]=argmaxj{0.10×0.2,0.16×0.5,0.28×0.3}=3\begin{aligned} \Psi_2(2)&=\arg\max_{1\leq j\leq3}[\delta_1(j)a_{j2}]\\&=\arg\max_j\{0.10×0.2,0.16×0.5,0.28×0.3\}\\&=3 \end{aligned}Ψ2(2)=arg1≤j≤3max[δ1(j)aj2]=argjmax{0.10×0.2,0.16×0.5,0.28×0.3}=3
δ2(3)=max1≤j≤3[δ1(j)aj3]b3(o2)=maxj{0.10×0.3,0.16×0.2,0.28×0.5}×0.3=0.042\begin{aligned} \delta_2(3)&=\max_{1\leq j\leq3}[\delta_1(j)a_{j3}]b_3(o_2)\\&=\max_j\{0.10×0.3,0.16×0.2,0.28×0.5\}×0.3\\&=0.042 \end{aligned}δ2(3)=1≤j≤3max[δ1(j)aj3]b3(o2)=jmax{0.10×0.3,0.16×0.2,0.28×0.5}×0.3=0.042
Ψ2(3)=argmax1≤j≤3[δ1(j)aj3]=argmaxj{0.10×0.3,0.16×0.2,0.28×0.5}=3\begin{aligned} \Psi_2(3)&=\arg\max_{1\leq j\leq3}[\delta_1(j)a_{j3}]\\&=\arg\max_j\{0.10×0.3,0.16×0.2,0.28×0.5\}\\&=3 \end{aligned}Ψ2(3)=arg1≤j≤3max[δ1(j)aj3]=argjmax{0.10×0.3,0.16×0.2,0.28×0.5}=3
同样,在t=3t=3t=3时,δ3(1)=0.00756,Ψ3(1)=2\delta_3(1)=0.00756,\ \ \ \Psi_3(1)=2δ3(1)=0.00756,Ψ3(1)=2
δ3(2)=0.01008,Ψ3(2)=2\delta_3(2)=0.01008,\ \ \ \Psi_3(2)=2δ3(2)=0.01008,Ψ3(2)=2
δ3(3)=0.01470,Ψ3(3)=3\delta_3(3)=0.01470,\ \ \ \Psi_3(3)=3δ3(3)=0.01470,Ψ3(3)=3
(3)以P∗P^*P∗表示最优路径的概率,则:P∗=max1≤i≤3δ3(i)=0.0147P^*=\max_{1\leq i\leq3}\delta_3(i)=0.0147P∗=1≤i≤3maxδ3(i)=0.0147
最优路径的终点是i3∗i_3^*i3∗:i3∗argmaxi[δ3(i)]=3i_3^*\arg\max_i[\delta_3(i)]=3i3∗argimax[δ3(i)]=3
(4)由最优路径的终点i3∗i_3^*i3∗,逆向找到i1∗i_1^*i1∗和i2∗i_2^*i2∗:在t=2时,i2∗=Ψ3(i3∗)=Ψ3(3)=3在t=2时,\ \ \ \ i_2^*=\Psi_3(i_3^*)=\Psi_3(3)=3在t=2时,i2∗=Ψ3(i3∗)=Ψ3(3)=3
在t=1时,i1∗=Ψ2(i2∗)=Ψ2(3)=3在t=1时,\ \ \ \ i_1^*=\Psi_2(i_2^*)=\Psi_2(3)=3在t=1时,i1∗=Ψ2(i2∗)=Ψ2(3)=3
于是求得最优路径,即最优状态序列I∗=(i1∗,i2∗,i3∗)=(3,3,3)I^*=(i_1^*,i_2^*,i_3^*)=(3,3,3)I∗=(i1∗,i2∗,i3∗)=(3,3,3)。
3. 隐马尔可夫模型总结
隐马尔可夫模型是关于时序的概率模型,描述一个隐藏的马尔可夫链随机生成不可观测的状态的序列,再由各个状态随机生成一个观测从而产生观测序列的过程。隐马尔可夫模型由初始状态向量π\piπ、状态转移概率矩阵AAA和观测概率矩阵BBB决定。因此,隐马尔可夫模型可以写成λ=(A,B,π)\lambda=(A,B,\pi)λ=(A,B,π)。