【大道至简】机器学习算法之隐马尔科夫模型(Hidden Markov Model, HMM)详解（3）---学习问题：Baum-Welch算法推导及Python代码实现

五点钟科技

已于 2024-01-23 16:01:01 修改

阅读量1.9k

点赞数 24

分类专栏：大道至简系列 # 机器学习算法系列文章标签：算法人工智能自然语言处理隐马尔科夫模型 Baum Welch

于 2023-02-07 11:09:47 首次发布

本文链接：https://blog.csdn.net/qq_36583400/article/details/128575653

版权

大道至简系列同时被 2 个专栏收录

11 篇文章 59 订阅

订阅专栏

机器学习算法系列

9 篇文章 2 订阅

订阅专栏

☕️ 本文系列文章汇总：

（1）HMM开篇：基本概念和几个要素

（2）HMM计算问题：前后向算法

（3）HMM学习问题：Baum-Welch算法
（4） HMM预测问题：维特比算法

☕️ 本文来自专栏：大道至简之机器学习系列专栏

❤️欢迎各位小伙伴们关注我的大道至简之机器学习系列专栏，一起学习各大机器学习算法

❤️还有更多精彩文章（NLP、热词挖掘、经验分享、技术实战等），持续更新中……欢迎关注我，个人主页：https://blog.csdn.net/qq_36583400，记得点赞+收藏哦！

📢个人GitHub地址：https://github.com/fujingnan

四、Baum-Welch算法的参数估计（重要核心）

写在前面

阳康后，赶紧捡起续篇，本篇是继上篇概率计算问题后，介绍的又一个重要的HMM经典问题---学习问题。学习问题主要算法就是Baum-Welch算法，在学习推导这个算法的过程中，简直要了我的命。特别是对于状态转移概率矩阵的参数估计，先后咨询了七月在线的学员、前同事（一个牛逼的数学专业的留美硕士），各大QQ群和微信群，百度搜索结果几乎遍历了个遍，甚至拉上我对象一块思考（她数学脑子比我好使）……总之，真的体会到了想要学好一个东西，就要绞尽脑汁千方百计不择手段……最终一定能成功！

在我搞懂算法推导后，我体会最深的就是学习一个算法，课本并非是唯一的选择，原因很简单，好几页的公式推导浓缩成几行，你看不懂很正常。网上搜出来的博客也并非是好的选择，那种和课本如出一辙的套路，能把你讲明白才怪。作者本身懂不懂我不清楚，反正肯定没法给读者讲懂。

所以总结一下我是怎么学习的：课本+所有可调动的资源（包括自己的脑子）。这篇文章，我要尽全力写成网上最好的Baum-welch详解博客，我会把我掌握的知识一步一步推给大家，让大家应懂尽懂。

一、何为学习问题

在概率计算问题中，我们都是已知了模型参数λ和观测结果，然后计算在已知模型λ的条件下，出现观测O的概率，因为该有的已知条件都有了，所以问题就简单了。但是，如果我们并不知道模型参数λ是啥，只观测到了结果是啥，这个时候如果我们想得到导致该结果出现的模型参数λ，该咋办？看上去有点无中生有了。那么对于这种无中生有的问题，只能靠已知的观测结果以及已有的数据来学习得到模型参数了，这就是学习问题需要做的事，这一点和深度学习中的模型训练思想特别类似。只不过这里我们用到的是无监督的学习方法，不需要人工进行样本的标注，隐马尔科夫模型中的学习方法用的就是Baum-welch算法，所以，本篇我们就来详细介绍一下Baum-welch算法的原理。

二、Baum-welch算法的概念

因为我们现在只有观测数据，而要仅通过观测数据来推导隐马尔科夫模型参数λ，怎么做呢？通过前几篇文章的介绍我们知道，每个观测值的背后都隐藏着一个隐状态，也就是说，一组观测值 $O_{1},O_{2},...,O_{T}$ 对应着一组状态序列 $i_{1},i_{2},...,i_{T}$ ，换句话说，学习问题中的隐马尔科夫模型其实就是一个含有隐变量I的隐马尔科夫模型。求解含有隐变量的模型，其实就是求解含有隐变量的模型参数估计问题。学过EM算法的同学一定对这句话特别熟，没错，Baum-welch算法本质上就是EM算法，EM算法就是要解决这种含有隐变量的参数估计问题。所以要学懂Baum-welch算法，还是需要先学好EM算法的基础，大家可以参考我写过关于EM算法原理介绍的博客。当然，Baum-welch算法提出的时候，还没有提出EM算法，但是很多教科书基本上都是把EM的学习放在前，Baum-welch算法放在后，因为EM更加通用更加基础。

总的来说：Baum-welch有如下描述

假设给定训练数据只包含S个长度为T的观测序列{ $O_{1},O_{2},...,O_{S}$ } 而没有对应的状态序列，比如上一篇的例子中，我们只观测到{红球，白球，红球}，目标是学习隐马尔可夫模型λ=(π，A，B)的参数π、A、B。如果将观测序列数据看做观测数据O，状态序列数据看做不可观测的隐数据I，那么隐马尔可夫模型事实上是一个含有隐变量的概率模型

$P(O;\lambda)=\sum_{I}P(O,I;\lambda)=\sum_{I}P(O|I;\lambda)P(I;\lambda) \qquad (10.32)$

为了等会推导方便，这里公式编号我直接和《统计学习方法》中的编号对应上。顺便插一嘴，公式里我用的分号，分号后的参数代表该参数参与公式计算，而|后的参数表示在该参数已知的条件下。书中不用分号，用的是|，这样容易和条件概率中的|混淆，所以我这边用分号代替。关于该公式为啥这么变换，我在我往期博客中说了很多遍了，就是对一个边缘概率转化成联合概率求和的形式。

三、重温EM算法

既然Baum-welch算法的本质就是EM算法，那我想还是稍微复习一下EM算法推导吧，大家也可以直接跳到这里详细阅读学习以下EM算法，为了高效复习，这一节EM算法的讲解我会按照《统计学习方法》中EM算法章节的内容做推导，而且比书上更详细易懂。磨刀不误砍柴工，走起~

EM算法的用武之地就是要通过公式10.32来得到概率模型的完全数据的对数似然函数（所谓完全数据，说白了就是可观测和不可观测变量都要体现在概率模型中），并通过极大化该似然函数来估计公式中的各项参数，从而得到模型。所以这里通过公式10.32得到的完全数据的对数似然函数就是

$\begin{aligned} L(\lambda)&=logP(O;\lambda) \\ &=log\sum_{I}P(O,I;\lambda) \\ &=log\sum_{I}P(O|I;\lambda)P(I;\lambda) \qquad (9.12) \end{aligned}$

由于EM算法的求解也是采用迭代的思想，因为我们要极大化似然函数，所以迭代的趋势就是每一次迭代得到的似然函数值都比上一次大，也就是说 $L(\lambda)>L(\overline{\lambda})$ 一定要成立。 $\overline{\lambda}$ 表示上一次迭代得到的参数值，是个已知量。现在我们来考虑一下两者的差：

$\begin{aligned} L(\lambda)-L(\overline{\lambda})&= \log\sum_{I}P(O|I;\lambda)P(I;\lambda)-\log P(O;\overline{\lambda}) \qquad (1) \\ &=\log \sum_{I}P(I|O;\overline{\lambda})\frac{P(O|I;\lambda)P(I;\lambda)}{P(I|O;\overline{\lambda})}-\log P(O;\overline{\lambda}) \qquad (2) \\ &\geqslant \sum_{I}\left [ P(I|O;\overline{\lambda})\log\frac{P(O|I;\lambda)P(I;\lambda)}{P(I|O;\overline{\lambda})} \right]-\sum_{I}P(I|O;\lambda)\log P(O;\overline{\lambda}) \qquad (3) \\ &=\sum_{I}P(I|O;\overline{\lambda})\log \frac{P(O|I;\lambda)P(I;\lambda)}{P(I|O;\overline{\lambda})P(O;\overline{\lambda})} \qquad (4) \end{aligned}$

上述公式是EM算法的核心，理解上面的公式非常之重要，现在我们来看一下每一步是怎么得到的：

（1）-（2）：这其实就是用到了一个技巧，分子分母同时乘以一个数，原式不变，这同时乘以了一个 $P(I|O;\overline{\lambda})$ ，至于为什么要乘以这个东西，是为了后面要用到的Jessen不等式，至于原创是怎么想到的，这个自己去问他~不过我推测原理大概就是因为，我们的EM算法要解决的是已知观测序列未知状态序列的问题，也就是在观测已知的条件下，其产生该观测的对应的状态的概率是多少，而 $P(I|O;\overline{\lambda})$ 恰好就是这个意思，而且随着我们后面的推导，该式子就是我们最终要计算的核心。

（2）-（3）：这一步用到了关键一步-Jessen不等式，我们这里用到的最重要的就是Jessen不等式的结论：

如果f是凸函数，X是随机变量，那么：E[f(X)]>=f(E[X])，通俗的说法是函数的期望大于等于期望的函数。

特别地，如果f是严格凸函数，当且仅当P(X = EX) = 1，即X是常量时，上式取等号，即E[f(X)] = f(EX)。

当f是（严格）凹函数，当且仅当-f是（严格）凸函数，不等号方向反向，也就是E[f(X)]<=f(E[X])

因为log是凹函数，所以对应公式中通俗点讲就是数学期望的log大于等于log的数学期望。

因为Jessen不等式是个纯数学知识，关于该不等式的定理，这里不再赘述了，可以参考通俗理解EM算法，博客里对对该公式在EM算法中的用法写的很清楚了哈。当然也可以自行搜索资料，网上一堆。

接着来看（3）后半部分，比（2）多了 $\sum_{I}P(I|O;\lambda)$ ，为啥多了这个呢？因为我们有（3）到（4）的变换；

（3）-（4）：在（4）中，我们把减号的前后两个部分合在了一起，归功于（3）中后半部分的 $\sum_{I}P(I|O;\lambda)$ ，这样一来，就可以根据log的计算法则将减法变成除法了，这里为啥可以这样？因为P(I|O;λ)本身是个条件概率，求的是在O已知条件下，I取各个值的概率，既然如此，对I取所有可能值的概率之和必然是1，即 $\sum_{I}P(I|O;\lambda)=1$ 。这里可能有人要问了，凭啥后头的 $\sum_{I}P(I|O;\lambda)$ 是1，前面的那个不也是1，这样的话不用在后面补上这个求和，只需要把前面的这个求和去掉就行了啊！注意哈，前面的 $\sum_{I}$ 还和log里的I关联的，所以这里是不可以简单等于1的。那么你又要问了，既然∑是对包含log的那一堆而言的，那 $\sum_{I}P(I|O;\lambda)$ 就不能和后半部分的 $\sum_{I}P(I|O;\lambda)$ 提取公因式啊！这里又要注意，∑的底数是I，I=i1,i2,...,iN，我们其实是对每个i进行计算的，也就是当I=i1时，有

$\left [ P(i_{1}|O;\overline{\lambda})\log\frac{P(O|i_{1};\lambda)P(i_{1};\lambda)}{P(i_{1}|O;\overline{\lambda})} \right]-P(i_{1}|O;\lambda)\log P(O;\overline{\lambda})$

单例情况下，我们可以提取公因式 $P(i_{1}|O;\overline{\lambda})$ 得到

$\begin{aligned} &P(i_{1}|O;\overline{\lambda})\left [ \log\frac{P(O|i_{1};\lambda)P(i_{1};\lambda)}{P(i_{1}|O;\overline{\lambda})} -\log P(O;\overline{\lambda})\right] \\ &=P(i_{1}|O;\overline{\lambda})\log\frac{P(O|i_{1};\lambda)P(i_{1};\lambda)}{P(i_{1}|O;\overline{\lambda})P(O;\overline{\lambda})} \end{aligned}$

同理，I=i2时，同样可以提取 $P(i_{2}|O;\overline{\lambda})$ ，以此类推，然后求和后，总的就可以写成（4）的形式了。我觉得我已经讲的够明白了==

至此，我们把EM算法最核心最难懂的部分讲完了。因为我们要求得是当前时刻的λ，即求解L(λ)，所以我们来看下它是多少。

通过上面的分析，我们得到这么一个式子：

$L(\lambda)\geqslant \sum_{I}P(I|O;\overline{\lambda})\log \frac{P(O|I;\lambda)P(I;\lambda)}{P(I|O;\overline{\lambda})P(O;\overline{\lambda})}+L(\overline{\lambda}) \qquad (9.14)$

对吧，而我们上面说了， $\overline{\lambda}$ 是个已知量，又观察9.14，现在我们的已知量有 $\overline{\lambda}$ 、O，未知量是I，而当O作为已知条件时，条件概率P(I|O; $\overline{\lambda}$ )可求，所以整个 P(I|O; $\overline{\lambda}$ )也是已知的。通过这个分析，我们可以把9.14再简化一下：

$L(\lambda)\geqslant \sum_{I}P(I|O;\overline{\lambda})\log \frac{P(O|I;\lambda)P(I;\lambda)}{P(I|O;\overline{\lambda})P(O;\overline{\lambda})}+L(\overline{\lambda}) \qquad (5)\\ \Rightarrow L(\lambda)\geqslant \sum_{I}P(I|O;\overline{\lambda})\log P(O|I;\lambda)P(I;\lambda)- \sum_{I}P(I|O;\overline{\lambda})\log P(I|O;\overline{\lambda})P(O;\overline{\lambda}) \qquad (6) \\ \Rightarrow L(\lambda)\geqslant \sum_{I}P(I|O;\overline{\lambda})\log P(O|I;\lambda)P(I;\lambda) \qquad (7)$

注意这里有个细节，我们不能简单的在（5）式中，直接将所有的 P(I|O; $\overline{\lambda}$ ) 抹去，虽然它是个已知量，但是这里它和log是相乘的关系，且log里的分子是未知量，且整体和前面的∑相关，就好比∑后面的那个P(I|O; $\overline{\lambda}$ )是log的系数，肯定是不可以简单去掉的。L>=2x+3y，你能简单的把它简化成L>=x+y吗？显然不行，所以从（5）到（6）我们做了一个log展开，展开后就明朗了，（6）中的后半部分全部都是已知量，所以直接消除，就剩下前半部分了，于是就得到了（7）。

高中的时候我们学过，要让 $L(\lambda)\geqslant L(\overline{\lambda})$ 恒成立，只需要让后者的最大值仍然小于前者就行，所以我们这里其实就需要求解

$\overline{\lambda}=\arg max \sum_{I}P(I|O;\overline{\lambda})\log P(O|I;\lambda)P(I;\lambda) \qquad (9.16)$

也就是说，通过求解极大化似然函数来估计促使该似然函数极大化的参数 $\overline{\lambda}$ 。令

$\begin{aligned} Q(\lambda,\overline{\lambda})&=\sum_{I}P(I|O;\overline{\lambda})\log P(O|I;\lambda)P(I;\lambda) \\ &=\sum_{I}\frac{P(I,O;\overline{\lambda})}{P(O;\overline{\lambda})}\log P(I,O;\lambda) \qquad (8) \\ &=\sum_{I}P(I,O;\overline{\lambda})\log P(I,O;\lambda) \qquad (10.33) \end{aligned}$

那么

$\begin{aligned} \lambda&=\arg max \sum_{I}P(I|O;\overline{\lambda})\log P(O|I;\lambda)P(I;\lambda) \\ &=\arg max \sum_{I}P(I,O;\overline{\lambda})\log P(I,O;\lambda) \\ &=\arg max\ Q(\lambda,\overline{\lambda}) \qquad (9.17) \end{aligned}$

上述的Q其实就是所谓的EM算法中的Q函数。另外，之所以有（8）式到10.33式的变换，实际就是条件概率和联合概率的转化公式，即P(A|B)=P(A,B)/P(B)，而P(O; $\overline{\lambda}$ )已知，故消去。至此，我们也详细推导了《统计学习方法》中10.33式是怎么来的，也顺手复习了一下EM算法。

总结一下，我们先由定义出发，写出了含有隐变量的完全数据的对数似然函数式，又从极大化似然函数的思想出发，写出了迭代式，即每次迭代都比上次计算结果大，从而一步步推导出Q函数。式9.17实际上是一次迭代，该次迭代用到了上次迭代计算得到的参数 $\overline{\lambda}$ ，所以我们每一次极大化Q得到的 $\overline{\lambda}$ 都能够反哺计算出更大的Q对应的λ，而首次迭代用到的λ需要我们初始化得到。

四、Baum-Welch算法的参数估计（重要核心）

通过EM算法的重温，我们发现了，Baum-Welch算法要解决的问题就是EM算法要解决的问题，就是对10.33式的参数的极大化似然估计。在隐马尔科夫模型中，我们的λ包括π、A、B，也就是我们需要通过10.33式来估计出这三个参数。从隐马尔科夫模型计算问题中我们得到：

$P(O,I|\lambda)=\pi_{i_{1}}b_{i_{1}}(o_{1})a_{i_{1}i_{2}}b_{i_{2}}(o_{2})...a_{i_{T-1}i_{T}}b_{i_{T}}(o_{T})$

所以，

$\begin{aligned} \arg max\ Q(\lambda,\overline{\lambda})&=\arg max \sum_{I}P(I,O;\overline{\lambda})\log P(I,O;\overline{\lambda}) \\ &=\arg max\sum_{I}P(I,O;\overline{\lambda})\log(\pi_{i_{1}}b_{i_{1}}(o_{1})a_{i_{i}i_{2}}b_{i_{2}}(o_{2})...a_{i_{T-1}i_{T}}b_{i_{T}}(o_{T})) \\ &=\arg max\sum_{I}P(I,O;\overline{\lambda})\left [ \log \pi_{i_{1}}+\sum_{t=1}^{T-1} \log a_{i_{t}i_{t+1}}+\sum_{t=1}^{T}\log b_{i_{t}}(o_{t})\right ] \qquad (10.34) \end{aligned}$

于是，极大化Q函数实际上就是分别极大化 $\sum_{I}P(I,O;\overline{\lambda})\log \pi_{i_{1}}$ 、 $\sum_{I}P(I,O;\overline{\lambda})\sum_{t=1}^{T-1} \log a_{i_{t}i_{t+1}}$ 、 $\sum_{I}P(I,O;\overline{\lambda})\sum_{t=1}^{T}\log b_{i_{t}}(o_{t})$ 这三项。接下来，我们就要分别对这三项求极大化来估计模型参数π、A、B，最痛苦最核心的部分来啦，跟上我的步伐。

1. π的估计

先来看下第一项：

$\sum_{I}P(I,O;\overline{\lambda})\log \pi_{i_{1}}$

其中， $\sum_{i=1}^{N}\pi_{i}=1$ ，我们回忆一下概率计算问题中的那个例子，例子中对于每一个盒子代表每一个状态，初始状态下，π表示初始状态的概率分布，即选取每一个盒子的概率分别为多少，那么这些概率的和即为1。现在有了 $\sum_{i=1}^{N}\pi_{i}=1$ 这个条件的限制，那么求解上式就不再是一个简单的无约束问题了，而是一个带约束的极大化问题，所以我们很自然的想到用拉格朗日乘子法来求。

首先，我们把I展开：

$\sum_{I}P(I,O;\overline{\lambda})\log \pi_{i_{1}}=\sum_{i=1}^{N}P(O,i_{1}=q_{i};\overline{\lambda})\log \pi_{q_{i}}$

这么变换的原理是因为我们的I表示i1, i2, ..., iT，而π只在初始状态有效，所以只取状态i1，又因为i1可取i=1到i=N这N种状态值，所以对I的求和需要变成从i=1到N的求和，所以我们的i1就可以写成i，注意it中的i和单独的那个i是不一样的，前者包含了所有N个状态值，后者是单独一个状态值的表示，本来为了避免歧义，我这里用qi表示每个状态，但为了后续公式敲打方便，下面qi我统一写成i：

$\sum_{I}P(I,O;\overline{\lambda})\log \pi_{i_{1}}=\sum_{i_{1}}P(O,i_{1};\overline{\lambda})\log \pi_{i_{1}}=\sum_{i=1}^{N}P(O,i_{1}=i;\overline{\lambda})\log \pi_{i}$

通过约束条件构建拉格朗日函数：

$L(\pi_{i}, \alpha)=\sum_{i=1}^{N}P(O,i_{1}=i;\overline{\lambda})\log \pi_{i}+\alpha[\sum_{i=1}^{N}\pi_{i}-1]$

常规操作，L对π求偏导，并令导数为0，求得α：

$\frac{\partial L}{\partial \pi_{i}}= \frac{1}{\pi_{i}}P(O,i_{1}=i;\overline{\lambda})+\alpha=0 \\ \\ \Rightarrow P(O,i_{1}=i;\overline{\lambda})+\alpha\pi_{i}=0 \qquad (9)$

注意这里的细节：之所以求偏导后求和符号不见了，这是因为在求和符号中，我们是对i取每个值而言的，如果将求和符号展开成多项式求和，那么我们求偏导就是得在i取不同值的情况下分别求偏导，当i取某个值的时候，多项式中其它项均为常数，求导为0，举个例子，对 $\sum_{i=0}^{1}\sum_{j=0}^{1}b_{ij}\log a_{ij}$ 求a的偏导，我们先把它展开：

$Y=\sum_{i=0}^{1}\sum_{j=0}^{1}b_{ij}\log a_{ij}=b_{00}\log a_{00}+b_{01}\log a_{01} +b_{10}\log a_{10}+b_{11}\log a_{11}$

所以，

$\begin{aligned} \frac{\partial Y}{\partial a_{00}}=b_{00}\frac{1}{a_{00}}\ , \\ \frac{\partial Y}{\partial a_{01}}=b_{01}\frac{1}{a_{01}}\ , \\ \frac{\partial Y}{\partial a_{10}}=b_{10}\frac{1}{a_{10}}\ , \\ \frac{\partial Y}{\partial a_{11}}=b_{11}\frac{1}{a_{11}}\ \ \end{aligned}$

所以 $\frac{\partial Y}{\partial a_{ij}}$ 可表示为：

$\frac{\partial Y}{\partial a_{ij}}=b_{ij}\frac{1}{a_{ij}}$

这就解释了为啥求偏导后∑没有了，对于这个细节，全网找不到解释，显得公式推导解释很不完整。虽然这是个基础，但是我觉得对于非数学科班的人来说，还是有必要解释一下的。

好了，让我们回到（9）式，现在要想办法消掉π，我们知道， $\sum_{i=1}^{N}\pi_{i}=1$ ，所以对（9）式两边求和，得到

$\sum_{i=1}^{N}P(O,i_{1}=i;\overline{\lambda})+\alpha\sum_{i=1}^{N}\pi_{i}=0 \\ \\ \\ \Rightarrow \sum_{I}P(O,I;\overline{\lambda})+\alpha=0 \\ \\ \\ \Rightarrow \alpha=-P(O;\overline{\lambda}) \qquad (10)$

（10）式的求和符号没有了是因为联合概率对I求和后会得到关于O的边缘概率。这种变化如果不懂的小伙伴参考这里。常规操作，我们把（10）式带回到（9），得到

$P(O,i_{1}=i;\overline{\lambda})-P(O;\overline{\lambda})\pi_{i}=0 \\ \\ \Rightarrow \pi_{i}=\frac{P(O,i_{1}=i;\overline{\lambda})}{P(O;\overline{\lambda})} \qquad (11)$

这样我们就得到了π的表达式（11）。

2. A的估计

现在我们来看一下第二项：

$\sum_{I}P(I,O;\overline{\lambda})\sum_{t=1}^{T-1} \log a_{i_{t}i_{t+1}}$

同样，我们有 $\sum_{j=1}^{N}a_{ij}=1$ ，即不管i是多少，所有的j的可能状态值的概率之和为1，在之前那个例子中，即A矩阵的每一行之和，就是对于当前时刻的某个盒子，转移到下一时刻所有盒子的可能性之和为1。同样是带约束问题，为了求导方便，我们需要把a中的t消去：

$\begin{aligned} \sum_{I}P(I,O;\overline{\lambda})\sum_{t=1}^{T-1} \log a_{i_{t}i_{t+1}}&=\sum_{I}\sum_{t=1}^{T-1}P(I,O;\overline{\lambda}) \log a_{i_{t}i_{t+1}} \\ &=\sum_{i_{1}}\sum_{i_{2}}...\sum_{i_{t}}\sum_{t=1}^{T-1}P(O,i_{1},i_{2},...,i_{t};\overline{\lambda})\log a_{i_{t}i_{t+1}} \end{aligned}$

因为我们的t从1到T-1中取任意一个时刻时，都对应一个状态it，每一个it又可以取i=1~N种状态值，所以原本我们除了最后一个∑外，前面的那几个∑都可以写成 $\sum_{i=1}^{N}$ ，但是我们注意到a的下标是两个，表示当前时刻的状态为it时转移到下一时刻的状态为 $i_{t+1}$ ，这两个时刻的状态都有N种可能的状态值，所以这里我们需要再添加一个∑符号表示转移到下一时刻状态的可能取值：

$\begin{aligned} &\sum_{i_{1}}\sum_{i_{2}}...\sum_{i_{t}}\sum_{t=1}^{T-1}P(O,i_{1},i_{2},...,i_{t};\overline{\lambda})\log a_{i_{t}i_{t+1}}\\ &=\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{t=1}^{T-1}P(O,i_{t}=i,i_{t+1}=j;\overline{\lambda})\log a_{ij} \qquad (12)\end{aligned}$

（12）就是我们待求解的形式，这个变换有些复杂，有限篇幅传达解释往往无法面面俱到，但是非常之重要，有看不懂的小伙伴请及时在评论区留言！

同求解π一样，我们构造拉格朗日函数：

$\begin{aligned} L(a_{ij},\beta)&=\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{t=1}^{T-1}P(O,i_{t}=i,i_{t+1}=j;\overline{\lambda})\log a_{ij}+\beta[\sum_{j=1}^{N}a_{ij}-1] \end{aligned}$

L对aij求偏导，并令偏导为0：

$\begin{aligned} \frac{\partial L}{\partial a_{ij}}=\sum_{t=1}^{T-1}\frac{1}{a_{ij}}P(O,i_{t}=i,i_{t+1}=j;\overline{\lambda})+\beta=0 \\ \Rightarrow \sum_{t=1}^{T-1}P(O,i_{t}=i,i_{t+1}=j;\overline{\lambda})+a_{ij}\beta=0 \qquad (13) \end{aligned}$

关于为啥求导后∑没了，这个我上面应该讲的很清楚了。因为我们的 $\sum_{j=1}^{N}a_{ij}=1$ ，所以对（13）式进行求和：

$\sum_{j=1}^{N}\sum_{t=1}^{T-1}P(O,i_{t}=i,i_{t+1}=j;\overline{\lambda})+\sum_{j=1}^{N}a_{ij}\beta=0 \\ \\ \\ \Rightarrow \sum_{t=1}^{T-1}P(O,i_{t}=i;\overline{\lambda})+\beta=0 \\ \\ \\ \Rightarrow \beta=- \sum_{t=1}^{T-1}P(O,i_{t}=i;\overline{\lambda}) \qquad(14)$

a的∑没了好理解，那么P关于j的∑咋也没了呢？因为关于j求和后，相当于P中的j已经取得所有可能性了，这个时候就和j没有关系了，所以可以把j消掉。常规操作，求得的β带回到（13）式中，得到：

$\sum_{t=1}^{T-1}P(O,i_{t}=i,i_{t+1}=j;\overline{\lambda})-a_{ij}\sum_{t=1}^{T-1}P(O,i_{t}=i;\overline{\lambda})=0 \\ \\ \\ \Rightarrow a_{ij}=\frac{\sum_{t=1}^{T-1}P(O,i_{t}=i,i_{t+1}=j;\overline{\lambda})}{\sum_{t=1}^{T-1}P(O,i_{t}=i;\overline{\lambda})} \qquad (15)$

这样我们就得到了aij的表达式（15）。

3. B的估计

我们来看最后一项，这一项比较麻烦，铁友们，我们再坚持一下：

$\sum_{I}P(I,O;\overline{\lambda})\sum_{t=1}^{T}\log b_{i_{t}}(o_{t})$

与上面相同，该项也是一个带约束的问题，其约束项为：

$\sum_{k=1}^{M}b_{i}(k)=1$

和前面有些不同，这里是k有M种可能得观测值，因为这里的b是观测概率矩阵，也就是每个盒子中取得白球和红球的概率，显然白球和红球的概率之和为1。

同样，为了求导方便，我们需要展开一下：

$\begin{aligned} \sum_{I}P(I,O;\overline{\lambda})\sum_{t=1}^{T}\log b_{i_{t}}(o_{t})&=\sum_{i_{1}}\sum_{i_{2}}...\sum_{i_{t}}\sum_{t=1}^{T}P(O,i_{1},i_{2},...,i_{t};\overline{\lambda})\log b_{i_{t}}(o_{t}) \end{aligned}$

这里在展开讲一下：

当t=1时，我们取的是 $\sum_{i_{1}}$ ， $i_{1}$ 取值为 $\sum_{i=1}^{N}$ ，此时 $o_{t}=o_{1}$ ， $o_{1}$ 取值为 $\sum_{k=1}^{M}$ ， $b_{i_{t}}=\sum_{i=1}^{N}b_{i}$ ；

当t=2时，我们取的是 $\sum_{i_{2}}$ ， $i_{2}$ 取值为 $\sum_{i=1}^{N}$ ，此时 $o_{t}=o_{2}$ ， $o_{2}$ 取值为 $\sum_{k=1}^{M}$ ， $b_{i_{t}}=\sum_{i=1}^{N}b_{i}$ ；

...

不管t取多少，都有对应的it取值为1到N种可能的状态值，注意看b的下标也为it，所以 $b_{i_{t}}$ 应写成 $b_{i}$ ，i取值为1到N种可能的状态值。所以上式可以转换成：

总之记住，每次只取一个状态 $i_{t}$ ，每个 $i_t$ 都有1到N种可能的状态值。同样的，因为某一状态下产生所有观测值的概率之和为1，所以有

$\sum_{k=1}^{M}b_{i}(k)=1$

这里b中的ot改写成k是因为对于每一个状态下，ot已经是对应状态下的一个观测变量，该观测变量是由M个可能的观测值组成的，所以没必要写成ot了，用k表示某一个观测值就行。那么好了，这又是一个带约束的问题，同样可以用拉格朗日函数来求解：

$\begin{aligned} L(b_{i}(k),\gamma)=\sum_{i=1}^{N}\sum_{t=1}^{T}P(O,i_{t}=i;\overline{\lambda})\log b_{i}(o_{t})+\gamma[\sum_{k=1}^{M}b_{i}(k)]\end{aligned}$

常规操作，L对b求偏导并令导数为0：

$\begin{aligned} \frac{\partial L}{\partial b_{i}(k)}=\sum_{t=1}^{T}P(O,i_{t}=i;\overline{\lambda})\frac{1}{b_{i}(o_{t})}+\gamma=0 \end{aligned}$

注意这里对b的参数估计有个非常不一样的地方，因为我们的观测值是个已知量，所以对于每个ot（ot=k=1，2，...，M）来讲，k取哪个值是确定的，未知的是那个i，这说明了对b求偏导的时候，只有在ot取得确定值vk的时候才有意义，对于ot取其它值时，上式固然是0。更近一步的解释，在之前的那个盒子取球的例子里，我们看下B矩阵：

上面那段话的意思就是说，对于每一个状态，我们取B的哪一列是已知的，例如，假设当前观测值为白球，那么在盒子2的状态下，我们一定取的是第2行底2列，其它列的可能性为0，在盒子3的状态下必然取第3行第2列，其它列的可能性为0，盒子1的状态亦如是。好了，经过上述分析，我们如何能够表示出 $\frac{\partial L}{\partial b_{i}(k)}$ 是对ot=vk而言的呢，即偏导不为0？很简单，只需要加个指示函数I即可：

$I=\left\{\begin{matrix} 1,o_{t}=v_{k} & \\ 0,o_{t}\neq v_{k} & \end{matrix}\right.$

所以，

$\begin{aligned} \frac{\partial L}{\partial b_{i}(k)}&=\sum_{t=1}^{T}P(O,i_{t}=i;\overline{\lambda})\frac{1}{b_{i}(k)}I+\gamma=0 \\ &\Rightarrow \sum_{t=1}^{T}P(O,i_{t}=i;\overline{\lambda})I+\gamma b_{i}(k)=0\qquad (17)\end{aligned}$

又因为 $\sum_{k=1}^{M}b_{i}(k)=1$ ，所以我们对（17）求和：

$\begin{aligned} &\sum_{k=1}^{M}\sum_{t=1}^{T}P(O,i_{t}=i;\overline{\lambda})I+\gamma \sum_{k=1}^{M}b_{i_{i}}(k)=0 \\ &\Rightarrow \sum_{t=1}^{T}P(O,i_{t}=i;\overline{\lambda})+\gamma =0\\ &\Rightarrow \gamma =-\sum_{t=1}^{T}P(O,i_{t}=i;\overline{\lambda}) \qquad (18) \end{aligned}$

中间那个式子中I不见了，是因为我们对指示函数I进行 $\sum_{k=1}^{M}$ ，而k从1到M中只有一个值才满足I=1，因此就可以把 $\sum_{k=1}^{M}$ 和I都去掉。常规操作，把（18）带回到（17）得：

$\sum_{t=1}^{T}P(O,i_{t}=i;\overline{\lambda})I- b_{i}(k)\sum_{t=1}^{T}P(O,i_{t}=i;\overline{\lambda})=0 \\ \\ \\ \Rightarrow b_{i}(k)=\frac{\sum_{t=1}^{T}P(O,i_{t}=i;\overline{\lambda})I}{\sum_{t=1}^{T}P(O,i_{t}=i;\overline{\lambda})} \qquad(19)$

至此，我们利用Baum-welch算法把模型参数π、A，B都估计出来了。注意上述（11）、（15）、（19）只是针对1次计算而言的，在我们的算法迭代中，是需要不断更新参数的，即我们每次迭代计算新的λ，都需要用到上次得到的λ，而最初始的λ是我们初始化赋值得到的。

那么我们推出了π、A、B的表达式后，具体怎么计算呢？毕竟上面一堆公式，都是用P来表示的，所以这些P应该等于啥呢？接下去看。

五、Baum-welch算法计算流程

通过第四节我们得到参数估计的公式如下：

$\pi_{i}=\frac{P(O,i_{1}=i;\overline{\lambda})}{P(O;\overline{\lambda})} \qquad (11)$

$a_{ij}=\frac{\sum_{t=1}^{T-1}P(O,i_{t}=i,i_{t+1}=j;\overline{\lambda})}{\sum_{t=1}^{T-1}P(O,i_{t}=i;\overline{\lambda})} \qquad (15)$

$b_{i}(k)=\frac{\sum_{t=1}^{T}P(O,i_{t}=i;\overline{\lambda})I}{\sum_{t=1}^{T}P(O,i_{t}=i;\overline{\lambda})} \qquad(19)$

回忆一下前向算法和后向算法文章中介绍的概念：

前向概率：

$\alpha _{t}(i)=P(o_{1},o_{2},...,o_{t},i_{t}=q_{i};\lambda)\qquad(20)$

后向概率：

$\beta_{t}(i)=P(o_{t+1},o_{t+2},...,o_{T}|i_{t}=q_{i};\lambda)\qquad (21)$

以时刻t为节点，t时刻往前的为到时刻t的这部分观序列的前向概率，t时刻往后的为t到T时刻这部分观测序列的后向概率，所以，它们联合起来就是：

$\alpha_{t}(i)\beta_{t}(i)=P(o_{1},o_{2},...,o_{t},i_{t}=i;\overline{\lambda})P(o_{t+1},o_{t+2},...,o_{T}|i_{t}=i;\overline{\lambda})$

这里补充插播一条很重要的知识，即隐马尔科夫模型成立的两个基本假设：

（1）其次马尔科夫假设：假设了隐马尔科夫链在任意时刻t的状态只依赖于前一时刻的状态，与其它时刻的状态及观测无关，也与时刻t无关；

（2）观测独立性假设：即假设了任意时刻的观测只依赖于该时刻的马尔科夫链的状态，与其它观测即状态无关。

根据假设（2），我们知道， $o_{t+1}$ 及其之后的观测和 $i_{t}$ 以及 $o_{t}$ 及其之前的观测无关，它们相互独立，于是，我们把 $o_{t+1},o_{t+2},...,o_{T}$ 看成A，把 $i_{t}=i$ 看成B，根据P(A|B)=P(A, B)/P(B)，有：

$\begin{aligned} P(o_{t+1},o_{t+2},...,o_{T}|i_{t}=i;\overline{\lambda})&=\frac{P(o_{t+1},o_{t+2},...,o_{T},i_{t}=i;\overline{\lambda})}{P(i_{t}=i;\overline{\lambda})}\\ &=\frac{P(o_{t+1},o_{t+2},...,o_{T};\overline{\lambda})P(i_{t}=i;\overline{\lambda})}{P(i_{t}=i;\overline{\lambda})}\\ &=P(o_{t+1},o_{t+2},...,o_{T};\overline{\lambda}) \end{aligned}$

所以，

$\begin{aligned} \alpha_{t}(i)\beta_{t}(i)&=P(o_{1},o_{2},...,o_{t},i_{t}=i;\overline{\lambda})P(o_{t+1},o_{t+2},...,o_{T}|i_{t}=i;\overline{\lambda}) \\ &=P(o_{1},o_{2},...,o_{t},i_{t}=i;\overline{\lambda})P(o_{t+1},o_{t+2},...,o_{T};\overline{\lambda})\\ &=P(o_{1},o_{2},...,o_{T},i_{t}=i;\overline{\lambda})\\ &=P(O,i_{t}=i;\overline{\lambda}) \qquad \qquad \qquad \qquad \qquad(22)\end{aligned}$

根据联合概率和边缘概率的换算公式，进一步得到：

$\sum_{i=1}^{N}P(O,i_{t}=i;\overline{\lambda})=\sum_{i=1}^{N}\alpha_{t}(i)\beta_{t}(i)=P(O;\overline{\lambda}) \qquad \qquad (23)$

将（22）和（23）代入到（11），我们就得到π的另一种形式：

$\pi_{i}=\frac{\alpha_{t}(i)\beta_{t}(i)}{\sum_{i=1}^{N}\alpha_{t}(i)\beta_{t}(i)},\ \ \ \ t=1 \qquad \qquad (24)$

我们再来看一下 $P(O,i_{t}=i,i_{t+1}=j;\overline{\lambda})$ 是个啥，看上去复杂些，其实也不难，这里把i分成了两个部分：t时刻 $i_{t}$ =i与 $o_{1},o_{2},...,o_{t}$ 的联合概率分布且t+1时刻 $i_{t+1}$ =j与 $o_{t+1},o_{t+2},...,o_{T}$ 的联合概率分布。回想一下前后向算法的表达式：

$P(O;\lambda)=\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{t}(i)a_{ij}b_{j}(o_{t+1})\beta_{t+1}(j),\quad t=1,2,...,T-1 \qquad (10.22)$

（统计学习方法公式10.22）

注意看 $P(O,i_{t}=i,i_{t+1}=j;\overline{\lambda})$ 和上式的区别，多了一个条件：在节点 $i_{t}$ =i和节点 $i_{t+1}$ =j处与O的联合概率分布，即在这两个特定节点上，已经不存在说有i=1~N和j=1~N这么一个条件了，因此对于 $P(O,i_{t}=i,i_{t+1}=j;\overline{\lambda})$ ，它等于：

$P(O,i_{t}=i,i_{t+1}=j;\overline{\lambda})=\alpha_{t}(i)a_{ij}b_{j}(o_{t+1})\beta_{t+1}(j) \qquad \qquad (25)$

其实我这么写应该会好看一点：

$P(O,i_{t}=q_{i},i_{t+1}=q_{j};\overline{\lambda})=\alpha_{t}(i)a_{ij}b_{j}(o_{t+1})\beta_{t+1}(j)$

将（25）和（22）代入到（15）得到a的另一种形式：

$a_{ij}=\frac{\sum_{t=1}^{T-1}\alpha_{t}(i)a_{ij}b_{j}(o_{t+1})\beta_{t+1}(j) }{\sum_{t=1}^{T-1}\alpha_{t}(i)\beta_{t}(i)} \qquad (26)$

好了，终于剩下最后一个了，最后一个b的改写就简单了不是，

$b_{i}(k)=\frac{\sum_{t=1}^{T}P(O,i_{t}=i;\overline{\lambda})I}{\sum_{t=1}^{T}P(O,i_{t}=i;\overline{\lambda})}$

根据我们上头的分析，得到：

$b_{i}(k)=\frac{\sum_{t=1,o_{t}=v_{k}}^{T}\alpha_{t}(i)\beta_{t}(i)}{\sum_{t=1}^{T}\alpha_{t}(i)\beta_{t}(i)}\qquad \qquad(27)$

以上过程中的α和β在前后向算法中都介绍过，都是可以通过a、b、π求得的，所以我们在迭代的时候：

已知观测数据： $O=(o_{1},o_{2},...,o_{T})$

1. 对于step=0步，随机初始化a、b、π，得到初始化模型参数λ(0)；

2. 对于step=1,2,...N步，每一步都利用公式（24）、（26）、（27）将上一步的a、b、π的计算结果代入，得到当前步的a、b、π；

3. 通过阈值或者最大迭代轮次设定终止条件，终止后，得到最终模型参数λ；

至此，隐马尔科夫模型的第二个基本问题的核心算法Baum-Welch算法介绍完毕，个人能力有限，文中若有错误不妥之处，请各位读者及时批评指正，不胜感激！代码我会在HMM的所有内容介绍完毕后补上。关注我，不迷路，创作不易，记得点赞哦！

六、代码实现（重要更新）

代码地址：ml_algorithm/baum_welch.py at master · fujingnan/ml_algorithm · GitHub

import random
import numpy as np

random.seed(1)  # 好像不起租用


class AttrDict(dict):
    # 一个小trick，将结果返回成一个字典格式
    def __init__(self, *args, **kwargs):
        super(AttrDict, self).__init__(*args, **kwargs)
        self.__dict__ = self


class Baum_Welch:

    def __init__(self, N, M, V):
        self.A = np.random.dirichlet(np.ones(N), size=N)  # 状态转移概率矩阵
        self.B = np.random.dirichlet(np.ones(M), size=N)  # 观测概率矩阵
        self.pi = np.array(np.random.dirichlet(np.ones(N), size=1))[0]  # 初始状态概率矩阵
        self.V = V # 所有可能的观测
        self.N = N # 所有可能的状态值
        self.M = M # 所有可能的观测值

    def forward(self):
        """
        前向算法，Baum welch算法需要用到
        :param O: 已知的观测序列
        :return: alpha_{i}
        """
        row, col = len(self.O), self.A.shape[0]
        alpha_t_plus_1 = np.zeros((row, col))
        obj_index = self.V.index(self.O[0])
        # 初值α 公式10.15
        alpha_t_plus_1[0][:] = self.pi * self.B[:].T[obj_index]
        for t, o in enumerate(self.O[1:]):
            t += 1
            # 递推 公式10.16
            obj_index = self.V.index(o)
            alpha_ji = alpha_t_plus_1[t - 1][:].T @ self.A
            alpha_t_plus_1[t][:] = alpha_ji * self.B[:].T[obj_index]

        self.alpha = alpha_t_plus_1

    def backward(self):
        """
        后向算法，Baum welch算法需要用到
        :param O: 已知的观测序列
        :return: beta_{i}
        """
        row, col = len(self.O), self.A.shape[0]
        betaT = np.zeros((row + 1, col))
        # 初值β 公式10.19
        betaT[0][:] = [1] * self.A.shape[0]
        for t, o in enumerate(self.O[::-1][1:]):
            t += 1
            # 反向递推 公式10.20
            obj_index = self.V.index(self.O[t - 1])
            beta_t = self.A * self.B[:].T[obj_index] @ betaT[t - 1][:].T
            betaT[t][:] = beta_t
        # betaT[-1][:] = [self.pi[i] * self.B[i][self.V.index(self.O[0])] * betaT[-2][i] for i in range(self.A.shape[0])]
        self.beta = betaT[:-1][::-1]

    def gamma(self, t, i):
        """
        根据课本公式【10.24】计算γ
        :param t: 当前时间点
        :param i: 当前状态节点
        :return: γ值
        """
        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 ksi(self, t, i, j):
        """
        根据公式【10.26】计算 ξ
        :param t: 当前时间点
        :param i: 当前状态节点
        :param j: 同i
        :return:
        """
        obj_index = self.V.index(self.O[t + 1])
        numerator = self.alpha[t][i] * self.A[i][j] * self.B[j][obj_index] * 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][obj_index] * self.beta[t + 1][j]

        return numerator / denominator

    def train(self, O, n):
        """
        根据算法【10.4】训练模型
        :param O: 已知观测序列
        :param n: 最大迭代步长
        :return: 模型参数λ=(π，A，B)
        """
        self.O = O
        self.T = len(O)
        maxIter = 0

        while maxIter < n:
            tempA = np.zeros((self.N, self.N))
            tempB = np.zeros((self.N, self.M))
            tempPi = np.array([0.] * self.N)

            # 根据前向算法和后向算法得到α和β
            self.forward()
            self.backward()

            maxIter += 1
            # a_{ij}，公式【10.39】
            for i in range(self.N):
                for j in range(self.N):
                    numerator = 0.
                    denominator = 0.
                    for t in range(self.T - 1):
                        numerator += self.ksi(t, i, j)
                        denominator += self.gamma(t, i)
                    tempA[i][j] = numerator / denominator

            # b_{i}{j}，公式【10.40】
            for j in range(self.N):
                for k in range(self.M):
                    numerator = 0.
                    denominator = 0.
                    for t in range(self.T):
                        if self.O[t] == self.V[k]:
                            numerator += self.gamma(t, j)
                        denominator += self.gamma(t, j)
                    tempB[j][k] = numerator / denominator

            # π_{i}，公式【10.41】
            for i in range(self.N):
                tempPi[i] = self.gamma(0, i)
            # 更新
            self.A = tempA
            self.B = tempB
            self.pi = tempPi

        return AttrDict(
            pi=self.pi,
            A=self.A,
            B=self.B
        )


if __name__ == '__main__':
    bm = Baum_Welch(N=3, M=2, V=['红', '白'])
    O = ['红', '白', '红']
    res = bm.train(O, 3)
    print('π：', res.pi)
    print('A：', res.A)
    print('B：', res.B)
    print(res.pi.shape)

参考文献：

[1]. 隐马尔可夫模型之Baum-Welch算法详解_Demon的黑与白的博客-CSDN博客_baum-welch

[2]. Baum-Welch算法的公式推导_cxy_19891109的博客-CSDN博客_baum welch算法公式推导

[3]. 《统计学习方法》李航

五点钟科技

关注

24
点赞
踩
26

收藏

觉得还不错? 一键收藏
打赏
3
评论
【大道至简】机器学习算法之隐马尔科夫模型(Hidden Markov Model, HMM)详解（3）---学习问题：Baum-Welch算法推导及Python代码实现

求解含有隐变量的模型，其实就是求解含有隐变量的模型参数估计问题。学过EM算法的同学一定对这句话特别熟，没错，Baum-welch算法本质上就是EM算法，EM算法就是要解决这种含有隐变量的参数估计问题。所以要学懂Baum-welch算法，还是需要先学好EM算法的基础，大家可以参考我写过关于EM算法原理介绍的博客。当然，Baum-welch算法提出的时候，还没有提出EM算法，但是很多教科书基本上都是把EM的学习放在前，Baum-welch算法放在后，因为EM更加通用更加基础。隐马尔科夫模型Python实现
复制链接

扫一扫