3月在线机器学习班第14课笔记--HMM
机器学习
紧随着贝叶斯网络的知识,我们来开启隐马尔科夫模型方面的内容。在之前一直认为与马尔科夫相关的知识点都很难,即使难也得上啊,跟着邹博老师一步一步来去学就没问题了。具体课程笔记如下:
简要回顾
在贝叶斯网络上,记录了一种特殊的贝叶斯网络--马尔科夫模型。
形如下图的贝叶斯网络我们称之为马尔科夫模型。
如上图所示,M个离散结点形成一条链(各个结点只有一个父结点),每一个结点有K个状态,那么根据之前的结论,可以知道描述这M个结点需要的参数个数为K−1+(M−1)K(K−1)(有一个结点不包含父结点,剩余M-1个都包含父结点)。可以知道,这是关于长度M的线性函数。
由有向分离可知,在xi给定的条件下,xi+1和x1,x2,...,xi−1条件独立。即xi+1的分布状态只和xi有关,和其他变量条件独立,这种顺序演变的随机过程模型,叫做马尔科夫模型,因为其满足马尔科夫性质:
P(Xn+1=x|X0,X1,...,Xn)=P(Xn+1=x|Xn)
HMM定义
隐马尔科夫模型(HMM, Hidden Markov Model)可用标注问题,在语音识别、NLP、生物信息、模式识别等领域被实践证明是有效的算法。HMM是关于时间序列的概率模型,下图即为HMM的贝叶斯网络:
HMM描述的是由一个隐藏的马尔科夫链随机生成不可观测的状态随机序列,再由各个状态生成一个观测而产生观测随机序列的过程。上图中的z1,z2,...,zn称为状态序列,它是不可观测的隐藏的马尔科夫链。每个状态生成一个观测,由此产生的观测随机序列,称为观测序列;x1,x2,...,xn即为观测序列。
所谓的"隐",表示的即是我们无法观测到状态序列,而只能通过观测序列来预测状态序列,有一种"听音辨人"的感觉。看到一个生动形象的HMM的例子:一对恋人中,女生有时候喜欢打男朋友,每天十二点可能发作一次。最奇怪是生气时候打,开心时候也打,有时候无缘无故平静还是打。(假设想当天十二点的心情需要参照前天十二点的心情。)能看到的就是妹子外在行动——打人,或者不打人。看不到的,就是隐含的变量:妹子的心情: 生气,平静,或者高兴。
就上面的图,请问,在z1,z2不可观察的前提下,x1和z2独立吗?x1和x2独立吗?
回答都是不独立的。结合贝叶斯网络中的"tail-to-tail"来判断即可。不仅仅x1和x2之间不独立,而且x1和xn之间也不独立,只不过此时x1和xn相关性相比于x1和x2要小很多,毕竟x1到xn已经经过了很多种状态。这种观测序列之间不独立的性质,在很大程度上是符合我们实际情况的。
HMM由初始概率分布π、状态转移概率分布A观测概率分布B确定,即HMM三要素。让我们来看看如何定义这三个要素。首先让我们定义:
Q是所有可能的状态的集合,N是所有可能的状态数,有:Q={q1,q2,...,qN}
V是所有可能的观测的集合,M是所有可能的观测数,有:V={v1,v2,...,vM}
I是长度为T的状态序列,O是对应的观测序列,有:
I={i1,i2,...,iT} O={o1,o2,...,oT}
A是状态转移概率矩阵,A=|aij|N×N,其中
aij=P(it+1=qj|it=qi)
aij是在时刻
t处于状态
qi的条件下,时刻
t+1转移到状态
qj的概率。
由定义可知,
A矩阵中的第
k行表示的是状态
qk在下一时刻转移到另一个状态的所有可能概率,每一行的概率之和为1。
B是观测概率矩阵,B=|aik|N×M,其中
bik=P(ot=vk|it=qi)
bik是在时刻
t处于状态
qi的条件下生成观测
vk的概率。
由定义可知,
B矩阵中的第
k行表示的是状态
qk在该时刻生成某一观测的所有可能概率,每一行的概率之和为1。
π是初始状态概率向量,π=(πi)N,其中πi=P(i1=qi)表示的是时刻t=1处于状态qi的概率。
由此可知,π和A决定状态序列,B决定观测序列。我们使用λ=(A,B,π)来表示HMM的三元素。
HMM的两个基本性质
1.齐次假设:
P(it|it−1,ot−1,it−2,ot−2,...,i1,o1)=P(it|it−1)
即状态序列具有马尔科夫性质,时刻
t的状态只与时刻
t−1的状态有关。
2.观测独立性假设:
P(ot|it,it−1,ot−1,...,i1,o1)=P(ot|it)
即时刻
t的观测结果只与时刻
t的状态结果有关,和其他无关。
HMM应用举例
假设有3个盒子,编号为1,2,3,每个盒子里面都装有红白两种颜色的小球,数目如下:
盒子号 |
1 |
2 |
3 |
红球数 |
5 |
4 |
7 |
白球数 |
5 |
6 |
3 |
按照下面的方法抽取小球,得到球颜色的观测序列:按照(0.2,0.4,0.4)的概率选择第1个盒子,从该盒子随机抽出1个球,记录颜色后放回盒子;按照A给定的概率选择新的盒子,重复上述过程;最终得到观测序列:“红红白白红”。
这个实例中的各个参数:
- 状态集合:Q={盒子1,盒子2,盒子3};
- 观测集合:V={红,白};
- 状态序列和观测序列的长度T=5;
- HMM三要素已知,分别如下表示:
我们可以针对3要素来解析其相应含义:
初始概率分布π:选择第一个盒子的时候,1,2,3号盒子被选到的概率分别是:(0.2,0.4,0.4);
状态转移概率分布A:先看第一行,t时刻选择到了1号盒子之时,那么t+1时刻1,2,3号盒子被选到的概率分别是:(0.5,0.2,0.3);再看第二行,t时刻选择到了2号盒子之时,那么t+1时刻1,2,3号盒子被选到的概率分别是:(0.3,0.5,0.2);第三行类似。
观测概率分布B:这个例子中该矩阵较为简单,每一行即为每个盒子中红白球数量所占的比例;
HMM的3个基本问题
当我们更加清楚了HMM的三要素之后,我们就有了HMM的3个基本问题。这三个问题的解答也是理解HMM的核心关键所在。
- 概率计算问题
给定模型 λ=(A,B,π)和观测序列O=o1,o2,...,oT ,计算模型λ下观测序列O出现的概率P(O|λ)。我们也可以将该问题看做是衡量一个给定模型对观测序列O的匹配好坏问题。
- 预测问题
即解码问题:已知模型λ=(A,B,π)和观测序列O=o1,o2,...,oT,求对给定观测序列条件概率P(I|O)最大的状态序列I。由已知的观测序列推测隐含的状态序列,相当于揭开该问题的神秘面纱。
- 学习问题
已知观测序列O=o1,o2,...,oT,估计模型λ=(A,B,π)的参数,使得在该模型下观测序列出现的概率P(O|λ)最大。这个问题会由是否知道状态序列而分两种情况:
1.如果有一组已知的状态序列以及与之相对应的观测序列,我们称之为训练序列,那么此时这些序列就是用于训练HMM的模型,使得模型参数在训练集中表现最好。这就是监督学习中的训练阶段,可以很快地通过极大似然估计得到相应参数。
2.如果只有一组已知的观测序列,而不知道状态序列,那么这个时候就是非监督学习,智能通过EM算法来解决。下面的相关解法针对的就是不知道状态序列的情况。
概率计算问题的解法
关于该问题的解答,有三种方法,包括直接算法(即暴力算法),前向算法和后向算法。后面两种算法是理解HMM的算法重点。此问题的解答过程如果能够很好理解了,那么后面的两个问题在理解上将会比较简单。可以说理解这个问题是基础。
直接计算法
直接计算法的思想是这样的:按照概率公式,列举所有可能的长度为T的状态序列I=i1,i2,...,iT,求各个状态序列I与观测序列的联合概率P(O,I|λ),然后对所有可能的P(O,I|λ)求和,从而得到P(O|λ)。
算法简单直接,具体步骤如下:
- 任意一个状态序列I=i1,i2,...,iT的概率是:
P(I|λ)=πi1ai1i2ai2i3...aiT−1iT
- 对固定的状态序列I,已知的观测序列O的概率是:
P(O|I,λ)=bi1o1bi2o2...biToiT
- 已知的观测序列O和I同时出现的联合概率是P(O,I|λ)=P(O|I,λ)P(I|λ) ,代入可得:
P(O,I|λ)=πi1bi1o1ai1i2bi2o2...aiT−1iTbiToiT
- 对所有可能的状态序列I求和,得到观测序列O的概率P(O|λ):
P(O|λ)=∑IP(O,I|λ)=∑Iπi1bi1o1ai1i2bi2o2...aiT−1iTbiToiT
针对最终的表达式,我们计算其时间复杂度:由于观测序列O是确定的,我们只需考虑状态序列I;枚举所有可能的状态序列I的时间复杂度为O(NT),再来看加和符号中有2T个相乘的因子,因此,计算最终表达式的时间复杂度为O(TNT)。这个复杂度真的太高了。我们需要更加有效的计算算法。
前向算法
该算法没有直接计算相关的概率,而是定义了一个新的概率:前向概率。具体定义:给定λ,定义到时刻t部分观测序列为o1,o2…ot且状态为qi的概率为前向概率,记做:
αt(i)=P(o1,o2…ot,it=qi|λ)
有了该定义,就可以通过递推的方法最终得到观测序列概率
P(O|λ)哦。挺屌的!让我们来看看具体流程:
1. 初值:
α1(i)=πibio1
2. 递推:对于
t=1,2,...,T−1,
αt+1(j)=(∑i=1Nαt(i)aij)bjot+1
3. 最终得到:
P(O|λ)=∑Ni=1αT(i)流程只有三步,但是真正理解起来却不简单。
首先,根据定义,初值α1(i)为在第一个时刻观测序列为o1且状态为qi的概率。结合观测序列是已知的,这理解起来不难。
递归的过程中较为重要的即为求和部分(∑Ni=1αt(i)aij),求和内部的含义为在t时刻状态为qi而t+1时刻状态为qj的概率,针对i求和即t+1时刻状态为qj的概率(即表示不管t时刻状态是什么,t+1时刻状态为qj的概率)。可以用下图来表示这个过程:
经过递推,我们可以得到时刻T观测序列为o1,o2…oT(即整个观测序列)且状态为qi的概率αT(i),对其求和即可求得最终我们需要的P(O|λ)。
可以验证,前向概率算法的时间复杂度为O(TN2)。
即学即用,考察盒子球模型,计算观测向量O=“红白红”的出现概率。模型参数不变,如下:
步骤一,计算初值(已知信息:观测序列O=“红白红”):
α1(1)=π1b1o1=0.2×0.5=0.1α1(2)=π2b2o1=0.4×0.4=0.16α1(3)=π3b3o1=0.4×0.7=0.28
步骤二:递推(t=1):
α2(1)=(∑j=1Nα1(j)aj1)b1o2=(0.1∗0.5+0.16∗0.3+0.28∗0.2)∗0.5=0.077α2(2)=(∑j=1Nα1(j)aj2)b2o2=(0.1∗0.2+0.16∗0.5+0.28∗0.3)∗0.6=0.1104α2(3)=(∑j=1Nα1(j)aj3)b3o2=(0.1∗0.3+0.16∗0.2+0.28∗0.5)∗0.3=0.0606
递推(t=2):
α3(1)=(∑j=1Nα2(j)aj1)b1o3=(0.077∗0.5+0.1104∗0.3+0.0606∗0.2)∗0.5=0.04187α3(2)=(∑j=1Nα2(j)aj2)b2o3=(0.077∗0.2+0.1104∗0.5+0.0606∗0.3)∗0.4=0.03551α3(3)=(∑j=1Nα2(j)aj3)b3o3=(0.077∗0.3+0.1104∗0.2+0.0606∗0.5)∗0.7=0.05284
步骤三:最终:
P(O|λ)=∑Ni=1αT(i)=0.04187+0.03551+0.05284=0.13022后向算法
实际上,使用前向算法已经能够很好解决概率计算问题了,但是为了更方便地解决接下来的问题,我们引入后向算法。和前向算法相似,该算法也定义了一个新的概率:后向概率。具体定义:给定模型λ,定义在时刻t的状态为qi的前提下,从t+1到T的部分观测序列为ot+1,ot+2,...,oT的概率为后向概率,记做:
βt(i)=P(ot+1,ot+2,...,oT|it=qi, λ)
和前向算法相似,经过递推,就可以得到观测序列概率
P(O|λ)。具体流程如下:
1. 初值:
βT(i)=1
2. 递推:对于
t=T−1,T−2,...,1,
βt(i)=∑j=1N(aijbjot+1βi+1(j))
3. 最终得到:
P(O|λ)=∑Ni=1πibio1β1(i)递归:假设已知t时刻的状态为qi,我需要得到ot+1,ot+2,...,oT出现的概率。首先,需要知道所有可能的在t+1时刻的状态qj(即需转化概率项aij);只有在t+1时刻的状态qj的前提下,才可能得到t+1时刻的观测ot+1(即需观测概率项bjot+1);得到了t+1时刻的观测ot+1,乘上状态qj时ot+2,ot+3,...,oT出现的概率(即后向概率βi+1(j))),即得到了我们需要的概率值。
可以用下图来表示这个过程:
预测问题的解法
和概率计算问题不同的是,预测问题不仅仅存在一个解,而是需要从众多的解中选取一个最优解:找到一个最优的状态序列,在给定一个观测序列和已知模型参数的前提下。该问题的难点在于,如何选取一个衡量是否是最优解的准则?一个可能的衡量准则是在时刻t选择最可能的状态it(individually most likely)。为了求解预测问题,我们定义变量:
γt(i)=P(it=qi|O,λ)
表示的是在模型
λ已知,观测序列
O给定的条件下,时刻
t状态为
qi的概率。
γt(i)可以使用前向-后项概率来表示:
γt(i)=P(it=qi,O|λ)P(O|λ)=αt(i)βt(i)∑Ni=1αt(i)βt(i)
因为根据前向概率与后向概率
αt(i)、βt(i)的定义,可以得到
P(it=qi,O|λ)=αt(i)βt(i)。
另外,从等式中可以看出:
∑i=1Nγt(i)=1
使用γt(i),我们可以求解得到在衡量准则--时刻t选择最可能的状态it(individually most likely)--的情况下,选择最优的解:
it=argmax1≤i≤N[γt(i)], 1≤t≤T
即每个时刻
t都选择使得
γt(i)最大的那个状态。
这样的衡量准则并不是没有缺点的,因为它只考虑了单个状态的最优,而并没有考虑与之相关的下一个状态的情况(我们知道隐藏的状态序列是一条马尔科夫链,相邻两个状态之间是有关联的)。当某个状态转移概率
aij=0时,该衡量准则得出的状态序列很有可能是无效的。
为了对这个衡量准则做出修正,我们可以选择考虑了相邻两个状态的准则,但是更常见的情况是我们考虑的是将整个状态序列(
T个状态同时考虑)都放进来考察,以此为基准建立衡量准则。即找到某一个最好的状态序列
I,使得
P(I|O,λ)最大。Viterbi算法即是一种满足要求的常用算法。
Viterbi算法
Viterbi算法实际上是利用动态规划的思想去解HMM,用动态规划求概率最大的路径(最优路径),这里的一条路径对应一个状态序列。在Viterbi算法中,定义变量:
δt(i)=maxi1,...,it−1P(it=i,it1,...,i1,ot,...,o1|λ)
表示沿着一条路径
i1,...,it−1,it,使得前
t个观测序列为
ot,...,o1,且
t时刻状态为
i的最大概率。
通过递归,我们可以得到:
δt+1(j)=[maxi(δt(i)αij)]bjot+1
为了能够跟踪状态每一步所取的值(最后需要求解的是状态序列),我们定义数组:
ψt+1(j)=[argmaxi(δt(i)αij)]
完整的Viterbi算法流程如下:
- 初始化:
δ1(i)=πibio1, 1≤i≤Nψ1(i)=0, 1≤i≤N
- 递归:
δt(j)=[maxi(δt−1(i)αij)]bjot, 2≤t≤T,1≤j≤Nψt(j)=argmaxi[δt−1(i)αij], 2≤t≤T,1≤j≤N
- 停止:
p∗=maxi[δT(i)]i∗T=argmaxi[δT(i)]
- 最优状态序列的查找:
i∗t=ψt+1(i∗t+1)
学习问题的解法
在观测序列已知、状态序列不可知的情况下,我们采用EM算法来解决学习问题,HMM中的EM算法我们称之为Baum-Welch算法。
定义变量:
εt(i,j)=P(it=qi,it+1=qj|O,λ)
表示的是在给定模型参数
λ与观测序列
O的条件下,时刻
t的状态为
qi且时刻
t+1的状态为
qj的概率。
而根据定义可以得到:
P(it=qi,it+1=qj,O|λ)=αt(i)aijbjot+1βt+1(j)
由此可得:
εt(i,j)=P(it=qi,it+1=qj,O|λ)P(O|λ)=αt(i)aijbjot+1βt+1(j)∑Ni=1∑Nj=1αt(i)aijbjot+1βt+1(j)
相关的变量在图表中的表示为:
之前我们就已经定义过
γt(i),表示的是在模型
λ已知,观测序列
O给定的条件下,时刻
t状态为
qi的概率。如今,我们可以将
γt(i)与
εt(i,j)关联在一起:
γt(i)=∑j=1Nεt(i,j)
如果我们将时刻
1≤t≤T−1的
γt(i)求和,可以得到一个新的变量,它用于衡量在
t=1到t=T−1这段时间内状态
qi在状态转换过程中被访问到的次数的均值。类似的,我们时刻
1≤t≤T−1的
εt(i)求和,表示的含义是在
t=1到t=T−1这段时间内从状态
qi被转换到状态
qj在状态转换过程中被访问到的次数的均值。即:
使用上面定义的公式,我们可以给出一种HMM对当前参数进行更新的方法。假设初始的模型为
λ=(A,B,π),那么更新的参数为:
π¯i=γ1(i), 1≤i≤Na¯ij=∑T−1t=1εt(i,j)∑T−1t=1γt(i) 1≤i,j≤Nb¯jk=∑T−1t=1,ot=vkγt(i)∑T−1t=1γt(i)
更新过后的模型参数
λ¯=(A¯,B¯,π¯)。Baum及其同事已经证明了,
P(O|λ¯)≥P(O|λ)。具体证明过程如下:
定义变量:
Q(λ,λ¯)=∑I[(lnP(O,I|λ))P(O,I|λ¯)]
根据之前得到的等式:
P(O,I|λ)=πi1bi1o1ai1i2bi2o2...aiT−1iTbiToiT
上式可写成:
Q(λ,λ¯)=∑I[(lnP(O,I|λ))P(O,I|λ¯)]=∑Ilnπi1P(O,I|λ¯)+∑I(∑t=1T−1lnaitit+1)P(O,I|λ¯)+∑I(∑t=1Tlnbitot)P(O,I|λ¯)
极大化
Q(λ,λ¯)所得的参数
λ¯,可推出:
maxλ¯[Q(λ,λ¯)]⟹P(O|λ¯)≥P(O|λ)
由于这三个参数分别位于三个项中,可分别极大化。领用拉格朗日乘子法即可分别求得三个参数的最优解表达式,具体表达式上面已经推导出,这里不再复述。值得注意的是,由Baum-Welch算法最终的到的最优参数结果并不是全局最优解,而只是局部最优解。
参考文献:
A Tutorial on Learning With Bayesian Networks,David Heckerman, 1996