隐马尔科夫模型(HMM)一前向与后向算法

在该篇文章中讲了隐马尔科夫模型(HMM)一基本模型与三个基本问题 - 知乎专栏隐马尔科夫链的基本概念和最经典的三个问题,这篇文章总结一下隐马尔科夫链(HMM)中的前向与后向算法,首先给出这俩个算法是为了解决HMM的第一个基本问题。先回忆一下第一个问题:
第一个问题是求,给定模型的情况下,求某种观测序列出现的概率。
比如,给定的HMM模型参数已知,求出三天观察是(Dizzy,Cold,Normal)的概率是多少?对应的HMM模型参数已知的意思,就是说的A(trainsition_probability),B(emission_probability),pi矩阵是已经知道的。
相关条件如下图所示:

由上图所示,也就是说,可以写成如下代码:

trainsition_probability = [[0.7,0.3],[0.4,0.6]]     
emission_probability = [[0.5,0.4,0.1],[0.1,0.3,0.6]]    
pi = [0.6,0.4]

在第一个问题中,我们需要求解出三天观察是(Dizzy,Cold,Normal)的概率是多少?
这里为了演示简单,我只求解出俩天观察为(Dizzy,Cold)的概率是多少!

这个问题太好求解了,最暴力的方法就是将路径全部遍历一遍。下面尽可能通俗易懂的说明一下:
首先画出时间序列状态图如下:

下面,我详细走一遍一条路径的暴力算法,这样既可以避开公式的晦涩,也不失正确性。其它路径完全类似
第一天为Healthy的概率为:0.6

在第一天为Healthy的基础上,观察为Dizzy的概率为:P(Dizzy|Healthy)=0.6P(Healthy->Dizzy)=0.6*0.1=0.06

然后求出在第一天为Healthy的基础上,并且第一天表现为Dizzy的前提下,第二天也为Healthy的概率为:
P(Healthy|Healthy,Dizzy)=0.06p(Healthy->Healthy) = P(Dizzy|healthy)07 = 0.06*0.7

上面求完的时候,代表下图中的红线已经转移完了。

好,那么当在前面基础上,第二天观察为Cold的概率为:
P(Cold|(Healthy,Dizzy),(Healthy)) = P(Healthy|Healthy,Dizzy)0.4 = 0.06*0.7*0.4
现在我们已经完成一条路径的完整结果了。

就是在第一天隐含状态为Healthy和第二天隐含状态为Healthy的基础上,观察序列为Dizzy,Cold的概率为
P(Dizzy,Cold|Healthy,Healthy) = 0.06*0.7*0.4=0.0168

那么同理,我们就可以求出其它三条路径。
(1)在第一天隐含状态为Healthy和第二天隐含状态为Fever的基础上,观察序列为Dizzy,Cold的概率
(2)在第一天隐含状态为Fever和第二天隐含状态为Healthy的基础上,观察序列为Dizzy,Cold的概率
(3)在第一天隐含状态为Fever和第二天隐含状态为Fever的基础上,观察序列为Dizzy,Cold的概率

然后最后的第一个问题的结果就是将这四个结果相加起来就可以了。是不是很简单,那么为了还需要前向后向算法来解决这个事呢?

其实这个算法在现实中是不可行的。我给的例子由于是为了讲解容易,状态值和观察值都很小,但是实际中的问题,隐状态的个数是非常大的。

那么我们的计算量是不可以忍受的。
我们可以稍微估计一下,加入状态值是N个,观察值是K个。总共观察长度为T。
那么我们的路径总个数就是N^{T} ,我的天,这个复杂度已经接受不了了,到达了每个隐含状态的时候,还需要算一遍观察值出现的概率(每个隐状态计算一遍到观察值的概率)。又要乘以NT(当然这已经对前面很大复杂度构成不了多少作用了)

所以我们得出结论,暴力法在这里并不实用,于是就引出了前向后向算法。它们都是基于动态规划思想求解。下面分别介绍一下:

前向算法

我们首先定义一下前向概率

定义:给定隐马科夫模型\lambda ,定义到时刻t为止的观测序列为o_{1} ,o_{2},....,o_{t} 且状态为q_{i} 的概率为前向概率,记作

可以递推地求得前向概率\alpha _{t} (i)及观测序列概率p(o|\alpha )

下面,我们可以整理一下前向算法的流程:

输入:隐马尔可夫模型\lambda ,观测序列O

输出:观测序列概率p(o|\alpha )

(1)初值

前向概率的定义中一共限定了两个条件,一是到当前为止的观测序列,另一个是当前的状态。所以初值的计算也有两项(观测和状态),一项是初始状态概率,另一项是发射到当前观测的概率。

(2)递推对t=1,2,3,.....,T-1


每次递推同样由两部分构成,大括号中是当前状态为i且观测序列的前t个符合要求的概率,括号外的是状态i发射观测t+1的概率。

下面稍微解释一下公式:

(3)终止

由于到了时间T,一共有N种状态发射了最后那个观测,所以最终的结果要将这些概率加起来(因为每个隐状态都可能产生我们需要的观测值,所以都要加起来)。

公式可以用下面的转移图表示,假设我要求第二层某个结点的前向概率,等于前一层所有结点到该结点的转移,如下图:


由于每次递推都是在前一次的基础上进行的,所以降低了复杂度(计算只存在于相邻的俩个时间点)。计算如下图所示:


下方标号表示时间节点,每个时间点都有N种状态,所以相邻两个时间之间的递推消耗N^2次计算。

而每次递推都是在前一次的基础上做的,所以只需累加O(T)次,所以总体复杂度是O(T)个N^2,即0(TN^2),这比起我们前面说的暴力法的复杂度已经好了太多了。

后向概率

后向概率与前向概率非常类似,也是基于动态规划的思想,下面介绍一下:

首先给出定义:

定义(后向概率)给定隐马尔可夫模型\lambda ,定义在时刻t状态为q_{i} 的条件下,从t+1到T的部分观测序列为o_{t+1},o_{t+2},...,o_{T} 的概率为后向概率,记作


可以用递推的方法求得后向概率,\beta _{t} (i)及观测序列概率p(o|\lambda ),下面给出后向算法的算法流程。

输入:隐马尔可夫模型\lambda ,观测序列o

输出:观测序列概率p(o|\lambda )

(1)初值


根据定义,从T+1到T的部分观测序列其实不存在,所以硬性规定这个值是1。(这个很重要)

(2)对t = T-1,T-2,...1


a_{ij} 表示状态i转移到j的概率,b_{j} 表示发射o_{t+1} \beta _{t+1} (j)表示j后面的序列对应的后向概率。

上面公式的求解可以用下图转移表示:假设我现在求第一层的某个结点的后向概率,等于第二层的所有结点转移到该层结点,如下:


(3)


最后的求和是因为,在第一个时间点上有N种后向概率都能输出从2到T的观测序列,所以乘上输出O1的概率后求和得到最终结果。(这些我在后面的代码中均有对应的部分)

前向算法和后向算法的python编程实现如下:(代码与上面的伪代码有着非常清楚的对应关系

# -*- coding: UTF-8 -*-
#代码是得到HMM的前向与后向算法,已经验证成功
import numpy as np

def Forward(trainsition_probability,emission_probability,pi,obs_seq):
    """
    :param trainsition_probability:trainsition_probability是状态转移矩阵
    :param emission_probability: emission_probability是发射矩阵
    :param pi: pi是初始状态概率
    :param obs_seq: obs_seq是观察状态序列
    :return: 返回结果
    """
    trainsition_probability = np.array(trainsition_probability)
    emission_probability  = np.array(emission_probability)
    print emission_probability[:,0]
    pi = np.array(pi)
    Row = np.array(trainsition_probability).shape[0]

    F = np.zeros((Row,Col))                      #最后要返回的就是F,就是我们公式中的alpha
    F[:,0] = pi * np.transpose(emission_probability[:,obs_seq[0]])  #这是初始化求第一列,就是初始的概率*各自的发射概率
    print F[:,0]
    for t in range(1,len(obs_seq)):              #这里相当于填矩阵的元素值
        for n in range(Row):                     #n是代表隐藏状态的
            F[n,t] = np.dot(F[:,t-1],trainsition_probability[:,n])*emission_probability[n,obs_seq[t]]   #对应于公式,前面是对应相乘


    return F

def Backward(trainsition_probability,emission_probability,pi,obs_seq):
    """
    :param trainsition_probability:trainsition_probability是状态转移矩阵
    :param emission_probability: emission_probability是发射矩阵
    :param pi: pi是初始状态概率
    :param obs_seq: obs_seq是观察状态序列
    :return: 返回结果
    """
    trainsition_probability = np.array(trainsition_probability)
    emission_probability = np.array(emission_probability)
    pi = np.array(pi)                 #要进行矩阵运算,先变为array类型

    Row = trainsition_probability.shape[0]
    Col = len(obs_seq)
    F = np.zeros((Row,Col))
    F[:,(Col-1):] = 1                  #最后的每一个元素赋值为1

    for t in reversed(range(Col-1)):
        for n in range(Row):
            F[n,t] = np.sum(F[:,t+1]*trainsition_probability[n,:]*emission_probability[:,obs_seq[t+1]])


    return F

if __name__ == '__main__':
    trainsition_probability = [[0.7,0.3],[0.4,0.6]]
    emission_probability = [[0.5,0.4,0.1],[0.1,0.3,0.6]]
    pi = [0.6,0.4]

    #然后下面先得到前向算法,在A,B,pi参数已知的前提下,求出特定观察序列的概率是多少?
    obs_seq = [0,1]
    Row = np.array(trainsition_probability).shape[0]
    Col = len(obs_seq)

    F = Forward(trainsition_probability,emission_probability,pi,obs_seq)                  #得到前向算法的结果
    F_backward = Backward(trainsition_probability,emission_probability,pi,obs_seq)        #得到后向算法的结果

    res_forward = 0
    for i in range(Row):                         #将最后一列相加就得到了我们最终的结果
        res_forward+=F[i][Col-1]                          #求和于最后一列
    emission_probability = np.array(emission_probability)
    #下面是得到后向算法的结果
    res_backword = 0
    res_backward = np.sum(pi*F_backward[:,0]*emission_probability[:,obs_seq[0]])           #一定要乘以发射的那一列
    print "res_backward = {}".format(res_backward)
    print "res_forward = {}".format(res_forward)

参考:

李航老师-《统计学习方法》

编辑于 2018-06-19

文章被以下专栏收录