李航《统计学习方法》第十章——用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. 训练如何停止?我现在的做法是设置一个训练次数上限。

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/wds2006sdo/article/details/75212599

【深度剖析HMM(附Python代码)】2.隐马尔科夫链HMM的EM训练过程

隐马尔科夫链HMM的参数θ的EM训练过程 现在回到前一节最后提出的参数θ的最大似然函数上来,先对其做个对数变换,做对数变换是考虑到序列X的概率计算公式中包含了连乘,为了方便计算同时避免序列X的概率过...
  • tostq
  • tostq
  • 2017-04-27 13:28:43
  • 2546

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

1.综述:隐马尔可夫模型是马尔可夫模型的进一步发展。马尔可夫模型是马尔可夫过程的模型化, 可以用图1 (a) 的框图形象表示。它把一个总随机过程看成一系列状态的不断转移。图1(b)是隐马尔可夫模型。1...
  • yywan1314520
  • yywan1314520
  • 2015-12-10 13:20:20
  • 5096

MCMC(Markov Chain Monte Carlo)的理解与实践(Python)

Markov Chain Monte Carlo (MCMC) methods are a class of algorithms for sampling from a probability di...
  • lanchunhui
  • lanchunhui
  • 2016-01-03 20:56:01
  • 5033

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

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

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

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

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

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

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

简介 隐马尔可夫模型(Hidden Markov Model,HMM)创立于20世纪70年代。主要用于行为识别,语音识别,文字识别等。 原理简述 隐马尔可夫模型由五个...
  • xiaotianlan
  • xiaotianlan
  • 2014-08-29 14:20:51
  • 2577

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

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

python自然语言处理-马尔科夫模型

本部分通过继续处理总统演讲稿的内容,提出了马尔科夫模型,并通过实例介绍了如何进行文本分析与写作。如何生成并使用一个马尔科夫链...
  • wanght89
  • wanght89
  • 2017-10-11 09:01:37
  • 274

<em>马尔</em>可<em>夫链</em>模型

这是数学模型中的<em>马尔</em>可<em>夫链</em>模型的PDF文档以及<em>Python</em>代码。欢迎学习数学建模以及机器学习的同行下载
  • 2018年04月18日 00:00
收藏助手
不良信息举报
您举报文章:李航《统计学习方法》第十章——用Python实现隐马尔科夫模型
举报原因:
原因补充:

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