马尔可夫模型及隐马尔可夫模型(HMM)

马尔可夫模型

马尔可夫模型是由Andrei A. Markov于1913年提出的

∙ ∙ 设  S S是一个由有限个状态组成的集合

S={1,2,3,…,n−1,n} S={1,2,3,…,n−1,n} 随机序列  X X 在  t t时刻所处的状态为  qt qt，其中  qt∈S qt∈S，若有：

P(qt=j|qt−1=i,qt−2=k,⋯)=P(qt=j|qt−1=i) P(qt=j|qt−1=i,qt−2=k,⋯)=P(qt=j|qt−1=i) aij≥0∑jnaij=1 aij≥0∑jnaij=1

则随机序列  X X构成一个一阶马尔科夫链。(Markov Chain)

∙ ∙ 令  P(qt=j|qt−1=i)=P(qs=j|qs−1=i) P(qt=j|qt−1=i)=P(qs=j|qs−1=i)(说明是齐次的，即转移概率不随时间改变)，则对于所有的  i,j i,j有下面的关系成立：

aij=P(qt=j|qt−1=i)1≤i,j≤n aij=P(qt=j|qt−1=i)1≤i,j≤n

 

∙ ∙ 一阶马尔科夫模型可以描述为一个二元组  (S,A)S (S,A)S是状态的集合，而 A A是所有状态转移概率组成的一个  n n行  n n列的矩阵，其中每一个元素  aij aij为从状态i转移到状态  j j 的概率。

∙ ∙ 同有限状态自动机类似，状态转移关系也可以用状态转换图来表示。

 

马尔科可模型举例

∙ ∙ 天气的变化，三种状态{1(阴天), 2(多云), 3(晴天)}

∙ ∙ 今天的天气情况仅和昨天的天气状况有关

∙ ∙ 根据对历史数据的观察得到下列状态转移关系

∙ ∙ 对于马尔科夫模型，给定了观察序列，同时也就确定了状态转换序列。例如有关天气状况的观察序列

(晴 晴 晴 阴 阴 晴 云 晴)

则状态转换序列为

(3, 3, 3, 1, 1, 3, 2, 3 )

∙ ∙ 如果把晴天称为状态3的输出，阴天称为状态1的输出，多云称为状态2的输出。根据观察到的输出序列就可以决定模型中的状态转换序列。(状态和输出是一对一的关系)

坛子与小球

∙ ∙ 在一个房间中，假定有  N N个坛子，每个坛子中都装有不同颜色的小球，并且假定总共有  M M种不同颜色的小球

∙ ∙ 一个精灵在房间中首先随机地选择一个坛子，再从这个坛子中随机选择一个小球，并把小球的颜色报告给房间外面的人员记录下来作为观察值

∙ ∙ 精灵然后把球放回到坛子中，以当前的坛子为条件再随机选择一个坛子，从中随机选择一个小球，并报告小球的颜色，然后放回小球，如此继续…，随着时间的推移，房间外的人会得到由这个过程产生的一个小球颜色的序列

∙ ∙ 如果令每一个坛子对应与一个状态，可以用一个一阶马尔科夫过程来描述坛子的选择过程

∙ ∙ 在马尔科夫过程中，每个状态只有一个输出，但在坛子和小球的问题中。可以从每个坛子中拿出不同颜色的小球。也就是每个状态能按照特定的概率分布产生多个输出

∙ ∙ 在坛子与小球问题中，如果给定一个观察序列(不同颜色的小球序列)，不能直接确定状态转换序列(坛子的序列)，因为状态转移过程被隐藏起来了。所以这类随机过程被称为隐马尔科夫过程

隐马尔可夫模型

∙ ∙ 隐马尔可夫模型  λ λ可以表示为一个五元组  (S,V,A,B,π) (S,V,A,B,π)

            S S 是一组状态的集合

                     S=1,2,3,…,N S=1,2,3,…,N (状态  n n对应坛子  n n)

            V V是一组输出符号组成的集合(取Visible之意？)

                     V=v1,v2,v3,⋯,vM V=v1,v2,v3,⋯,vM ( v1 v1对应红色小球)

            A A是状态转移矩阵， N N行  N N列

                     A=[aij],aij=P(qt+1=j|qt=i),1≤i,j≤N A=[aij],aij=P(qt+1=j|qt=i),1≤i,j≤N

            B B是输出符号的概率分布

                     B=bj(k)bj(k) B=bj(k)bj(k) 表示在状态  j j时输出符号  vk vk的概率

                     bj(k)=P(vk|j),1≤k≤M,1≤j≤N bj(k)=P(vk|j),1≤k≤M,1≤j≤N

            π π是初始状态概率分布  π={πi} π={πi}

                    πi=P(q1=i) πi=P(q1=i)表示时刻1选择某个状态的概率

∙ ∙ 隐马尔可夫过程是一个双重随机过程，其中一重随机过程不能直接观察到，通过状态转移概率矩阵描述。另一重随机过程输出可以观察到的观察符号，这由输出概率来定义

利用隐马尔可夫模型生成观察序列

∙ ∙ 可以把隐马尔可夫模型看做符号序列的生成装置，按照一定的步骤，隐马尔可夫模型可以生成下面的符号序列：

O=(o1,o2,⋯,oT) O=(o1,o2,⋯,oT)

\1. \1.  令  t=1 t=1，按照初始状态概率分布  π π选择一个初始状态  q1=i q1=i

\2. \2.  按照状态i输出符号的概率分布  bi(k) bi(k)选择一个输出值  ot=vk ot=vk

\3. \3.  按照状态转移概率分布  aij aij选择一个后继状态  qt+1=j qt+1=j

\4. \4.  若  t<T t<T，令  t=t+1 t=t+1，并且转移到算法第2步继续执行，否则结束

抛掷硬币

∙ ∙ 三枚硬币，随机选择一枚，进行抛掷，记录抛掷结果。可以描述为一个三个状态的隐马尔科夫模型  λ λ

                  λ=(S,V,A,B,π) λ=(S,V,A,B,π)，其中  S={1,2,3}V={H,T} S={1,2,3}V={H,T}

                  A A 如下表所示                                        B B 如下表所示

                       

                 π={1/3,1/3,1/3} π={1/3,1/3,1/3}

∙ ∙  问题一： 给定上述模型，观察到下列抛掷结果的概率是多少？      O=(HHHHTHTTTT) O=(HHHHTHTTTT)

∙ ∙  问题二： 给定上述模型，若观察到上述抛掷结果，最可能的硬币选择序列(状态转换序列)是什么？

∙ ∙  问题三： 若上述模型中的状态转移矩阵  A A、状态输出概率  B B和初始状态分布  π π 均未知，如何根据观察序列得到它们？

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

∙ ∙ 给定HMM   λ=(A,B,π) λ=(A,B,π) 给定观察序列  O=(o1o2o3⋯oT) O=(o1o2o3⋯oT) 如何有效地计算出观察序列的概率，即  P(O|λ) P(O|λ)? (估算问题) (另一种语言模型)

∙ ∙ 给定 HMM  λ=(A,B,π) λ=(A,B,π) 给定观察序列  O=(o1o2o3⋯oT) O=(o1o2o3⋯oT) 如何寻找一个状态转换序列  q=(q1q2q3⋯qT) q=(q1q2q3⋯qT)，使得该状态转换序列最有可能产生上述观察序列？ (解码问题)

∙ ∙ 在模型参数未知或不准确的情况下，如何根据观察序列 O=(o1o2o3⋯oT) O=(o1o2o3⋯oT)求得模型参数或调整模型参数，即如何确定一组模型参数，使得  P(O|λ) P(O|λ)最大？ (学习问题 或 训练问题)

估算观察序列概率

∙ ∙ 对隐马尔可夫模型而言，状态转换序列是隐藏的，一个观察序列可能由任何一种状态转换序列产生。因此要计算一个观察序列的概率值，就必须考虑所有可能的状态转换序列。

∙ ∙ 上图表示了产生观察序列  O=(o1o2o3⋯oT) O=(o1o2o3⋯oT)的所有可能的状态转换序列。

∙ ∙ 给定  λ λ, 以及状态转换序列  q=(q1q2q3⋯qT) q=(q1q2q3⋯qT)产生观察序列  O=(o1o2o3⋯oT) O=(o1o2o3⋯oT)的概率可以通过下面的公式计算：            P(O|q,λ)=bq1(o1)bq2(o2)bq3(o3)⋯bqT(oT) P(O|q,λ)=bq1(o1)bq2(o2)bq3(o3)⋯bqT(oT) 给定  λ λ, 状态转换序列  q=(q1q2q3⋯qT) q=(q1q2q3⋯qT)的概率可以通过下面的公式计算：            P(q|λ)=πq1aq1q2aq2q3⋯aqT−1qT P(q|λ)=πq1aq1q2aq2q3⋯aqT−1qT 则  O O和  q q的联合概率为：            P(O,q|λ)=P(O|q,λ)P(q|λ) P(O,q|λ)=P(O|q,λ)P(q|λ)

∙ ∙ 考虑所有的状态转换序列，则

P(O|λ)=∑qP(O,q|λ)=∑q1q2,⋯,qTπq1bq1(o1)aq1q2bq2(o2)⋯aqT−1qTbqT(oT) P(O|λ)=∑qP(O,q|λ)=∑q1q2,⋯,qTπq1bq1(o1)aq1q2bq2(o2)⋯aqT−1qTbqT(oT)

∙ ∙ 理论上，可以通过穷举所有状态转换序列的办法计算观察序列  O O的概率

∙ ∙ 实际上，这样做并不现实。     可能的状态转换序列共有  NT NT个     需要做  (2T–1)NT (2T–1)NT次乘法运算， NT–1 NT–1次加法运算     若  N=5，T=100 N=5，T=100，则  (2×100−1)×5100≈1072 (2×100−1)×5100≈1072

∙ ∙ 需要寻找更为有效的计算方法

向前算法

∙ ∙ 向前变量  αt(i) αt(i) αt(i)=P(o1o2o3⋯ot,qt=i|λ) αt(i)=P(o1o2o3⋯ot,qt=i|λ)

∙ ∙  αt(i) αt(i)的含义是，给定模型  λ λ ，时刻  t t，处在状态 i i，并且部分观察序列为  o1o2o3⋯ot o1o2o3⋯ot的概率。

∙ ∙ 显然有  α1(i)=πibi(o1)(1≤i≤N) α1(i)=πibi(o1)(1≤i≤N)

∙ ∙ 若  αt(i)(1≤i≤N) αt(i)(1≤i≤N) 已知，如何计算  αt+1(i) αt+1(i)?

         αt+1(j)=[∑Ni=1αt(i)aij]bj(Ot+1)1≤t≤T−1,1≤j≤N αt+1(j)=[∑i=1Nαt(i)aij]bj(Ot+1)1≤t≤T−1,1≤j≤N

向前算法计算步骤

\1. \1.  初始化      α1(i)=πibi(o1)(1≤i≤N) α1(i)=πibi(o1)(1≤i≤N)

\2. \2.  迭代计算      αt+1(j)=[∑Ni=1αt(i)aij]bj(Ot+1)1≤t≤T−1,1≤j≤N αt+1(j)=[∑i=1Nαt(i)aij]bj(Ot+1)1≤t≤T−1,1≤j≤N

\3. \3.  终止      P(O|λ)=∑Ni=1αT(i) P(O|λ)=∑i=1NαT(i)

向后算法

∙ ∙ 向后变量  βt(i) βt(i)      βt(i)=P(ot+1ot+2⋯oT|qt=i,λ) βt(i)=P(ot+1ot+2⋯oT|qt=i,λ)

∙ ∙  βt(i) βt(i)的含义是，在给定模型 λ λ，时刻 t t，处在状态 i i，并且部分观察序列为 ot+1ot+2⋯oT ot+1ot+2⋯oT的概率

∙ ∙  βT(i)=1(1≤i≤N) βT(i)=1(1≤i≤N)

∙ ∙ 若  βt+1(j)(1≤j≤N) βt+1(j)(1≤j≤N)已知，如何计算  βt(i) βt(i) βt(i)=∑Nj=1aijbj(ot+1)βt+1(j)1≤t≤T−1,1≤j≤N βt(i)=∑j=1Naijbj(ot+1)βt+1(j)1≤t≤T−1,1≤j≤N

向后算法计算步骤

\1. \1.  初始化      βT(i)=1(1≤i≤N) βT(i)=1(1≤i≤N)

\2. \2.  迭代计算      βt(i)=∑Nj=1aijbj(ot+1)βt+1(j)1≤t≤T−1,1≤j≤N βt(i)=∑j=1Naijbj(ot+1)βt+1(j)1≤t≤T−1,1≤j≤N

\3. \3.  终止      P(O|λ)=∑Ni=1πibi(o1)β1(i) P(O|λ)=∑i=1Nπibi(o1)β1(i)

求解最佳状态转换序列

∙ ∙ 隐马尔可夫模型的第二个问题是计算出一个能最好解释观察序列的状态转换序列

∙ ∙ 理论上，可以通过枚举所有的状态转换序列，并对每一个状态转换序列  q q计算  P(O,q|λ) P(O,q|λ)，能使  P(O,q|λ) P(O,q|λ)取最大值的状态转换序列  q∗ q∗就是能最好解释观察序列的状态转换序列，即：

                                        

q∗=argmaxqP(O,q|λ) q∗=arg⁡maxqP(O,q|λ)

∙ ∙ 同样，这不是一个有效的计算方法，需要寻找更好的计算方法

维特比算法

∙ ∙ 维特比变量  δt(i) δt(i)          δt(i)=maxq1,q2,\cots,qt−1P(q1q2⋯qt−1,qt=i,o1o2⋯ot|λ) δt(i)=maxq1,q2,\cots,qt−1P(q1q2⋯qt−1,qt=i,o1o2⋯ot|λ)

∙ ∙  δt(i) δt(i)的含义是，给定模型  λ λ，在时刻  t t处于状态  i i，观察到  o1o2o3⋯ot o1o2o3⋯ot的最佳状态转换序列为  q1q2⋯qt q1q2⋯qt的概率。

∙ ∙  δ1(i)=πibi(o1),1≤i≤N δ1(i)=πibi(o1),1≤i≤N

∙ ∙ 若  δt(i)(1≤i≤N) δt(i)(1≤i≤N)已知，如何计算  δt+1(i) δt+1(i)?         \delta_{t+1}(j) = \left[{\mathop {\max}\limits_i \;\delta_t(i)a_{ij}\right]b_j(o_{t+1})\delta_{t+1}(j) = \left[{\mathop {\max}\limits_i \;\delta_t(i)a_{ij}\right]b_j(o_{t+1})

∙ ∙ 如何记录路径? 设定  T T个数组  ψ1(N),ψ2(N),⋯,ψT(N) ψ1(N),ψ2(N),⋯,ψT(N)          ψt(i) ψt(i)记录在时刻  t t到达状态i的最佳状态转换序列  t−1 t−1时刻的最佳状态

维特比算法步骤

\1. \1.  初始化

         δ1(i)=πibi(o1),1≤i≤N δ1(i)=πibi(o1),1≤i≤N

         ψ1(i)=0 ψ1(i)=0

\2. \2.  迭代计算

         δt(j)=max1≤i≤N[δt−1(i)aij]bj(ot)ψt(j)=argmax1≤i≤N[δt−1aij] δt(j)=max1≤i≤N⁡[δt−1(i)aij]bj(ot)ψt(j)=arg⁡max1≤i≤N⁡[δt−1aij]

\3. \3.  终止

        P^* = \mathop {\max}\limits_{1\le i \le N}}\ \left[delta_T(i)\right] \;\;\;q_T^* = \mathop {\arg \max}\limits_{1 \le i \le N}}\ \left[delta_T(i)\right]P^* = \mathop {\max}\limits_{1\le i \le N}}\ \left[delta_T(i)\right] \;\;\;q_T^* = \mathop {\arg \max}\limits_{1 \le i \le N}}\ \left[delta_T(i)\right]

\4. \4.  求解最佳路径

         q∗t=ψt+1(q∗t+1),t=T−1,T−2,\codts,1 qt∗=ψt+1(qt+1∗),t=T−1,T−2,\codts,1

参数学习

∙ ∙ 隐马尔科夫模型的第三个问题是如何根据观察序列  O=(o1o2o3⋯oT) O=(o1o2o3⋯oT)求得模型参数或调整模型参数，即如何确定一组模型参数使得  P(O|λ) P(O|λ)最大？

∙ ∙ 隐马尔科夫模型的前两个问题均假设模型参数已知，第三个问题是模型参数未知，求最佳模型的问题，是三个问题中最为困难的问题

有指导的参数学习

∙ ∙ 在模型  λ λ未知的情况下，如果给定观察序列的同时，也给定了状态转换序列，此时可以通过有指导的学习方法学习模型参数。例如给定下面的训练数据，可以通过最大似然估计法估计模型参数：

H/1H/1T/1T/2H/3T/5… H/1H/1T/1T/2H/3T/5…

T/2H/1T/2H/3H/3H/1… T/2H/1T/2H/3H/3H/1…

∙ ∙ 参数学习非常简单，在训练数据足够大的前提下，效果不错

∙ ∙ 缺点，状态信息未知时无法使用。或者要由人工标注状态信息，代价高

∙ ∙ 在NLP中，在无指导学习效果不佳时，需要采用有指导学习

无指导的参数学习

∙ ∙ 在模型  λ λ未知的情况下，如果仅仅给定了观察序列，此时学习模型的方法被称做无指导的学习方法

∙ ∙ 对于隐马尔科夫模型而言，采用无指导学习方法，没有解析方法。通常要首先给定一组不准确的参数，再通过反复迭代逐步求精的方式调整模型参数，最终使参数稳定在一个可以接受的精度

∙ ∙ 利用无指导的学习方法估计隐马尔科夫模型参数时，并不能一定保证求得最优模型，一般能得到一个局部最优模型

直观的想法

∙ ∙ 给定一组初始参数  (ABπ) (ABπ)

∙ ∙ 由于没有给定状态转换序列，无法计算状态转移的频率、状态输出的频率以及初始状态频率

∙ ∙ 假定任何一种状态转换序列都可能

∙ ∙ 对每种状态转换序列中的频次加权处理，计算状态转移、状态输出、以及初始状态的期望频次

∙ ∙ 利用计算出的期望频次更新  (ABπ) (ABπ)

∙ ∙ 权值如何选择         对状态转换序列  q q而言，选择  P(q|O,λ) P(q|O,λ)

∙ ∙ 理论上可行，现实不可行,要考虑所有的状态转移路径需要多次迭代，问题更为严重

∙ ∙ 需要更为有效的算法，即Baum-Welch算法

Baum-Welch Algorithm

∙ ∙ 定义变量  ξt(i,j) ξt(i,j)      ξt(i,j)=P(qt=i,qt+1=j|O,λ) ξt(i,j)=P(qt=i,qt+1=j|O,λ)

∙ ∙  ξt(i,j)=P(qt=i,qt+1=j|O,λ) ξt(i,j)=P(qt=i,qt+1=j|O,λ)含义是，给定模型  λ λ和观察序列  O O，在时刻  t t处在状态  i i，时刻  t+1 t+1处在状态  j j的期望概率

∙ ∙  ξt(i,j) ξt(i,j)可以进一步写成：

ξt(i,j)=P(qt=i,qt+1=j,O|λ)P(O|λ)=αt(i)aijbj(ot+1)βt+1(j)P(O|λ)=αt(i)aijbj(ot+1)βt+1(j)∑i=1N∑j=1Nαt(i)aijbj(ot+1)βt+1(j) ξt(i,j)=P(qt=i,qt+1=j,O|λ)P(O|λ)=αt(i)aijbj(ot+1)βt+1(j)P(O|λ)=αt(i)aijbj(ot+1)βt+1(j)∑i=1N∑j=1Nαt(i)aijbj(ot+1)βt+1(j)

∙ ∙ 定义变量  γt(i) γt(i)，令其表示在给定模型以及观察序列的情况下， t t时刻处在状态  i i的概率，则有：

γt(i)=∑j=1Nξt(i,j) γt(i)=∑j=1Nξt(i,j)

∑T−1t=1γt(i) ∑t=1T−1γt(i)观察序列  O O中从状态  i i出发的转换的期望概率

∑T−1t=1ξt(i,j) ∑t=1T−1ξt(i,j)观察序列  O O中从状态  i i到状态  j j的转换的期望概率

∙ ∙ 关于  π,A,B π,A,B，一种合理的估计方法如下

                

∙ ∙ 利用上述结论，即可进行模型估算

∙ ∙ 选择模型参数初始值，初始值应满足隐马尔科夫模型的要求，即：

∑i=1Nπi=1∑j=1N\aij=1,1≤i≤N∑k=1Mbj(k)=1,1≤j≤N ∑i=1Nπi=1∑j=1N\aij=1,1≤i≤N∑k=1Mbj(k)=1,1≤j≤N

∙ ∙ 将初始值代入前面的公式中，计算一组新的参数  π¯¯¯,A¯¯¯¯,B¯¯¯¯ π¯,A¯,B¯

∙ ∙ 再将新的参数代入，再次计算更新的参数

∙ ∙ 如此反复，直到参数收敛

∙ ∙ Baum-Welch算法是一种EM算法

∙ ∙ E-step:计算  ξt(i,j) ξt(i,j)和  γt(i) γt(i)

∙ ∙ M-step:估计模型  λ¯¯¯ λ¯

∙ ∙ 终止条件  |log(P(O|λi+1))−log(P(O|λi))|<ε |log⁡(P(O|λi+1))−log⁡(P(O|λi))|<ε

∙ ∙ Baum等人证明要么估算值 \overline{\lambda}}\overline{\lambda}}和估算前的参数值  λ λ相等，要么估算值  \overline{\lambda}}\overline{\lambda}}比估算前的参数值  λ λ更好的解释了观察序列  O O

∙ ∙ 参数最终的收敛点并不一定是一个全局最优值，但一定是一个局部最优值

 

L.R.Rabiner, A Tutorial on Hidden Markov Models and Selected Applications in Speech recognition, Proc. IEEE, 77(2): 257-286, 1989

隐马尔科夫模型的实现

∙ ∙ 浮点溢出问题

对于韦特比算法，采用取对数的方式对于Baum-Welch算法，采用放大因子对于向前算法采用放大因子以及取对数的方式 http://www.cnblogs.com/gemstone/archive/2012/09/05/2671577.html