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()  

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

您可能感兴趣的文章:

  • python使用tensorflow保存、加载和使用模型的方法
  • Python实现感知器模型、两层神经网络
  • VTK与Python实现机械臂三维模型可视化详解
  • 浅谈Python基础之I/O模型
  • Python自定义进程池实例分析【生产者、消费者模型问题】
  • 总结网络IO模型与select模型的Python实例讲解
  • 理解生产者消费者模型及在Python编程中的运用实例
  • python基于隐马尔可夫模型实现中文拼音输入
  • 基于python yield机制的异步操作同步化编程模型
  • 用Python给文本创立向量空间模型的教程
(0)

相关推荐

  • VTK与Python实现机械臂三维模型可视化详解

    三维可视化系统的建立依赖于三维图形平台, 如 OpenGL.VTK.OGRE.OSG等, 传统的方法多采用OpenGL进行底层编程,即对其特有的函数进行定量操作, 需要开发人员熟悉相关函数, 从而造成了开发难度大. 周期长等问题.VTK. ORGE.OSG等平台使用封装更好的函数简化了开发过程.下面将使用Python与VTK进行机器人上位机监控界面的快速原型开发. 完整的上位机程序需要有三维显示模块.机器人信息监测模块(位置/角度/速度/电量/温度/错误信息...).通信模块(串口/USB/WI

  • 基于python yield机制的异步操作同步化编程模型

    本文总结下如何在编写python代码时对异步操作进行同步化模拟,从而提高代码的可读性和可扩展性. 游戏引擎一般都采用分布式框架,通过一定的策略来均衡服务器集群的资源负载,从而保证服务器运算的高并发性和CPU高利用率,最终提高游戏的性能和负载.由于引擎的逻辑层调用是非抢占式的,服务器之间都是通过异步调用来进行通讯,导致游戏逻辑无法同步执行,所以在代码层不得不人为地添加很多回调函数,使一个原本完整的功能碎片化地分布在各个回调函数中. 异步逻辑 以游戏中的副本评分逻辑为例,在副本结束时副本管理进程需要

  • python基于隐马尔可夫模型实现中文拼音输入

    在网上看到一篇关于隐马尔科夫模型的介绍,觉得简直不能再神奇,又在网上找到大神的一篇关于如何用隐马尔可夫模型实现中文拼音输入的博客,无奈大神没给可以运行的代码,只能纯手动网上找到了结巴分词的词库,根据此训练得出隐马尔科夫模型,用维特比算法实现了一个简单的拼音输入法.githuh地址:https://github.com/LiuRoy/Pinyin_Demo 原理简介隐马尔科夫模型 抄一段网上的定义: 隐马尔可夫模型 (Hidden Markov Model) 是一种统计模型,用来描述一个含有隐含未

  • Python自定义进程池实例分析【生产者、消费者模型问题】

    本文实例分析了Python自定义进程池.分享给大家供大家参考,具体如下: 代码说明一切: #encoding=utf-8 #author: walker #date: 2014-05-21 #function: 自定义进程池遍历目录下文件 from multiprocessing import Process, Queue, Lock import time, os #消费者 class Consumer(Process): def __init__(self, queue, ioLock):

  • 用Python给文本创立向量空间模型的教程

    我们需要开始思考如何将文本集合转化为可量化的东西.最简单的方法是考虑词频. 我将尽量尝试不使用NLTK和Scikits-Learn包.我们首先使用Python讲解一些基本概念. 基本词频 首先,我们回顾一下如何得到每篇文档中的词的个数:一个词频向量. #examples taken from here: http://stackoverflow.com/a/1750187 mydoclist = ['Julie loves me more than Linda loves me', 'Jane

  • 浅谈Python基础之I/O模型

    一.I/O模型 IO在计算机中指Input/Output,也就是输入和输出.由于程序和运行时数据是在内存中驻留,由CPU这个超快的计算核心来执行,涉及到数据交换的地方,通常是磁盘.网络等,就需要IO接口. 同步(synchronous) IO和异步(asynchronous) IO,阻塞(blocking) IO和非阻塞(non-blocking)IO分别是什么,到底有什么区别? 这个问题其实不同的人给出的答案都可能不同,比如wiki,就认为asynchronous IO和non-blockin

  • python使用tensorflow保存、加载和使用模型的方法

    使用Tensorflow进行深度学习训练的时候,需要对训练好的网络模型和各种参数进行保存,以便在此基础上继续训练或者使用.介绍这方面的博客有很多,我发现写的最好的是这一篇官方英文介绍: http://cv-tricks.com/tensorflow-tutorial/save-restore-tensorflow-models-quick-complete-tutorial/ 我对这篇文章进行了整理和汇总. 首先是模型的保存.直接上代码: #!/usr/bin/env python #-*- c

  • Python实现感知器模型、两层神经网络

    本文实例为大家分享了Python实现感知器模型.两层神经网络,供大家参考,具体内容如下 python 3.4 因为使用了 numpy 这里我们首先实现一个感知器模型来实现下面的对应关系 [[0,0,1], --- 0 [0,1,1], --- 1 [1,0,1], --- 0 [1,1,1]] --- 1 从上面的数据可以看出:输入是三通道,输出是单通道. 这里的激活函数我们使用 sigmoid 函数 f(x)=1/(1+exp(-x)) 其导数推导如下所示: L0=W*X; z=f(L0);

  • 总结网络IO模型与select模型的Python实例讲解

    网络I/O模型 人多了,就会有问题.web刚出现的时候,光顾的人很少.近年来网络应用规模逐渐扩大,应用的架构也需要随之改变.C10k的问题,让工程师们需要思考服务的性能与应用的并发能力. 网络应用需要处理的无非就是两大类问题,网络I/O,数据计算.相对于后者,网络I/O的延迟,给应用带来的性能瓶颈大于后者.网络I/O的模型大致有如下几种: 同步模型(synchronous I/O) 阻塞I/O(bloking I/O) 非阻塞I/O(non-blocking I/O) 多路复用I/O(multi

  • 理解生产者消费者模型及在Python编程中的运用实例

    什么是生产者消费者模型 在 工作中,大家可能会碰到这样一种情况:某个模块负责产生数据,这些数据由另一个模块来负责处理(此处的模块是广义的,可以是类.函数.线程.进程等).产 生数据的模块,就形象地称为生产者:而处理数据的模块,就称为消费者.在生产者与消费者之间在加个缓冲区,我们形象的称之为仓库,生产者负责往仓库了进商 品,而消费者负责从仓库里拿商品,这就构成了生产者消费者模型.结构图如下: 生产者消费者模型的优点: 1.解耦 假设生产者和消费者分别是两个类.如果让生产者直接调用消费者的某个方法,

随机推荐