隐马尔可夫模型HMM
文章目录
一.简介
我们之前学的很多模型,比如说线性回归模型,甚至包括与HMM有很强联系的朴素贝叶斯算法都有这样一个前提:样本之间默认是相互独立的。
对于样本之间并不独立的情况,就比如说LDA模型应用的词分布,主题分布当中。一个文章的一个单词,和与这个单词紧挨着的词往往具有很强的关联性。在这个情况下,HMM的效果往往要更好一些。
那么,这个HMM到底怎么理解呢?首先得先知道什么是马尔可夫模型
马尔可夫模型(Markov Model)是通过寻找事物状态的规律对未来事物状态进行预测的概率模型,在马尔可夫模型中假设当前事物的状态只与之前的n个状态有关。这个,我们与贝叶斯网络联系在一起。马尔可夫模型其实就是一种特殊的贝叶斯网络:
我们看到,所有结点,都和其所有的父亲结点有着密切的联系。
那么什么是隐马尔可夫模型呢?
隐马尔可夫模型在马尔可夫模型的基础上,加入了新结点,如下(图一):
这当中,z1, z2等都是不可观察的。因此不同于马尔可夫模型,如果说马尔可夫模型是:根据在现在时间之前,出现雨天的概率来预测明天是不是雨天;那么隐马尔可夫模型就是,我并不知道之前下雨的概率,但我可以通过现在观察到的现象,比如说,现在是不是乌云密布,气压很低,来预测一下下雨的概率。
二.HMM原理
1. 隐马尔可夫图的说明
我们就从上面这个隐马尔可夫图入手,来解释一下这幅图是什么意思。
首先,隐马尔科夫模型随机生成的状态随机序列(就是图中的z1, z2等等,每一个结点都有若干个状态属性),称为状态序列;每个状态生成一个观测,由此产生的观测随机序列,称为观测序列(就是下方的x1, x2等等)。
在贝叶斯网络当中,两个都没有入度的结点是独立的。那么显然,在这幅图的情况下,x1和z2并不是独立的,x1,x2也是不独立的。
2. HMM原理
HMM由初始概率分布π、状态转移概率分布A以及观测概率分布B确定,表达式如下:
λ
=
(
A
,
B
,
π
)
\lambda = (A,B,\pi)
λ=(A,B,π)
A,B π 到底是啥意思呢?我们回图一
我们首先看从z1到z2,到z3……一直到zn,每一个z都有不同的隐状态,我们直接看z1到z2,这两个结点之间的转换就会有不同的状态转移,我们假设z1有n个隐状态,z2也有n个隐状态。设:从z1的i状态转移到z2的j状态的概率为aij,于是我们就有了一个z1,z2的联合密度:
z1 z2 | 1 | 2 | …… | n |
---|---|---|---|---|
1 | a11 | a12 | a1n | |
2 | a21 | a22 | a2n | |
…… | …… | …… | …… | …… |
n | an1 | an2 | …… | ann |
如此一来,组成了一个n*n的方阵,我们记作:A,称为:状态转移概率矩阵。
从z到x也是可以转移的,我们看z1到x1:
z1有n个隐状态,x1有m个观测值。通过z1的i状态观测到x1的j号观测值的概率,我们记作bij,这样一来,也可以组成一个n*m的矩阵,我们记作:B,称为:观测概率分布矩阵
现在从z(n-1)到z(n)都满足A矩阵。从zn到xn都满足矩阵B。这当中z1又比较特殊,因为他是所有结点共同的父节点,我们知道,z1有n个隐性状态(i1,i2,……in),取一号状态的概率为π1, 取二号状态的概率是π2……,这些概率在一起,就会组成一个向量π={π1, π2, π3, ……},我们记作:π。
这就是A, B, π的含义。
对于z的隐状态,要求:一定是离散的,这个是隐性马尔科夫的硬性要求。如果是连续的,那是别的模型(卡尔曼滤波模型)。
对于X中的可观测值,这个没有硬性要求,若(Xi|Zi)服从高斯分布(~N(ui, ai^2)),那么整体上就是:高斯隐性马尔科夫模型
3. HMM当中的概率计算
这个涉及的问题如下:
给定模型 : λ = ( A , B , π ) \lambda=(A,B,\pi) λ=(A,B,π)和观测序列 : O = { o 1 , o 2 , … … , o T } O=\{o_{1},o_{2},……,o_{T}\} O={o1,o2,……,oT} ,计算模型λ下观测序列O出现的概率P(O| λ)
这个问题,有不同的计算方法。
3.1 直接计算(不做重点)
这种方法,其实应该称为“暴力法”,比较简单直接:
按照概率公式,列举所有可能的长度为T的状态序列 ,求各个状态序列:
I
=
{
i
1
,
i
2
,
…
…
,
i
T
}
I=\{i_{1},i_{2},……,i_{T}\}
I={i1,i2,……,iT}与观测序列 :
O
=
{
o
1
,
o
2
,
…
…
,
o
T
}
O=\{o_{1},o_{2},……,o_{T}\}
O={o1,o2,……,oT} 的联合概率P(O,I|λ),然后对所有可能的状态序列求和,从而得到P(O|λ)。如果把这段文字写成数学表达式,则如下所示:
P
(
O
∣
λ
)
=
∑
I
P
(
O
,
I
∣
λ
)
=
∑
I
P
(
O
∣
I
,
λ
)
P
(
I
∣
λ
)
P(O|\lambda) = \sum_{I}P(O,I|\lambda)=\sum_{I}P(O|I,\lambda)P(I|\lambda)
P(O∣λ)=I∑P(O,I∣λ)=I∑P(O∣I,λ)P(I∣λ)
那么实际上,你只要求出后面这两项,就可以求出结果了,那么最后这两个又怎么计算呢?其实很简单:
当模型确定时,状态序列选择的概率,其实就是(结合图一),你先选中z1,然后一层一层的传递下去
P
(
I
∣
λ
)
=
π
i
1
a
i
1
i
2
a
i
2
i
3
…
…
a
i
T
−
1
i
T
P(I|\lambda)= \pi_{i_{1}}a_{i_{1}i_{2}}a_{i_{2}i_{3}}……a_{i_{T-1}i_{T}}
P(I∣λ)=πi1ai1i2ai2i3……aiT−1iT
如果状态已经固定了,那么计算观察到某一个状态的概率,其实就是从z转移到x的过程:
P
(
O
∣
I
,
λ
)
=
b
i
1
o
1
b
i
2
o
2
…
…
b
i
T
o
T
P(O|I,\lambda) = b_{i_{1}o_{1}}b_{i_{2}o_{2}}……b_{i_{T}o_{T}}
P(O∣I,λ)=bi1o1bi2o2……biToT
然后把上面这些代入到原式即可。
3.2 前向-后向算法(**)
这个是需要重点掌握的。
alpha是前向概率,beta是后向概率,t是分隔结点
为了和3.1部分统一,我们现在把y部分全都改成o,状态则全用I表示。
3.2.1 前向概率
前向概率:给定λ,定义到时刻t部分观测序列为o1,o2…ot且状态为qi的概率称为前向概率,记作:
α
t
(
i
)
=
P
(
o
1
,
o
2
,
.
.
.
,
o
t
,
i
t
=
q
i
∣
λ
)
\alpha_{t}(i) = P(o_{1},o_{2},...,o_{t},i_{t}=q_{i}|\lambda)
αt(i)=P(o1,o2,...,ot,it=qi∣λ)
根据公式,我们很容易得到前向概率的初值:
α
1
(
i
)
=
P
(
i
1
=
q
i
∣
λ
)
=
π
i
b
i
o
1
\alpha_{1}(i) =P(i_{1}=q_{i}|\lambda) = \pi_{i}b_{io_{1}}
α1(i)=P(i1=qi∣λ)=πibio1
能否通过当前时间的前向概率,来向后预测呢?
这当中,需要注意的是:alpha_(t+1)要求我们观测o1,o2,……,o_(t+1)这些隐属性, alpha_t要求我们观测o1, o2, o3, ……ot。这是有区别的。其实,我们只要从alpha_{t} 转移到alpha_{t+1},然后,通过alpha_{t+1}往o_{t+1}转移即可。由于,具体是从alpha_{t}的哪个属性转移到alpha_{t+1}的i属性并不知道,所以,我们应该把alpha_{t}的所有属性都考虑在内(因此我们要做一个累加)。最后再乘以从alpha_{t+1}转移到o_{t+1}的概率。于是我们就得到了如下公式:
α
t
+
1
(
i
)
=
(
∑
i
=
1
n
α
t
(
j
)
a
j
i
)
⋅
b
i
o
t
+
1
\alpha_{t+1}(i) = (\sum_{i=1}^n \alpha_{t}(j)a_{ji})\cdot b_{io_{t+1}}
αt+1(i)=(i=1∑nαt(j)aji)⋅biot+1
这个,就是前向算法的递推公式。
最后的结果就是:
P
(
O
∣
λ
)
=
∑
i
=
1
N
α
T
(
i
)
P(O|\lambda) = \sum_{i=1}^{N}\alpha_{T}(i)
P(O∣λ)=i=1∑NαT(i)
3.2.2 后向概率
有了前向概率的介绍,后向概率其实就不难了,我们还是从定义入手:
**后向概率:**给定λ,定义到时刻t状态为qi的前提下,从t+1到T的部分观测序列为o t+1 ,o t+2 …o T的概率为后向概率,记做:
β
t
(
i
)
=
P
(
o
t
+
1
,
o
t
+
2
,
.
.
.
,
o
T
∣
i
t
=
q
i
,
λ
)
\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
)
=
1
\beta_{T}(i) = 1
βT(i)=1
递推公式(通过现在的时间,向前推前一个),具体思路是:为了计算在时刻t状态为qi条件下时刻t+1之后的观测序列为o t+1 ,o t+2 …o T 的后向概率βt (i),只需要考虑在时刻t+1所有可能的N个状态qj的转移概率(aij ),以及在此状态下的观测o_{t+1} 的观测概率(B矩阵元素) ,然后考虑状态qj之后的观测序列的后向概率βt+1 (j)
β
t
(
i
)
=
∑
j
=
1
N
(
a
i
j
b
j
o
t
+
1
β
t
+
1
(
j
)
)
\beta_{t}(i) = \sum_{j=1}^{N}(a_{ij}b_{jo_{t+1}}\beta_{t+1}(j))
βt(i)=j=1∑N(aijbjot+1βt+1(j))
最终:
P
(
O
∣
λ
)
=
∑
i
=
1
N
π
i
b
i
o
1
β
1
(
i
)
P(O|\lambda) = \sum_{i=1}^{N}\pi_{i}b_{io_{1}}\beta_{1}(i)
P(O∣λ)=i=1∑Nπibio1β1(i)
3.2.3 前向概率与后向概率的关系
直接看公式推导:
此时,我们给出单个状态概率的定义:当给定模型λ和观测序列O的时候,在时刻t处于状态qi的概率,记作:
γ
t
(
i
)
=
P
(
i
t
=
q
i
∣
O
,
λ
)
\gamma_{t}(i) = P(i_{t}=q_{i}|O,\lambda)
γt(i)=P(it=qi∣O,λ)
我们根据条件概率公式:
假如,在给定模型λ和观测序列O的情况下,我们想求在时刻t处于状态qi并且时刻t+1处于状态qj的概率呢?
其中:
P
(
i
t
=
q
i
,
i
t
+
1
=
q
j
,
O
∣
λ
)
=
α
t
(
i
)
a
i
j
b
j
o
t
+
1
β
t
+
1
(
j
)
P(i_{t}=q_{i},i_{t+1}=q_{j},O|\lambda) = \alpha_{t}(i)a_{ij}b_{jo_{t+1}}\beta_{t+1}(j)
P(it=qi,it+1=qj,O∣λ)=αt(i)aijbjot+1βt+1(j)
4. HMM三个问题之二:学习问题(Baum-Welch算法)
不同于概率计算。这个是另外一个情景:若训练数据只有观测序列,则HMM的学习需要使用EM算法,是非监督学习
既然这个算法,是基于EM算法的,那么不要忘记EM算法是啥:
学习问题是什么呢?已知观测序列 O={o1, o2, …, oT} ,估计模型 的参数 λ=(A , B, Pi),使得在该模型下观测序列P(O| λ)最大
在HMM 的学习问题当中,我们通常用的是Baum-Welch算法。
首先,我们已经知道了观测序列O,我们把所有的隐状态写成:I={i1, i2, ……, iT}。这样,完全的数据就是:(O,I)=(o1,o2,…,oT,i1,i2,…kT)。我们看看学习问题的内容。完全数据的对数似然函数是lnP(O,I|λ)。假设
λ
ˉ
\bar{\lambda}
λˉ
是HMM参数的当前估计值,λ为待求的参数
我们按照EM算法所给出的函数来构造一个函数:
我们在3.1当中,已经推出:
P
(
O
∣
I
,
λ
)
=
P
(
O
∣
I
,
λ
)
P
(
I
∣
λ
)
=
b
i
1
o
1
b
i
2
o
2
…
…
b
i
T
o
T
P(O|I,\lambda)=P(O|I,\lambda)P(I|\lambda) = b_{i_{1}o_{1}}b_{i_{2}o_{2}}……b_{i_{T}o_{T}}
P(O∣I,λ)=P(O∣I,λ)P(I∣λ)=bi1o1bi2o2……biToT
把这个式子代入到我们构造的函数当中:
注意:我们在求解这个问题的时候,别忘了一个隐藏条件:Pi向量当中所有的元素加和肯定是1,即:
π
1
+
π
2
+
.
.
.
+
π
T
=
1
\pi_{1} + \pi_{2} + ... + \pi_{T}=1
π1+π2+...+πT=1
所以,整个就可以转变成为一个条件极值问题:A, B, Pi取什么值,能够使得上面等式成立的同时,保证Q最大。很显然,我们应该用:拉格朗日乘数法。
上面的Q已经拆成了三项,因此,我们可以把这三项分别构造拉格朗日函数,去求。
我们先看第一项:
然后对pi求导,得到:
P
(
O
,
i
1
=
i
∣
λ
ˉ
)
+
γ
π
i
=
0
P(O,i_{1}=i|\bar{\lambda}) + \gamma\pi_{i} = 0
P(O,i1=i∣λˉ)+γπi=0
对于所有的pi_i,我们都求一次导数,那么就会得到T个上面这个形式的式子,我们把这T个式子进行加和就会得到:
γ
=
−
P
(
O
∣
λ
ˉ
)
\gamma = -P(O|\bar{\lambda})
γ=−P(O∣λˉ)
然后在代回到上面这个式子:
π
i
=
P
(
O
,
i
1
=
i
∣
λ
ˉ
)
P
(
O
∣
λ
ˉ
)
=
P
(
O
,
i
1
=
i
∣
λ
ˉ
)
∑
i
=
1
N
P
(
O
,
i
1
=
i
∣
λ
ˉ
)
=
γ
1
(
i
)
\pi_{i}=\frac{P(O,i_{1}=i|\bar{\lambda})}{P(O|\bar{\lambda})} = \frac{P(O,i_{1}=i|\bar{\lambda})}{\sum_{i=1}^{N} P(O,i_{1}=i|\bar{\lambda})} = \gamma_{1}(i)
πi=P(O∣λˉ)P(O,i1=i∣λˉ)=∑i=1NP(O,i1=i∣λˉ)P(O,i1=i∣λˉ)=γ1(i)
然后我们看第二项:
然后我们继续使用拉格朗日乘数法,只不过这次,我们对A求偏导(对aij求偏导):
然后,我们用展开式的第三项,用相同的方法,对B求偏导,就可以得到:
5. HMM三个问题之三:预测问题(Viterbi算法)
预测问题要干什么呢?已知模型λ=(A , B, Pi)和观测序列 O={o1, o2, …, oT},求给定观测序列条件概率P(I|O, λ )最大的状态序列I
所以,不难看出,这个其实是一个路径问题。预测问题的原理很简单:动态规划。这个时候,用到了另一个算法:Viterbi算法
Viterbi算法实际是用动态规划解HMM预测问题,用DP求概率最大的路径(最优路径),这是一条路径对应一个状态序列。
我们定义一个变量:
δ
t
(
i
)
\delta_{t}(i)
δt(i),用来表示在时刻t状态为i的所有路径中,概率的最大值,即:
δ
t
(
i
)
=
m
a
x
i
1
,
i
2
,
.
.
,
i
t
−
1
P
(
i
t
=
i
,
i
t
−
1
,
.
.
.
,
i
1
,
o
t
,
.
.
.
o
1
∣
λ
)
\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)
δt(i)=maxi1,i2,..,it−1P(it=i,it−1,...,i1,ot,...o1∣λ)
然后,我们就可以用动态规划的方式去递推:
终止条件:
P
∗
=
m
a
x
1
<
=
i
<
=
N
δ
T
(
i
)
P^{*} = max_{1<=i<=N}\delta_{T}(i)
P∗=max1<=i<=NδT(i)
6. HMM简单举例
给了这么多理论,比较抽象,实际上,HMM一个最简单的应用,就是盒子拿小球的例子。
假设有3个盒子,编号为1、2、3,每个盒子都装有红白两种颜色的小球,数目如下:
盒子编号 | 1 | 2 | 3 |
---|---|---|---|
红球数量 | 5 | 4 | 7 |
白球数量 | 5 | 6 | 3 |
假设,盒子大小并不相同,这就导致,拿盒子的概率并不相等。按照π=(0.2,0.4,0.4)的概率选择1个盒子,从盒子随机抽出
1个球,记录颜色后放回盒子。
于是,我们就可以建立一个初始概率分布:
这里很显然pi就是拿到盒子的概率
我们看A矩阵,第一行第一列的那个0.5,意思就是:先拿到1号盒子的前提下,下次又拿到了一号盒子的概率是0.5(这个A矩阵当中的数字是随意给的,并不是通过某个方式计算出来的,因此不必纠结)
再看B矩阵,这个就很明显了,第一行的意思是:拿到1号盒子的前提下,取到红球,白球的概率……
在这个模型之下,我们求:观测向量O="红白红"的出现概率。
首先,我们需要计算初值(alpha_{t}(i)表示:拿了i号盒子,并取到红球的概率(我们用1表示红球,2表示白球)),由于硬性要求O=“红白红”,因此,我们肯定要求第一次一定要拿到红球,但是具体从哪个盒子取到红球就不一定了。因此初值计算,就是要计算:第一次拿到红球的概率。
然后,我们按照前向概率的递推公式,一层层的往下递推:其中第二次,要求拿到白球了,因此在后面要乘以b_io_2,才可以保证拿到白球。相应的,在迭代alpha3的时候,我们后面要乘以b_i_o_1才能保证拿到红球。经过迭代之后,我们得到下面结果:
最终的结果就是:
那么,得到:o="红白红"的概率最大路径是什么呢?这个时候就需要Viterbi算法了:
首先,初始化:在t=1时,对于每一个状态i,求状态为i观测到o1=红的概率
δ
1
(
i
)
=
π
i
b
i
o
1
=
π
i
b
i
红
\delta_{1}(i) = \pi_{i}b_{io_{1}} =\pi_{i}b_{i红}
δ1(i)=πibio1=πibi红
根据这个公式很容易就得到:
δ
1
(
1
)
=
0.1
,
δ
1
(
2
)
=
0.16
,
δ
1
(
3
)
=
0.28
\delta_{1}(1) = 0.1, \delta_{1}(2)=0.16,\delta_{1}(3) = 0.28
δ1(1)=0.1,δ1(2)=0.16,δ1(3)=0.28
然后,我们看在t=2时:
对每个状态i,求在t=1时状态为j,且观测到为红的前提下,在t=2时状态为i,且观测为白的路径的最大概率(注意,这里的状态就是拿到哪个盒子),记作:
求得:
同理:
δ
2
(
2
)
=
0.0504
,
δ
2
(
3
)
=
0.042
\delta_{2}(2)=0.0504,\delta_{2}(3) = 0.042
δ2(2)=0.0504,δ2(3)=0.042
同样的方式,我们还可以得到:
δ
3
(
1
)
=
0.00756
,
δ
3
(
2
)
=
0.01008
,
δ
3
(
3
)
=
0.0147
\delta_{3}(1) = 0.00756,\delta_{3}(2) = 0.01008,\delta_{3}(3)= 0.0147
δ3(1)=0.00756,δ3(2)=0.01008,δ3(3)=0.0147
根据Viterbi算法的最终取值:
P
∗
=
m
a
x
1
<
=
i
<
=
N
δ
T
(
i
)
P^{*} = max_{1<=i<=N}\delta_{T}(i)
P∗=max1<=i<=NδT(i)
我们知道最大是0.0147,每一步都取最大,不难发现序列是I=(3,2,3)