统计学习方法笔记-隐马尔可夫模型（内含Python代码实现)

一马尔可夫模型

我们通过一个具体的例子来介绍一下什么是马尔可夫模型

我们假设天气有3种情况，阴天，雨天，晴天，它们之间的转换关系如下:
(稍微解释一下这个图，我们可以这样认为，已知第一天是阴天，那第二天是阴天的概率是0.1，第二天是晴天的概率是0.3，第二天是雨天的概率是0.6)

每一个状态转换到另一个状态或者自身的状态都有一定的概率。

马尔可夫模型就是用来表述上述问题的一个模型。

有这样一个状态链，第一天是阴天，第二天是晴天，第三天是雨天。这样一个状态链就是马尔可夫链。

在上述例子中，如果我们并不知道今天天气属于什么状况，我们只知道今明后三天的水藻的干燥湿润状态，因为水藻的状态与天气有关，我们用水藻状态来推测出这三天的天气，就需要用到隐马尔科夫模型。

下面我们将介绍一下隐马尔可夫模型

二隐马尔可夫模型

隐马尔可夫模型定义

隐马尔可夫模型是关于时序的概率模型，描述由一个隐藏的马尔可夫链随机生成不可观测的状态随机序列，再由各个状态生成一个观测从而产生观测随机序列的过程。隐藏的马尔可夫链随机生成的状态的序列，称为状态序列 ;每个状态生成一个观测，而由此产生的观测的随机序列，称为观测序列。序列的每一个位置又可以看作是一个时刻。

隐马尔可夫模型的组成
(这一块看的比较懵没关系，下面会结合例子解释哦)

设 Q是所有可能的状态集合 Q ={q₁,q₂,…q_N} V是所有可能的观测的集合 V={v₁,v₂,…v_M}

N是可能的状态数，M是可能的观测数

I 是长度为T的状态序列 I={i₁,i₂,…i_T}

O 是状态序列对应的观测序列 O={o₁,o₂,…o_T}

隐马尔可夫模型由以下三部分组成

1 初始状态概率向量 π

π=（π_i)

其中 π_i=P(i₁ = q_i) 表示在时刻 t=1 处于状态q_i的概率

2 状态转移概率矩阵 A

A=[a_ij] _N*N

其中a_ij = P( i_t+1 = q_j | i_t = q_i) 表示在t时刻处于状态q_i 的情况下在t+1时刻转移到状态q_j 的概率。

3 观测概率矩阵 B

B = [ b_j(k) ] _N*M

b_j(k) = P(o_t = v_k| i_t = q_j) 表示在时刻 t 处于状态q_j 的条件下生成观测v_k 的概率。

隐马尔可夫模型一般表示为

我们来看一个具体的例子体会一下这些符号都是什么意思

我们有三种骰子分别是六面体骰子（它的值为1，2，3，4，5，6）四面体骰子(它的值为1，2，3，4) 八面体骰子 (它的值为1，2，3，4，5，6，7，8) 把这三种骰子放到一个盒子中有放回的去抽一个骰子，拿到这个骰子之后，抛起这个骰子记录它的数字。重复这个动作多次，得到一个数字序列。
假设我们做了10次这个动作得到的序列为 1，6，3，5，2，7，3，5，2，4

假设状态转换图是这样的

(解释一下: 假如某一次的骰子是D₄ 那下一次抽到的骰子是D₄的概率是0.1 ，下一次抽到的骰子是D₈的概率是 0.4 下一次抽到的骰子是D₆ 的概率是0.5.）

题目介绍完了，我们来分析一下

在本例中，我们的状态就是不同种类的骰子
所以 Q={D₄,D₆,D₈} N=3

我们的观测集合就是骰子上的数字可能的情况
V={1,2,3,4,5,6,7,8} M=8

我们的观测序列是多次掷骰子得到的一个数字序列 1，6，3，5，2，7，3，5，2，4

初始时刻我们从盒子里选择一个骰子，因为是随机选择，抽到每一个骰子的概率都是 1/3.

则初始状态向量 π =（1/3，1/3，1/3）

状态转移概率矩阵 A (根据上方给出的状态转换图得到)
观测概率矩阵 B (在某一状态下，生成某一观测的概率)

隐马尔可夫模型的两个假设

齐次马尔可夫性假设，即假设隐藏的马尔可夫链在任意时刻 t 的状态只依赖于其前一时刻的状态

观测独立性假设，即假设任意时刻的观测只依赖于该时刻的马尔可夫链的状态

隐马尔可夫模型有三个需要解决的问题分别是概率计算问题，学习问题，预测问题。在下文将一一介绍。

三概率计算问题

概率计算问题指的是给定一个隐马尔可夫模型 λ =（A,B,π) ,和观测序列O = (o₁， o₂， … ， o_T), 在模型 λ 下观测序列O 出现的概率P(OIλ)。

直接计算法

这是一种理论上可行，但是由于计算量太大，导致计算上不可行的方法，了解即可。

前向算法

先介绍一个概念: 前向概率

给定隐马尔可夫模型 λ，定义到时刻 t 部分观测序列为 o₁， o₂， … ， o_t 且状态为 q_i 的概率为前向概率，记作 α_t(i) = P(o₁,o₂,…,o_t, i_t=q_i| λ)

前向算法可以递推地求得前向概率α_t(i) 从而得到观测序列概率 P(OIλ)。

观测序列概率的前向算法如下:

前向算法的一些推导

前向算法计算例子

后向算法

后向概率

给定隐马尔可夫模型 λ，定义在时刻 t 状态为 q_i 的条件下，从 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， λ)

观测序列概率的后向算法

后向算法相关推导

利用前向概率和后向概率的定义可以将观测序列概率 P(OIλ) 统一写成

我们可以根据前向概率和后向概率来表示一些概率和期望值的计算。
根据以上则

在观测O下状态 i 出现的期望值:
在观测 O 下由状态 i 转移的期望值:
在观测O下由状态 i 转移到状态 j 的期望值:

四隐马尔可夫模型的学习问题

假设给定训练数据只包含 S 个长度为 T 的观测序列 {O₁，O₂，…， O_s} 而没有对应的状态序列，隐马尔可夫模型的学习问题的目标是学习隐马尔可夫模型 λ = (A， B， π )的参数。我们将观测序列数据看作观测数据 O，状态序列数据看作不可观测的隐数据 I，那么隐马尔可夫模型事实上是一个含有隐变量的概率模型.

它的参数学习可以由 EM 算法实现，该算法命名为Baum-Welch.

下面给出用EM算法求解隐马尔可夫模型参数的详细推导

参数估计为

五隐马尔可夫模型的预测问题

隐马尔可夫模型的预测问题，也称为解码问题。已知模型 λ = (A， B， π ) 和观测序列O=(o₁,o₂,…o_T),求对给定观测序列条件概率P(I |O) 最大的状态序列I=(i₁,i₂,…i_T). 即给定观测序列，求最有可能的对应的状态序列。

我们用的方法是维特比算法，下面介绍一下维特比算法。

维特比算法是一个基于动态规划思想的方法，它把求解状态序列看作寻找一条最优路径。

维特比算法是这样的

维特比算法示例
该示例代码实现如下

class Viterbi:def __init__(self,o,A,B,PI,index):self.o=oself.A=Aself.B=Bself.index = indexself.N=self.A.shape[0] #状态集合有多少元素self.M = len(self.index)  # 观测集合有多少元素self.T =len(self.o) #观测序列一共有多少值self.PI=PI # 初始状态概率self.delte=[[0]*self.N]*self.Tself.I=[] #得到的状态序列是self.keci=[[0]*self.N]*self.Tdef cal_delte(self):# 在书中时刻t的取值是 1到 T 但是写代码数组是从0 开始的 方便起见 我们讲t也从0开始o1=self.o[0]#第一个观测变量是o1_index=self.index[o1] #第一个观测变量的下标是for i in range(self.N):self.delte[0][i] = self.PI[i]*self.B[i][o1_index]for t in range(1,self.T):#从时刻t=1 开始 到T-1ot=self.o[t]ot_index=self.index[ot]for i in range(self.N):max=0maxj=0for j in range(self.N):a = self.delte[t-1][j] *self.A[j][i]*self.B[i][ot_index]if a>max:max=amaxj=jself.delte[t][i]= maxself.keci[t][i]=maxjdef cal_max_path(self):max=0maxi=0path=[]for i in range(self.N):a=self.delte[self.T-1][i]if a>max:max=amaxi=ipath.append(maxi+1)for t in range(self.T-1,0,-1):maxi=self.keci[t][maxi]path.append(maxi+1)for i in range(len(path)-1,-1,-1):self.I.append(path[i])print(self.I)
A=np.array([[0.5,0.2,0.3],[0.3,0.5,0.2],[0.2,0.3,0.5]])
B=np.array([[0.5,0.5],[0.4,0.6],[0.7,0.3]])
PI=np.array([[0.2],[0.4],[0.4]])
o=['红','白','红']
index={'红':0,'白':1}
hmm=Viterbi(o,A,B,PI,index)
hmm.cal_delte()
hmm.cal_max_path()

六隐马尔可夫模型代码实现

class HMM:def __init__(self,o,status,observe,n):self.o=o #观测数据self.status= status #状态集合self.observe=observe# 观测集合self.N = len(self.status)  # 状态集合有多少元素self.M = len(observe)  # 观测集合有多少元素self.A = [[1 / self.N] * self.N] * self.Nself.B = [[1/self.M]*self.M]* self.Nself.T = len(self.o)  # 观测序列一共有多少值self.PI = [1/self.N]*self.N  # 初始状态概率 N个状态 每个状态的概率是1/Nself.delte = [[0] * self.N] * self.Tself.I = []  # 得到的状态序列是self.psi = [[0] * self.N] * self.Tself.a=self.cal_a()self.b=self.cal_b()self.n=ndef cal_a(self):#计算前向概率o1=self.o[0]o1_index=self.observe[o1]a=[[0]*self.N]*self.Tfor i in range(self.N):a[0][i]=self.PI[i]*self.B[i][o1_index]for t in range(1, self.T):  # 从时刻t=1 开始 到T-1ot=self.o[t]ot_index = self.observe[ot]for i in range(self.N):sum=0for j in range(self.N):sum += a[t-1][j]*self.A[j][i]a[t][i]=sum*self.B[i][ot_index]return adef cal_b(self):#计算后向概率b = [[0] * self.N] * self.Tfor i in range(self.N):b[self.T-1][i] = 1for t in range(self.T-2,-1,-1):ot_add_1 = self.o[t+1]ot_ot_add_1_index = self.observe[ot_add_1]for i in range(self.N):sum=0for j in range(self.N):sum+=self.A[i][j]*self.B[j][ot_ot_add_1_index]*b[t+1][j]b[t][i]=sumreturn bdef cal_gamma(self, t, i):# 计算李航《统计学习方法》 p202 公式10.24sum = 0for j in range(self.N):sum += self.a[t][j] * self.b[t][j]# print(self.a)# print(self.b)return self.a[t][i] * self.b[t][i] / sumdef cal_xi(self, t, i1, j1):# 计算李航《统计学习方法》 p203 公式10.26sum = 0ot_add_1 = self.o[t + 1]ot_ot_add_1_index = self.observe[ot_add_1]for i in range(self.N):for j in range(self.N):sum += self.a[t][i] * self.A[i][j] * self.B[j][ot_ot_add_1_index] * self.b[t + 1][j]p = self.a[t][i1] * self.A[i1][j1] * self.B[j1][ot_ot_add_1_index] * self.b[t + 1][j1]return p / sumdef update_A(self):for i in range(self.N):for j in range(self.N):sum1=0sum2=0for t in range(self.T - 1):sum1+=self.cal_xi(t,i,j)sum2+=self.cal_gamma(t,i)self.A[i][j]=sum1/sum2def update_B(self):for j in range(self.N):for k in range(self.M):sum1=0sum2=0for t in range(self.T):ot=self.o[t]ot_index=self.observe[ot]if ot_index == k:sum1+=self.cal_gamma(t,j)sum2+=self.cal_gamma(t,j)self.B[j][k]=sum1/sum2def update_pi(self):for i in range(self.N):self.PI[i]=self.cal_gamma(0,i)def fit(self):for i in range(self.n):print(i)self.update_A()self.update_B()self.update_pi()self.a = self.cal_a()self.b = self.cal_b()def cal_delte(self):# 在书中时刻t的取值是 1到 T 但是写代码数组是从0 开始的 方便起见 我们讲t也从0开始o1=self.o[0]#第一个观测变量是o1_index=self.observe[o1] #第一个观测变量的下标是for i in range(self.N):self.delte[0][i] = self.PI[i]*self.B[i][o1_index]for t in range(1,self.T):#从时刻t=1 开始 到T-1ot=self.o[t]ot_index=self.observe[ot]for i in range(self.N):max=0maxj=0for j in range(self.N):a = self.delte[t-1][j] *self.A[j][i]*self.B[i][ot_index]if a>max:max=amaxj=jself.delte[t][i]= maxself.psi[t][i]=maxjdef cal_max_path(self):max=0maxi=0path=[]for i in range(self.N):a=self.delte[self.T-1][i]if a>max:max=amaxi=ipath.append(maxi+1)for t in range(self.T-1,0,-1):maxi=self.psi[t][maxi]path.append(maxi+1)for i in range(len(path)-1,-1,-1):self.I.append(path[i])print(self.I)o=['红','白','红']
observe={'红':0,'白':1}
status=[1,2,3]
hmm=HMM(o,status,observe,200)
hmm.fit()
hmm.cal_delte()
hmm.cal_max_path()

参考

李航《统计学习方法》

b站视频机器学习-白板推导系列(十四)-隐马尔可夫模型HMM