|
隐马尔科夫作为最经典的机器学习方法之一，是动态规划的经典工具，衍生出了Ngram等诸多模型，对其方法原理的理解非常重要。本文可以说是最详尽的一篇关于隐马尔科夫方法的解读了，就是比较长，可以收藏下来慢慢看。
目录 1.基本的概念 1.1 知道骰子是什么（知道参数是什么），每种骰子是什么，每次投的是什么骰子，根据骰子的结果求产生这个结果的概率 1.2 知道骰子有几种（隐含状态数量），也就是有几种骰子，在上面的例子中就有三种，每种骰子是什么（转换概率），根据扔骰子投出的结果（观测状态链），我想知道每次投的是哪种骰子（隐含状态链） 1.3 知道骰子有几种（隐含状态的数量），每种骰子是什么（转换概率），观测到很多次骰子的结果（观测序列），想知道投出这个结果的概率 1.4 知道骰子有几种（隐含状态的数量），不知道每种骰子是什么，检测到很多次骰子的结果（观测序列），我现在想知道每种骰子是什么 2.公式推导
2.1 HMM模型的确定 2.2 三个问题 2.3 暴力求解 2.4 前向算法 2.5 后向算法 2.6 前向概率和后向概率的关系 2.7 学习算法 2.8 baum welch算法 2.9 预测算法 3.代码实现：中文分词案例 3.1 有监督学习案例 3.2 无监督学习案例 首先举一个比较经典的例子，扔骰子：有三个骰子，第一个是我们最常见的6个面，（1,2,3,4，5,6）每一个面概率1/6；第二个只有4个面，（1,2,3,4）每一个面概率1/4；第三个有8个面，（1,2,3,4,5,6,7,8），每一个面出现的概率1/8。 
假如我们先挑一个骰子，挑到每一个骰子的概率就是1/3，然后投掷，可以得到数字中1,2,3,4,5,6,7,8中的一个。不停重复上述的过程，最后会得到一串数字，例如我们可能得到一串这样的数字，13725857356。这种数字叫可见状态序列，但在隐马尔科夫模型中，我们不仅仅有这么一串可见序列，还有一串不可见的隐含状态链。比如这个隐含状态链可能是：D6 D8 D8 D6 D4 D8 D6 D6 D4 D8
一般来说，HMM说到的马尔科夫链其实就是隐含的状态序列，而上面的13725857356这些就叫做观测序列。而隐含状态之间其实是存在有转换概率的。在我们这个例子里面，D6的下一个状态是D4，D6，D8的概率都是1/3。这样的设定是为了最开始容易说清楚，但是我们其实是可以随意设定转换概率的。比如，D6后面不能接D4，D4到D6概率是0.1等等。这样就是又是一个新的HMM了。 同样的，尽管可见的状态之间没有转换概率，但是隐含状态和可见状态之间有一个概率叫做输出概率，也叫发射概率，就这个例子来说，D6输出1的概率就是6了，产生23456的概率都是1/6。 示意图： 
隐含状态关系图： 
其实如果知道转换概率和所有隐含的概率做模拟是相当容易的。但有时候总是会缺失一些信息，比如知道了观测序列和参数，但是不知道这个观测序列的隐含序列是什么，或者是已知观测序列，想要估计参数，又或者是已知参数和观测序列，想要知道这个观测序列出现的概率是什么。其实这些基本问题就是对应后面我们要解决的HMM3个基本问题，后面都是使用python代码进行实现。 三个基本问题：
1.知道骰子有几种（隐含状态数量），也就是有几种骰子，在上面的例子中就有三种，每种骰子是什么（转换概率），根据扔骰子投出的结果（观测状态链），我想知道每次投的是哪种骰子（隐含状态链）。 这个问题在语音识别领域叫解码问题。这个问题其实有两种解法。 第一种是求最大似然的状态路径，其实就是找一串骰子序列，使得这串骰子产生观测结果的概率最大。第二种就是求出每个序列属于某种骰子的概率。 2.知道骰子有几种（隐含状态的数量），每种骰子是什么（转换概率），观测到很多次骰子的结果（观测序列），想知道投出这个结果的概率。 这个问题看起来其实意义不太大，因为一般结果都会对应着最大的概率，但是这个问题的目的其实是为了确认这个模型和结果是不是吻合的，这样就可以知道模型是否出错了。 3.知道骰子有几种（隐含状态的数量），不知道每种骰子是什么，检测到很多次骰子的结果（观测序列），我现在想知道每种骰子是什么。 这个问题是最重要的问题了，后面的分词会用到解决这个方法。 要解决上面的问题，其实还有一个0号问题要解决： 0.知道骰子是什么（知道参数是什么），每种骰子是什么，每次投的是什么骰子，根据骰子的结果求产生这个结果的概率。 这个问题的解决其实就是概率相乘： 
回到第一个问题：
1.求骰子的序列，已知参数，求隐含序列。 这里我说的是第一种解法，解最大似然路径问题。
举例来说，我知道我有三个骰子，六面骰，四面骰，八面骰。我也知道我掷了十次的结果（1 6 3 5 2 7 3 5 2 4），我不知道每次用了那种骰子，我想知道最有可能的骰子序列。 其实最简单的方法就是穷举所有的方法，从第一个开始求概率，把最大的概率算出来即可，如果马尔科夫链不长还行，要是长的就很难算出来了。 这里就要使用到Viterbi algorithm算法，其实就是动态规划的思想： 首先从第一个开始看： 
看到的的结果是1，那么很大的概率就是D4了，因为他投出的概率是1/4，是最大的。
把这个情况拓展，扔两次骰子： 
这时结果就是1,6了，这个时候还是和上面的一样要计算他们之间的最大概率。第一个骰子为D4，第二个骰子为D6的概率为： 
同样的，我们可以计算第二个骰子是D4或D8时的最大概率。我们发现，第二个骰子取到D6的概率最大。而使这个概率最大时，第一个骰子为D4。所以最大概率骰子序列就是D4 D6。
继续拓展，投第三次骰子： 同样，我们要计算的三个骰子分别是D6,D4,D8的最大概率，要取到最大的概率，第二个骰子一定是D6，所以取到最大概率是： 
同上，我们可以计算第三个骰子是D6或D8时的最大概率。我们发现，第三个骰子取到D4的概率最大。而使这个概率最大时，第二个骰子为D6，第一个骰子为D4。所以最大概率骰子序列就是D4 D6 D4。
写到这里大概就知道怎么回事了，我们发现我们要求最大概率骰子序列要做那么几件事情。首先，不管序列多长，要从序列长度为1算一遍在这个长度为1时取到每一个骰子的最大概率，因为上一个长度下的取到每一个骰子的最大概率都算过了，重新计算的话其实不难。当算到最后一位的时候，就已经可以知道最后一位是哪一个骰子的概率最大了，然后再把对应的鸟哥最大概率的序列从后面往前反推就好了。 2.谁动了我的骰子 比如你怀疑模型的不正确或者是因为骰子被动了手脚，可以这个骰子投到1的概率更加大。那么你就可以算一算正常的三个骰子投出的一段序列的概率，再算一下不正常的序列的概率，如果前者比后者小，那么就应该有问题了。 比如投骰子的结果是： 
要算用正常的三个骰子投出这个结果的概率，其实就是将所有可能的情况的概率进行加和。同样，简单暴力的方法就是穷举所有骰子的序列，还是计算每一个骰子对应的概率，但是这回我们不挑最大值了，而是找到所有的概率进行相加和，得到的总概率就是我们要求的结果了。这个方法依然不能用于太长的马尔科夫链。 我们会应用到和前面一个问题类似的解决方法。前一个问题关心的是最大值，而这个问题关心的是总和，也就是前向算法了。 首先，如果要只投一个骰子： 
每一个骰子都可以投出1这个值，所以有三种：
再拓展一下情况，投掷两次骰子：
看到结果为1,6,。产生这个结果的概率可以如下计算：
上一个状态转移到当前这个状态的概率。因为这里都是1/3，所以没有什么不同，之后再乘上发射概率。所以一共就0.05。 继续，这次就要投三个骰子了： 同样的，这样往下算，多长的都能算出来。比较一下就知道了。 3.知道扔出的是是什么数字（观测序列），想要知道是什么骰子。 这个要用到baum_welch算法，在后面的公式推导中会详细讲解。 =========================================================================================================================================================================== ②下面就是公式的推导了： HMM马尔科夫模型的确定：
初始分布pi，其实就是刚刚最开始第一次选中骰子的概率1/3。状态转移概率其实就是从第一个骰子到选中第二个骰子的概率，观测概率其实就是这个骰子投到1,2,3,4,5这些的概率分布。一起就叫做是lamda参数。而还要给出隐含状态的集合，在刚刚的例子中就是三了，因为有三个骰子；和所有可能观测到的集合，就是1,2,3,4,5,6,7,8了。我们记为Q和V：
这里的I是隐含状态序列，O是观测序列，别和上面的Q和V集合搞混了。A就是从i状态到j状态的概率了。A[1][2]就是从1状态转移到2状态的概率了。观测概率矩阵。
隐马尔科夫其实还有两个条件：
其次假设，当前的隐含状态只与前一次的隐含状态相关，观测独立性：当前时刻的观测状态只与当前时刻相关。 然后就是三个问题了：
给定参数，试计算出现观测序列O的概率是多少。这个问题有三种做法：暴力求解，前向算法，后向算法。 暴力求解： 这个算法其实没有什么卵用，只是用于理解整个模型而已：
所以，由上面的公式就可以退出：
这就是暴力求解的过程。但是上面也说过了，这里如果马尔科夫链过长的话，完不了的，所以出现了一个前向算法和后向算法，本质上就是动态规划。
前向算法：
求在给定t时刻，出现有观测序列O1,O2,O3...Ot的概率是多少。
这个算法其实很好理解，刚刚也在骰子那个例子讲过的了。
这里就不过多阐述了。
后向算法：
给定T时刻，出现Qt+1,Qt+2...的概率。 这个我直观理解不了，但是我用公式推了一下：
接下来讲一些其他的内容：
前向概率和后向概率的关系： 这个内容会在后面用到。
求给定参数和观测序列，t时刻处于qi，t+1时刻处于qj的概率。
要注意是已知了所有的观测序列，从上图其实就可以知道了：
自然也可以得到：
学习算法：
这里要用到EM算---Baum-welch算法，无监督学习。
首先回顾一下EM算法： 上一篇博客其实是讲到了EM算法，讲的其实是通俗角度和正式角度的一种结合，由于我没有看统计学习方法，我并不知道那个是通俗的理解。所以今天重新讲一下：
通俗的角度：
求极大，肯定就是极大似然估计了，都是在给定数据的情况下，求解使得似然函数最大的参数的取值。
通常的做法是对似然函数求偏导，然后令偏导等于零，参数取得的数值就是近似最优值。但是，有些含有隐变量的模型没办法直接进行似然函数的偏导，但是如果假设已经知道隐变量的值，就可以将似然函数简化进行下一步的求偏导。
框架如下：
E:求观测数据下隐变量的期望。
M:求经过隐变量改写的极大似然函数的极大。
正式的角度：
在统计学习方法里面，有一个函数叫Q函数，EM算法其实就是通过求解下界极大值来逼近的，所以会存在局部极大值和对初值敏感。
从形式上看，Q函数是完全数据的对数似然函数关于在给定观测数据和当前参数下对未观测数据的条件概率分布的期望。
框架：
E:求Q函数 M:求给定Q函数取极大值的时候的参数
在使用Jensen不等式的时候，需要假设隐变量服从某种形式的概率分布，才可以将推导过程的一部分看成是期望的表达形式从而应用Jensen不等式。然而这个分布不是随便指定的。我们令Jensen不等式取等号的时候，可以计算出这个分布其实就是：已知观测数据的隐变量的后验概率分布。由于求Q函数需要先求出隐变量的后验概率的期望，因此，这就可以解释为什么EM算法的“通俗”理解角度的E步骤是求隐变量的期望了。
有时候在用EM算法解决某个具体问题的时候，会发现M步骤极大化的居然是完全数据的对数似然函数。这是因为，Q函数虽然是完全数据的对数似然函数的某种期望，但是求这个期望的过程有时其实就是将隐变量的后验概率的期望代入就可以了。因此，本质上我们其实还是在求Q函数的极大。
所以，我们在求GMM混合高斯模型的时候，第一步其实就是求隐变量的期望，为的就是代入似然函数使其变成Q函数。
最后其实可以变成这个步骤。
baum welch算法
最后可以写成这样
可以看到，其实三个要求的参数位于三个不同的部分中。
求pi
求A
这里是求和aij = 1
同样拉格朗日有：
求B
同样的求法：
这样就出结果了。 预测算法：
Viterbi算法本质上就是动态规划思想。
在上面骰子的例子中已经详细解释了，下面代码再详细解答了。
===========================================================================================================================================================================
③代码实现
使用的例子是中文分词
忘了说了，HMM还有一个监督学习算法，非常简单：
直接计数算概率就好了。这里使用分词的参数预测就是这个算法。
隐状态，B开始，M中间，E结束，S单个。
首先是一些工具的准备：
import matplotlib.pyplot import math import codecs import random as r infinite = float(-2**31) def log_normalize(a): s = 0 for x in a: s += x if s == 0: print ('log标准化出错') return s = math.log(s) for i in range(len(a)): if a[i] == 0: a[i] = infinite else: a[i] = math.log(a[i]) - s
上面是一个log的归一化的过程，我们使用log来表示各个概率。首先传进来的数字都是一些随机数，而不是概率，一般来说归一化操作是求和然后各个数字除总和即可。
比如这里有a,b,c三个数字，s = math.log(s)这步其实就是得到log(a + b + c)。而下面的减号，其实就是除法，
得到log(a/(a+b+c))。这样完成归一化操作的。 监督算法得到参数： 首先要有一个人工已经标注好的数据：
这个数据是用jieba分词器分词后写入的文本里面。
完整的监督学习算法。
def train(): pi = [0] * 4 a = [ [0] * 4 for x in range(4) ] b = [ [0] * 65535 for x in range(4) ] f = open('.\\pku_training.utf8' , encoding='UTF-8') data = f.read()[3:] f.close() tokens = data.split(' ') #开始训练 设置第一个状态为结束状态 last_q = 2 old_process = 0 print('进度：') for k , token in enumerate(tokens): process = float(k) / float(len(tokens)) if process > old_process + 0.1: print ('%.3f%%' % (process * 100)) old_process = process token = token.strip() n = len(token) if n <=>0: continue #如果只是一个词 if n == 1: pi[3] += 1 #当前的结束状态有一个 a[last_q][3] += 1 #上一个状态到这里的有一个 b[3][ord(token[0])] += 1 #这个状态到当前词的又有一个 last_q = 3 continue #初始向量 #如果不是，则证明上一次的全部都是了 pi[0] += 1 pi[2] += 1 pi[1] += (n-2) #转移矩阵 a[last_q][0] += 1 last_q = 2 if n == 2: a[0][2] += 1 else: a[0][1] += 1 a[1][1] += (n-3) a[1][2] += 1 #发射矩阵 b[0][ord(token[0])] += 1 b[2][ord(token[n-1])] += 1 for i in range(1, n-1): b[1][ord(token[i])] += 1 # 正则化 log_normalize(pi) for i in range(4): log_normalize(a[i]) log_normalize(b[i]) return pi, a, b
只讲主要部分：
if n <=>0: continue #如果只是一个词 if n == 1: pi[3] += 1 #当前的结束状态有一个 a[last_q][3] += 1 #上一个状态到这里的有一个 b[3][ord(token[0])] += 1 #这个状态到当前词的又有一个 last_q = 3 continue
单个词，说明是S，单个状态，是第三号状态，所以pi[3]号首先加一，后面的状态转移语句也要加1，last_q是上一个状态的，一开始设为2，代表是结束E,2状态。完了之后要把状态设置为当前的单个状态。
#初始向量 #如果不是，则证明上一次的全部都是了 pi[0] += 1 pi[2] += 1 pi[1] += (n-2) #转移矩阵 a[last_q][0] += 1 last_q = 2 if n == 2: a[0][2] += 1 else: a[0][1] += 1 a[1][1] += (n-3) a[1][2] += 1 #发射矩阵 b[0][ord(token[0])] += 1 b[2][ord(token[n-1])] += 1 for i in range(1, n-1): b[1][ord(token[i])] += 1
如果不是呢？说明这个词是多个了，大于2个了，首先结束，开始一个，剩下的都是中间M了，所以pi是这样。转移矩阵开始是肯定有一个的，因为大于1了。上一个状态开始有一个，加一，这一次结束之后，肯定是E，结束状态，可以再最后设置，但是这里提前了而已，然后如果是两个，那么直接第二个结束，而如果是两个，那么肯定第一个就是开始0状态了，所以是a[0][2]，下面的也是一样了，如果不是开始到中间加一，中间到中间不知道多少个，中间到结束加一。下面发射矩阵的也是一样。之后是log的归一化了。 得到了pi,A,B之后就是维特比算法来预测了： def viterbi(pi , A , B , o): T = len(o) delta = [ [0 for i in range(4)] for t in range(T) ] pre = [ [0 for i in range(4)] for t in range(T)] for i in range(4): delta[0][i] = pi[i] + B[i][ord(o[0])] for t in range(1 , T): for i in range(4): delta[t][i] = delta[t-1][0] + A[0][i] for j in range(1 , 4): vj = delta[t-1][j] + A[j][i] if delta[t][i] <> delta[t][i] = vj pre[t][i] = j delta[t][i] += B[i][ord(o[t])] decode = [-1 for t in range(T)] # 解码：回溯查找最大路径 q = 0 for i in range(1, 4): if delta[T-1][i] > delta[T-1][q]: q = i decode[T-1] = q for t in range(T-2, -1, -1): q = pre[t+1][q] decode[t] = q return decode pass
都是按照上面公式来的，不过要记住在log里，乘是加，减是除法。 for i in range(4): delta[0][i] = pi[i] + B[i][ord(o[0])]
得到第一个delta，用于和其他的比较： for t in range(1 , T): for i in range(4): delta[t][i] = delta[t-1][0] + A[0][i] for j in range(1 , 4): vj = delta[t-1][j] + A[j][i] if delta[t][i] <> delta[t][i] = vj pre[t][i] = j delta[t][i] += B[i][ord(o[t])]
得到接下来的1到T个
计算最大的一个，另外要保存当前最大的一个，因为要回溯寻找最大的。最后返回一个序列。
之后就是解码显示了： def segment(sentence , decode): N = len(sentence) i = 0 while i <> if decode[i] == 0 or decode[i] == 1: j = i+1 while j <> if decode[j] == 2: break j += 1 print (sentence[i:j+1], '|', end=' ') i = j+1 elif decode[i] == 3 or decode[i] == 2: # single print (sentence[i:i + 1], '|', end=' ') i += 1 else: print ('Error:', i, decode[i] , end=' ') i += 1
这是用于显示分词效果的。 运行： if __name__ == '__main__': pi, A, B = train() save_parameter(pi, A, B) #A , B , pi = read_lamda() f = open('novel.txt' , encoding='UTF-8') data = f.read()[3:] f.close() decode = viterbi(pi, A, B, data) segment(data, decode) print('训练完成...') read_lamda()
最后的结果：
这个效果其实可以，但这只是监督学习的。
维特比和监督的实现讲完了。 ================================================================= 无监督学习的效果不好，甚至可以说是惨不忍睹，不想看的可以不看啦！ ================================================================= 然后就是无监督学习的实现了： 工具的准备其实和上面的差不多，但是要多一个相加的工具： def log_sum(a): if not a: # a为空 return infinite m = max(a) s = 0 for t in a: s += math.exp(t-m) return m + math.log(s)
log求和，这个应该不难看懂。 计算alpha的工具： def calc_alpha(pi , A , B , o , alpha): print('开始计算alpha') for i in range(4): alpha[0][i] = pi[i] + B[i][ord(o[0])] T = len(o) temp = [ 0 for i in range(4) ] del i for t in range(1 , T): for i in range(4): for j in range(4): temp[j] = (alpha[t-1][j] + A[j][i]) alpha[t][i] = log_sum(temp) alpha[t][i] += B[i][ord(o[t])] print('结束计算alpha') pass
根据这个公式
为什么要计算呢？因为：
为什么要计算呢？因为：
计算beta：
def calc_beta(pi , A , B , o , beta): print('开始计算beta') T = len(o) for i in range(4): beta[T-1][i] = 1 temp = [ 0 for i in range(4) ] del i for t in range(T-2, -1, -1): for i in range(4): beta[t][i] = 0 for j in range(4): temp[j] = A[i][j] + B[j][ord(o[t+1])] + beta[t+1][j] beta[t][i] += log_sum(temp) print('结束计算beta') pass
根据公式：
计算gamma： def calc_gamma(alpha , beta , gamma): print('开始计算gamma') for t in range(len(alpha)): for i in range(4): gamma[t][i] = alpha[t][i] + beta[t][i] s = log_sum(gamma[t]) for i in range(4): gamma[t][i] -= s print('结束计算gamma')
根据计算公式：
计算ksi工具：
def calc_ksi(alpha , beta , A , B , o , ksi): print('开始计算ksi') T = len(alpha) temp = [0 for x in range(16)] for t in range(T - 1): k = 0 for i in range(4): for j in range(4): ksi[t][i][j] = alpha[t][i] + A[i][j] + B[j][ord(o[t+1])] + beta[t+1][j] temp[k] = ksi[t][i][j] k += 1 s = log_sum(temp) for i in range(4): for j in range(4): ksi[t][i][j] -= s print('结束计算ksi')
根据公式：
baum_welch算法更新：
def bw(pi, A, B, alpha, beta, gamma, ksi, o): print('开始更新') T = len(alpha) for i in range(4): pi[i] = gamma[0][i] s1 = [0 for x in range(T-1)] s2 = [0 for x in range(T-1)] for i in range(4): for j in range(4): for t in range(T-1): s1[t] = ksi[t][i][j] s2[t] = gamma[t][i] A[i][j] = log_sum(s1) - log_sum(s2) s1 = [0 for x in range(T)] s2 = [0 for x in range(T)] for i in range(4): for k in range(65536): if (k % 5000 == 0): print(i , k) valid = 0 for t in range(T): if ord(o[t]) == k: s1[valid] = gamma[t][i] valid += 1 s2[t] = gamma[t][i] if valid == 0: B[i][k] = -log_sum(s2) # 平滑 else: B[i][k] = log_sum(s1[:valid]) - log_sum(s2) print('结束更新')
都是根据公式来的。。。根据没有什么好说的。
包装一下：
def baum_welch(pi, A, B): f = open('novel.txt' , encoding='UTF-8') sentence = f.read()[3:] f.close() T = len(sentence) # 观测序列 alpha = [[0 for i in range(4)] for t in range(T)] beta = [[0 for i in range(4)] for t in range(T)] gamma = [[0 for i in range(4)] for t in range(T)] ksi = [[[0 for j in range(4)] for i in range(4)] for t in range(T-1)] for time in range(100): calc_alpha(pi, A, B, sentence, alpha) # alpha(t,i):给定lamda，在时刻t的状态为i且观测到o(1),o(2)...o(t)的概率 calc_beta(pi, A, B, sentence, beta) # beta(t,i)：给定lamda和时刻t的状态i，观测到o(t+1),o(t+2)...oT的概率 calc_gamma(alpha, beta, gamma) # gamma(t,i)：给定lamda和O，在时刻t状态位于i的概率 calc_ksi(alpha, beta, A, B, sentence, ksi) # ksi(t,i,j)：给定lamda和O，在时刻t状态位于i且在时刻i+1，状态位于j的概率 bw(pi, A, B, alpha, beta, gamma, ksi, sentence) #baum_welch算法 save_parameter(pi, A, B) print('保存完毕')
可以直接运行的总代码： import math import matplotlib.pyplot as plt import numpy as np import codecs import random infinite = float(-2**31) def log_normalize(a): s = 0.0 for x in a: s += x s = math.log(s) for i in range(len(a)): if a[i] == 0: a[i] = infinite else: a[i] = math.log(a[i]) - s def log_sum(a): if not a: # a为空 return infinite m = max(a) s = 0 for t in a: s += math.exp(t-m) return m + math.log(s) def calc_alpha(pi , A , B , o , alpha): print('开始计算alpha') for i in range(4): alpha[0][i] = pi[i] + B[i][ord(o[0])] T = len(o) temp = [ 0 for i in range(4) ] del i for t in range(1 , T): for i in range(4): for j in range(4): temp[j] = (alpha[t-1][j] + A[j][i]) alpha[t][i] = log_sum(temp) alpha[t][i] += B[i][ord(o[t])] print('结束计算alpha') pass def calc_beta(pi , A , B , o , beta): print('开始计算beta') T = len(o) for i in range(4): beta[T-1][i] = 1 temp = [ 0 for i in range(4) ] del i for t in range(T-2, -1, -1): for i in range(4): beta[t][i] = 0 for j in range(4): temp[j] = A[i][j] + B[j][ord(o[t+1])] + beta[t+1][j] beta[t][i] += log_sum(temp) print('结束计算beta') pass def calc_gamma(alpha , beta , gamma): print('开始计算gamma') for t in range(len(alpha)): for i in range(4): gamma[t][i] = alpha[t][i] + beta[t][i] s = log_sum(gamma[t]) for i in range(4): gamma[t][i] -= s print('结束计算gamma') def calc_ksi(alpha , beta , A , B , o , ksi): print('开始计算ksi') T = len(alpha) temp = [0 for x in range(16)] for t in range(T - 1): k = 0 for i in range(4): for j in range(4): ksi[t][i][j] = alpha[t][i] + A[i][j] + B[j][ord(o[t+1])] + beta[t+1][j] temp[k] = ksi[t][i][j] k += 1 s = log_sum(temp) for i in range(4): for j in range(4): ksi[t][i][j] -= s print('结束计算ksi') def bw(pi, A, B, alpha, beta, gamma, ksi, o): print('开始更新') T = len(alpha) for i in range(4): pi[i] = gamma[0][i] s1 = [0 for x in range(T-1)] s2 = [0 for x in range(T-1)] for i in range(4): for j in range(4): for t in range(T-1): s1[t] = ksi[t][i][j] s2[t] = gamma[t][i] A[i][j] = log_sum(s1) - log_sum(s2) s1 = [0 for x in range(T)] s2 = [0 for x in range(T)] for i in range(4): for k in range(65536): if (k % 5000 == 0): print(i , k) valid = 0 for t in range(T): if ord(o[t]) == k: s1[valid] = gamma[t][i] valid += 1 s2[t] = gamma[t][i] if valid == 0: B[i][k] = -log_sum(s2) # 平滑 else: B[i][k] = log_sum(s1[:valid]) - log_sum(s2) print('结束更新') def save_parameter(pi, A, B): f_pi = open('pi.txt', 'w') list_write(f_pi, pi) f_pi.close() f_A = open('A.txt', 'w') for a in A: list_write(f_A, a) f_A.close() f_B = open('B.txt', 'w') for b in B: list_write(f_B, b) f_B.close() def list_write(f, v): for a in v: f.write(str(a)) f.write(' ') f.write('\n') def baum_welch(pi, A, B): f = open('novel.txt' , encoding='UTF-8') sentence = f.read()[3:] f.close() T = len(sentence) # 观测序列 alpha = [[0 for i in range(4)] for t in range(T)] beta = [[0 for i in range(4)] for t in range(T)] gamma = [[0 for i in range(4)] for t in range(T)] ksi = [[[0 for j in range(4)] for i in range(4)] for t in range(T-1)] for time in range(100): calc_alpha(pi, A, B, sentence, alpha) # alpha(t,i):给定lamda，在时刻t的状态为i且观测到o(1),o(2)...o(t)的概率 calc_beta(pi, A, B, sentence, beta) # beta(t,i)：给定lamda和时刻t的状态i，观测到o(t+1),o(t+2)...oT的概率 calc_gamma(alpha, beta, gamma) # gamma(t,i)：给定lamda和O，在时刻t状态位于i的概率 calc_ksi(alpha, beta, A, B, sentence, ksi) # ksi(t,i,j)：给定lamda和O，在时刻t状态位于i且在时刻i+1，状态位于j的概率 bw(pi, A, B, alpha, beta, gamma, ksi, sentence) #baum_welch算法 save_parameter(pi, A, B) print('保存完毕') def inite(): # 初始化pi,A,B pi = [random.random() for x in range(4)] # 初始分布 log_normalize(pi) A = [[random.random() for y in range(4)] for x in range(4)] # 转移矩阵：B/M/E/S A[0][0] = A[0][3] = A[1][0] = A[1][3]\ = A[2][1] = A[2][2] = A[3][1] = A[3][2] = 0 # 不可能事件 B = [[random.random() for y in range(65536)] for x in range(4)] for i in range(4): log_normalize(A[i]) log_normalize(B[i]) return pi , A , B def run(): pi , A , B = inite() print('初始化完成') baum_welch(pi , A , B) pass if __name__ == '__main__': run()
当然，数据要自己找了。 最后把保存的A,B,pi用进来，使用维特比算法预测，效果很惨。
跑了40次跌代了，已经4个小时了，电脑不行。
|
