当前位置: 首页 > 编程笔记 >

python实现隐马尔科夫模型HMM

江航
2023-03-14
本文向大家介绍python实现隐马尔科夫模型HMM,包括了python实现隐马尔科夫模型HMM的使用技巧和注意事项,需要的朋友参考一下

一份完全按照李航<<统计学习方法>>介绍的HMM代码,供大家参考,具体内容如下

#coding=utf8 
''''' 
Created on 2017-8-5 
里面的代码许多地方可以精简,但为了百分百还原公式,就没有精简了。 
@author: adzhua 
''' 
 
import numpy as np 
 
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]  # 总共观察值个数   
    
   
  # 输出HMM的参数信息 
  def printHMM(self): 
    print ("==================================================") 
    print ("HMM content: N =",self.N,",M =",self.M) 
    for i in range(self.N): 
      if i==0: 
        print ("hmm.A ",self.A[i,:]," hmm.B ",self.B[i,:]) 
      else: 
        print ("   ",self.A[i,:],"    ",self.B[i,:]) 
    print ("hmm.pi",self.pi) 
    print ("==================================================") 
           
   
  # 前向算法  
  def forwar(self, T, O, alpha, prob): 
    ''''' 
    T: 观察序列的长度 
    O: 观察序列 
    alpha: 运算中用到的临时数组 
    prob: 返回值所要求的概率 
    '''   
     
    # 初始化 
    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]]     
     
    # 终止 
    sum = 0.0 
    for i in range(self.N): 
      sum += alpha[T-1, i] 
     
    prob[0] *= sum   
 
   
  # 带修正的前向算法 
  def forwardWithScale(self, T, O, alpha, scale, prob): 
    scale[0] = 0.0 
     
    # 初始化 
    for i in range(self.N): 
      alpha[0, i] = self.pi[i] * self.B[i, O[0]] 
      scale[0] += alpha[0, i] 
       
    for i in range(self.N): 
      alpha[0, i] /= scale[0] 
     
    # 递归 
    for t in range(T-1): 
      scale[t+1] = 0.0 
      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]] 
        scale[t+1] += alpha[t+1, j] 
       
      for j in range(self.N): 
        alpha[t+1, j] /= scale[t+1] 
      
    # 终止 
    for t in range(T): 
      prob[0] += np.log(scale[t])     
       
       
  def back(self, T, O, beta, prob):  
    ''''' 
    T: 观察序列的长度  len(O) 
    O: 观察序列 
    beta: 计算时用到的临时数组 
    prob: 返回值;所要求的概率 
    '''  
     
    # 初始化         
    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 
     
    # 终止 
    sum = 0.0 
    for i in range(self.N): 
      sum += self.pi[i]*self.B[i,O[0]]*beta[0,i] 
     
    prob[0] = sum   
     
     
  # 带修正的后向算法 
  def backwardWithScale(self, T, O, beta, scale): 
    ''''' 
    T: 观察序列的长度 len(O) 
    O: 观察序列 
    beta: 计算时用到的临时数组 
    ''' 
    # 初始化 
    for i in range(self.N): 
      beta[T-1, i] = 1.0 
     
    # 递归         
    for t in range(T-2, -1, -1): 
      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 / scale[t+1]     
         
   
  # viterbi算法       
  def viterbi(self, O): 
    ''''' 
    O: 观察序列 
    ''' 
    T = len(O) 
    # 初始化 
    delta = np.zeros((T, self.N), np.float) 
    phi = np.zeros((T, self.N), np.float) 
    I = np.zeros(T) 
     
    for i in range(self.N): 
      delta[0, i] = self.pi[i] * self.B[i, O[0]] 
      phi[0, i] = 0.0 
     
    # 递归 
    for t in range(1, T): 
      for i in range(self.N): 
        delta[t, i] = self.B[i, O[t]] * np.array([delta[t-1, j] * self.A[j, i] for j in range(self.N)] ).max() 
        phi = np.array([delta[t-1, j] * self.A[j, i] for j in range(self.N)]).argmax() 
       
    # 终止 
    prob = delta[T-1, :].max() 
    I[T-1] = delta[T-1, :].argmax() 
     
    for t in range(T-2, -1, -1): 
      I[t] = phi[I[t+1]] 
       
     
    return prob, I 
   
   
  # 计算gamma(计算A所需的分母;详情见李航的统计学习) : 时刻t时马尔可夫链处于状态Si的概率 
  def computeGamma(self, T, alpha, beta, gamma): 
    '''''''' 
    for t in range(T): 
      for i in range(self.N): 
        sum = 0.0 
        for j in range(self.N): 
          sum += alpha[t, j] * beta[t, j] 
         
        gamma[t, i] = (alpha[t, i] * beta[t, i]) / sum   
   
  # 计算sai(i,j)(计算A所需的分子) 为给定训练序列O和模型lambda时 
  def computeXi(self, T, O, alpha, beta, Xi): 
     
    for t in range(T-1): 
      sum = 0.0 
      for i in range(self.N): 
        for j in range(self.N): 
          Xi[t, i, j] = alpha[t, i] * self.A[i, j] * self.B[j, O[t+1]] * beta[t+1, j] 
          sum += Xi[t, i, j] 
       
      for i in range(self.N): 
        for j in range(self.N): 
          Xi[t, i, j] /= sum 
   
   
  # 输入 L个观察序列O,初始模型:HMM={A,B,pi,N,M} 
  def BaumWelch(self, L, T, O, alpha, beta, gamma):                   
    DELTA = 0.01 ; round = 0 ; flag = 1 ; probf = [0.0] 
    delta = 0.0; probprev = 0.0 ; ratio = 0.0 ; deltaprev = 10e-70 
     
    xi = np.zeros((T, self.N, self.N)) # 计算A的分子 
    pi = np.zeros((T), np.float)  # 状态初始化概率 
     
    denominatorA = np.zeros((self.N), np.float) # 辅助计算A的分母的变量 
    denominatorB = np.zeros((self.N), np.float) 
    numeratorA = np.zeros((self.N, self.N), np.float)  # 辅助计算A的分子的变量 
    numeratorB = np.zeros((self.N, self.M), np.float)  # 针对输出观察概率矩阵 
    scale = np.zeros((T), np.float) 
     
    while True: 
      probf[0] =0 
       
      # E_step 
      for l in range(L): 
        self.forwardWithScale(T, O[l], alpha, scale, probf) 
        self.backwardWithScale(T, O[l], beta, scale) 
        self.computeGamma(T, alpha, beta, gamma)  # (t, i) 
        self.computeXi(T, O[l], alpha, beta, xi)  #(t, i, j) 
         
        for i in range(self.N): 
          pi[i] += gamma[0, i] 
          for t in range(T-1): 
            denominatorA[i] += gamma[t, i] 
            denominatorB[i] += gamma[t, i] 
          denominatorB[i] += gamma[T-1, i] 
         
          for j in range(self.N): 
            for t in range(T-1): 
              numeratorA[i, j] += xi[t, i, j] 
             
          for k in range(self.M): # M为观察状态取值个数 
            for t in range(T): 
              if O[l][t] == k: 
                numeratorB[i, k] += gamma[t, i]   
                 
       
      # M_step。 计算pi, A, B 
      for i in range(self.N): # 这个for循环也可以放到for l in range(L)里面 
        self.pi[i] = 0.001 / self.N + 0.999 * pi[i] / L 
         
        for j in range(self.N): 
          self.A[i, j] = 0.001 / self.N + 0.999 * numeratorA[i, j] / denominatorA[i]           
          numeratorA[i, j] = 0.0 
         
        for k in range(self.M): 
          self.B[i, k] = 0.001 / self.N + 0.999 * numeratorB[i, k] / denominatorB[i] 
          numeratorB[i, k] = 0.0   
         
        #重置 
        pi[i] = denominatorA[i] = denominatorB[i] = 0.0 
         
      if flag == 1: 
        flag = 0 
        probprev = probf[0] 
        ratio = 1 
        continue 
       
      delta = probf[0] - probprev  
      ratio = delta / deltaprev   
      probprev = probf[0] 
      deltaprev = delta 
      round += 1 
       
      if ratio <= DELTA : 
        print('num iteration: ', round)   
        break 
     
 
if __name__ == '__main__': 
  print ("python my HMM") 
   
  # 初始的状态概率矩阵pi;状态转移矩阵A;输出观察概率矩阵B; 观察序列 
  pi = [0.5,0.5] 
  A = [[0.8125,0.1875],[0.2,0.8]] 
  B = [[0.875,0.125],[0.25,0.75]] 
  O = [ 
     [1,0,0,1,1,0,0,0,0], 
     [1,1,0,1,0,0,1,1,0], 
     [0,0,1,1,0,0,1,1,1] 
    ] 
  L = len(O) 
  T = len(O[0])  # T等于最长序列的长度就好了 
   
  hmm = HMM(A, B, pi) 
  alpha = np.zeros((T,hmm.N),np.float) 
  beta = np.zeros((T,hmm.N),np.float) 
  gamma = np.zeros((T,hmm.N),np.float) 
   
  # 训练 
  hmm.BaumWelch(L,T,O,alpha,beta,gamma) 
   
  # 输出HMM参数信息 
  hmm.printHMM()  

以上就是本文的全部内容,希望对大家的学习有所帮助,也希望大家多多支持小牛知识库。

 类似资料:
  • 1. hmmlearn概述 hmmlearn安装很简单,"pip install hmmlearn"即可完成。 hmmlearn实现了三种HMM模型类,按照观测状态是连续状态还是离散状态,可以分为两类。GaussianHMM和GMMHMM是连续观测状态的HMM模型,而MultinomialHMM是离散观测状态的模型,也是我们在HMM原理系列篇里面使用的模型。 对于MultinomialHMM的模型

  • 首先我们来看看什么样的问题解决可以用HMM模型。使用HMM模型时我们的问题一般有这两个特征:1)我们的问题是基于序列的,比如时间序列,或者状态序列。2)我们的问题中有两类数据,一类序列数据是可以观测到的,即观测序列;而另一类数据是不能观察到的,即隐藏状态序列,简称状态序列。 有了这两个特征,那么这个问题一般可以用HMM模型来尝试解决。这样的问题在实际生活中是很多的。比如:我现在在打字写博客,我在键

  • 隐马尔可夫模型基础 摘要 我们如何将机器学习应用于随时间变化观察到的一系列数据中来?例如,我们可能对根据一个人讲话的录音来发现他所说的话的顺序感兴趣。或者,我们可能对用词性标记来注释单词序列感兴趣。本小节的内容对马尔可夫模型的概念进行了全面的数学介绍,该模型是一种关于状态随时间变化的推理一种学习形式。并且使用隐马尔可夫模型,我们希望从一系列观察数据中恢复这一系列模型的初始状态。最后一节包含一些特定

  • 本文向大家介绍python基于隐马尔可夫模型实现中文拼音输入,包括了python基于隐马尔可夫模型实现中文拼音输入的使用技巧和注意事项,需要的朋友参考一下 在网上看到一篇关于隐马尔科夫模型的介绍,觉得简直不能再神奇,又在网上找到大神的一篇关于如何用隐马尔可夫模型实现中文拼音输入的博客,无奈大神没给可以运行的代码,只能纯手动网上找到了结巴分词的词库,根据此训练得出隐马尔科夫模型,用维特比算法实现了一

  • 1. 马尔科夫链概述 马尔科夫链定义本身比较简单,它假设某一时刻状态转移的概率只依赖于它的前一个状态。举个形象的比喻,假如每天的天气是一个状态的话,那个今天是不是晴天只依赖于昨天的天气,而和前天的天气没有任何关系。当然这么说可能有些武断,但是这样做可以大大简化模型的复杂度,因此马尔科夫链在很多时间序列模型中得到广泛的应用,比如循环神经网络RNN,隐式马尔科夫模型HMM等,当然MCMC也需要它。 如