马尔科夫系列——三、隐马尔可夫模型 - 学习问题 - Baum-Welch算法
转载的过程中发现，原文有些地方不太理解，就阅读了其他的文章，然后把代码的实现也引进来了。之前并没有验证代码的准确性，后面有人说，代码可能有问题，我尝试了修改。把简单修改版本的也放上来。
目录
一、训练数据包含观测序列和状态序列
1、初始概率的计算
2、转移概率的计算
3、发射概率
二、训练数据中只有观测序列
Baum-Welch算法 - π求解
Baum-Welch算法 - A求解
Baum-Welch算法 - B求解
代码实现
若训练数据包含观测序列和状态序列，则HMM的学习问题非常简单，是监督学习算法。
若训练数据只包含观测序列，则HMM的学习问题需要使用EM算法求解，是非监督学习算法。
一、训练数据包含观测序列和状态序列
直接利用大数定理的结论“频率的极限是概率”，直接给出HMM的参数估计；
下面一一给出解释和说明。
1、初始概率的计算
公式中|Si| : Si的个数，表示状态i的总个数；∑|Si|：表示所有时间点上的状态个数。举例如下图：
2、转移概率的计算
公式中aij 表示i号状态转移到j号状态的概率；|Sij| 表示i号状态转移到j号状态的个数；∑|Sij| 表示所有时间点上i号状态转移到j号状态的个数。举例如下图所示：
注意上图中的公式有点小问题，分母中的求和符号和sij中的j应该修改为其他的字母——如K，这样才没有歧义。
3、发射概率
bij 从当前状态转移到某个观测值的个数。
比如从1号盒子这个状态，转移到取出白球这个观测值的可能性 = b1白；
b1白 = (1->白)的个数 / [(1->白)的个数+(1->黑)的个数]
同样的它的公式中分母中的求和符号和qij中的J应该修改为其他的字母——如k，这样才更好理解。
二、训练数据中只有观测序列
只有观测序列，没有状态序列，这个时候就得采用Baum-Welch算法。本文中的Baum-Welch算法解释写的不是很详细，具体的推导和实现需要参考别的博客——隐马尔可夫模型之Baum-Welch算法详解和”相亲记“之从EM算法到Baum-Welch算法。在Baum-Welch算法的推导过程中需要用到EM算法的原理和HMM概率问题中的前向和后向算法结论。
根据EM算法的E步：得到Q函数，然后对Q函数极大化，
Baum-Welch算法 - π求解
极大化L，使用拉格朗日乘子法，求解π的值：
Baum-Welch算法 - A求解
极大化L，使用拉格朗日乘子法，求解aij的值：
Baum-Welch算法 - B求解
极大化L，使用拉格朗日乘子法，求解bij的值：
Baum-Welch算法 - 极大化L函数，分别可以求得π、a、b的值
其中:
代码实现
来自这篇博文：”相亲记“之从EM算法到Baum-Welch算法，理解代码可以配合博文大学食堂之HMM模型（三）——Baum-Walch算法来进行——我还没有来得及去验证这个代码，以后有机会在验证，先把HMM的理论理解了进入CRF中，CRF才是最重要的，目前对我来说！
#!/usr/bin/env python
# -*- coding:utf-8 -*-
import numpy'''
Created on 2017年2月9日@author: 薛沛雷
'''class HMM:def __init__(self,A,B,Pi):self.A=Aself.B=Bself.Pi=Pi#前向算法def forward(self,O):row=self.A.shape[0]col=len(O)alpha=numpy.zeros((row,col))#初值alpha[:,0]=self.Pi*self.B[:,O[0]]#递推for t in range(1,col):for i in range(row):alpha[i,t]=numpy.dot(alpha[:,t-1],self.A[:,i])*self.B[i,O[t]]#终止return alpha#后向算法def backward(self,O):row=self.A.shape[0]col=len(O)beta=numpy.zeros((row,col))#初值beta[:,-1:]=1#递推for t in reversed(range(col-1)):for i in range(row):beta[i,t]=numpy.sum(self.A[i,:]*self.B[:,O[t+1]]*beta[:,t+1])#终止return beta#前向-后向算法(Baum-Welch算法):由 EM算法 & HMM 结合形成def baum_welch(self,O,e=0.05):row=self.A.shape[0]col=len(O)done=Falsewhile not done:zeta=numpy.zeros((row,row,col-1))alpha=self.forward(O)beta=self.backward(O)#EM算法：由 E-步骤 和 M-步骤 组成#E-步骤：计算期望值zeta和gammafor t in range(col-1):#分母部分denominator=numpy.dot(numpy.dot(alpha[:,t],self.A)*self.B[:,O[t+1]],beta[:,t+1])for i in range(row):#分子部分以及zeta的值numerator=alpha[i,t]*self.A[i,:]*self.B[:,O[t+1]]*beta[:,t+1]zeta[i,:,t]=numerator/denominatorgamma=numpy.sum(zeta,axis=1)final_numerator=(alpha[:,col-1]*beta[:,col-1]).reshape(-1,1)final=final_numerator/numpy.sum(final_numerator)gamma=numpy.hstack((gamma,final))#M-步骤：重新估计参数Pi,A,BnewPi=gamma[:,0]newA=numpy.sum(zeta,axis=2)/numpy.sum(gamma[:,:-1],axis=1)newB=numpy.copy(self.B)b_denominator=numpy.sum(gamma,axis=1)temp_matrix=numpy.zeros((1,len(O)))for k in range(self.B.shape[1]):for t in range(len(O)):if O[t]==k:temp_matrix[0][t]=1newB[:,k]=numpy.sum(gamma*temp_matrix,axis=1)/b_denominator#终止阀值if numpy.max(abs(self.Pi-newPi))<e and numpy.max(abs(self.A-newA))<e and numpy.max(abs(self.B-newB))<e:done=True self.A=newAself.B=newBself.Pi=newPireturn self.Pi#将字典转化为矩阵
def matrix(X,index1,index2):#初始化为0矩阵m = numpy.zeros((len(index1),len(index2)))for row in X:for col in X[row]:#转化m[index1.index(row)][index2.index(col)]=X[row][col]return mif __name__ == "__main__": #初始化,随机的给参数A,B,Pi赋值status=["相处","拜拜"]observations=["撒娇","低头玩儿手机","眼神很友好","主动留下联系方式"] #撒娇:小拳拳捶你胸口A={"相处":{"相处":0.5,"拜拜":0.5},"拜拜":{"相处":0.5,"拜拜":0.5}}B={"相处":{"撒娇":0.4,"低头玩儿手机":0.1,"眼神很友好":0.3,"主动留下联系方式":0.2},"拜拜":{"撒娇":0.1,"低头玩儿手机":0.5,"眼神很友好":0.2,"主动留下联系方式":0.2}}Pi=[0.5,0.5]O=[1,2,0,2,3,0]A=matrix(A,status,status)B=matrix(B,status,observations)hmm=HMM(A,B,Pi)print(hmm.baum_welch(O))
以上代码，最后的结果中P的概率之和为1.A的和B的每行概率之和明显不是1，这里可能就有点问题了。我做了简单的修改，只是修改了一句代码，PAB的各自概率之和差不多都近似为1了。
#!/usr/bin/env python
# -*- coding:utf-8 -*-
import numpy
import numpy as np'''
Created on 2020-04-18
@author: HY
'''class HMM:def __init__(self, A, B, Pi):self.A = Aself.B = Bself.Pi = Pi# 前向算法def forward(self, O):row = self.A.shape[0]col = len(O)alpha = numpy.zeros((row, col))# 初值alpha[:, 0] = self.Pi * self.B[:, O[0]]# 递推for t in range(1, col):for i in range(row):alpha[i, t] = numpy.dot(alpha[:, t - 1], self.A[:, i]) * self.B[i, O[t]]# 终止return alpha# 后向算法def backward(self, O):row = self.A.shape[0]col = len(O)beta = numpy.zeros((row, col))# 初值beta[:, -1:] = 1# 递推for t in reversed(range(col - 1)):for i in range(row):# beta[i, t] = numpy.sum(self.A[i, :] * self.B[:, O[t + 1]] * beta[:, t + 1])beta[i, t] = numpy.dot(beta[:, t + 1],self.A[i,:])* self.B[i,O[t + 1]]# 终止return beta# 前向-后向算法(Baum-Welch算法):由 EM算法 & HMM 结合形成def baum_welch(self, observations, criterion=0.05):n_states = self.A.shape[0]# 观察序列的长度Tn_samples = len(observations)done = Falsewhile not done:# alpha_t(i) = P(o_1,o_2,...,o_t,q_t = s_i | hmm)# Initialize alpha# 获得所有前向传播节点值 alpha_t(i)alpha = self.forward(observations)# beta_t(i) = P(o_t+1,o_t+2,...,o_T | q_t = s_i , hmm)# Initialize beta# 获得所有后向传播节点值 beta_t(i)beta = self.backward(observations)# 计算 xi_t(i,j) -> xi(i,j,t)xi = np.zeros((n_states, n_states, n_samples - 1))# 在每个时刻for t in range(n_samples - 1):# 计算P(O | hmm)denom = sum(alpha[:, -1])for i in range(n_states):# numer[1,:] = 行向量，alpha[i,t]=实数，slef.A[i,:] = 行向量# self.B[:,observations[t+1]].T = 行向量,beta[:,t+1].T = 行向量numer = alpha[i, t] * self.A[i, :] * self.B[:, observations[t + 1]].T * beta[:, t + 1].Txi[i, :, t] = numer / denom# 计算gamma_t(i) 就是对j进行求和gamma = np.sum(xi, axis=1)# need final gamma elements for new Bprod = (alpha[:, n_samples - 1] * beta[:, n_samples - 1]).reshape((-1, 1))# 合并T时刻的节点gamma = np.hstack((gamma, prod / np.sum(prod)))# 列向量newpi = gamma[:, 0]newA = np.sum(xi, 2) / np.sum(gamma[:, :-1], axis=1).reshape((-1, 1))newB = np.copy(self.B)# 观测状态num_levels = self.B.shape[1]sumgamma = np.sum(gamma, axis=1)temp_matrix = numpy.zeros((1, len(observations)))for lev in range(num_levels):for t in range(len(observations)):if observations[t]== lev:temp_matrix[0][t] = 1newB[:, lev] = np.sum(gamma*temp_matrix, axis=1) / sumgammatemp_matrix = numpy.zeros((1, len(observations)))if np.max(abs(self.Pi - newpi)) < criterion and np.max(abs(self.A - newA)) < criterion and np.max(abs(self.B - newB)) < criterion:done = 1self.A[:], self.B[:], self.Pi = newA, newB, newpireturn self.Pi, self.A, self.B# 将字典转化为矩阵
def matrix(X, index1, index2):# 初始化为0矩阵m = numpy.zeros((len(index1), len(index2)))for row in X:for col in X[row]:# 转化m[index1.index(row)][index2.index(col)] = X[row][col]return mif __name__ == "__main__":# 初始化,随机的给参数A,B,Pi赋值status = ["相处", "拜拜"]observations = ["撒娇", "低头玩儿手机", "眼神很友好", "主动留下联系方式"] # 撒娇:小拳拳捶你胸口A = {"相处": {"相处": 0.5, "拜拜": 0.5}, "拜拜": {"相处": 0.5, "拜拜": 0.5}}B = {"相处": {"撒娇": 0.4, "低头玩儿手机": 0.1, "眼神很友好": 0.3, "主动留下联系方式": 0.2},"拜拜": {"撒娇": 0.1, "低头玩儿手机": 0.5, "眼神很友好": 0.2, "主动留下联系方式": 0.2}}Pi = [0.5, 0.5]Ob = [1,2,0,2,3,0]A = matrix(A, status, status)B = matrix(B, status, observations)hmm = HMM(A, B, Pi)p = hmm.baum_welch(Ob)[0]a = hmm.baum_welch(Ob)[1]b = hmm.baum_welch(Ob)[2]print('p:',p)print('a:',a)print('b:',b)最后的结果：
马尔科夫系列——三、隐马尔可夫模型 - 学习问题 - Baum-Welch算法相关推荐
- HMM（马尔科夫过程及隐马尔科夫过程）
转载地址(http://blog.csdn.net/xinzhangyanxiang/article/details/8522078) 学习概率的时候,大家一定都学过马尔科夫模型吧,当时就觉得很有意思 ...
- 隐马尔可夫模型（Baum Welch算法与Viterbi算法）
1.Baum Welch算法就是EM算法,用于求解隐马尔可夫模型的Learing问题 2.隐马尔可夫模型的Decoding问题是指给定X与λ,求使得概率最大的隐状态序列 3.Decoding问题采取V ...
- ONNX系列三 --- 使用ONNX使PyTorch AI模型可移植
目录 PyTorch简介 导入转换器 快速浏览模型 将PyTorch模型转换为ONNX 摘要和后续步骤 参考文献 下载源547.1 KB 系列文章列表如下: ONNX系列一 --- 带有ONNX的便携 ...
- 马氏系列（马尔科夫过程、马尔科夫决策、隐马尔科夫模型）
马尔科夫过程 马尔可夫链就是这样一个任性的过程,它将来的状态分布只取决于现在,跟过去无关! 具体内容参见:https://zhuanlan.zhihu.com/p/26453269 马尔科夫决策 也是 ...
- Python隐形马尔科夫实战_隐马尔科夫模型的Python3实现代码
下面给出计算隐马尔科夫模型的编程代码:from hmmlearn.hmm importGaussianHMMimportdatetimeimportnumpy as npfrom matplotlib ...
- 数学之美 系列三 -- 隐含马尔可夫模型在语言处理中的应用
发表者:吴军,Google 研究员 前言:隐含马尔可夫模型是一个数学模型,到目前为之,它一直被认为是实现快速精确的语音识别系统的最成功的方法.复杂的语音识别问题通过隐含马尔可夫模型能非常简单地被表述. ...
- 学术报告系列(四) - 端云协同下分布式模型学习与进化
- 20200918：【广发金融工程】2018年重磅专题系列之九：基于隐马尔科夫模型的选股策略研究
参考链接:https://www.sohu.com/a/252454782_465470 [广发金融工程]2018年重磅专题系列之九:基于隐马尔科夫模型的选股策略研究 2018-09-07 11:26 ...
- 【深度学习】隐马尔科夫
文章目录 一 隐马尔可夫 1. 自动机 2. 马尔可夫链和马尔可夫假设 3. 隐马尔可夫模型 3.1 马尔科夫模型中的三个问题: 3.2 似然度问题: 3.3 解码问题 案例 如何对句子进行分词 统计 ...
最新文章
- 江苏“超牛”女博导：16岁考入北大，跨专业读博，成为全球第4个获此奖项的学者！...
- objective-c对NSArray的学习
- 用Visio制作周行事历
- oracle 存储过程(分析理解)
- 语音怎么进入滤波器matlab,基于Matlab的语音信号滤波器的设计与实现
- 数据可视化【七】 更新模式
- JUnit 4和JUnit 5区别
- linux路由内核实现分析(二)---FIB相关数据结构(4)
- （转）OpenLayers3基础教程——OL3 介绍interaction
- 音视频开发（30）---新一代视频编码标准：HEVC、AVS2和AV1性能对比报告
- leetcode题库1--两数之和
- MySql常用函数大全（详细）
- 免费分享9本经典的MySQL书籍。
- 软考初级程序员真题资料（2009年上半年——2019年上半年）
- 计算机开放电子书 2016 归档
- Android 不得不知Tag的妙用
- day 19 C# 窗体金额换算
- Docker入门实战(三)-Docker容器镜像
- 递归算法的原理(js)
- ACM 海贼王之伟大航路(深搜剪枝)