隐马尔可夫模型(HMM)——从理论证明、算法实现到实际应用

目录

历史简述

建模

HMM一般表述

观测序列生成机理

经典问题一之前向算法

经典问题二

经典问题三之viterbi算法

viterbi算法求解案例

viterbi算法代码实现

HMM模型应用


历史简述

        对于马尔可夫链,比较多的说法是:由俄国数学家安德雷·马尔可夫(Андрей Андреевич Марков)在1906-1907年间发表的一篇研究中而来,研究中为了证明随机变量间的独立性不是弱大数定律(weak law of large numbers)和中心极限定理(central limit theorem)成立的必要条件,构造了一个按条件概率相互依赖的随机过程,并证明其在一定条件下收敛于一组向量。在这个研究被提出之后,真的是一发不可收拾,后人在此基础上相继提出了各种模型,比如保罗·埃伦费斯特(Paul Ehrenfest)和Tatiana Afanasyeva在1907年使用马尔可夫链建立了Ehrenfest扩散模型(Ehrenfest model of diffusion) 。1912年亨利·庞加莱(Jules Henri Poincaré)研究了有限群上的马尔可夫链并得到了庞加莱不等式(Poincaré inequality)  。1931年,安德雷·柯尔莫哥洛夫(Андрей Николаевич Колмогоров)在对扩散问题的研究中将马尔可夫链推广至连续指数集得到了连续时间马尔可夫链,并推出了其联合分布函数的计算公式  。独立于柯尔莫哥洛夫,1926年,Sydney Chapman在研究布朗运动时也得到了该计算公式,即后来的Chapman-Kolmogorov等式  。二十世纪50年代,前苏联数学家Eugene Borisovich Dynkin完善了柯尔莫哥洛夫的理论并通过Dynkin公式(Dynkin formula)将平稳马尔可夫过程与鞅过程(martingale process)相联系。

        不幸地是,今天,马尔可夫模型被应用在了机器学习中,困扰着我…… ,对于历史真实性,不做考证,随它遗留在历史角落中……


建模

        在机器学习中,马尔可夫链(Markov chain)得到了极大的应用。它究竟是什么?又能解决什么样问题呢?先从一个最简单的生活案例谈起(以下简称“案例“)。

        小宝有3天假期,每天可以从“郊游”、“逛街”和“打游戏”中选择一种活动。如果是晴天,小宝更可能选择去郊游,而如果下雨,小宝就可能会待在家打游戏了。小宝知道如下信息:

(1)当天的天气只与前一天的天气有关;

(2)第一天下雨的概率为0.6,晴天的概率为0.4;

(3)天气变化概率如下:如果前一天下雨,那么第二天仍然下雨的概率为0.7,而转为晴天的概率为0.3;但是如果前一天是晴天,则第二天有0.6的概率是晴天,而有0.4的概率将会下雨;

小宝根据天气情况如下选择活动:

(a)如果下雨,3个活动的被选择的概率分别为:0.1,  0.4,  0.5;

(b)如果是晴天,3个活动的选择概率分别是:0.6,  0.3,  0.1;

现在如果已经知道了小宝3天假期安排的活动按顺序分别是:郊游、逛街、游戏,我们会面临如下三个问题:

(i)小宝这样安排活动的可能性多大?

(ii)上面(1)~(3)假设的这些概率值是不是导致小宝最后选择这3种活动的最佳概率呢?也就是说有没有其他的这些概率值更可能导致小宝选择这3种活动?

(iii)这3天的天气情况分别是?

 

        下面以案例为基础进行建模来解决这三个问题。这个模型就是马尔可夫模型(Hidden Markov Model,HMM)。下面开始建模。

(1)天气——隐状态

        使用y_1,y_2,...,y_T表示连续T天的天气情况,每个y_t可以看成一个随机变量,可能的取值为S=\{s_1,s_2,...,s_M\},称为隐状态空间,因为通常天气情况是无法观察的。对于案例来说,天气情况序列为y_1,y_2,y_3,状态空间为S={ 下雨, 晴天},如果第一天是晴天,则y_1=晴天。

(2)活动——观测状态

        使用x_1,x_2,...,x_T表示连续T天小宝选择的活动,每个x_t也是一个随机变量,可能的取值为O=\{o_1,o_2,...,o_N\},称为观测状态空间。对于案例来说,活动序列为x_1,x_2,x_3,观测状态空间为O={ 郊游,逛街,游戏 }。

(3)齐次马尔科夫链假设

       天气变化的概率可以使用条件概率P(y_t|y_1,y_2,...,y_{t-1})表示,它表示在第1天到第t-1天的天气情况之后,第t天出现何种天气的概率。根据案例中的(1),当天的天气只与前一天的天气有关,我们有

P(y_t|y_1,y_2,...,y_{t-1})=P(y_t|y_{t-1})

这就是齐次马尔科夫链假设。根据隐状态空间的各种状态的转移的可能,构成如下M\times M阶矩阵:

Q=\begin{bmatrix} q_{11} &q_{12} &... &q_{1M} \\ q_{21}&q_{22} &... &q_{2M} \\ \vdots& \vdots &\ddots & \vdots\\ q_{M1} &q_{M2} &... &q_{MM} \end{bmatrix}

q_{ij}表示由隐状态s_i转移到s_j的概率,Q称之为隐状态转移矩阵。对案例来说,转移矩阵Q如下:

Q=\begin{bmatrix} 0.7&0.3 \\ 0.4&0.6 \end{bmatrix}   

注意到隐状态转移矩阵和时刻t无关。可以使用如下图来表示:

(4)观测独立性假设

        即任意时刻t的观测状态x_t只仅仅依赖于当前时刻的隐藏状态y_t,这是为了简化模型。如果在时刻t的隐藏状态y_t=s_i, 而对应的观察状态为x_t=o_j, 则该时刻在隐藏状态y_t下观察值到x_t的的概率记为g_{ij},如下:

P(x_t=o_j|y_t=s_i)

考虑所有可能的情况,构成一个M\times N阶矩阵G,如下:

G=\begin{bmatrix} g_{11} &g_{12} &... &g_{1N} \\ g_{21}&g_{22} &... &g_{2N} \\ \vdots& \vdots &\ddots & \vdots\\ g_{M1} &g_{M2} &... &g_{MN} \end{bmatrix}

g_{ij}表示在状态为s_i的情况下,观测到o_j的概率,G称之为观测状态矩阵。对于案例来说,G如下:

G=\begin{bmatrix} 0.1 &0.4 &0.5 \\ 0.6&0.3 & 0.1 \end{bmatrix}.

注意到观测状态转移矩阵和t也无关。观测状态矩阵可以使用如下图表示:

(5)初始隐状态分布

        案例中的信息(3)是天气的初始状态,也就是y_1,概率分布F(y_1)如下:

y_1下雨晴天
P(y_1)0.6 0.4

称之为隐状态初始分布。

        有了上面这些之后,马尔科夫模型就可以如下形式地表示:

\lambda=(F,Q,G)

其中F是隐状态初始分布,Q是隐状态转移矩阵,G是观测状态矩阵。

案例的马尔可夫模型可以直观地用下面这个图来表示:

图中表示三天的隐状态和观察状态之间的关系。箭头表示依赖关系。熟悉概率图的知道,这就是HMM的概率图。现在对案例中三个问题进行建模,分别如下:

(1) 估计观察序列X概率

        即给定模型\lambda=(F,Q,G)和观测序列X=\{​{x_1,x_2,...,x_T}\},计算在模型λ下观测序列X出现的概率P(X|λ)。

(2)模型参数的学习问题

        即给定观测序列X=\{​{x_1,x_2,...,x_T}\},在满足马尔科夫条件和独立性假设前提下,求解模型\lambda=(F,Q,G),使该模型下观测到序列X的条件概率P(X|λ)最大,即:\hat{\lambda}=argmax_{\lambda }P(X|\lambda)

(3)预测问题,也称为解码问题

        即给定模型\lambda=(F,Q,G)和观测序列X=\{​{x_1,x_2,...,x_T}\},求给定观测序列条件下,最可能出现的对应的状态序列Y=\{​{y_1,y_2,...,y_T}\}

这三个问题是HMM模型的三个经典问题。好在我们都可以解决。


HMM一般表述

        下面我们给出HMM模型的一般表述。首先我们假设S是所有可能的隐藏状态的集合,O是所有可能的观测状态的集合,即:

S=\{s_1,s_2,...,s_M\}O=\{o_1,o_2,...,o_N\}

其中,M是隐状态的个数,N是观测状态的个数。

  对于一个长度为T的序列,Y是对应的状态序列, X是对应的观测序列,即:

X=\{x_1,x_2,...,x_T\}Y=\{y_1,y_2,...,y_T\}

其中x_t\in Oy_t\in S

        我们对这里提到的序列做些解释,在t时刻和t+1时刻的状态序列和观测序列分别是:

t 时刻:Y=\{​{y_1,y_2,...,y_t}\},X=\{​{x_1,x_2,...,x_t}\}

t+1时刻:Y=\{​{y_1,y_2,...,y_{t+1}\},X=\{​{x_1,x_2,...,x_{t+1}\}

它们之间有一个关系,也就是在t+1时刻只是在t时刻的基础上增加了状态y_{t+1}和对应的观测值x_{t+1}。因为只有经历了t时刻才能到达t+1时刻,然后才出现状态y_{t+1}和对应的观察值x_{t+1}。理解这一点对后面的理论证明会很有帮助。


观测序列生成机理

        观测序列生成机理指的是HMM在已知模型\lambda=(F,Q,G)后,是如何一步步生成观测序列X=\{​{x_1,x_2,...,x_T}\}的。它有助于我们理解HMM为什么可以表示成\lambda=(F,Q,G)形式。以上面的HMM模型表示图为例:

假设现在你站在y_1的位置上,因为隐状态的初始分布式是已知的,而且隐状态到观测状态的概率也是已知的,所以此时你可沿着箭头方向向下走一步生成观测值x_1,此时序列的第一个值已经生成,因为y_1y_2的状态转移概率是已知的,此时可以得到y_2的隐状态的概率,也就是你可以从y_1走到y_2,达到了y_2之后,那么现在从y_2走到x_2就与从y_1走到x_1是类似的了,此时会生成第2个观测值,依次走下去,就可以生成整个观测序列X=\{​{x_1,x_2,...,x_T}\}了。下面使用案例按照这个思路具体来生成。

模型\lambda=(F,Q,G)对应的隐状态初始分布,隐状态转移矩阵,观测状态矩阵分别如下:

y_1下雨晴天
P(y_1)0.6 0.4

Q=\begin{bmatrix} 0.7&0.3 \\ 0.4&0.6 \end{bmatrix}

G=\begin{bmatrix} 0.1 &0.4 &0.5 \\ 0.6&0.3 & 0.1 \end{bmatrix}

小宝应该如下做出选择:

(1)第1天

根据如下公式来计算

P(x_1)=\sum_{y_1}P(x_1,y_1)=\sum_{y_1}P(x_1|y_1)P(y_1)

最后一项中的项都是已知的,P(x_1|y_1),P(y_1)分别是观测状态概率和隐状态初始分布,如下:

P(x1=郊游)=P(郊游|晴天)P(晴天)+P(郊游|下雨)P(下雨)=0.6*0.4+0.1*0.6=0.3

同理,

P(x1=逛街)=0.3*0.4+0.4+0.6=0.36,

P(x1=游戏)=0.1*0.4+0.5+0.6=0.34

因为P(x1=逛街)概率最大,所以小宝第一天很可能会选择逛街。

(2)第2天

根据如下两步计算:

P(y_2)=\sum_{y_1}P(y_1,y_2)=\sum_{y_1}P(y_2|y_1)P(y_1)

P(x_2)=\sum_{y_2}P(x_2,y_2)=\sum_{y_2}P(x_2|y_2)P(y_2)

P(y2=晴天)=P(y2=晴天|y1=晴天)P(y1=晴天)+P(y2=晴天|y1=雨天)P(y1=雨天)

同理计算出P(y2=下雨)。

P(x2=郊游)=P(x2=郊游|y2=晴天)P(y2=晴天)+P(x2=郊游|y2=下雨)P(y2=下雨)

同理得到P(x2=逛街)和P(x2=游戏),最后根据计算结果,做出活动选择。

(3)第3天

计算与第二天相同。

        以上HMM观测序列的生成思路,可以一般地如下叙述:

(1)生成x1

由初始分布F和观测矩阵G,得到第一个观测值x1,计算公式如下:

P(x_1)=\sum_{y_1}P(x_1,y_1)=\sum_{y_1}P(x_1|y_1)P(y_1)

(2)生成x_tt = 2,...,T

分为两步:

a、由转移矩阵Q和y_{t-1}的状态概率分布Q生成y_t的状态概率分布;

P(y_t)=\sum_{y_{t-1}}P(y_{t-1},y_t)=\sum_{y_{t-1}}P(y_t|y_{t-1})P(y_{t-1})

b、将y_i的状态概率分布看成是初始分布,重复步骤(1)生成xi;

P(x_t)=\sum_{y_t}P(x_t,y_t)=\sum_{y_t}P(x_t|y_t)P(y_t)


        接下来,我们来解决HMM的三大经典问题。

经典问题一之前向算法

问题:给定模型\lambda=(F,Q,G)和观测序列X=\{​{x_1,x_2,...,x_T}\},求P(X|\lambda )

:为了方便起见,将P(X|\lambda )简记为P(X),在本经典问题求解过程中关于λ的条件概率也都如此简记。

        对于任意一个状态序列Y=\{​{y_1,y_2,...,y_T}\},可能出现的状态序列总共有M^T个。这里先假设Y=\{​{y_1,y_2,...,y_T}\}是其中一个。那么

P(Y)=P(y_1,y_2,...,y_T)=P(y_T|y_1,..,y_{T-1})P(y_1,..,y_{T-1})=P(y_T|y_{T-1})P(y_1,..,y_{T-1})

上式最后一步是由观测独立性假设得到,即y_T只与y_1,..,y_{T-1}中的y_{T-1}有关。同理将P(y_1,..,y_{T-1})递推下去,最终得到

P(Y)=P(y_T|y_{T-1})P(y_{T-1}|y_{T-2})...P(y_2|y_1)P(y1)

由隐状态转移矩阵和隐状态初始分布可以计算得到P(Y)

我们再来计算在隐状态序列Y下生成观测序列X的条件概率P(X|Y)

P(X|Y)\\ =P(x_1,x_2,...,x_T|y_1,y_2,...,y_T)\\ =P(x_1|y_1,y_2,...,y_T)P(x_2|y_1,y_2,...,y_T)...P(x_T|y_1,y_2,...,y_T)\\ =P(x_1|y_1)P(x_2|y_2)...P(x_T|y_T)

第2步由观测独立性假设,即观测值x_t只与当前的状态y_t有关。由观测矩阵可以得到P(X|Y)

有了P(Y)P(X|Y),我们就能计算X和Y的联合概率P(X,Y)了,如下: 

P(X,Y)=P(X|Y)P(Y)

那么最终P(X)就可以如下计算:

P(X)=\sum_{Y}P(X,Y)

等式右边表示在各种可能的状态下(N^T个),所以

P(X)=\sum_{Y}P(X,Y)=\sum_{Y}P(X|Y)P(Y)\\ =\sum_{Y}P(x_1|y_1)P(x_2|y_2)...P(x_T|y_T)\cdot P(y_T|y_{T-1})P(y_{T-1}|y_{T-2})...P(y_2|y_1)P(y1)

到此就完成了P(X)的计算。

        下面对以上方法计算的复杂度做个简单说明。根据P(X)的计算公式,完成一个观测值X=\{​{x_1,x_2,...,x_T}\}的概率P(X)计算,需要计算2TM^T个运算(计算每一个Y有2T个乘积,总共有M^T的Y),使用大O表示就是复O(TM^T)。这种暴力计算手段,计算复杂度随着T和N的增长成指数级增长,实际中根本无法计算。实际上,我们在<观测序列生成机理>一节使用的就是这种暴力方法。下面介绍另一种更简单的计算方法,称之为“前向算法”。

        假设在时刻t,状态为y_t=s_i,观察值为X_t=\{x_1,x_2,...,x_t\}的概率记为P(X_t,y_t),称之为前向概率。现在假设我们已经找到了t时刻,所有可能状态下,观察值为X的前向概率,现在我们使用递推法计算出在t+1时刻各个状态的前向概率P(X_{t+1},y_{t+1}),其中X_{t+1}=\{x_1,x_2,...,x_t,x_{t+1}\}即:

P(X_{t+1},y_t,y_{t+1})表示t时刻状态为y_j,且t+1时刻状态为y_i,观察序列为X的概率,所以

P(X_{t+1},y_{t+1})

=\sum_{y_t}P(X_{t+1},y_t,y_{t+1})

=\sum_{y_t}P(X_{t+1}|y_t,y_{t+1})P(y_t,y_{t+1})

=\sum_{y_t}P(X_{t},x_{t+1}|y_t,y_{t+1})P(y_t,y_{t+1})

=\sum_{y_t}P(X_{t}|y_t,y_{t+1})P(x_{t+1}|y_t,y_{t+1})P(y_t,y_{t+1})       (由观测独立性假设得到)

=\sum_{y_t}P(X_{t}|y_t)P(x_{t+1}|y_{t+1})P(y_t,y_{t+1})                                   

=\sum_{y_t}P(X_{t}|y_t)P(x_{t+1}|y_{t+1})P(y_{t+1}|y_t)P(y_t)

=\sum_{y_t}P(X_{t},y_t)P(x_{t+1}|y_{t+1})P(y_{t+1}|y_t)

=P(x_{t+1}|y_{t+1})\sum_{y_t}P(X_{t},y_t)P(y_{t+1}|y_t)

这样我们就得到了P(X_{t+1},y_{t+1})P(X_t,y_t)的递推关系。另外,递推关系的初始值,也就是t=1时刻的前向概率为:

P(X_1,y_1)=P(x_1,y_1)=P(x_1|y_1)P(y_1)

这样有了前向概率的初始值和递推关系,就可以计算任意时刻t的前向概率P(X_t,y_t)了。

经典问题一中需要求的P(X),其实就是在时刻T,观察到X_T=\{x_1,x_2,...,x_T\}的概率,即:

P(X)=P(X_T)=\sum_{y_T}P(X_T,y_T)       (1)

        下面对该算法的计算复杂度说明,我们先计算从P(X_t,y_t)P(X_{t+1},y_{t+1})的计算量,即:

P(X_{t+1},y_{t+1})=P(x_{t+1}|y_{t+1})\sum_{y_t}P(X_{t},y_t)P(y_{t+1}|y_t)

总共做2M+1个乘积和加法运算,那么从P(X_{1},y_{1})P(X_{T},y_{T})就需要计算(T-1)(2M+1)个乘积和加法运算,

再根据(1)式,总共需要计算(T-1)(2M+1)M,大O表示为O(TM^2)。从复杂度看,比之前的暴力方法提升了很多。实际上,暴力方法是穷举法,前向概率是递推法。

 

经典问题二

略。

 

经典问题三之viterbi算法

问题:给定模型\lambda=(F,Q,G)和观测序列X=\{​{x_1,x_2,...,x_T}\},预测最可能出现的隐状态序列Y=\{​{y_1,y_2,...,y_T}\}。数学表达如下:

\hat{Y}=arg\max_YP(Y|X)

解:因为P(Y|X)=P(X,Y)/P(X)X已经给定,所以P(X)是常数,问题等价于找到Y使得联合概率P(X,Y)最大即可。

下面结合图介绍viterbi算法,如图:

图中每个实心圆表示某个时刻可能出现的状态,每个方块表示对应的观测值。实心圆之间的箭头,表示隐状态转移,比如图中的红色箭头,表示t-1时刻状态Y1,在t时刻转移为状态Y2。

         我们注意到,任意一个状态序列Y=\{​{y_1,y_2,...,y_T}\}正好与图中的一条路径一一对应,这条路径从t=1时刻状态为y_1的实心圆开始,然后沿着箭头连接到t=2时刻状态为y_2的实心圆,一直连接下去,直到t=T时刻状态为y_T的实心圆。对应的这条路径记为p=(y_1,y_2,...,y_T)。所以求解Y,其实就是从图中寻找一条最优路径p使得联合概率P(X,Y)最大。这是viterbi算法思想一。下面我们来寻找这样的路径。

        我们知道从时刻1到达时刻t总共有M^t条不同路径,将它们根据时刻t的不同状态分为M类,也就是在时刻t隐状态为s_1的所有路径分为一类,时刻t到达状态s_2的路径分为一类,等等,记类别分别为C_i(i=1,2,...,M)。现在从所有路径中寻找一条使得P(X,Y)最大的最优路径问题就转化为从M个类中先分别找出一条使得P(X,Y)最大的路径,然后再从这M条路径中找到使P(X,Y)最大的路径。这样我们就将一个问题分解成了M个小问题。这是viterbi算法思想二。现在固定t时刻的状态,假设为s_i,在类别C_i中使P(X,Y)达到最大的路径记为\hat{p}=(\hat{y_1},\hat{y_2},...,\hat{y_{t-1}},s_i),最大概率记为\delta_t(i)。另外,我们还需要记录下此时\hat{y}_{t-1}}的状态,也就是倒数第二个节点的状态,记为\varphi_t(i),即\hat{y}_{t-1}}=\varphi_t(i)。为了书写方便,引用之前的记号

X_t=\{x_1,x_2,...,x_t\}Y_t=\{y_1,y_2,...,y_t\}

根据\delta_t(i)的定义,它可以如下表示:

\delta_t(i)=\max_{Y_{t-1}}P(Y_{t-1},y_t=s_i,X_t)

同理t+1时刻对应的\delta_{t+1}(i),如下表示:

\delta_{t+1}(i)=\max_{Y_{t}}P(Y_{t},y_{t+1}=s_i,X_{t+1})

下面推导\delta_t(i)\delta_{t+1}(i)两者的关系,注意到x_{t+1}x_1,...x_t,y_1,...,y_t无关,y_{t+1}只与y_{t}有关,所以

P(Y_{t},y_{t+1}=s_i,X_{t+1})

=P(Y_{t},y_{t+1}=s_i,X_{t},x_{t+1})

=P(y_{t+1}=s_i,x_{t+1}|X_t,Y_{t})P(X_t,Y_{t})

=P(x_{t+1}|X_t,Y_{t},y_{t+1}=s_i)P(y_{t+1}=s_i|X_t,Y_{t})P(X_t,Y_{t})

=P(x_{t+1}|y_{t+1}=s_i)P(y_{t+1}=s_i|y_{t})P(X_t,Y_{t})   (2)

等式(2)两边先对变量Y_{t-1}求最大值,然后对变量y_t求最大值:

\delta_{t+1}(i)=\max_{Y_{t}}P(Y_{t},y_{t+1}=s_i,X_{t+1})

=\max_{y_{t}}\max_{Y_{t-1}}P(Y_{t},y_{t+1}=s_i,X_{t+1})

=\max_{y_{t}}\max_{Y_{t-1}}P(x_{t+1}|y_{t+1}=s_i)P(y_{t+1}=s_i|y_{t})P(X_t,Y_{t})

=P(x_{t+1}|y_{t+1}=s_i)\max_{y_{t}}\max_{Y_{t-1}}P(y_{t+1}=s_i|y_{t})P(X_t,Y_{t})

=P(x_{t+1}|y_{t+1}=s_i)\max_{y_{t}}[P(y_{t+1}=s_i|y_{t})\max_{Y_{t-1}}P(X_t,Y_{t})]

=P(x_{t+1}|y_{t+1}=s_i)\max_{y_{t}}[P(y_{t+1}=s_i|y_{t})\delta_t(y_t)]                            (3)

最后一项最大值是从下面这些值中选取:

\\P(y_{t+1}=s_i|y_{t}=s_1)\delta_t(1),\\P(y_{t+1}=s_i|y_{t}=s_2)\delta_t(2),\\...,\\P(y_{t+1}=s_i|y_{t}=s_M)\delta_t(M)

这样\delta_t(i)\delta_{t+1}(i)递归就建立起来了。

这个递归说明:如果我们已经找到了从时刻1到时刻t的所有可能状态s_j对应的最大概率\delta_t(j),j=1,2,...,M,这个概率乘以转移概率\\P(y_{t+1}=s_i|y_{t}=s_j),然后求得它们乘积中的最大值,最后再乘以观测概率=P(x_{t+1}|y_{t+1}=s_i),也就找到了从时刻1到时刻t+1状态为s_i的最大概率\delta_{t+1}(i)。可以用如下图表示:

最优路径就是概率乘积值最大的那条路径。假设它们中的最大者是\delta_{t+1}(\hat{y}_{t+1}),那么现在我们至少能够确定最优的路径的最后两个节点,分别是\hat{y}_t=\varphi (\hat{y}_{t+1})\hat{y}_{t+1}。最优路径中前面这些节点如何找到呢?事实上,上面的递推公式(3)隐含了一个性质:

\delta_{t+1}(i)的最优路径p=(\hat{y_1},...,\hat{y}_{t-1},\hat{y_t}=j,\hat{y}_{t+1}=i)的前t个节点组成的路径q=(\hat{y_1},...,\hat{y}_{t-1},\hat{y_t}=j)\delta_{t}(j)对应的最优路径。因为p\delta_{t+1}(i)的最优路径,所以公式(3)简化为

\delta_{t+1}(i)=P(x_{t+1}|y_{t+1}=i)P(y_{t+1}=i|y_t=j)\delta_t(j)

当状态i和j固定后,\delta_{t+1}(i)\delta_{t}(j)成正比,那么q就是\delta_{t}(j)的最优路径,如果不是,假设\delta_{t}(j)的最优路径是m=(\hat{k_1},...,\hat{k}_{t-1},\hat{y_t}=j)(m\neq q),那么路径n=(\hat{k_1},...,\hat{k}_{t-1},\hat{y_t}=j,\hat{y}_{t+1}=i)(n\neq p)将是\delta_{t+1}(i)的最优路径,矛盾。

 

基于上面这些介绍,我们们使用\varphi_{t+1}(i)记录的状态往回追溯就可以找出最优路径。具体如下:

(1)\delta_{t+1}(\hat{y}_{t+1})使得P(Y^{(t+1)},X^{(t+1)})最大,所以\hat{y}_{t+1}就是t+1时刻的状态;

(2)\varphi_{t+1}(i)记录了在t+1时刻时,第t个节点的状态,\varphi_{t+1}(\hat{y}_{t+1})就是这个节点的状态,记作\hat{y}_{t}=\varphi_{t+1}(\hat{y}_{t+1})

(3)\varphi_{t}(\hat{y}_{t})的值记录的就是在t时刻,t-1时刻的状态,记作\hat{y}_{t-1}=\varphi_{t}(\hat{y}_{t})

...

按照递归\hat{y}_{t}=\varphi_{t+1}(\hat{y}_{t+1})可以追溯到t=1时刻的状态,也就能得到

\hat{Y}=\{​{\hat{y_1},\hat{y_2},...,\hat{y}_{t+1}}\}

 

viterbi算法基本思路:

(1)由初始时的\delta_1(i)和递推公式,计算出任意时刻t的\delta_t(i),同时记录前一个节点的状态\varphi_t (i);

(2)找到\delta_t(i)中的最大者,它对应的路径就是我们要找的,此时可以得到最后两个节点状态i和\varphi_t (i)

(3)使用\varphi_t(i)记录的前一个状态,回溯出所有时刻的状态,得到\hat{Y}=\{​{\hat{y_1},\hat{y_2},...,\hat{y}_{t+1}}\}

 

viterbi算法求解案例

下面我们以“小宝出行”案例来熟悉下vertibi算法。

(1)对于时刻t=1,也就是第一天,可能的状态为{下雨,晴天},小宝第一天选择的是郊游,另外,从Start开始只有一条路径到达晴天,所以

\delta_1(晴天)=P(x1=郊游,y1=晴天)=P(x1=郊游|y1=晴天)P(y1=晴天)=0.6*04=0.24

\delta_1(下雨)=P(x1=郊游,y1=下雨)=P(x1=郊游|y1=下雨)P(y1=下雨)=0.1*0.6=0.06

因为第一天的前一天没有状态,所以\varphi_1(i)不存在。

 

(2)对于时刻t=2,也就是第二天,小宝选择的是逛街,根据递推有:

\delta_2(i)=P(x_{2}|y_{2}=i)max_{j\in YY}P(y_{2}=i|y_1=j)\delta_1(j),所以,

\delta_2(晴天)=P(x2=逛街|y2=晴天)*max{P(y2=晴天|y1=晴天)\delta_1(晴天),P(y2=晴天|y1=下雨)\delta_1(下雨)}

            =0.3*max{0.6*0.24,0.3*0.06}=0.3*max{0.144,0.018}=0.0432

\varphi_2(晴天)=晴天

同理,

\delta_2(下雨)=P(x2=逛街|y2=下雨)*max{P(y2=下雨|y1=晴天)\delta_1(晴天),P(y2=下雨|y1=下雨)\delta_1(下雨)}

             =0.4*max{0.4*0.24,0.7*0.06}=0.4*max{0.096,0.042}=0.0384

\varphi_2(下雨)=晴天

 

(3)对于时刻t=3,也就是第三天,小宝选择游戏,根据递推:

\delta_3(i)=P(x_{3}|y_{3}=i)max_{j\in YY}P(y_{3}=i|y_2=j)\delta_2(j),所以,

\delta_3(晴天)=P(x3=游戏|y3=晴天)*max{P(y3=晴天|y2=晴天)\delta_2(晴天),P(y3=晴天|y2=下雨)\delta_2(下雨)}

             =0.1*max{0.6*0.0432,0.3*0.0384}=0.1*max{0.02592,0.01152}=0.002592

\varphi_3(晴天)=晴天

\delta_3(下雨)=P(x3=游戏|y3=下雨)*max{P(y3=下雨|y2=晴天)\delta_2(晴天),P(y3=下雨|y2=下雨)\delta_2(下雨)}

             =0.5*max{0.4*0.0432,0.7*0.0384}=0.5*max{0.01728,0.02688}=0.01344

\varphi_3(下雨)=下雨

 

下面开始回溯:

\delta_3(下雨)>\delta_3(晴天),所以第3天下雨,而\varphi_3(下雨)=下雨,所以第二天也是下雨,根据\varphi_2(下雨)=晴天,所以第一天是晴天。

推算出3天天气是:晴天,下雨,下雨。

 

viterbi算法代码实现

import numpy as np


def viterbi(trans_prob, emit_prob, init_prob, views, states, obs):
    """
    viterbi算法

    :param trans_prob: 状态转移概率矩阵
    :param emit_prob: 观察概率矩阵,也称为发射概率矩阵
    :param init_prob: 初始状态分布
    :param views: 所有可能的观测值集合
    :param states: 所有可能的状态集合
    :param obs: 实际观测值序列
    :return:
    """

    state_num, obs_len = len(states), len(obs)
    delta = np.array([[0] * state_num] * obs_len, dtype=np.float64)
    phi = np.array([[0] * state_num] * obs_len, dtype=np.int64)
    print('state_num=', state_num, 'obs_len=', obs_len)
    print('delta=', delta)
    print('phi=', phi)
    # 初始化
    for i in range(state_num):
        delta[0, i] = init_prob[i] * emit_prob[i][views.index(obs[0])]
        phi[0, i] = 0
    print('初始化后delta=', delta)
    print('初始化后phi=', phi)

    # 递归计算
    for i in range(1, obs_len):
        for j in range(state_num):
            tmp = [delta[i - 1, k] * trans_prob[k][j] for k in range(state_num)]
            delta[i, j] = max(tmp) * emit_prob[j][views.index(obs[i])]
            phi[i, j] = tmp.index(max(tmp))

    # 最终的概率及节点
    max_prob = max(delta[obs_len - 1, :])
    last_state = int(np.argmax(delta[obs_len - 1, :]))

    # 最优路径path
    path = [last_state]
    for i in reversed(range(1, obs_len)):
        end = path[-1]
        path.append(phi[i, end])

    hidden_states = [states[i] for i in reversed(path)]

    return max_prob, hidden_states


def main():
    # 所有可能的状态集合
    states = ('晴天', '下雨')
    # 观测集合
    views = ['郊游', '逛街', '游戏']
    # 转移概率: Q -> Q
    trans_prob = [[0.6, 0.4],
                  [0.3, 0.7]]
    # 观测概率, Q -> V
    emit_prob = [[0.6, 0.3, 0.1],
                 [0.1, 0.4, 0.5]]
    # 初始概率
    init_prob = [0.4, 0.6]
    # 观测序列
    obs = ['郊游', '逛街', '游戏']

    max_prob, hidden_states = viterbi(trans_prob, emit_prob, init_prob, views, states, obs)
    print('最大的概率为: %.5f.' % max_prob)
    print('隐藏序列为:%s.' % hidden_states)


if __name__ == '__main__':
    main()

 

HMM模型应用

中文分词

 

股市分析

 

  • 14
    点赞
  • 37
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 3
    评论
HMM是一种常见的统计模型,用于描述随机生成的序列。Viterbi算法是一种在HMM中进行解码的动态规划算法,用于寻找最可能的隐藏状态序列。下面是HMM的MATLAB实现,包括Viterbi算法。 首先,我们需要定义一个HMM模型。我们假设这个模型有3个隐藏状态和2个可见状态。我们使用矩阵A表示隐藏状态之间的转移概率,矩阵B表示每个隐藏状态生成每个可见状态的概率,向量pi表示初始隐藏状态的概率分布。 ```matlab % 定义HMM模型参数 A = [0.5 0.2 0.3; 0.3 0.5 0.2; 0.2 0.3 0.5]; B = [0.5 0.5; 0.4 0.6; 0.7 0.3]; pi = [0.2 0.4 0.4]; ``` 接下来,我们需要生成一个可见序列。我们使用HMM模型中的随机过程生成一个长度为10的可见序列。 ```matlab % 生成可见序列 T = 10; q = zeros(1, T); o = zeros(1, T); q(1) = randsrc(1, 1, [1:3; pi]); o(1) = randsrc(1, 1, [1:2; B(q(1), :)]); for t = 2:T q(t) = randsrc(1, 1, [1:3; A(q(t-1), :)]); o(t) = randsrc(1, 1, [1:2; B(q(t), :)]); end ``` 接下来,我们使用Viterbi算法解码这个可见序列,得到最可能的隐藏状态序列。我们定义一个矩阵V表示每个时间步的最大概率,以及一个矩阵path表示每个时间步的最大概率对应的前一个状态。 ```matlab % Viterbi算法解码 V = zeros(3, T); path = zeros(3, T); V(:, 1) = pi' .* B(:, o(1)); for t = 2:T for j = 1:3 [V(j, t), path(j, t)] = max(V(:, t-1) .* A(:, j)); V(j, t) = V(j, t) * B(j, o(t)); end end ``` 最后,我们找到最可能的隐藏状态序列。我们首先找到最后一个时间步的最大概率对应的隐藏状态,然后从后往前依次寻找每个时间步的最大概率对应的隐藏状态,最终得到整个隐藏状态序列。 ```matlab % 找到最可能的隐藏状态序列 [~, q(T)] = max(V(:, T)); for t = T-1:-1:1 q(t) = path(q(t+1), t+1); end ``` 完整代码如下: ```matlab % 定义HMM模型参数 A = [0.5 0.2 0.3; 0.3 0.5 0.2; 0.2 0.3 0.5]; B = [0.5 0.5; 0.4 0.6; 0.7 0.3]; pi = [0.2 0.4 0.4]; % 生成可见序列 T = 10; q = zeros(1, T); o = zeros(1, T); q(1) = randsrc(1, 1, [1:3; pi]); o(1) = randsrc(1, 1, [1:2; B(q(1), :)]); for t = 2:T q(t) = randsrc(1, 1, [1:3; A(q(t-1), :)]); o(t) = randsrc(1, 1, [1:2; B(q(t), :)]); end % Viterbi算法解码 V = zeros(3, T); path = zeros(3, T); V(:, 1) = pi' .* B(:, o(1)); for t = 2:T for j = 1:3 [V(j, t), path(j, t)] = max(V(:, t-1) .* A(:, j)); V(j, t) = V(j, t) * B(j, o(t)); end end % 找到最可能的隐藏状态序列 [~, q(T)] = max(V(:, T)); for t = T-1:-1:1 q(t) = path(q(t+1), t+1); end ```

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

leboop-L

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值