#HMM隐马尔科夫模型 ###①通俗的理解 首先举一个例子，扔骰子，有三种骰子，第一个是比较常见的6个面，每一个面的概率都是1/6。第二个只有4个面，，每一个面的概率是1/4。第三个有8个面，，每一个面的概率是1/8。

首先先选择一个骰子，挑到每一个骰子的概率1/3，然后投掷，可以得到。最后会得到一堆的序列，比如会得到等等，这种序列叫可见状态序列，但在HMM里面，还存在一个隐含状态链，比如这个状态链可能是。从8个面出来的，6个面投出来的。所以隐马尔科夫模型里面有两串序列，可观测序列和隐含状态序列。 一般来讲，说到隐马尔科夫链其实就是隐含的状态序列了，markov链各个元素直接是存在了转化关系的，比如一开始是D6下一个状态是D4，D6，D8的概率都是1/3，这样设置平均值只是为了初始化而已，一般这种值其实是可以随意改变的，比如这里下一个状态序列是D8，那么就可以设置转化到D6的概率是0.1，D4的概率是0.1，D8的概率是0.8。这样就又是一个很新的HMM了。 同样的，尽管可见的状态之间没有转换概率，但是隐含状态和可见状态之间有一个概率叫做输出概率，也叫发射概率，就这个例子来说，D6输出1的概率就是6了，产生23456的概率都是1/6。 所以，HMM首先是有两个序列，可观测序列，隐含序列。每一个隐含序列里面的元素直接是存在着转化概率的，也就是用概率来描述转化到下一个状态的可能性；而隐含序列和观测序列之间也有一个概率，发射概率，也叫输出概率。 隐含状态关系图： 每一个箭头都会对应着一个概率，转换概率。 如果知道了转换概率和隐含概率那么求解观测序列会很简单，但如果是这样那么这个模型就没有什么可以研究的了，所以肯定是会缺失一点值的，比如知道观测序列和参数，但是不知道隐含序列是什么，或者是知道了 观测序列想求参数，又或者是知道了这个参数和隐含序列求这个观测序列出现的参数是什么。其实这些就对应了后面的三个问题。HMM的核心其实也就是这三个问题。 ###②三个基本问题通俗解释 1.知道骰子有几种（知道隐含序列有多少种），也就是几种骰子，上面就有三种，每种骰子是什么（转换概率），根据骰子扔出的结果（可观测状态链），想知道每次投的是哪种骰子（隐含状态序列）？ 这问题其实就是解码问题了，知道观测序列求隐含序列。其实有两种解法，一种就是暴力求解：其实就是maximum likelihood全部乘起来一起算即可。另一种就厉害带了，递推的方法，前向算法，后向算法，前向-后向算法。 2.还是知道骰子有几种隐含状态（隐含状态的数量），知道骰子是什么（知道转换概率），根据骰子投出的结果我想知道投出这个结果的概率是多少？ 这个问题看起来貌似没有什么太大的作用，但是实际上是用于对于模型的检测，如果我们投出来的结果都对应了很小的概率，那么模型可能就是错误的，有人换了我们的骰子。 3.知道骰子有几种（知道隐含序列的种类），不知道每种骰子是什么（不知道转换概率），实验得到多次骰子的结果（可观测序列），现在想知道每种骰子是什么（转换概率） 这很重要，后面会用到来求解分词的状态。 到这里先引出一个比较简单的问题，0号问题： 0.知道骰子的种类，骰子是什么，每次扔什么骰子，根据这个结果，求投出这个结果的概率是多少。 其实就是该知道的都知道了，直接求概率。 求解无非就是概率相乘了：

###③三个问题的简单解答 ####1.看见了观测序列，求隐藏序列 这里解释的就是第一种解法了，最大似然估计，也就是暴力解法。首先我知道我有3个骰子，六面骰子，八面骰子，四面骰子，同时结果我也知道了（1 6 3 5 2 7 3 5 2 4），现在想知道我每一次投的哪一个骰子，是六面的还是四面的还是八面的呢？ 最简单的方法就是穷举了，从零个一直到第最后一个一次把每一个概率算出来，链不长还行，链要是长了就算不了了，穷举的数量太大。 接下来是要讨论另一个比较牛逼的算法，viterbi algorithm。 首先只扔一次骰子：

结果是1的话那么概率最大的就是四面骰子了，1/4，其他的都是1/6和1/8。接着再扔一次：这个时候我们就要计算三个值了，分别是D6，D8，D4的概率。要取到最大的概率，那么第一次肯定就是取到D4了，所以取到D6的最大概率：$$P_2(D6) = P(D4)P(D4->1)P(D6)P(D6->6) = \frac{1}{3}\frac{1}{4}\frac{1}{3}\frac{1}{6}$$ 取到D8的最大概率：$$P_2(D8) = P(D4)P(D4->1)P(D8)P(D8->6) = \frac{1}{3}\frac{1}{4}\frac{1}{3}\frac{1}{8}$$取到D4这辈子都不可能了，四面骰子没有6。所以最大的序列就是D4，D6 继续拓展一下，投多一个：这个时候又要计算三种情况了，分别计算为D4，D6，D8的概率是多少：

可以看到，概率最大的就是D4了，所以三次投骰子概率最大的隐含序列就是D4，D6，D4。所以无论这个序列多长，直接递推计算可以算出来的，算到最后一位的时候，直接反着找就好了。 ####2.求可观测序列的概率 如果你的骰子有问题，要怎么检测出来？首先，可以先算一下正常骰子投出一段序列的概率，再算算不正常的投出来序列的概率即可。如果小于正常的，那么就有可能是错的了。

比如结果如上，我们这个时候就不是要求隐含序列了，而是求出现这个观测序列的总概率是多少。首先如果是只投一个骰子：这时候三种骰子都有可能：再投一次：同样是要计算总的分数。再投一次： 就是这样按照套路计算一遍再比较结果即可。 ####3.只是知道观测序列，想知道骰子是怎么样的？ 这个算法在后面会细讲，Baum-welch算法。 ###④HMM的公式角度 下面正式开始讲解，上面只是过一个印象。 ####1.马尔科夫模型 要看隐马尔科夫自然先动马尔科夫是什么。已知现在有N个有序的随机变量，根据贝叶斯他们的联合分布可以写成条件连乘：

再继续拆：$$p(x_1,x_2,x_3...x_N) = \prod_{i=1}^{n}p(x_i|x_{i-1},...,x_1)$$ 马尔科夫链就是指，序列中的任何一个随机变量在给定它的前一个变量的分布于更早的变量无关。也就是说当前的变量只与前一个变量有关，与更早的变量是没有关系的。用公式表示：

代进去：$$Q(\lambda,\hat\lambda) = \sum_Iln\pi_{i_1}P(O,I|\hat\lambda) + \sum_I(\sum_{t=1}^{T-1}lna_{i_ti_{t+1}})P(O,I|\hat\lambda)+\sum_I(\sum_{t=1}^{T}lnb_{i_to_t})P(O,I|\hat\lambda)$$ 为什么要分成三个部分呢？因为求导一下就没一边了，开心。 求： 有一个条件：。既然是有条件了，拉格朗日乘子法就用上了。化简一下上式：$$ \sum_Iln\pi_{i_1}P(O,I|\hat\lambda) = \sum_Iln\pi_{i_1}P(O,i_1=i|\hat\lambda) $$拉格朗日乘子法：

求导： 上面那几条式子全部加起来： 所以求A： 这里的条件就是，还是延用上面的方法，拉格朗日求导，其实就是改一下上面的公式就好了。所以有：

求B： 这里的条件就是，于是：

这样就是整个更新算法了。 ######有监督学习 如果已经有了一堆标注好的数据，那就没有必要再去玩baum-welch算法了。这个算法就很简单了，多的不说，直接上图：

#####问题3：如果有了参数和观测序列，求的概率最大的隐含序列是什么。 上面骰子的例子已经举过了，再讲一个下雨天的例子：

转移概率如上

发射概率入上 初始概率已知模型参数和观测三天行为(walk,shop,clean)，求天气。 ①首先是初始化：
②第二天： ③最后一天： 接着就是回溯了，查看最后一天哪个最大，倒着找，最后就是2,1,1了，也就是(sun,rain,rain)。 ###⑤Hidden Markov Model代码实现 ####1.工具类的实现

class Tool(object):infinite = float(-2**31)@staticmethoddef log_normalize(a):s = 0for x in a:s += xif s == 0:print('Normalize error,value equal zero')returns = np.log(s)for i in range(len(a)):if a[i] == 0:a[i] = Tool.infiniteelse:a[i] = np.log(a[i]) - sreturn a@staticmethoddef log_sum(a):if not a:return Tool.infinitem = max(a)s = 0for t in a:s += np.exp(t - m)return m + np.log(s)@staticmethoddef saveParameter(pi,A, B, catalog):np.savetxt(catalog + '/' + 'pi.txt', pi)np.savetxt(catalog + '/' + 'A.txt', A)np.savetxt(catalog + '/' + 'B.txt', B)
复制代码

准备了几个静态方法，一个就是log标准化，也就是把所有的数组都归一化操作，变成一个概率，log算加和，保存矩阵。所有的操作都是使用log表示，如果单单是只用数值表示的话，因为涉及到很多的概率相乘，很多时候就变很小，根本检测不出。而用log之后下限会大很多。 ####2.有监督算法的实现 其实就是根据上面的公式统计即可。

def supervised(filename):'''The number of types of lastest variable is four,0B(begin)|1M(meddle)|2E(end)|3S(sigle):param filename:learning fron this file:return: pi A B matrix'''pi = [0]*4a = [[0] * 4 for x in range(4)]b = [[0] * 65535 for x in range(4)]f = open(filename,encoding='UTF-8')data = f.read()[3:]f.close()tokens = data.split('  ')#start traininglast_q = 2old_process = 0allToken = len(tokens)print('schedule : ')for k, token in enumerate(tokens):process = float(k) /float(allToken)if process > old_process + 0.1:print('%.3f%%' % (process * 100))old_process = processtoken = token.strip()n = len(token)#empty we will choose anotherif n <= 0:continue#if just only oneif n == 1:pi[3] += 1a[last_q][3] += 1b[3][ord(token[0])] += 1last_q = 3continue#if notpi[0] += 1pi[2] += 1pi[1] += (n-2)#transfer matrixa[last_q][0] += 1last_q = 2if n == 2:a[0][2] += 1else:a[0][1] += 1a[1][1] += (n-3)a[1][2] += 1#launch matrixb[0][ord(token[0])] += 1b[2][ord(token[n-1])] += 1for i in range(1, n-1):b[1][ord(token[i])] += 1pi = Tool.log_normalize(pi)for i in range(4):a[i] = Tool.log_normalize(a[i])b[i] = Tool.log_normalize(b[i])return pi, a, b
复制代码

按照公式来，一个一个单词判断，如果是单个的那自然就是single，所以当的时候直接就是。初始状态是因为一开始肯定是结束之后才接着开始的，所以自然就是2，end。之后都是比较常规了。 ####3.viterbi算法的实现

def viterbi(pi, A, B, o):'''viterbi algorithm:param pi:initial matrix:param A:transfer matrox:param B:launch matrix:param o:observation sequence:return:I'''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):#first iterationdelta[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[0][j]if delta[t][i] < vj:delta[t][i] = vjpre[t][i] = jdelta[t][i] += B[i][ord(o[t])]decode = [-1 for t in range(T)]q = 0for i in range(1, 4):if delta[T-1][i] > delta[T-1][q]:q = idecode[T-1] = qfor t in range(T-2, -1, -1):q = pre[t+1][q]decode[t] = qreturn decode
复制代码

根据上面的例子来就好，先找到转移概率最大的，迭代到最后找到概率最大的序列之后倒着回来找即可。最后就得到一串编码，然后使用这段编码来进行划分。

def segment(sentence, decode):N = len(sentence)i = 0while i < N:if decode[i] == 0 or decode[i] == 1:j = i+1while j < N:if decode[j] == 2:breakj += 1print(sentence[i:j+1],"|",end=' ')i = j+1elif decode[i] == 3 or decode[i] == 2:  # singleprint (sentence[i:i + 1], "|", end=' ')i += 1else:print ('Error:', i, decode[i] , end=' ')i += 1
复制代码

这个就是根据编码划分句子的函数了。首先是通过有监督的学习得到参数，然后用viterbi算法得到编码序列，再用segment函数做划分即可。

if __name__ == '__main__':pi = np.loadtxt('supervisedParam/pi.txt')A = np.loadtxt('supervisedParam/A.txt')B = np.loadtxt('supervisedParam/B.txt')f = open("../Data/novel.txt" , encoding='UTF-8')data = f.read()[3:]f.close()decode = viterbi(pi, A, B, data)segment(data, decode)
复制代码

执行代码。

效果当然可以了，毕竟是有监督。无监督的就没有这么秀了。 ####4.baum-welch算法的实现 参考上面三个公式，除了B有点难更新之外其他的都很简单。

def cal_alpha(pi, A, B, o, alpha):print('start calculating 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 ifor 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] = Tool.log_sum(temp)alpha[t][i] += B[i][ord(o[t])]print('The calculation of alpha have been finished......')def cal_beta(pi, A, B, o, beta):print('start calculating beta......')T = len(o)for i in range(4):beta[T-1][i] = 1temp = [0 for i in range(4)]del ifor t in range(T-2, -1, -1):for i in range(4):beta[t][i] = 0for j in range(4):temp[j] = A[i][j] + B[j][ord(o[t + 1])] + beta[t + 1][j]beta[t][i] += Tool.log_sum(temp)print('The calculation of beta have been finished......')def cal_gamma(alpha, beta, gamma):print('start calculating gamma......')for t in range(len(alpha)):for i in range(4):gamma[t][i] = alpha[t][i] + beta[t][i]s = Tool.log_sum(gamma[t])for i in range(4):gamma[t][i] -= sprint('The calculation of gamma have been finished......')def cal_kesi(alpha, beta, A, B, o, ksi):print('start calculating ksi......')T = len(o)temp = [0 for i in range(16)]for t in range(T - 1):k = 0for 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 += 1s = Tool.log_sum(temp)for i in range(4):for j in range(4):ksi[t][i][j] -= sprint('The calculation of kesi have been finished......')复制代码

的更新按照公式即可。

def update(pi, A, B, alpha, beta, gamma, ksi, o):print('start updating......')T = len(o)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] = Tool.log_sum(s1) - Tool.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 = 0for t in range(T):if ord(o[t]) == k:s1[valid] = gamma[t][i]valid += 1s2[t] = gamma[t][i]if valid == 0:B[i][k] = -Tool.log_sum(s2)else:B[i][k] = Tool.log_sum(s1[:valid]) - Tool.log_sum(s2)print('baum-welch algorithm have been finished......')
复制代码

最后再封装一把：

def baum_welch(pi, A, B, filename):f = open(filename , encoding='UTF-8')sentences = f.read()[3:]f.close()T = len(sentences)   # 观测序列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(1000):print('Time : ', time)sentence = sentencescal_alpha(pi, A, B, sentence, alpha)cal_beta(pi, A, B, sentence, beta)cal_gamma(alpha, beta, gamma)cal_kesi(alpha, beta, A, B, sentence, ksi)update(pi, A, B, alpha, beta, gamma, ksi, sentence)Tool.saveParameter(pi, A, B, 'unsupervisedParam')print('Save matrix successfully!')
复制代码

初始化矩阵：

def inite():pi = [random.random() for x in range(4)]Tool.log_normalize(pi)A = [[random.random() for y in range(4)] for x in range(4)]A[0][0] = A[0][3] = A[1][0] = A[1][3] = A[2][1] = A[2][2] = A[3][1] = A[3][2] = 0B = [[random.random() for y in range(65536)] for x in range(4)]for i in range(4):A[i] = Tool.log_normalize(A[i])B[i] = Tool.log_normalize(B[i])return pi , A , B
复制代码

最后跑一遍就OK了。 无监督算法使用的是EM框架，肯定会存在局部最大的问题和初始值敏感，跑了56次，用的谷歌GPU跑，代码没有经过优化，跑的贼慢。

最后的效果也是惨不忍睹。 ###总结 概率图模型大多都是围绕着三个问题展开的，求观测序列的概率最大，求隐含序列的概率最大，求参数。MEMM,RCF大多都会围绕这几个问题展开。求观测序列的概率，暴力求解是为了理解模型，前向后向算法才是真正有用的；概率最大的隐含序列viterbi算法，动态规划的思想最后回溯回来查找；求参数就是套用EM框架求解。 HMM是属于生成模型，首先求出，转移概率和表现概率直接建模，也就是对样本密度建模，然后才进行推理预测。事实上有时候很多NLP问题是和前后相关，而不是只是和前面的相关，HMM这里明显是只和前面的隐变量有关，所以还是存在局限性的。对于优化和模型优缺点等写了MEMM和RCF再一起总结。

###最后附上GitHub代码： github.com/GreenArrow2…