李航《统计学习方法》第十章——用Python实现隐马尔科夫模型

原创 2017年07月16日 22:52:45

相关文章:

隐马尔科夫模型有3个基本问题:
1. 概率计算问题(前向算法,后向算法)
2. 学习问题 (Baum-Welch算法)
3. 预测问题 (维特比算法)

我实现了Baum-Welch算法,且该算法也包含了前向算法与后向算法。

Baum-Welch算法

这里先贴上书中的算法

这里写图片描述
这里写图片描述

数据集

本来打算试一下用自己写的HMM跑一下中文分词,但很可惜,代码运行的比较慢。
所以改成 模拟 三角波 以及 正弦波

代码

代码已放到Github上,这边也贴出来

# encoding=utf8

import numpy as np
import csv

class HMM(object):
    def __init__(self,N,M):
        self.A = np.zeros((N,N))        # 状态转移概率矩阵
        self.B = np.zeros((N,M))        # 观测概率矩阵
        self.Pi = np.array([1.0/N]*N)   # 初始状态概率矩阵

        self.N = N                      # 可能的状态数
        self.M = M                      # 可能的观测数

    def cal_probality(self, O):
        self.T = len(O)
        self.O = O

        self.forward()
        return sum(self.alpha[self.T-1])

    def forward(self):
        """
        前向算法
        """
        self.alpha = np.zeros((self.T,self.N))

        # 公式 10.15
        for i in range(self.N):
            self.alpha[0][i] = self.Pi[i]*self.B[i][self.O[0]]

        # 公式10.16
        for t in range(1,self.T):
            for i in range(self.N):
                sum = 0
                for j in range(self.N):
                    sum += self.alpha[t-1][j]*self.A[j][i]
                self.alpha[t][i] = sum * self.B[i][self.O[t]]

    def backward(self):
        """
        后向算法
        """
        self.beta = np.zeros((self.T,self.N))

        # 公式10.19
        for i in range(self.N):
            self.beta[self.T-1][i] = 1

        # 公式10.20
        for t in range(self.T-2,-1,-1):
            for i in range(self.N):
                for j in range(self.N):
                    self.beta[t][i] += self.A[i][j]*self.B[j][self.O[t+1]]*self.beta[t+1][j]

    def cal_gamma(self, i, t):
        """
        公式 10.24
        """
        numerator = self.alpha[t][i]*self.beta[t][i]
        denominator = 0

        for j in range(self.N):
            denominator += self.alpha[t][j]*self.beta[t][j]

        return numerator/denominator

    def cal_ksi(self, i, j, t):
        """
        公式 10.26
        """

        numerator = self.alpha[t][i]*self.A[i][j]*self.B[j][self.O[t+1]]*self.beta[t+1][j]
        denominator = 0

        for i in range(self.N):
            for j in range(self.N):
                denominator += self.alpha[t][i]*self.A[i][j]*self.B[j][self.O[t+1]]*self.beta[t+1][j]

        return numerator/denominator

    def init(self):
        """
        随机生成 A,B,Pi
        并保证每行相加等于 1
        """
        import random
        for i in range(self.N):
            randomlist = [random.randint(0,100) for t in range(self.N)]
            Sum = sum(randomlist)
            for j in range(self.N):
                self.A[i][j] = randomlist[j]/Sum

        for i in range(self.N):
            randomlist = [random.randint(0,100) for t in range(self.M)]
            Sum = sum(randomlist)
            for j in range(self.M):
                self.B[i][j] = randomlist[j]/Sum

    def train(self, O, MaxSteps = 100):
        self.T = len(O)
        self.O = O

        # 初始化
        self.init()

        step = 0
        # 递推
        while step<MaxSteps:
            step+=1
            print(step)
            tmp_A = np.zeros((self.N,self.N))
            tmp_B = np.zeros((self.N,self.M))
            tmp_pi = np.array([0.0]*self.N)

            self.forward()
            self.backward()

            # a_{ij}
            for i in range(self.N):
                for j in range(self.N):
                    numerator=0.0
                    denominator=0.0
                    for t in range(self.T-1):
                        numerator += self.cal_ksi(i,j,t)
                        denominator += self.cal_gamma(i,t)
                    tmp_A[i][j] = numerator/denominator

            # b_{jk}
            for j in range(self.N):
                for k in range(self.M):
                    numerator = 0.0
                    denominator = 0.0
                    for t in range(self.T):
                        if k == self.O[t]:
                            numerator += self.cal_gamma(j,t)
                        denominator += self.cal_gamma(j,t)
                    tmp_B[j][k] = numerator / denominator

            # pi_i
            for i in range(self.N):
                tmp_pi[i] = self.cal_gamma(i,0)

            self.A = tmp_A
            self.B = tmp_B
            self.Pi = tmp_pi

    def generate(self, length):
        import random
        I = []

        # start
        ran = random.randint(0,1000)/1000.0
        i = 0
        while self.Pi[i]<ran or self.Pi[i]<0.0001:
            ran -= self.Pi[i]
            i += 1
        I.append(i)

        # 生成状态序列
        for i in range(1,length):
            last = I[-1]
            ran = random.randint(0, 1000) / 1000.0
            i = 0
            while self.A[last][i] < ran or self.A[last][i]<0.0001:
                ran -= self.A[last][i]
                i += 1
            I.append(i)

        # 生成观测序列
        Y = []
        for i in range(length):
            k = 0
            ran = random.randint(0, 1000) / 1000.0
            while self.B[I[i]][k] < ran or self.B[I[i]][k]<0.0001:
                ran -= self.B[I[i]][k]
                k += 1
            Y.append(k)

        return Y



def triangle(length):
    '''
    三角波
    '''
    X = [i for i in range(length)]
    Y = []

    for x in X:
        x = x % 6
        if x <= 3:
            Y.append(x)
        else:
            Y.append(6-x)
    return X,Y



def show_data(x,y):
    import matplotlib.pyplot as plt
    plt.plot(x, y, 'g')
    plt.show()

    return y


if __name__ == '__main__':
    hmm = HMM(10,4)
    tri_x, tri_y = triangle(20)

    hmm.train(tri_y)
    y = hmm.generate(100)
    print(y)
    x = [i for i in range(100)]
    show_data(x,y)

运行结果

三角波

三角波比较简单,我设置N=10,扔进去长度为20的序列,训练100次,下图是其生成的长度为100的序列。
可以看出效果还是很不错的。

这里写图片描述

正弦波

正弦波有些麻烦,因为观测序列不能太大,所以我设置N=15,M=100,扔进去长度为40的序列,训练100次,训练的非常慢,下图是其生成的长度为100的序列。
可以看出效果一般,如果改成C的代码,并增加N应该是能模拟的。
这里写图片描述

balabala

隐马尔科夫模型的初值真的非常重要,我曾尝试过每行平均和单位矩阵两种方式,结果100轮迭代后啥都没变,所以最好还是随机初始化。

还有几点细节不是很理解
1. 如果有多个序列如何同时训练?
2. 训练如何停止?我现在的做法是设置一个训练次数上限。

版权声明:本文为博主原创文章,未经博主允许不得转载。

Python实现HMM(隐马尔可夫模型)

HMM的代码,这次用Python写了一遍,依据仍然是李航博士的《统计学习方法》 由于第一次用python,所以代码可能会有许多缺陷,但是所有代码都用书中的例题进行了测试,结果正确。 这...
  • liuyukuan
  • liuyukuan
  • 2017年05月22日 14:03
  • 1698

Python实现HMM(隐马尔可夫模型)

前几天用MATLAB实现了HMM的代码,这次用python写了一遍,依据仍然是李航博士的《统计学习方法》 由于第一次用python,所以代码可能会有许多缺陷,但是所有代码都用书中的例题进行了测试,结果...
  • sinat_36005594
  • sinat_36005594
  • 2017年04月07日 16:13
  • 2837

隐马尔可夫模型 HMM 的python实现

摘要: 本文主要介绍隐马尔可夫模型HMM的python实现,参考的文献主要是: [1]. Lawrence R. Rabiner, ‘A Tutorial on Hidden Markov...
  • qq_27824601
  • qq_27824601
  • 2016年08月12日 15:40
  • 552

【深度剖析HMM(附Python代码)】3.隐马尔科夫链所解决的问题

通过隐马尔科夫链,有以下几方面功能: 1. 解码问题  已知某一序列,找到最可能的隐藏状态序列(即所谓的解码问题,利用维比特算法来解决) 解码过程的相关pytho...
  • tostq
  • tostq
  • 2017年04月27日 13:51
  • 1124

【深度剖析HMM(附Python代码)】1.前言及隐马尔科夫链HMM的背景

1. 前言 隐马尔科夫HMM模型是一类重要的机器学习方法,其主要用于序列数据的分析,广泛应用于语音识别、文本翻译、序列预测、中文分词等多个领域。虽然近年来,由于RNN等深度学习方法的发展,HMM模型...
  • tostq
  • tostq
  • 2017年04月27日 12:24
  • 2484

HMM:隐马尔可夫模型 - 表示

http://blog.csdn.net/pipisorry/article/details/50722178隐马尔可夫模型 Hidden Markov Models隐马尔可夫模型(Hidden Ma...
  • pipisorry
  • pipisorry
  • 2016年02月23日 12:37
  • 4744

隐马尔可夫模型(HMM)实现分词

隐马尔可夫模型(HMM)实现分词最近在调研分词的一些工作,又重新看了一遍HMM,发现原来之前看过的,竟然给忘的干干净净!之前看的,真是喂了狗了,于是决定把这次的理解记录下来,俗话说,好记性不如赖笔头嘛...
  • zkl33a
  • zkl33a
  • 2017年03月31日 17:00
  • 168

隐马尔科夫模型(HMM)及其实现

马尔科夫模型 马尔科夫模型是单重随机过程,是一个2元组:(S,A)。 其中S是状态集合,A是状态转移矩阵。 只用状态转移来描述随机过程。   马尔科夫模型的2个假设 有限历史性假设:t+l时刻系统状态...
  • u012428391
  • u012428391
  • 2014年08月27日 10:01
  • 3374

ML--HMM(隐马尔可夫模型及python的实现2)

1.HMM的应用1,这个代码不知道出处了,若有侵权请联系本文作者删除,注释为本人所加。2.对基本的HMM需要进一步了解的,请戳这里3.下面是HMM代码的解释之一# _*_ coding:utf-8 _...
  • yywan1314520
  • yywan1314520
  • 2016年01月17日 22:26
  • 4293

隐马尔可夫模型 HMM 原理及实现

简介 隐马尔可夫模型(Hidden Markov Model,HMM)创立于20世纪70年代。主要用于行为识别,语音识别,文字识别等。 原理简述 隐马尔可夫模型由五个...
  • xiaotianlan
  • xiaotianlan
  • 2014年08月29日 14:20
  • 2276
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:李航《统计学习方法》第十章——用Python实现隐马尔科夫模型
举报原因:
原因补充:

(最多只允许输入30个字)