1. 前言
隐马尔科夫 HMM 模型是一类重要的机器学习方法, 其主要用于序列数据的分析, 广泛应用于语音识别, 文本翻译, 序列预测, 中文分词等多个领域. 虽然近年来, 由于 RNN 等深度学习方法的发展, HMM 模型逐渐变得不怎么流行了, 但并不意味着完全退出应用领域, 甚至在一些轻量级的任务中仍有应用. 本系列博客将详细剖析隐马尔科夫链 HMM 模型, 同以往网络上绝大多数教程不同, 本系列博客将更深入地分析 HMM, 不仅包括估计序列隐状态的维特比算法(HMM 解码问题), 前向后向算法等, 而且还着重的分析 HMM 的 EM 训练过程, 并将所有的过程都通过数学公式进行推导.
由于隐马尔科夫 HMM 模型是一类非常复杂的模型, 其中包含了大量概率统计的数学知识, 因此网络上多数博客一般都采用举例等比较通俗的方式来介绍 HMM, 这么做会让初学者很快明白 HMM 的原理, 但要丢失了大量细节, 让初学者处于一种似懂非懂的状态. 而本文并没有考虑用非常通俗的文字描述 HMM, 还是考虑通过详细的数学公式来一步步引导初学者掌握 HMM 的思想. 另外, 本文重点分析了 HMM 的 EM 训练过程, 这是网络上其他教程所没有的, 而个人认为相比于维特比算法, 前向后向算法, HMM 的 EM 训练过程虽然更为复杂, 但是一旦掌握这个训练过程, 那么对于通用的链状图结构的推导, EM 算法和网络训练的理解都会非常大的帮助. 另外通过总结 HMM 的数学原理, 也能非常方便将数学公式改写成代码.
最后, 本文提供了一个简单版本的隐马尔科夫链 HMM 的 Python 代码, 包含了高斯模型的 HMM 和离散 HMM 两种情况, 代码中包含了 HMM 的训练, 预测, 解码等全部过程, 核心代码总共只有 200~300 行代码, 非常简单! 个人代码水平比较渣 =_=||, 大家按照我的教程, 应该都可以写出更鲁棒性更有高效的代码, 附上 Github 地址: https://github.com/tostq/Easy_HMM
觉得好, 就点星哦! 觉得好, 就点星哦!
觉得好, 就点星哦!
重要的事要说三遍!!!!
为了方便大家学习, 我将整个 HMM 代码完善成整个学习项目, 其中包括
hmm.py:HMM 核心代码, 包含了一个 HMM 基类, 一个高斯 HMM 模型及一个离散 HMM 模型
DiscreteHMM_test.py 及 GaussianHMM_test.py: 利用 unnitest 来测试我们的 HMM, 同时引入了一个经典 HMM 库 hmmlearn 作为对照组
Dice_01.py: 利用离散 HMM 模型来解决丢色子问题的例子
Wordseg_02.py: 解决中文分词问题的例子
Stock_03.py: 解决股票预测问题的例子
2. 隐马尔科夫链 HMM 的背景
首先, 已知一组序列 :
我们从这组序列中推导出产生这组序列的函数, 假设函数参数为 , 其表示为
即使得序列 X 发生概率最大的函数参数, 要解决上式, 最简单的考虑是将序列 的每个数据都视为独立的, 比如建立一个神经网络. 然后这种考虑会随着序列增长, 而导致参数爆炸式增长. 因此可以假设当前序列数据只与其前一数据值相关, 即所谓的一阶马尔科夫链:
有一阶马尔科夫链, 也会有二阶马尔科夫链(即当前数据值取决于其前两个数据值)
当前本文不对二阶马尔科夫链进行深入分析了, 着重考虑一阶马尔科夫链, 现在根据一阶马尔科夫链的假设, 我们有:
因此要解一阶马尔科夫链, 其关键在于求数据 (以下称观测值) 之间转换函数 , 如果假设转换函数 同序列中位置 (时间)无关, 我们就能根据转换函数 而求出整个序列的概率:
然而, 如果观测值 x 的状态非常多(特别极端的情况是连续数据), 转换函数 会变成一个非常大的矩阵, 如果 x 的状态有 K 个, 那么转换函数 就会是一个 K*(K-1) 个参数, 而且对于连续变量观测值更是困难.
为了降低马尔科夫链的转换函数的参数量, 我们引入了一个包含较少状态的隐状态值, 将观测值的马尔科夫链转换为隐状态的马尔科夫链(即为隐马尔科夫链 HMM)
其包含了一个重要假设: 当前观测值只由当前隐状态所决定. 这么做的一个重要好处是, 隐状态值的状态远小于观测值的状态, 因此隐藏状态的转换函数 的参数更少.
此时我们要决定的问题是:
即在所有可能隐藏状态序列情况下, 求使得序列 发生概率最大的函数参数 .
这里我们再总结下:
隐马尔科夫链 HMM 三个重要假设:
1. 当前观测值只由当前隐藏状态确定, 而与其他隐藏状态或观测值无关(隐藏状态假设)
2. 当前隐藏状态由其前一个隐藏状态决定(一阶马尔科夫假设)
3. 隐藏状态之间的转换函数概率不随时间变化(转换函数稳定性假设)
隐马尔科夫链 HMM 所要解决的问题:
在所有可能隐藏状态序列情况下, 求使得当前序列 X 产生概率最大的函数参数θ.
代码
http://blog.csdn.net/sinat_36005594/article/details/69568538
前几天用 MATLAB 实现了 HMM 的代码, 这次用 python 写了一遍, 依据仍然是李航博士的统计学习方法
由于第一次用 python, 所以代码可能会有许多缺陷, 但是所有代码都用书中的例题进行了测试, 结果正确.
这里想说一下 python, 在编写 HMM 过程中参看了之前写的 MATLAB 程序, 发现他们有太多相似的地方, 用到了 numpy 库, 在 python 代码过程中最让我头疼的是数组角标, 和 MATLAB 矩阵角标从 1 开始不同, numpy 库数组角标都是从 0 开始, 而且数组的维数也需要谨慎, 一不小心就会出现 too many indices for array 的错误. 程序中最后是维特比算法, 在运行过程中出现了__main__:190: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future 的警告, 还没有去掉这个警告, 查了一下说不影响结果, 后面会去解决这个问题, 下面贴出我的代码
# -*- coding: utf-8 -*-
"""
Created on Thu Feb 16 19:28:39 2017
2017-4-2
ForwardBackwardAlg 函数功能: 实现前向算法
理论依据: 李航统计学习方法
2017-4-5
修改了 ForwardBackwardAlg 函数名称为 ForwardAlgo 以及输出的 alpha 数组形式
完成了 BackwardAlgo 函数功能: 后向算法
以及函数 FBAlgoAppli: 计算在观测序列和模型参数确定的情况下,
某一个隐含状态对应相应的观测状态的概率
2017-4-6
完成 BaumWelchAlgo 函数一次迭代
2017-4-7
实现维特比算法
@author: sgp """
import numpy as np
# 输入格式如下:
#A = np.array([[.5,.2,.3],[.3,.5,.2],[.2,.3,.5]])
#B = np.array([[.5,.5],[.4,.6],[.7,.3]])
#Pi = np.array([[.2,.4,.4]])
#O = np.array([[1,2,1]])
# 应用 ndarray 在数组之间进行相互运算时, 一定要确保数组维数相同!
# 比如:
#In[93]:m = np.array([1,2,3,4])
#In[94]:m
#Out[94]: array([1, 2, 3, 4])
#In[95]:m.shape
#Out[95]: (4,)
# 这里表示的是一维数组
#In[96]:m = np.array([[1,2,3,4]])
#In[97]:m
#Out[97]: array([[1, 2, 3, 4]])
#In[98]:m.shape
#Out[98]: (1, 4)
# 而这里表示的就是二维数组
# 注意 In[93]和 In[96]的区别, 多一对中括号!!
#N = A.shape[0]为数组 A 的行数, H = O.shape[1]为数组 O 的列数
# 在下列各函数中, alpha 数组和 beta 数组均为 N*H 二维数组, 也就是横向坐标是时间, 纵向是状态
def ForwardAlgo(A,B,Pi,O):
N = A.shape[0]# 数组 A 的行数
M = A.shape[1]# 数组 A 的列数
H = O.shape[1]# 数组 O 的列数
sum_alpha_1 = np.zeros((M,N))
alpha = np.zeros((N,H))
r = np.zeros((1,N))
alpha_1 = np.multiply(Pi[0,:], B[:,O[0,0]-1])
alpha[:,0] = np.array(alpha_1).reshape(1,N)#alpha_1 是一维数组, 在使用 np.multiply 的时候需要升级到二维数组.# 错误是 IndexError: too many indices for array
for h in range(1,H):
for i in range(N):
for j in range(M):
sum_alpha_1[i,j] = alpha[j,h-1] * A[j,i]
r = sum_alpha_1.sum(1).reshape(1,N)# 同理, 将数组升级为二维数组
alpha[i,h] = r[0,i] * B[i,O[0,h]-1]
#print("alpha矩阵: \n % r " % alpha)
p = alpha.sum(0).reshape(1,H)
P = p[0,H-1]
#print("观测概率: \n % r " % P)
#return alpha
return alpha, P
def BackwardAlgo(A,B,Pi,O):
N = A.shape[0]# 数组 A 的行数
M = A.shape[1]# 数组 A 的列数
H = O.shape[1]# 数组 O 的列数
#beta = np.zeros((N,H))
sum_beta = np.zeros((1,N))
beta = np.zeros((N,H))
beta[:,H-1] = 1
p_beta = np.zeros((1,N))
for h in range(H-1,0,-1):
for i in range(N):
for j in range(M):
sum_beta[0,j] = A[i,j] * B[j,O[0,h]-1] * beta[j,h]
beta[i,h-1] = sum_beta.sum(1)
#print("beta矩阵: \n % r " % beta)
for i in range(N):
p_beta[0,i] = Pi[0,i] * B[i,O[0,0]-1] * beta[i,0]
p = p_beta.sum(1).reshape(1,1)
#print("观测概率: \n % r " % p[0,0])
return beta, p[0,0]
def FBAlgoAppli(A,B,Pi,O,I):
#计算在观测序列和模型参数确定的情况下, 某一个隐含状态对应相应的观测状态的概率
#例题参考李航统计学习方法 P189 习题 10.2
#输入格式:
#I 为二维数组, 存放所求概率 P(it = qi,O|lambda)中 it 和 qi 的角标 t 和 i, 即 P=[t,i]
alpha,p1 = ForwardAlgo(A,B,Pi,O)
beta,p2 = BackwardAlgo(A,B,Pi,O)
p = alpha[I[0,1]-1,I[0,0]-1] * beta[I[0,1]-1,I[0,0]-1] / p1
return p
def GetGamma(A,B,Pi,O):
N = A.shape[0]# 数组 A 的行数
H = O.shape[1]# 数组 O 的列数
Gamma = np.zeros((N,H))
alpha,p1 = ForwardAlgo(A,B,Pi,O)
beta,p2 = BackwardAlgo(A,B,Pi,O)
for h in range(H):
for i in range(N):
Gamma[i,h] = alpha[i,h] * beta[i,h] / p1
return Gamma
def GetXi(A,B,Pi,O):
N = A.shape[0]# 数组 A 的行数
M = A.shape[1]# 数组 A 的列数
H = O.shape[1]# 数组 O 的列数
Xi = np.zeros((H-1,N,M))
alpha,p1 = ForwardAlgo(A,B,Pi,O)
beta,p2 = BackwardAlgo(A,B,Pi,O)
for h in range(H-1):
for i in range(N):
for j in range(M):
Xi[h,i,j] = alpha[i,h] * A[i,j] * B[j,O[0,h+1]-1] * beta[j,h+1] / p1
#print("Xi矩阵: \n % r " % Xi)
return Xi
def BaumWelchAlgo(A,B,Pi,O):
N = A.shape[0]# 数组 A 的行数
M = A.shape[1]# 数组 A 的列数
Y = B.shape[1]# 数组 B 的列数
H = O.shape[1]# 数组 O 的列数
c = 0
Gamma = GetGamma(A,B,Pi,O)
Xi = GetXi(A,B,Pi,O)
Xi_1 = Xi.sum(0)
a = np.zeros((N,M))
b = np.zeros((M,Y))
pi = np.zeros((1,N))
a_1 = np.subtract(Gamma.sum(1),Gamma[:,H-1]).reshape(1,N)
for i in range(N):
for j in range(M):
a[i,j] = Xi_1[i,j] / a_1[0,i]
#print(a)
for y in range(Y):
for j in range(M):
for h in range(H):
if O[0,h]-1 == y:
c = c + Gamma[j,h]
gamma = Gamma.sum(1).reshape(1,N)
b[j,y] = c / gamma[0,j]
c = 0
#print(b)
for i in range(N):
pi[0,i] = Gamma[i,0]
#print(pi)
return a,b,pi
def BaumWelchAlgo_n(A,B,Pi,O,n):# 计算迭代次数为 n 的 BaumWelch 算法
for i in range(n):
A,B,Pi = BaumWelchAlgo(A,B,Pi,O)
return A,B,Pi
def viterbi(A,B,Pi,O):
N = A.shape[0]# 数组 A 的行数
M = A.shape[1]# 数组 A 的列数
H = O.shape[1]# 数组 O 的列数
Delta = np.zeros((M,H))
Psi = np.zeros((M,H))
Delta_1 = np.zeros((N,1))
I = np.zeros((1,H))
for i in range(N):
Delta[i,0] = Pi[0,i] * B[i,O[0,0]-1]
for h in range(1,H):
for j in range(M):
for i in range(N):
Delta_1[i,0] = Delta[i,h-1] * A[i,j] * B[j,O[0,h]-1]
Delta[j,h] = np.amax(Delta_1)
Psi[j,h] = np.argmax(Delta_1)+1
print("Delta矩阵: \n % r " % Delta)
print("Psi矩阵: \n % r " % Psi)
P_best = np.amax(Delta[:,H-1])
psi = np.argmax(Delta[:,H-1])
I[0,H-1] = psi + 1
for h in range(H-1,0,-1):
I[0,h-1] = Psi[I[0,h]-1,h]
print("最优路径概率: \n % r " % P_best)
print("最优路径: \n % r " % I)"
来源: http://www.bubuko.com/infodetail-2485555.html