第一句子网 > 【机器学习】【隐马尔可夫模型-3】后向算法：算法详解+示例讲解+Python实现

【机器学习】【隐马尔可夫模型-3】后向算法：算法详解+示例讲解+Python实现

时间：2019-07-29 06:17:53

0.前排提示

csdn有些数学公式编辑不出来，所以本博用容易书写的表达式来表示专业数学公式，如：

（1）在本博客中用α<T>(i)来表示

（2）在本博客中用[i=1, N]∑来表示

注：这是为了表示一些无法编辑出来的数学公式而本博自己想出来的表示方法，不是专业使用，在离开此博客时请忘记他们~~

1.隐马尔可夫模型简介

隐马尔可夫模型的详细讲解，请详见：【机器学习】【隐马尔可夫模型-1】基本概念+观测序列的生成算法+示例讲解

2.后向算法简介

2.1后向算法的需求背景

首先简要介绍下后向算法的需求背景。

隐马尔可夫模型有3个问题：概率计算问题，学习问题，预测问题。后向算法就是为了解决概率计算问题的。

概率计算问题描述如下：已知一个隐马尔可夫模型λ=（A，B，λ），已知一个观测序列O=（o1, o2,……，oT）问题：请计算在模型λ下观测序列O出现的概率P(O|λ)

HMM的这个概率计算问题有前向算法和后向算法，本博客介绍后向算法，

前向算法在前面博客已经有介绍，详见：【机器学习】【隐马尔可夫模型-2】前向算法：算法详解+示例讲解+Python实现

到此后向算法的需求背景讲解完了，下面开始详解后向算法。

2.2后向概率

已知隐马尔可夫模型λ=（A, B, π），定义到时刻t(t<T)状态为q的条件下，从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=1,2,...,N

2.3观测序列概率的后向算法

可以递推得后向概率以及观测序列概率P(O|λ)

输入：隐马尔卡夫模型λ=(A, B, π)，观测序列O

输出：观测序列概率P(O | λ)

核心是通过模型λ=(A, B, π)以及时刻t+1的β<t+1>(i)来递推得到前一个时刻t的β<t>(i).直到递推得到β<0>(i)

上面步骤比较容易理解，如果看不懂，可以借助下面的示例来理解就变得很容易啦。

话不多说，直接上例子。

3.观测序列概率的后向算法示例

3.1模型λ

从前面介绍的内容可知，模型λ由三个对象确定表示：λ=(A, B, π)=(状态转移概率分布，观测概率分布，初始状态概率分布)

考虑盒子和球模型λ=(A, B, π)，此模型相关信息我们假设给定为：

1）状态集合Q={盒子1，盒子2，盒子3}，N=3

2）观测集合V={红，白}，M=2

3）状态转移概率分布 A = [0.5，0.2，0.3

0.3，0.5，0.2

0.2，0.3，0.5]

4）观测概率分布B=[0.5，0.5

0.4，0.6

0.7，0.3]

5）初始状态概率分布π=[ [0.2],[0.4],[0.4]] #注意：shape=(3,1)

3.2问题

请用后向算法求观测序列概率：已知观测序列O=(红，白，红)，T=3，请用后向算法求P(O|λ)

3.3简单介绍

在进行求解之前，做下简单介绍：

π表示状态为i的初始概率，i是模型中状态集合I={1， 2，3}中的元素状态值，则i可以取1、2或者3，举例：

π<1>表示初始状态为盒1的概率，即π<1>=0.2π<2>表示初始状态为盒1的概率，即π<2>=0.4π<3>表示初始状态和盒3的概率，即π<3>=0.4

b表示状态为盒i时观测概率分布，平民说法就是从盒i取出球，球颜色的概率分布，举例：

b<1>表示状态为盒1时观测概率分布，即b<1> = [0.5，0.5] b<2>表示状态为盒2时观测概率分布，即b<2> = [0.4，0.6]b<3>表示状态为盒3时观测概率分布，即b<3>=[0.7，0.3]

b(红)表示状态为盒i时观测为红的概率，平民说法就是从盒i中取出红球的概率~，

b(白)表示状态为盒i时观测为白的概率，平民说法就是从盒i中取出白球的概率~，举例：

b<1>(红)=0.5，b<1>(白)=0.5b<2>(红)=0.4，b<2>(白)=0.6b<3>(红)=0.7，b<3>(白)=0.3

o表示观测序列O=(o<1>, o<2>，……，o<T>)中的第i个观测，举例：

示例中给出了观测序列O=(红，白，红)，则o<1>=红，o<2>=白，o<3>=红

a<ij>表示从状态i转移到状态j的概率，举例：

示例中给出了状态转移概率分布 A = [0.5，0.2，0.3 #a<11>=0.5，a<12>=0.2，a<13>=0.3

0.3，0.5，0.2 #a<21>=0.3，a<22>=0.5，a<23>=0.2

0.2，0.3，0.5] #a<31>=0.2，a<32>=0.3，a<33>=0.5

a<11>表示当前盒子1取球后，下次取球的盒子是盒1的概率，a<13>表示当前盒子1取球后，下次取球的盒子是盒3的概率

3.4后向算法求解观测序列的概率

下面开始后向算法的求解过程（从T-1到t=0依次计算后向概率，注：前向算法是从t=0到T-1时刻计算概率~）：

Step1：计算初值

后向概率的初值为：

beta0 = [1, 1, 1] #状态0,1,2的后向概率的初值都为1，时刻T

Step2：递推计算

下面用初始的后向概率beta0 递推计算时刻T-1的后向概率：

时刻t=T-1=2, 状态i=0的后向概率：beta[T-1, 0] = A[0,0] * B[0, O[T-1]] * beta0[0]+A[0,1] * B[1, O[T-1]] * beta0[1]+A[0,2] * B[2, O[T-1]] * beta0[2]=0.5*0.5*1 + 0.2 *0.4*1 +0.3*0.7*1 = 0.25 + 0.08 + 0.21 = 0.54如果用Python思想写，就是：a=A[i0] [ 0.5 0.2 0.3]b=B[:,O[t2]]: [ 0.5 0.4 0.7]backward=beta0: [[ 1. 1. 1.]]beta[T-2, 0] = sum(np.multiply(np.multiply(a, b), backward))=0.5*0.5*1 + 0.2 *0.4*1 +0.3*0.7*1 = 0.25 + 0.08 + 0.21 = 0.54

时刻t=T-1=2, 状态i=1的后向概率：beta[T-1, 1] = A[1,0] * B[0, O[T-1]] * beta0[0]+A[1,1] * B[1, O[T-1]] * beta0[1]+A[1,2] * B[2, O[T-1]] * beta0[2]=0.3*0.5*1 + 0.5*0.4*1 +0.2*0.7*1 = 0.15 + 0.2 + 0.14 = 0.49如果用Python思想写，就是：a=A[i1] [ 0.3 0.5 0.2] b=B[:,O[t2]]: [ 0.5 0.4 0.3]backward=beta0: [[ 1. 1. 1.]]beta[T-2, 1] = sum(np.multiply(np.multiply(a, b), backward))=0.3*0.5*1 + 0.5*0.4*1 +0.2*0.7*1 = 0.15 + 0.2 + 0.14 = 0.49

时刻t=T-1=2, 状态i=2的后向概率：beta[T-1, 2] = A[2,0] * B[0, O[T-1]] * beta0[0]+A[2,1] * B[1, O[T-1]] * beta0[1]+A[2,2] * B[2, O[T-1]] * beta0[2]=0.2*0.5*1 + 0.3*0.4*1 +0.5*0.7*1 = 0.1 + 0.12 + 0.35 = 0.57如果用Python思想写，就是：a=A[i2]=[ 0.2 0.3 0.5] b=B[:,O[t2]]=[ 0.5 0.4 0.3]backward=beta0=[[ 1. 1. 1.]]beta[T-1, 2] = sum(np.multiply(np.multiply(a, b), backward))=0.2*0.5*1 + 0.3*0.4*1 +0.5*0.7*1 = 0.1 + 0.12 + 0.35 = 0.57

则求得时刻t=T-1=2每个状态的后向概率为:

beta[2][0]=0.54是时刻2状态0的后向概率，

beta[2][1]=0.49是时刻2状态1的后向概率，

beta[2][2]=0.57是时刻2状态2的后向概率

则时刻2每个状态的后向概率为：beta[2] = [0.54, 0.49, 0.57]，但是时刻2的后向概率=sum(beta[2]) = 1.6

下面用时刻2的后向概率递推计算时刻1的后向概率：

时刻t=T-2=1, 状态i=0的后向概率：beta[T-2, 0] = A[0,0] * B[0, O[T-2]] * beta[T-1, 0]+A[0,1] * B[1, O[T-2]] * beta[T-1, 1]+A[0,2] * B[2, O[T-2]] * beta[T-1, 2]=0.5*0.5*0.54 + 0.2*0.3*0.49 +0.3*0.6*0.57 = 0.135 +0.0588 +0.0513 = 0.2451如果用Python思想写，就是：a=A[i0] [ 0.5 0.2 0.3]b=B[:,O[t2]]: [ 0.5 0.3 0.6]backward=beta0: [[ 0.54 0.49 0.57]]beta[T-2, 0] = sum(np.multiply(np.multiply(a, b), backward))=0.2451

时刻t=T-2=1, 状态i=1的后向概率：beta[T-2, 1] = A[1,0] * B[0, O[T-2]] * beta[T-1, 0]+A[1,1] * B[1, O[T-2]] * beta[T-1, 1]+A[1,2] * B[2, O[T-2]] * beta[T-1, 2]=0.3*0.5*0.54 + 0.5*0.3*0.49 +0.2*0.6*0.57 = 0.2622如果用Python思想写，就是：a=A[i1] [ 0.3 0.5 0.2] b=B[:,O[t2]]: [ 0.5 0.3 0.6]backward=beta0: [[ 0.54 0.49 0.57]]beta[T-2, 1] = sum(np.multiply(np.multiply(a, b), backward))=0.2622

时刻t=T-2=1, 状态i=2的后向概率：beta[T-1, 2] = A[2,0] * B[0, O[T-2]] * beta[T-1, 0]+A[2,1] * B[1, O[T-2]] * beta[T-1, 1]+A[2,2] * B[2, O[T-2]] * beta[T-1, 2]=0.2*0.5*0.54 + 0.3*0.3*0.49 +0.5*0.6*0.57 = 0.2277如果用Python思想写，就是：a=A[i2]=[ 0.2 0.3 0.5] b=B[:,O[t2]]=[ 0.5 0.3 0.6]backward=beta0=[[ 0.54 0.49 0.57]]beta[T-2, 2] = sum(np.multiply(np.multiply(a, b), backward))= 0.2277

则求得时刻t=T-2=1每个状态的后向概率为:

beta[1][0]=0.2451是时刻1状态0的后向概率，

beta[1][1]=0.2622是时刻1状态1的后向概率，

beta[1][2]=0.2277是时刻1状态2的后向概率

则时刻1每个状态的后向概率为：beta[1] = [0.2451, 0.2622, 0.2277]，但是时刻1的后向概率=sum(beta[1]) = 0.735

(3)终止

下面用时刻1的后向概率求时刻0的后向概率：

时刻t=T-3=0, 状态i=0的后向概率：beta[T-3, 0] = pai[0] * B[0, O[T-3]] * beta[T-2, 0] = 0.2 * 0.5 *0.2451 = 0.02451

时刻t=T-3=0, 状态i=1的后向概率：beta[T-3, 1] = pai[1] * B[1, O[T-3]] * beta[T-2, 1] = 0.4 * 0.4 * 0.2622 = 0.041952

时刻t=T-3=0, 状态i=2的后向概率：beta[T-3, 2] = pai[2] * B[2, O[T-3]] * beta[T-2, 2] = 0.4 * 0.7 * 0.2277 = 0..062756

则求得时刻t=T-3=0每个状态的后向概率为:

beta[0][0]=0.02451是时刻0状态0的后向概率，

beta[0][1]=0.041952是时刻0状态1的后向概率，

beta[0][2]=0.062756是时刻0状态2的后向概率

则时刻0每个状态的后向概率为：beta[0] = [0.02451, 0.041952, 0.062756]，时刻0的后向概率=sum(beta[0]) =0.130218

而时刻0的后向概率，就是给定的观测序列O=(红，白，红)出现的概率~

4.Python实现后向算法

4.1后向算法的python代码

完全人肉，代码如下，为便于理解后向思想而加注释的代码，如果想看精简的后向算法代码可知直接看5章节。

# -*- coding: utf-8 -*-"""@author: 蔚蓝的天空tomAim:实现HMM的后向算法(Backward Algorithm)Aim:已知一个马尔可夫模型λ=(A, B, π)，给定一个观测序列O，请用后向算法求此观测序列出现的概率"""import numpy as npclass CHMMBackward(object):'''实现隐马尔科夫模型的后向算法(Backward Algorithm)'''def __init__(self, A, B, pai):self.A = Aself.B = Bself.pai = paidef backward(self, A, B, pai, O):N = np.shape(A)[0] #隐马尔科夫模型的状态个数T = np.shape(O)[0] #给定的状态序列的总长度(即时刻t的个数)#(1)计算初值beta0 = np.ones((1,N)) #后向概率的初值#(2)递推计算每个时刻的序列观测概率, beta_t, t=[T-1, T-2, ..., 0]beta = np.zeros((T,N)) #最终shape=(T,N)#[[0.54,0.49, 0.57],#t=T-1=2, 记作beta_2=[beta_2_0, beta_2_1, beta_2_2]# [0.2451, 0.2604, 0.2277], #t=T-2=1, 记作beta_1=[beta_1_0, beta_1_1, beta_1_2]# [0.02451, 0.0417, 0.0638]] #t=t-3=0, 记作beta_0=[beta_0_0, beta_0_1, beta_0_2]t_range = -1 * np.array(sorted(-1 * np.arange(T))) #[2 1 0]for t in t_range: #t=[T-1, T-2, ..., 0]print('\n==========================time t=%d'%t)beta_t = np.zeros((1, N)) #递推计算每个时刻的序列观测概率beta_tif 0 == t:#终止递推，计算观测序列出现的概率a, b, backward = pai.T, B[:, O[t]], beta[t+1]print('a=π', a, 'b=B[:,O[t0]]:'%t, b, 'backward=beta[t%d+1]:'%t, backward)beta_t = np.multiply(a, b)beta_t = np.multiply(beta_t, backward)else: #递推计算时刻T-1,T-2,...1观测序列出现的概率for i in range(N):print('\n******state i=%d'%i)a,b,backward = A[i],B[:, O[t]],[]if T-1 == t:backward = beta0print('a=A[i%d]'%i, a, 'b=B[:,O[t%d]]:'%t, b, 'backward=beta0:'%t, backward)else:backward = beta[t+1]print('a=A[i%d]'%i, a, 'b=B[:,O[t%d]]:'%t, b, 'backward=beta[t%d+1]:'%t, backward)tmp = np.multiply(a, b)tmp = np.multiply(tmp, backward)beta_t_i = tmp.sum()print('beta_t%d_i%d math:'%(t,i), tmp, 'beta_t%d_i%d value:'%(t,i), beta_t_i)beta_t[0,i] = beta_t_iprint('beta_t%d after insert i%d:'%(t,i), beta_t)beta[t] = beta_tprint('beta_t:',beta_t)print('beta:\n', beta)#后向算法计算出来的O出现的概率是t=0时刻每个状态i对应的序列状态出现的概率总和prob = beta[0].sum()return probdef Prob(self, O):return self.backward(self.A, self.B, self.pai, O)def CHMMBackward_manual():#隐马尔可夫模型λ=(A, B, pai)#A是状态转移概率分布，状态集合Q的大小N=np.shape(A)[0]#从下给定A可知：Q={盒1, 盒2, 盒3}, N=3A = np.array([[0.5, 0.2, 0.3],[0.3, 0.5, 0.2],[0.2, 0.3, 0.5]])#B是观测概率分布，观测集合V的大小T=np.shape(B)[1]#从下面给定的B可知：V={红，白}，T=2B = np.array([[0.5, 0.5],[0.4, 0.6],[0.7, 0.3]])#pai是初始状态概率分布，初始状态个数=np.shape(pai)[0]pai = np.array([[0.2],[0.4],[0.4]])O = [0, 1, 0] #0表示红色，1表示白，就是(红，白，红)观测序列hmm = CHMMBackward(A, B, pai)prob = hmm.Prob(O)print('\n通过HMM的后向算法计算得到：观测序列',O,'出现的概率为:', prob)returnif __name__=='__main__':CHMMBackward_manual()

4.2运行结果

runfile('C:/Users/l13277/test_multiply.py', wdir='C:/Users/l13277')==========================time t=2******state i=0a=A[i0] [ 0.5 0.2 0.3] b=B[:,O[t2]]: [ 0.5 0.4 0.7] backward=beta0: [[ 1. 1. 1.]]beta_t2_i0 math: [[ 0.25 0.08 0.21]] beta_t2_i0 value: 0.54beta_t2 after insert i0: [[ 0.54 0. 0. ]]******state i=1a=A[i1] [ 0.3 0.5 0.2] b=B[:,O[t2]]: [ 0.5 0.4 0.7] backward=beta0: [[ 1. 1. 1.]]beta_t2_i1 math: [[ 0.15 0.2 0.14]] beta_t2_i1 value: 0.49beta_t2 after insert i1: [[ 0.54 0.49 0. ]]******state i=2a=A[i2] [ 0.2 0.3 0.5] b=B[:,O[t2]]: [ 0.5 0.4 0.7] backward=beta0: [[ 1. 1. 1.]]beta_t2_i2 math: [[ 0.1 0.12 0.35]] beta_t2_i2 value: 0.57beta_t2 after insert i2: [[ 0.54 0.49 0.57]]beta_t: [[ 0.54 0.49 0.57]]beta:[[ 0. 0. 0. ][ 0. 0. 0. ][ 0.54 0.49 0.57]]==========================time t=1******state i=0a=A[i0] [ 0.5 0.2 0.3] b=B[:,O[t1]]: [ 0.5 0.6 0.3] backward=beta[t1+1]: [ 0.54 0.49 0.57]beta_t1_i0 math: [ 0.135 0.0588 0.0513] beta_t1_i0 value: 0.2451beta_t1 after insert i0: [[ 0.2451 0.0. ]]******state i=1a=A[i1] [ 0.3 0.5 0.2] b=B[:,O[t1]]: [ 0.5 0.6 0.3] backward=beta[t1+1]: [ 0.54 0.49 0.57]beta_t1_i1 math: [ 0.081 0.147 0.0342] beta_t1_i1 value: 0.2622beta_t1 after insert i1: [[ 0.2451 0.2622 0. ]]******state i=2a=A[i2] [ 0.2 0.3 0.5] b=B[:,O[t1]]: [ 0.5 0.6 0.3] backward=beta[t1+1]: [ 0.54 0.49 0.57]beta_t1_i2 math: [ 0.054 0.0882 0.0855] beta_t1_i2 value: 0.2277beta_t1 after insert i2: [[ 0.2451 0.2622 0.2277]]beta_t: [[ 0.2451 0.2622 0.2277]]beta:[[ 0.0.0. ][ 0.2451 0.2622 0.2277][ 0.54 0.49 0.57 ]]==========================time t=0a=π [[ 0.2 0.4 0.4]] b=B[:,O[t0]]: [ 0.5 0.4 0.7] backward=beta[t0+1]: [ 0.2451 0.2622 0.2277]beta_t: [[ 0.02451 0.041952 0.063756]]beta:[[ 0.02451 0.041952 0.063756][ 0.2451 0.2622 0.2277 ][ 0.540.490.57 ]]通过HMM的后向算法计算得到：观测序列 [0, 1, 0] 出现的概率为: 0.130218