马尔可夫链与隐马尔可夫模型
马尔可夫链状态转移概率有限状态机隐马尔可夫模型HMM概率计算问题:前向算法后向算法HMM常用概率的计算序列问题维特比算法流程学习问题(训练算法)监督学习无监督学习EM算法简介Baum-Welch算法过程:马尔可夫链
考察一个随机变量序列,这些随机变量并不是相互独立的,每个随机变量的值依赖于这个序列前面的状态。如果一个系统有NNN个有限状态S=(s1,s2,s3,...sN)S=(s_1,s_2,s_3,...s_N)S=(s1,s2,s3,...sN),那么随着时间推移,该系统将从某一状态转移到另一状态。设Q=(q1,q2,q3,...,qT)Q=(q_1,q_2,q_3,...,q_T)Q=(q1,q2,q3,...,qT)为一个随机变量序列,取值为状态集SSS中的某个状态,对该系统的描述通常需要给出当前时刻ttt的状态和其前面所有状态的关系:系统在时刻ttt处于状态sjs_jsj的概率为:P(qt=sj∣qt−1=si,qt−2=sk,...)P(q_t=s_j|q_{t-1}=s_i,q_{t-2}=s_k,...)P(qt=sj∣qt−1=si,qt−2=sk,...)
如果特定条件下,系统在时间ttt的状态只与其在时间t−1t-1t−1的状态相关,即:
P(qt=sj∣qt−1=si,qt−2=sk,...)=P(qt=sj∣qt−1=si)P(q_t=s_j|q_{t-1}=s_i,q_{t-2}=s_k,...)=P(q_t=s_j|q_{t-1}=s_i)P(qt=sj∣qt−1=si,qt−2=sk,...)=P(qt=sj∣qt−1=si)
则该系统构成一个离散的一阶马尔可夫链。
状态转移概率
考虑独立于时间ttt的一阶马尔可夫过程:
P(qt=sj∣qt−1=si)=aij,1≤i,j≤NP(q_t=s_j|q_{t-1}=s_i)=a_{ij}, 1\le i,j\le NP(qt=sj∣qt−1=si)=aij,1≤i,j≤N
满足:aij≥0,∑j=1Naij=1a_{ij} \ge 0,\sum_{j=1}^Na_{ij}=1aij≥0,∑j=1Naij=1
有NNN个状态的一阶马尔可夫过程有N2N^2N2次状态转移,其状态转移概率可以表示成一个状态转移矩阵A=[aij]A=[a_{ij}]A=[aij]。
有限状态机
如图,马尔可夫模型又可视为随机的有限状态机,圆圈表示状态,状态之间的转移用带箭头的弧表示,弧上的数字为状态转移的概率。初始状态用标记为start的输入箭头表示,任何一个状态都可作为终止状态。
一个马尔可夫链的状态序列的概率可以通过计算形成该状态序列的所有状态之间转移弧上的概率乘积而得出,即:
P(q1,q2,...,qT)=P(q1)P(q2∣q1)P(q3∣q1,q2)...P(qT∣q1,q2,...,qT−1)=P(q1)P(q2∣q1)P(q3∣q2)...P(qT∣qT−1)=πq1∏t=1T−1aqtqt+1\begin{aligned} P(q_1,q_2,...,q_T) =& P(q_1)P(q_2|q_1)P(q_3|q_1,q_2)...P(q_T|q_1,q_2,...,q_{T-1})\\ =& P(q_1)P(q_2|q_1)P(q_3|q_2)...P(q_T|q_{T-1})\\ =& \pi_{q_1}\prod_{t=1}^{T-1}a_{q_tq_{t+1}} \end{aligned} P(q1,q2,...,qT)===P(q1)P(q2∣q1)P(q3∣q1,q2)...P(qT∣q1,q2,...,qT−1)P(q1)P(q2∣q1)P(q3∣q2)...P(qT∣qT−1)πq1t=1∏T−1aqtqt+1
其中,πq1=P(q1)\pi_{q_1}=P(q_1)πq1=P(q1)。
隐马尔可夫模型HMM
马尔可夫模型中,每一个状态代表了一个可观察的事件,所以又称作可视马尔可夫模型。对于有些问题,我们并不知道所经过的状态序列,只知道状态序列的概率函数(观测序列,可见的)。也就是说模型的状态转换过程是不可观察的,即隐蔽的。这种情况下,需要使用隐马尔可夫模型进行建模,它可以看作是马尔可夫链的一个扩展。
详细介绍推荐/pinard/p/6945257.html
一般的,一个HMM记为一个五元组μ=(S,K,A,B,π)\mu=(S,K,A,B,\pi)μ=(S,K,A,B,π),SSS为状态的集合,KKK为输出符号的集合,A,B,πA,B,\piA,B,π分别表示状态转移概率、符号发射概率(也叫生成概率,来自观测独立性假设,即任意时刻的观察状态只仅仅依赖于当前时刻的隐藏状态,bj(k)=P(ot=vk∣it=qj)b_j(k)=P(o_t=v_k|i_t=q_j)bj(k)=P(ot=vk∣it=qj))和初始状态的概率分布。为了简单,有时也将五元组记为三元组λ=(A,B,π)\lambda=(A,B,\pi)λ=(A,B,π)。
HMM中有三个基本问题:
概率计算问题:
给定HMM模型三元组λ\lambdaλ和观察序列O=(o1,o2,o3...)O=(o_1,o_2,o_3...)O=(o1,o2,o3...),计算观察序列的概率P(O∣λ)P(O|\lambda)P(O∣λ)。
如果用穷举所有可能的隐藏状态序列I=(q1,q2,...,qT)I=(q_1,q_2,...,q_T)I=(q1,q2,...,qT),分别求观测序列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∣λ)。但该算法的时间复杂度将是O(TNT)O(TN^T)O(TNT),因此利用下列简洁算法——前向后向算法
前向算法
前向和后向算法本质都属于动态规划算法。
定义时刻ttt时隐藏状态为qtq_tqt, 观测状态的序列为o1,o2,...,oto_1,o_2,...,o_to1,o2,...,ot的概率为前向概率,即:
α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∣λ)
ajia_{ji}aji为状态转移概率,bi(ot+1)b_i(o_{t+1})bi(ot+1)为生成概率,则有概率关系的递推公式:
αt+1(i)=[∑j=1Nαt(j)aji]bi(ot+1)\alpha_{t+1}(i)=[\sum_{j=1}^N \alpha_t(j)a_{ji}]b_i(o_{t+1}) αt+1(i)=[j=1∑Nαt(j)aji]bi(ot+1)
总结前向算法:
输入:HMM模型λ=(A,B,π)\lambda=(A,B,\pi)λ=(A,B,π),观测序列O=(o1,o2,...,oT)O=(o_1,o_2,...,o_T)O=(o1,o2,...,oT)
输出:观测序列概率P(O∣λ)P(O|\lambda)P(O∣λ)
计算时刻1各隐藏状态的前向概率:
α1(i)=πibi(o1),i=1,2,...N\alpha_1(i)=\pi_ib_i(o_1),i=1,2,...N α1(i)=πibi(o1),i=1,2,...N递推时刻2,3,…T的前向概率:
αt+1(i)=[∑j=1Nαt(j)aji]bi(ot+1),i=1,2,...,N\alpha_{t+1}(i)=[\sum_{j=1}^N \alpha_t(j)a_{ji}]b_i(o_{t+1}), i=1,2,...,N αt+1(i)=[j=1∑Nαt(j)aji]bi(ot+1),i=1,2,...,N计算最终结果:
P(O∣λ)=∑i=1NαT(i)P(O|\lambda)=\sum_{i=1}^N \alpha_T(i) P(O∣λ)=i=1∑NαT(i)
可以看到算法复杂度是O(TN2)O(TN^2)O(TN2)
后向算法
相对于前向概率,定义:时刻ttt时刻隐藏状态为qiq_iqi的条件下, 从时刻t+1t+1t+1到最后时刻TTT的观测序列为ot+1,ot+2,...oTo_{t+1},o_{t+2},...o_Tot+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)=∑j=1Naijbj(ot+1)βt+1(j)\beta_t(i)=\sum_{j=1}^N a_{ij}b_j(o_{t+1})\beta_{t+1}(j) βt(i)=j=1∑Naijbj(ot+1)βt+1(j)
推导过程:在时间ttt状态为qiq_iqi的条件下,HMM输出观察序列的过程可以分解为如下两个步骤:
从时间ttt到时间t+1t+1t+1,隐藏状态变为qjq_jqj,并输出观察值Ot+1O_{t+1}Ot+1,概率为aijbj(ot+1)a_{ij}b_j(o_{t+1})aijbj(ot+1)在时间t+1t+1t+1,状态为qjq_jqj的条件下,输出观察序列ot+2,...oTo_{t+2},...o_Tot+2,...oT的概率,根据定义就是后向概率βt+1(j)\beta_{t+1}(j)βt+1(j)
所以得到上式的递推关系。
总结后向算法:
初始化时刻TTT各个隐藏状态的后向概率:
βT(i)=1,i=1,2,...,N\beta_T(i)=1,i=1,2,...,N βT(i)=1,i=1,2,...,N
2) 递推时刻T−1T-1T−1、T−2T-2T−2、T−3T-3T−3…1时刻的后向概率:
βt(i)=∑j=1Naijbj(ot+1)βt+1(j),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(i)=j=1∑Naijbj(ot+1)βt+1(j),i=1,2,...,N
3) 计算最终结果:
P(O∣λ)=∑i=1Nπibi(o1)β1(i)P(O|\lambda)=\sum_{i=1}^N\pi_ib_i(o_1)\beta_1(i) P(O∣λ)=i=1∑Nπibi(o1)β1(i)
算法复杂度依然为O(TN2)O(TN^2)O(TN2)
HMM常用概率的计算
根据前向和后向概率的定义有:
P(O,it=qi∣λ)=P(o1,...,ot,it=qi,ot+1,...,oT∣λ)=P(o1,...,ot,it=qi∣λ)P(ot+1,...,oT∣o1,...ot,it=qi,λ)=P(o1,...,ot,it=qi∣λ)P(ot+1,...,oT∣it=qi,λ)=αt(i)βt(i)\begin{aligned} P(O,i_t=q_i|\lambda) =& P(o_1,...,o_t,i_t=q_i,o_{t+1},...,o_T|\lambda)\\ =& P(o_1,...,o_t,i_t=q_i|\lambda)P(o_{t+1},...,o_T|o_1,...o_t,i_t=q_i,\lambda)\\ =& P(o_1,...,o_t,i_t=q_i|\lambda)P(o_{t+1},...,o_T|i_t=q_i,\lambda)\\ =&\alpha_t(i)\beta_t(i) \end{aligned} P(O,it=qi∣λ)====P(o1,...,ot,it=qi,ot+1,...,oT∣λ)P(o1,...,ot,it=qi∣λ)P(ot+1,...,oT∣o1,...ot,it=qi,λ)P(o1,...,ot,it=qi∣λ)P(ot+1,...,oT∣it=qi,λ)αt(i)βt(i)
P(O∣λ)=∑i=1Nαt(i)βt(i),1≤t≤TP(O|\lambda)=\sum_{i=1}^N\alpha_t(i)\beta_t(i), 1\le t \le T P(O∣λ)=i=1∑Nαt(i)βt(i),1≤t≤T
给定模型λ\lambdaλ和观测序列OOO,在时刻ttt处于状态qiq_iqi的概率记为:
γt(i)=P(it=qi∣O,λ)=P(it=qi,O∣λ)P(O∣λ)=αt(i)βt(i)∑j=1Nαt(j)βt(j)\begin{aligned} \gamma_t(i)=& P(i_t=q_i|O,\lambda)\\ =& {P(i_t=q_i,O|\lambda) \over P(O|\lambda)}\\ =& {\alpha_t(i) \beta_t(i) \over \sum_{j=1}^N\alpha_t(j)\beta_t(j)} \end{aligned} γt(i)===P(it=qi∣O,λ)P(O∣λ)P(it=qi,O∣λ)∑j=1Nαt(j)βt(j)αt(i)βt(i)给定模型λ\lambdaλ和观测序列OOO,在时刻ttt处于状态qiq_iqi,且时刻t+1t+1t+1处于状态qjq_jqj的概率记为:
ξt(i,j)=P(it=qi,it+1=qj∣O,λ)=Pit=qi,it+1=qj,O∣λ)P(O∣λ)=αt(i)aijbj(ot+1)βt+1(j)∑r=1N∑s=1Nαt(r)arsbs(ot+1)βt+1(s)\begin{aligned} \xi_t(i,j)=& P(i_t=q_i,i_{t+1}=q_j|O,\lambda)\\ =& {Pi_t=q_i,i_{t+1}=q_j,O|\lambda) \over P(O|\lambda)}\\ =& \alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j) \over \sum_{r=1}^N \sum_{s=1}^N \alpha_t(r)a_{rs}b_s(o_{t+1})\beta_{t+1}(s) \end{aligned} ξt(i,j)===P(it=qi,it+1=qj∣O,λ)P(O∣λ)Pit=qi,it+1=qj,O∣λ)∑r=1N∑s=1Nαt(r)arsbs(ot+1)βt+1(s)αt(i)aijbj(ot+1)βt+1(j)将γt(i)\gamma_t(i)γt(i)和ξt(i,j)\xi_t(i,j)ξt(i,j)在各个时刻ttt求和,可以得到:
在观测序列OOO下状态iii出现的(概率)期望值:
∑t=1Tγt(i)\sum_{t=1}^T\gamma_t(i) t=1∑Tγt(i)
在观测序列OOO下状态iii转移的期望值:
∑t=1T−1γt(i)\sum_{t=1}^{T-1}\gamma_t(i) t=1∑T−1γt(i)
在观测序列OOO下状态iii转移到状态j的期望值:
∑t=1T−1ξt(i,j)\sum_{t=1}^{T-1}\xi_t(i,j) t=1∑T−1ξt(i,j)
序列问题
序列问题提出:给定模型λ=(A,B,Π)λ=(A,B,Π)λ=(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∗},即P(I∗∣O)P(I^∗|O)P(I∗∣O)要最大化。可以利用维特比算法进行问题的求解。维特比算法是一个通用的解码算法,是基于动态规划的求序列最短路径的方法。维特比算法定义了两个变量,用于描述局部状态:
1.定义维特比变量δt(i)\delta_t(i)δt(i),是在时间ttt时,HMM沿着某一条路径到达状态iii,并输出观察序列o1,o2,o3,...,oto_1,o_2,o_3,...,o_to1,o2,o3,...,ot的最大概率:
δt(i)=maxi1,i2,...,itP(it=i,i1,i2,...,it−1,ot,ot−1,...,o1∣λ),i=1,2,...,N\delta_t(i)=max_{i_1,i_2,...,i_t}P(i_t=i,i_1,i_2,...,i_{t-1},o_t,o_{t-1},...,o_1|\lambda),i=1,2,...,N δt(i)=maxi1,i2,...,itP(it=i,i1,i2,...,it−1,ot,ot−1,...,o1∣λ),i=1,2,...,N
2.定义变量Ψt(i)\Psi_t(i)Ψt(i)用于路径记忆,记录某一条概率最大路径上状态it=ii_t=iit=i的前一个状态:
Ψt(i)=argmax1≤j≤N[δt−1(j)aji]\Psi_t(i)=argmax_{1\le j\le N}[\delta_{t-1}(j)a_{ji}] Ψt(i)=argmax1≤j≤N[δt−1(j)aji]
通过这两个局部状态,我们就可以从时刻000一直递推到时刻TTT,然后利用Ψt(i)Ψ_t(i)Ψt(i)记录的前一个最可能的状态节点回溯,直到找到最优的隐藏状态序列。
维特比算法流程
算法流程:
初始化局部状态:
δ1(i)=πibi(o1),i=1,2,...,NΨ1(i)=0,i=1,2,...,N\delta_1(i)=\pi_ib_i(o_1),i=1,2,...,N\\ \Psi_1(i)=0,i=1,2,...,N δ1(i)=πibi(o1),i=1,2,...,NΨ1(i)=0,i=1,2,...,N进行动态规划,递推时刻t=2,3,...,Tt=2,3,...,Tt=2,3,...,T的局部状态:
δt(i)=maxj[δt−1(j)aji]bi(ot),i=1,2,...,NΨt(i)=argmaxj[δt−1(j)aji],i=1,2,...,N\delta_t(i)=max_j[\delta_{t-1}(j)a_{ji}]b_i(o_{t}),i=1,2,...,N\\ \Psi_t(i)=argmax_j[\delta_{t-1}(j)a_{ji}],i=1,2,...,N δt(i)=maxj[δt−1(j)aji]bi(ot),i=1,2,...,NΨt(i)=argmaxj[δt−1(j)aji],i=1,2,...,N计算时刻TTT最大的δT(i)\delta_T(i)δT(i),即为最可能隐藏状态序列出现的概率。计算时刻TTT最大的Ψt(i)Ψ_t(i)Ψt(i),即为时刻TTT最可能的隐藏状态:
P∗=max1≤j≤NδT(j)iT∗=argmax1≤j≤N[δT(j)]P*=max_{1 \le j \le N}\delta_T(j)\\ i_T^*=argmax_{1 \le j \le N}[\delta_T(j)] P∗=max1≤j≤NδT(j)iT∗=argmax1≤j≤N[δT(j)]利用局部状态Ψ(i)Ψ(i)Ψ(i)开始回溯。对于t=T−1,T−2,...,1t=T−1,T−2,...,1t=T−1,T−2,...,1:
it∗=Ψt+1(it+1∗)i_t^*=Ψ_{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∗}。
学习问题(训练算法)
已知观测序列,估计模型λ\lambdaλ参数,使得在该模型下观测序列概率P(O∣λ)P(O|\lambda)P(O∣λ)最大。
隐马尔可夫模型的学习,根据训练数据包括观测序列和对应的状态序列,还是只有观测序列,可以分别由监督学习和非监督学习算法(Baum-Welch算法,也就是EM算法)实现。
监督学习
已知SSS个长度为同为TTT的观测序列和对应的隐藏状态序列,即{(O1,I1),(O2,I2),...(OS,IS)}\{(O_1,I_1),(O_2,I_2),...(O_S,I_S)\}{(O1,I1),(O2,I2),...(OS,IS)}是已知的,此时我们使用最大似然估计(MLE)方法来求解模型参数。
转移概率aija_{ij}aij的估计
假设样本中ttt时刻处于状态 qiq_iqi,t+1t+1t+1转移到状态qjq_jqj的频数是AijA_{ij}Aij,那么状态转移概率估计为:
a^ij=Aij∑j=1NAij\hat{a}_{ij}={A_{ij} \over \sum_{j=1}^NA_{ij}} a^ij=∑j=1NAijAij生成概率bj(k)b_j(k)bj(k)的估计
假设样本隐藏状态为qjq_jqj且观测值为vkv_kvk的频数是BjkB_{jk}Bjk
b^j(k)=Bjk∑s=1MBjs\hat{b}_j(k)={B_{jk} \over \sum_{s=1}^MB_{js}} b^j(k)=∑s=1MBjsBjk初始状态π(i)\pi(i)π(i)的概率估计
假设所有样本中初始隐藏状态为qiq_iqi的频数为C(i)C(i)C(i),那么初始概率分布为:
π^(i)=C(i)∑s=1NC(s)\hat{\pi}(i)={C(i) \over \sum_{s=1}^NC(s)} π^(i)=∑s=1NC(s)C(i)
无监督学习
很多时候,我们的训练数据中只有观察序列,而无隐藏的状态序列,此时模型训练问题/学习问题就是:给定SSS个长度为同为TTT的观测序列,即{O1,O2,...OS}\{O_1,O_2,...O_S\}{O1,O2,...OS}是已知的,求合适的模型参数λ\lambdaλ,使得训练数据发生的概率最大,即:
λ=argmaxλP(Otrain∣λ)\lambda=argmax_{\lambda}P(O_{train}|\lambda) λ=argmaxλP(Otrain∣λ)
该问题的求解最常用的是鲍姆-韦尔奇(Baum-Welch)算法,其实就是基于EM算法的求解。此时隐马尔可夫模型实际上是含有一个隐变量的概率模型:
P(O∣λ)=∑tP(O∣I,λ)P(I∣λ)P(O|\lambda)=\sum_tP(O|I,\lambda)P(I|\lambda) P(O∣λ)=t∑P(O∣I,λ)P(I∣λ)
EM算法简介
EM算法也称期望最大化(Expectation-Maximum,简称EM)算法,它是一个基础算法。每次迭代分为E步和M步两步。
算法的引入:概率模型有时既含有观测变量,又含有隐含变量。如果概率模型的变量都是观测变量,那么给定数据,可用利用最大似然估计方法,或者贝叶斯方法去估计模型参数。但当模型含有隐含变量时,前面的方法就不能直接应用,此时需要引入EM算法进行解决。EM算法是一种迭代算法,用于含有隐含变量的概率模型参数的最大似然估计。
算法的基本思想:初始随机给模型参数赋值并满足参数限制;给参数赋初值之后,得到模型μ0\mu_0μ0,根据模型μ0\mu_0μ0可以得到隐变量的期望值。然后利用隐变量的期望值代替隐变量的真实值,得到模型参数新的估计值μ1\mu_1μ1。执行迭代这个过程,直到参数收敛于最大似然估计值。
一个最直观了解EM算法思路的是K-Means算法。在K-Means聚类时,每个聚类簇的质心是隐含数据。我们会假设K个初始化质心,即EM算法的E步;然后计算得到每个样本最近的质心,并把样本聚类到最近的这个质心,即EM算法的M步。重复这个E步和M步,直到质心不再变化为止,这样就完成了K-Means聚类。
对于样本的观察数据x=(x(1),x(2),...,x(m))x=(x^{(1)},x^{(2)},...,x^{(m)})x=(x(1),x(2),...,x(m)),模型参数为θ\thetaθ,未观察到的隐含数据为z=(z(1),z(2),...,z(m))z=(z^{(1)},z^{(2)},...,z^{(m)})z=(z(1),z(2),...,z(m))。极大化模型分布的对数似然函数,则θ\thetaθ的估计值满足:
θ=argmaxθ∑i=1mlogP(x(i)∣θ)=argmaxθ∑i=1mlog∑z(i)P(x(i),z(i)∣θ)\theta=argmax_\theta\sum_{i=1}^mlogP(x^{(i)}|\theta)=argmax_\theta\sum_{i=1}^mlog\sum_{z(i)}P(x^{(i)},z^{(i)}|\theta) θ=argmaxθi=1∑mlogP(x(i)∣θ)=argmaxθi=1∑mlogz(i)∑P(x(i),z(i)∣θ)
这个式子是无法直接求出θ\thetaθ的,需要一些特殊技巧。对上式进行缩放:
∑i=1mlog∑z(i)P(x(i),z(i)∣θ)=∑i=1mlog∑z(i)Q(z(i))P(x(i),z(i)∣θ)Q(z(i))......(1)≥∑i=1m∑z(i)Q(z(i))logP(x(i),z(i)∣θ)Q(z(i))......(2)\begin{aligned} &\sum_{i=1}^mlog\sum_{z(i)}P(x^{(i)},z^{(i)}|\theta)\\ =& \sum_{i=1}^mlog\sum_{z(i)}Q(z^{(i)}){P(x^{(i)},z^{(i)}|\theta) \over Q(z^{(i)})}......(1)\\ \ge& \sum_{i=1}^m\sum_{z(i)}Q(z^{(i)})log{P(x^{(i)},z^{(i)}|\theta) \over Q(z^{(i)})}......(2) \end{aligned} =≥i=1∑mlogz(i)∑P(x(i),z(i)∣θ)i=1∑mlogz(i)∑Q(z(i))Q(z(i))P(x(i),z(i)∣θ)......(1)i=1∑mz(i)∑Q(z(i))logQ(z(i))P(x(i),z(i)∣θ)......(2)
上式引入未知分布Q(z(i))Q(z{(i)})Q(z(i)),满足∑z(i)Q(z(i))=1\sum_{z(i)}Q(z{(i)})=1∑z(i)Q(z(i))=1不等号用到jensen不等式,即对于log函数,有:
log(E(x))≥E(log(x))log(E(x)) \ge E(log(x)) log(E(x))≥E(log(x))
当且仅当xxx为常数时取等号。
因此,当jensen不等式取等号时,有:
P(x(i),z(i)∣θ)Q(z(i))=c,c为常数{P(x^{(i)},z^{(i)}|\theta) \over Q(z^{(i)})}=c,c为常数 Q(z(i))P(x(i),z(i)∣θ)=c,c为常数
得到:
∑z(i)Q(z(i))=∑z(i)P(x(i),z(i)∣θ)c=1\sum_{z(i)}Q(z{(i)})=\sum_{z(i)}{P(x^{(i)},z^{(i)}|\theta) \over c}=1 z(i)∑Q(z(i))=z(i)∑cP(x(i),z(i)∣θ)=1
Q(z(i))=P(x(i),z(i)∣θ)c=P(x(i),z(i)∣θ)∑z(i)P(x(i),z(i)∣θ)=P(x(i),z(i)∣θ)P(x(i)∣θ)=P(z(i)∣x(i),θ)\begin{aligned} Q(z{(i)})=&{P(x^{(i)},z^{(i)}|\theta) \over c} \\ =& {P(x^{(i)},z^{(i)}|\theta) \over {\sum_{z(i)}P(x^{(i)},z^{(i)}|\theta)}}\\ =& {P(x^{(i)},z^{(i)}|\theta) \over {P(x^{(i)}|\theta)}}\\ =& P(z^{(i)}|x^{(i)},\theta) \end{aligned} Q(z(i))====cP(x(i),z(i)∣θ)∑z(i)P(x(i),z(i)∣θ)P(x(i),z(i)∣θ)P(x(i)∣θ)P(x(i),z(i)∣θ)P(z(i)∣x(i),θ)
如果Q(z(i))=P(z(i)∣x(i),θ)Q(z{(i)})=P(z^{(i)}|x^{(i)},\theta)Q(z(i))=P(z(i)∣x(i),θ),则上面(2)的等号是对数似然函数的下界。如果最大化这个下界,则也在尝试最大化我们的对数似然。即我们需最大化下式:
argmaxθ∑i=1m∑z(i)Q(z(i))logP(x(i),z(i)∣θ)Q(z(i))argmax_{\theta}\sum_{i=1}^m\sum_{z(i)}Q(z^{(i)})log{P(x^{(i)},z^{(i)}|\theta) \over Q(z^{(i)})} argmaxθi=1∑mz(i)∑Q(z(i))logQ(z(i))P(x(i),z(i)∣θ)
EM算法流程:
输入:观察数据x=(x(1),x(2),...x(m))x=(x^{(1)},x^{(2)},...x^{(m)})x=(x(1),x(2),...x(m)),联合分布p(x,z∣θ)p(x,z|\theta)p(x,z∣θ), 条件分布p(z∣x,θ)p(z|x,\theta)p(z∣x,θ), 最大迭代次数JJJ。
随机初始化模型参数θ\thetaθ的初值θ0\theta_0θ0。for JJJ from 1 to JJJ,开始EM算法迭代:
a) E步(固定θθθ):计算联合分布的条件概率期望:
Q(z(i))=P(z(i)∣x(i),θj))Q(z^{(i)})=P(z^{(i)}|x^{(i)},θ_j)) Q(z(i))=P(z(i)∣x(i),θj))
L(θ,θj)=∑i=1m∑z(i)Q(z(i))logP(x(i),z(i)∣θ)L(θ,θ_j)=\sum_{i=1}^m∑_{z^{(i)}}Q(z^{(i)})logP(x^{(i)},z^{(i)}|θ) L(θ,θj)=i=1∑mz(i)∑Q(z(i))logP(x(i),z(i)∣θ)
b) M步[Q(z(i))Q(z^{(i)})Q(z(i))固定,是个常数]:极大化L(θ,θj)L(θ,θ_j)L(θ,θj),得到θj+1θ_{j+1}θj+1:
θj+1=argmaxθL(θ,θj)θ_{j+1}=argmax_θL(θ,θ_j) θj+1=argmaxθL(θ,θj)
c) 如果θj+1θ_{j+1}θj+1已收敛,则算法结束。否则继续回到步骤a)进行E步迭代。
输出:模型参数θ\thetaθ。
EM算法也是一个逐步求解局部极值,并且结果与初值选择有关。
Baum-Welch算法过程:
可以当作是EM算法的一个具体实现,步骤如下:
输入:DDD个观测序列样本{(O1),(O2),...(OD)}\{(O_{1}),(O_{2}),...(O_{D})\}{(O1),(O2),...(OD)}
输出:HMM模型参数
随机初始化所有的πi,aij,bj(k)π_i,a_{ij},b_j(k)πi,aij,bj(k)对于每个样本d=1,2,...Dd=1,2,...Dd=1,2,...D,用前向后向算法计算γt(d)(i),ξt(d)(i,j)\gamma_t^{(d)}(i),\xi_t^{(d)}(i,j)γt(d)(i),ξt(d)(i,j),见常见概率计算部分。更新模型参数:
πi=∑d=1Dγt(d)(i)D\pi_i={\sum_{d=1}^D\gamma_t^{(d)}(i) \over D}\\ πi=D∑d=1Dγt(d)(i)
aij=∑d=1D∑t=1T−1ξt(d)(i,j)∑d=1D∑t=1T−1γt(d)(i)a_{ij}={\sum_{d=1}^D\sum_{t=1}^{T−1}\xi_t^{(d)}(i,j) \over \sum_{d=1}^D\sum_{t=1}^{T−1}\gamma_t^{(d)}(i)} aij=∑d=1D∑t=1T−1γt(d)(i)∑d=1D∑t=1T−1ξt(d)(i,j)
bj(k)=∑d=1D∑t=1,ot(d)=vkTγt(d)(j)∑d=1D∑t=1Tγt(d)(j)b_j(k)={\sum_{d=1}^D\sum_{{t=1},o^{(d)}_t=v_k}^T\gamma_t^{(d)}(j) \over \sum_{d=1}^D\sum_{t=1}^T\gamma_t^{(d)}(j)} bj(k)=∑d=1D∑t=1Tγt(d)(j)∑d=1D∑t=1,ot(d)=vkTγt(d)(j)如果πi,aij,bj(k)π_i,a_{ij},b_j(k)πi,aij,bj(k)的值已经收敛,则算法结束,否则回到第(2)步继续迭代。
以上就是鲍姆-韦尔奇算法的整个过程。