马尔科夫链是离散的马尔科夫过程,而马尔科夫过程是一类平稳的随机过程。如果一个时间序列是马尔科夫链,设时间序列在每个时刻的状态有n种状态S={s1,s2,s3,..sn},在m时刻为si,则在m+k时刻的状态sj只与m时刻相关,与m-1,m-2,m-3...,1时刻无关:
上面的条件概率说明了马尔科夫链未来状态(m+k)仅与当前状态有关(时刻m)有关,而平稳过程表明时间序列状态的转移与起始点m也无关,仅与k有关,关于平稳随机过程的定义请参见时间序列分析:AR(p),MA(q),满足平稳性的马尔科夫过程也称齐次马尔科夫过程。常研究马尔科夫链k=1时各个状态转移情况,用aij表示从状态si到下一时刻状态为sj的概率,如果隐藏变量有n种状态,则马尔科夫链的状态转移矩阵大小是。
隐马尔科夫链HMM是含有隐变量马尔科夫过程,关于隐变量在EM算法中曾提过。举一个例子,有三个箱子分别标记为1号箱、2号箱、3号箱,每个箱子各有一定比例的白色球和红色球,现将箱子标记遮住,每次随机选择一个箱子、抽取一个球,得的的一组按次序摆放的白、红球称为观察变量,每次选择的箱子是隐藏变量,隐马尔科夫链暗含了两个随机过程,一个是选择箱子的序列,另一个是观察球序列,箱子之间状态的转移称为转移概率矩阵(A),每个箱子含有白色球、红色球的比例称为发射矩阵(B),另外,初始时选择各个箱子的概率称为初始概率向量(Pi),三个变量构成一个参数集合可用表示。HMM常见的有三类问题,分别是:
1、评估问题:已知λ,求出现某种观察结果的概率。
2、解码问题:已知λ,通过观察结果求出隐藏状态序列。
3、学习问题:λ未知,通过观察结果求解λ。
一、评估问题
1.1 前向算法、后向算法
可以用穷举法计算观察结果的概率P(O|λ),但当观察结果序列很长时穷举法的复杂度很高,前向、后向算法是利用递归方式实现评估问题求解,可大大降低算法运算量,前向算法中首先引入前向概率:
(1.1)
前向概率是联合概率,表示当前面t个观察结果已知时,在t时刻隐藏状态为i观察结果为ot的概率,利用加法原理可知在t+1时刻的前向概率为:
(1.2)
(1.2)中n为隐藏状态的数量,aji是转移概率矩阵中相对应的元素,bi(ot+1)是发射概率矩阵中相对应的元素,如果观察结果的长度为T,很容易得到观察结果概率P(O|λ):
(1.3)
前向算法的初值需要利用初始概率矩阵:
与前向算法相反的是后向算法,后向算法中后向概率的定义为:
(1.4)
后向概率从观察结果后端向前计算,另外观察(1.4)可知后向概率并不知道时刻t的观察结果ot,此时只有t+1,t+2,..T时刻的观察结果是已知的,根据(1.4)可得到形式:
假设t-1时刻隐藏状态为si、t 时刻隐藏状态为sj,可以看成由两个部分组成,第一部分是出现观察结果ot的概率,即为:,另一个是t时刻隐藏状态为sj、观察结果ot+1,ot+2,ot+3,oT的概率:,利用乘法定理可知,此时的概率为:
而t时刻隐藏状态可能性有n个,利用加法原理可得:
(1.5)
由上面推导可递归得到
在t=1时刻,从隐藏状态si得到o1的概率为,由加法原理可知后向算法P(O|λ)为:
(1.6)
后向概率的初始概率为,可以表示为:
实际上oT+1是观察不到的,可以认为在T时刻隐藏状态每种可能性都有,可设,利用上面公式来实现开篇箱子抽取红白球的评估问题:
from hmmlearn import hmm import numpy as np #隐变量,3个箱子 states = ['box1', 'box2', 'box3'] #观察结果状态:红色球、白色球 observations = ["red", "white"] n_states = len(states) start_probability = np.array([0.2, 0.5, 0.3],np.float64) transition_probability = np.array([ [0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5] ],np.float64) emission_probability = np.array([ [0.5, 0.5], [0.4, 0.6], [0.7, 0.3] ],np.float64) class HMM(object): def __init__(self, A, B, pi): ''''' A: 状态转移概率矩阵 B: 输出观察概率矩阵 pi: 初始化状态向量 ''' self.A = np.array(A) self.B = np.array(B) self.pi = np.array(pi) self.N = self.A.shape[0] # 总共状态个数 self.M = self.B.shape[1] # 总共观察值个数 # 前向算法 def forward(self, O ): ''''' T: 观察序列的长度 O: 观察序列 alpha: 前向概率数组 ''' T=len(O) alpha = np.zeros((T, self.N), np.float64) # 初始化 for i in range(self.N): alpha[0, i] = self.pi[i] * self.B[i, O[0]] # 递归 for t in range(T-1): for j in range(self.N): sum = 0.0 for i in range(self.N): sum += alpha[t, i] * self.A[i, j] alpha[t+1, j] = sum * self.B[j, O[t+1]] # 终止 公式(1.3) sum = 0.0 for i in range(self.N): sum += alpha[T-1, i] return sum def back(self, O): ''''' T: 观察序列的长度 len(O) O: 观察序列 beta: 后向概率数组 ''' T=len(O) beta = np.zeros((T, self.N), np.float64) # 初始化 for i in range(self.N): beta[T-1, i] = 1.0 # 递归 for t in range(T-2, -1, -1): # 从T-2开始递减;即T-2, T-3, T-4, ..., 0 for i in range(self.N): sum = 0.0 for j in range(self.N): sum += self.A[i, j] * self.B[j, O[t+1]] * beta[t+1, j] beta[t, i] = sum # 终止 公式(1.6) sum = 0.0 for i in range(self.N): sum += self.pi[i]*self.B[i,O[0]]*beta[0,i] return sum #评估问题 def assessment(): obs = [1,1, 0 ,0 ,1,0,1,0,1,0,0] model = hmm.MultinomialHMM(n_components=n_states, n_iter=100, tol=1e-3) model.startprob_ = start_probability model.transmat_ = transition_probability model.emissionprob_ = emission_probability score = model.score(np.array([obs]).T) #利用hmm库计算观察结果概率 print("MultinomialHMM P(O|λ):%f" % (np.exp(score))) _hmm = HMM(transition_probability, emission_probability, start_probability) print("前向概率 P(O|λ):%f"%(_hmm.forward( obs )) ) print("后向概率 P(O|λ):%f" % (_hmm.back( obs ))) if __name__ == '__main__': assessment()
1.2 单个状态概率,两个状态概率
由前向概率和后向概率可得到另外两个有用的概率:单个状态概率,两个状态概率,这两个概率都是条件概率,表示固定一个或两个时刻,计算特定隐藏状态转至观察结果的概率,在后面的'学习问题'中将要使用这两个概率。单个状态概率定义为:
(1.7)
上式分子的联合分布概率表示在t时刻,隐藏状态为i而观察结果为ot的概率。ot出现在整个观察结果的中间部分,整个观察结果在时刻t可划分为两部分:o1,o2,o3,...,ot;ot+1,ot+2,ot+1,...,oT。前半段可用前向概率表示,后半段可用后向概率表示,利用乘法原理可知联合概率分布为:
(1.7.1)
分母P(ot)表示t时刻观察结果为ot的概率,由加法原理知:
(1.7.2)
(1.7.1),(1.7.2)代入(1.7)中可得单个状态概率:
(1.8)
两个时刻状态概率定义为在时刻t隐藏变量状态为i,t+1时刻隐藏变量为j:
(1.9)
可利用前后概率计算出联合概率,此时观察结果仍可视为两个部分:o1,o2,o3,...,ot;ot+1,ot+2,ot+1,...,oT。前半段用表示,而可表示为:
,表示观察结果序列少了在t+1时刻观察值ot+1的概率,而在t+1时刻隐藏变量值是j,所以ot+1的概率为,由乘法原理知联合概率结果为:
(1.9.1)
(1.9.1)并结合加法原理可得分母公式形式:
(1.9.2)
(1.9.1)、(1.9.2)代入公式(1.9)得两时刻概率公式为:
(1.10)
二、解码问题
2.1 维特比算法
解码问题是在已知参数集合λ(A,B,Pi)前提下,通过观察结果得到隐藏变量的序列,如开篇介绍的抽取小球问题中,通过红白球的结果得到每次选择的箱子编号。解码问题可用穷举法实现,实际应用使用基于动态规划的维特比算法,关于动态规划知识点请参见动态规划解决背包问题。如有一个隐马尔科夫链隐藏状态为3个,并得到含有3个数值的观察结果,如下图所示:
在t=1时刻称为第1阶段,分别列出3个隐藏状态转为观察结果o1概率:
称为指标函数值,k表示处在动态规划的哪个阶段,si表示在该阶段结束时的状态,同样的可以列出在第2阶段指标函数:
f2(s1)的函数中有三项,代表着从阶段1到阶段2状态为s1时有三种可能,如下图所示:
根据最大似然法,利用max函数将三种可能性最大值作为指标函数f2(s1)的值,按此规律可以得到f2(s2)、f2(s3):
从S端到状态1或状态2、3称为一个决策,由于f1(s2)、f1(s2)、f1(s3)都是已知值,所以f2(s2)、f2(s2)、f2(s3)也可以求出具体数值,可列出第3阶段各个指标函数:
到E代表动态规划终止,此时有:
f4(E)为观察结果的最大概率,从E点逆向追溯到开始状态S,取各个阶段指标函数的最大值所对应的状态就可得到解码问题的解,同时也是隐藏状态最优解,维特比算法是逆向推导的动态规划,下面代码实现了维特比算法:
from hmmlearn import hmm import numpy as np #隐变量,3个箱子 states = ['box1', 'box2', 'box3'] #观察结果状态:红色球、白色球 observations = ["red", "white"] n_states = len(states) start_probability = np.array([0.2, 0.5, 0.3],np.float64) transition_probability = np.array([ [0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5] ],np.float64) emission_probability = np.array([ [0.5, 0.5], [0.4, 0.6], [0.7, 0.3] ],np.float64) class HMM(object): def __init__(self, A, B, pi): ''''' A: 状态转移概率矩阵 B: 输出观察概率矩阵 pi: 初始化状态向量 ''' self.A = np.array(A) self.B = np.array(B) self.pi = np.array(pi) self.N = self.A.shape[0] # 总共状态个数 self.M = self.B.shape[1] # 总共观察值个数 #维特比算法 def viterbi(self, O): self.O = O T = len(self.O) I = np.zeros(T, int) #f是每个阶段的指标函数 f = np.zeros((T, self.N), np.float64) #strategies 保存指标函数中max{...}最大值项对应的状态,方便在追溯时得到隐藏状态序列 strategies = np.zeros((T, self.N), np.float64) for i in range(self.N): f[0, i] = self.pi[i] * self.B[i, self.O[0]] strategies[0, i] = 0 for t in range(1, T): for i in range(self.N): #计算指标函数 f[t, i] = self.B[i, self.O[t]] * np.array([f[t - 1, j] * self.A[j, i] for j in range(self.N)]).max() #得到max{...}最大值项对应的状态 strategies[t, i] = np.array([f[t - 1, j] * self.A[j, i] for j in range(self.N)]).argmax() I[T - 1] = f[T - 1, :].argmax() # 逆向回溯 for t in range(T - 2, -1, -1): I[t] = strategies [t + 1, I[t + 1]] return I #解码问题 def decode(): obs = [1,1, 0 ,0 ,1,0,1,0,1,0,0] model = hmm.MultinomialHMM(n_components=n_states, n_iter=100, tol=1e-3) model.startprob_ = start_probability model.transmat_ = transition_probability model.emissionprob_ = emission_probability logprob, box = model.decode( np.array([obs]).T, algorithm="viterbi") print('MultinomialHMM viterbi算法') print([states[i] for i in box]) print('实现viterbi算法') _hmm = HMM(transition_probability, emission_probability, start_probability) hiddens=_hmm.viterbi(obs) print([ states[i] for i in hiddens]) if __name__ == '__main__': decode()
用hmmlearn库验证我们实现的维特比方法,得到如下结果:
三、学习问题
学习问题仅知道隐藏状态的个数,不妨为N个,通过观察结果求出参数集合:
Baum-Welch算法解决学习问题,Baum-Welch算法本质是EM算法(期望最大化算法):通过E、M步骤迭代得到最优解,其中E步求出目标函数的下界函数,M步求出下界函数本轮最大值,当误差小于设定值后得到最优解,详细算法参见EM算法详解。当每个时刻的隐藏状态已知、观察结果为T个,HMM的观察结果概率可表达为:
代表每个时刻隐藏状态的序号,下标数字代表时刻,显然有:
学习问题中未知、且每个时刻隐藏状态也是未知的,用I表示所有隐藏状态可能性的集合:
由加法原理知,此时观察结果的概率为:
上式取对数后得到函数f(λ):
根据Jesson不等式,可得函数f(λ)的下界函数Q(λ):
根据数学期望公式,下界函数Q(λ)可以进一步展开:
(3.1)
以上是EM算法的E步骤,得到最大似然法函数f(λ)的下界函数Q(λ)。接下来通过M步骤,对公式(3.1)求导得到本轮参数集合λ的最优解。EM算法随机会初始化λ值,不妨设初始化参数集合为,引入后公式(3.1)中以下三个元素都是常数:
此时下界函数Q(λ)有以下形式:
(3.1.1)
函数中分为三部分,分别是初始概率矩阵、发射概率矩阵、转移概率矩阵,这三者之间是相互独立的,分别求出这三部分最值,也就得到了函数最值。
3.1 初始状态概率Pi
先求初始概率矩阵Pi,目标函数为:
另有约束条件:
利用拉格朗日乘数法得拉格朗日函数:
上一篇 动态规划解决背包问题 | 下一篇 利用均匀分布生成其他分布随机数、蒙特卡罗方法原理 |
评论区 |