隐马尔科夫链HMM详解

马尔科夫链是离散的马尔科夫过程,而马尔科夫过程是一类平稳的随机过程。如果一个时间序列马尔科夫链,设时间序列在每个时刻的状态有n种状态S={s1,s2,s3,..sn},在m时刻为si,则在m+k时刻的状态sj只与m时刻相关,与m-1,m-2,m-3...,1时刻无关:

png.png

上面的条件概率说明了马尔科夫链未来状态(m+k)仅与当前状态有关(时刻m)有关,而平稳过程表明时间序列状态的转移与起始点m也无关,仅与k有关,关于平稳随机过程的定义请参见时间序列分析:AR(p),MA(q),满足平稳性的马尔科夫过程也称齐次马尔科夫过程。常研究马尔科夫链k=1时各个状态转移情况,用aij表示从状态si到下一时刻状态为sj的概率,如果隐藏变量有n种状态,则马尔科夫链的状态转移矩阵大小是png (1).png

隐马尔科夫链HMM是含有隐变量马尔科夫过程,关于隐变量在EM算法中曾提过。举一个例子,有三个箱子分别标记为1号箱、2号箱、3号箱,每个箱子各有一定比例的白色球和红色球,现将箱子标记遮住,每次随机选择一个箱子、抽取一个球,得的的一组按次序摆放的白、红球称为观察变量,每次选择的箱子是隐藏变量,隐马尔科夫链暗含了两个随机过程,一个是选择箱子的序列,另一个是观察球序列,箱子之间状态的转移称为转移概率矩阵(A),每个箱子含有白色球、红色球的比例称为发射矩阵(B),另外,初始时选择各个箱子的概率称为初始概率向量(Pi),三个变量构成一个参数集合可用png (2).png表示。HMM常见的有三类问题,分别是:

1、评估问题:已知λ,求出现某种观察结果的概率。

2、解码问题:已知λ,通过观察结果求出隐藏状态序列。

3、学习问题:λ未知,通过观察结果求解λ。

一、评估问题

1.1 前向算法、后向算法

可以用穷举法计算观察结果的概率P(O|λ),但当观察结果序列很长时穷举法的复杂度很高,前向、后向算法是利用递归方式实现评估问题求解,可大大降低算法运算量,前向算法中首先引入前向概率:

alpha.png  (1.1)

前向概率是联合概率,表示当前面t个观察结果已知时,在t时刻隐藏状态为i观察结果为ot的概率,利用加法原理可知在t+1时刻的前向概率为:

alpha2.png    (1.2)

(1.2)中n为隐藏状态的数量,aji是转移概率矩阵中相对应的元素,bi(ot+1)是发射概率矩阵中相对应的元素,如果观察结果的长度为T,很容易得到观察结果概率P(O|λ):

png (3).png       (1.3)

前向算法的初值png (4).png需要利用初始概率矩阵:

png (5).png

与前向算法相反的是后向算法,后向算法中后向概率的定义为:

png (6).png   (1.4)

后向概率从观察结果后端向前计算,另外观察(1.4)可知后向概率bata.png并不知道时刻t的观察结果ot,此时只有t+1,t+2,..T时刻的观察结果是已知的,根据(1.4)可得到2112.png形式:

png (7).png

假设t-1时刻隐藏状态为sit 时刻隐藏状态为sj,2112.png可以看成由两个部分组成,第一部分是出现观察结果ot的概率,即为:png (8).png,另一个是t时刻隐藏状态为sj、观察结果ot+1ot+2ot+3oT的概率:bata.png,利用乘法定理可知,此时的概率为:

png (12).png

而t时刻隐藏状态可能性有n个,利用加法原理可得2112.png

png (11).png  (1.5)

由上面推导可递归得到test.image.latex (1).png

test.image.latex.png

在t=1时刻,从隐藏状态si得到o1的概率为back1.png,由加法原理可知后向算法P(O|λ)为:

back.png  (1.6)

后向概率的初始概率为bi.png,可以表示为:

png (13).png

实际上oT+1是观察不到的,可以认为在T时刻隐藏状态每种可能性都有,可设bi1.png,利用上面公式来实现开篇箱子抽取红白球的评估问题:

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  单个状态概率,两个状态概率

由前向概率和后向概率可得到另外两个有用的概率:单个状态概率gamma.png,两个状态概率xi.png,这两个概率都是条件概率,表示固定一个或两个时刻,计算特定隐藏状态转至观察结果的概率,在后面的'学习问题'中将要使用这两个概率。单个状态概率gamma.png定义为:

gamma2.png (1.7)

上式分子的联合分布概率表示在t时刻,隐藏状态为i而观察结果为ot的概率。ot出现在整个观察结果的中间部分,整个观察结果在时刻t可划分为两部分:o1,o2,o3,...,ot;ot+1,ot+2,ot+1...,oT。前半段可用前向概率test.image.latex (3).png表示,后半段可用后向概率bata.png表示,利用乘法原理可知联合概率分布为:

test.image.latex (4).png   (1.7.1)

分母P(ot)表示t时刻观察结果为ot的概率,由加法原理知:

test.image.latex (5).png(1.7.2)

(1.7.1),(1.7.2)代入(1.7)中可得单个状态概率gamma.png

test.image.latex (6).png     (1.8)

两个时刻状态概率xi.png定义为在时刻t隐藏变量状态为i,t+1时刻隐藏变量为j:

test.image.latex.png(1.9)

可利用前后概率计算出联合概率test.image.latex (1).png,此时观察结果仍可视为两个部分:o1,o2,o3,...,ot;ot+1,ot+2,ot+1,...,oT。前半段用test.image.latex (3).png表示,而test.image.latex (3).png可表示为:

test.image.latex (2).png

test.image.latex (3).pngtest.image.latex (3).png表示观察结果序列少了在t+1时刻观察值ot+1的概率,而在t+1时刻隐藏变量值是j,所以ot+1的概率为ot.png,由乘法原理知联合概率test.image.latex (1).png结果为:

test.image.latex.png   (1.9.1)

(1.9.1)并结合加法原理可得分母test.image.latex (1).png公式形式:

test.image.latex (2).png        (1.9.2)

(1.9.1)、(1.9.2)代入公式(1.9)得两时刻概率xi.png公式为:

test.image.latex (3).png   (1.10)

二、解码问题

2.1 维特比算法

解码问题是在已知参数集合λ(A,B,Pi)前提下,通过观察结果得到隐藏变量的序列,如开篇介绍的抽取小球问题中,通过红白球的结果得到每次选择的箱子编号。解码问题可用穷举法实现,实际应用使用基于动态规划的维特比算法,关于动态规划知识点请参见动态规划解决背包问题。如有一个隐马尔科夫链隐藏状态为3个,并得到含有3个数值的观察结果,如下图所示:

v2-8a0e2bc64b7305a04c3aef847fa62bdb_720w.jpg

在t=1时刻称为第1阶段,分别列出3个隐藏状态转为观察结果o1概率:

test.image.latex (7).png

test.image.latex (8).png称为指标函数值,k表示处在动态规划的哪个阶段,si表示在该阶段结束时的状态,同样的可以列出在第2阶段指标函数:

test.image.latex (9).png

f2(s1)的函数中有三项,代表着从阶段1到阶段2状态为s1时有三种可能,如下图所示:

v2-ab128a2158510b6885daba65c7de50e0_720w.jpg

根据最大似然法,利用max函数将三种可能性最大值作为指标函数f2(s1)的值,按此规律可以得到f2(s2)、f2(s3):

test.image.latex (10).png

从S端到状态1或状态2、3称为一个决策,由于f1(s2)、f1(s2)、f1(s3)都是已知值,所以f2(s2)、f2(s2)、f2(s3)也可以求出具体数值,可列出第3阶段各个指标函数:

test.image.latex (11).png

到E代表动态规划终止,此时有:

test.image.latex (12).png

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库验证我们实现的维特比方法,得到如下结果:

无标题.png

三、学习问题

学习问题仅知道隐藏状态的个数,不妨为N个,通过观察结果求出参数集合:

png (2).png     test.image.latex (14).png

Baum-Welch算法解决学习问题,Baum-Welch算法本质是EM算法(期望最大化算法):通过E、M步骤迭代得到最优解,其中E步求出目标函数的下界函数,M步求出下界函数本轮最大值,当误差小于设定值后得到最优解,详细算法参见EM算法详解当每个时刻的隐藏状态已知、观察结果为T个,HMM的观察结果概率可表达为:

pl.png

test.image.latex (5).png代表每个时刻隐藏状态的序号,下标数字代表时刻,显然有:

test.image.latex (4).png

学习问题中png (2).png未知、且每个时刻隐藏状态test.image.latex (5).png也是未知的,用I表示所有隐藏状态可能性的集合:

test.image.latex (6).png

由加法原理知,此时观察结果的概率为:

PL23.png

上式取对数后得到函数f(λ):

test.image.latex (7).png

根据Jesson不等式,可得函数f(λ)的下界函数Q(λ):

test.image.latex (8).png

根据数学期望公式,下界函数Q(λ)可以进一步展开:

test.image.latex (10).png  (3.1)

以上是EM算法的E步骤,得到最大似然法函数f(λ)的下界函数Q(λ)。接下来通过M步骤,对公式(3.1)求导得到本轮参数集合λ的最优解。EM算法随机会初始化λ值,不妨设初始化参数集合test.image.latex (11).png,引入test.image.latex (11).png公式(3.1)中以下三个元素都是常数:

test.image.latex (12).png

此时下界函数Q(λ)有以下形式:

test.image.latex (13).png    (3.1.1)

test.image.latex (15).png函数中分为三部分,分别是初始概率矩阵、发射概率矩阵、转移概率矩阵,这三者之间是相互独立的,分别求出这三部分最值,也就得到了函数test.image.latex (15).png最值

3.1 初始状态概率Pi

先求初始概率矩阵Pi,目标函数为:

test.image.latex.png

另有约束条件:

test.image.latex (1).png

利用拉格朗日乘数法得拉格朗日函数:

-免费试读结束-
登录|注册后打赏作者吧! 1.0元