本文接上一篇，继续探讨隐马尔科夫模型三个基本问题中最后一个问题，也是最复杂的一个问题。在此之前还是欢迎主角HMM登场！

3. 隐马尔科夫模型的三个基本问题

Problem 3. 学习问题 (Learning Problem)

已知观测序列O=(o₀, o₁, o₂, o₃, …, o_T-1)，需通过调整HMM的参数λ= (A, B, π)使是概率P(O|λ)最大。用气温与年轮的例子来讲就是，已知感兴趣的这十年年轮大小分别为S-M-S-L-M-L-M-S-M-M，求最可能出现这种情况的初始概率、状态转移矩阵和观测矩阵。

对于这个问题枚举法显然是行不通了。因为不可控的元素太多了，不可能像第一类问题(Evaluation Problem)或第二类问题(Decoding Problem)那样得到一个唯一解，不同的初始值和不同的初始矩阵会得到不同的结果，学习类问题通常是这样。同时对于这个问题，通常需要观察链比较长（即需要输入更多信息）。

3.1 算法推导

目前有一种专门的算法来处理这个问题，即Baum-Welch算法，它是EM算法在解决HMM学习问题时的一个特例。

PS. Baum-Welch算法的两位作者分别是Leonard Baum和Lloyd Welsh，其中Leonard Baum是一位传奇人物，他是著名金融公司——文艺复兴科技公司的核心创始人，咱们将要讲到的Baum-Welch算法帮助该公司在对冲基金中连续27年回报率打败巴菲特。有兴趣的朋友可以搜索这段故事来看看，相信看完对学习Baum-Welch算法会充满动力！

另外，有关EM算法(Expectation Maximization Algorithm)欢迎参考我的另一篇博文：EM算法简介

大家在看完EM算法后，对EM算法的套路已经有了一个大致的了解，那么具体到Baum-Welch算法，推导如下：

在推导之前，我们回顾一下在第一类问题(Evaluation Problem)中有介绍到前向算法(Forward Algorithm)，其中前向概率

同理，还有一个后向算法(Backward Algorithm)，其中后向概率

(1) E-step

我们先定义两个变量ξ_t(i, j)和γ_t(i)：

                         (1)

-------------------------------------------------------------------------------------------------------

                          (2)

这两个变量看起来好像很复杂，但是它们想要表达的概念并不复杂。ξ_t(i, j)表示在t时刻状态为i，在t+1时刻状态转为状态j的概率。再具体一些来讲，例如ξ₂(H, L)表示第二年是高温(H)，而第三年为低温(L)的概率。

而γ_t(i)表示t时刻为状态i的概率，例如γ₂(H)表示第二年是高温(H)，而第三年为低温(L)或高温(H)的概率之和，也即仅考虑第二年为高温(H)的概率。

如果还是不太明白，可以看看下图

于是有

(2) M-step

有了上述的两个式子后，我们可以轻松地写出当前状态转移矩阵中A的元素a_ij和观测矩阵B中的元素b_j(k)：

如何理解这个式子？

打个比方，如果用a_HL表示当前参数条件下整个链中由高温转低温的概率，那么式子中的分子表示所有（t从第一年至倒数第二年，因为是有转移，这里是写为到倒数第二年）高温(H)转到低温(L)时刻的概率之和；分母t时刻为高温(H)的概率(无论t+1时刻的状态是什么)。两者相除即由i状态转移到j状态的概率。

式子中，分子是所有满足观测为k的这些时刻，状态为j的概率之和，例如所有观测到大年轮(L)的那些年份状态为高温(H)的期望值；分母为所有t时刻(第一年至倒数第二年)状态为j的概率之和，例如在所有年份中状态为高温(H)的期望值。两者相除即观测为k时状态为j的概率，即由状态j转为观测k的概率。

另外，还需指定初始状态分布。如果有经验值，可按照经验；如果没有可随机取值：

这样就获取到当前状态(t=1)的转移矩阵A₁和当前观测矩阵B₁，以及初始状态π，即当前λ₁=(A₁, B₁, π)。

(3) 用λ1代替λ0，不断地迭代M-step和E-step，直到收敛

在计算的过程中通常会用log函数来简化计算，因此也即：

或者直接指定迭代多少次后就退出迭代过程，终止计算

3.2 实例分析

提到HMM当然要提到晴雨预测例子，假设观测了三天，观测结果为O={o₁=soggy, o₂=dry, o₃=dryish}。再次强调，实际应用中，必须要有足够多的信息才能得到较好的结果，也就是说观测链通常会比较长。这里为了简化计算，观测链长度仅为3。

先初始化一组参数λ₁=(A₁, B₁, π)，其中的值是随机的(当然如果有经验的人可以找到一些初始化里的细微套路，使得算出的局部最优解为全局最优解的可能性更大)，本例中设置的值如下：

=> 状态转移矩阵A₁

=> 初始状态分布π

=> 观测矩阵B₁

将状态转移矩阵A₁和观测矩阵B₁的所有元素标记在一张图中，如下：

(1) 基于以上数据，利用向前向后算法计算出所有的前向概率和后向概率

再利用前向算法计算t=2和t=3时的前向概率，如下：

在计算完前向概率后，我们可以算出在当前参数条件下，观测链O出现的概率：

P(O|λ₁)=α₃(1) +α₃(2)+α₃(3)=0.013

同理我们也计算出所有的后向概率。注意：后向概率是从后往前算的(即先计算t=3，再算t=2，最后算t=1)，t=3时：

再利用后向算法求出所有的后向概率：

(2) 计算中间变量ξ_t(i, j)和γ_t(i)

(3) 计算当前的参数λ₂=(A₂, B₂, π)，并用当前的参数代替λ₁，完成一次迭代

(4) 不断迭代直到收敛

这里将t=1和t=2时计算结果列出来，有兴趣的读者可以试着计算一下

我们可以看到表中P(O|Δ)变化较大，表明未收敛。

注：这里写的P(O|Δ)就是上文中的P(O|λ)，只是写法不同，内容一致

当迭代至第9次后结果收敛，计算终止。如下表

因此，最终结果如下表：

参考材料：

[1] Dmitry Shemetov. Bayesian structural inference and the Baum-Welch algorithm.

[2] Leonard E. Baum, Ted Petrie, George Soules and Norman Weiss. A Maximization Technique Occurring in the Statistical Analysis of Probabilistic Functions of Markov Chains.

[3] Loc Nguyen. Tutorial on Hidden Markov Mode