自然语言处理从小白到大白系列（3）聊聊序列模型中的隐马模型

如标题一样，隐马尔可夫模型谜一样的推导和应用，一直是机器学习入门朋友们的一个拦路虎。就是那种，提起来大致知道：噢！隐马模型啊，就是那个转移来转移去的一个模型，要解决三个基本问题，哪三个来着？对了，还有那啥维特比算法，前向后向算法，好了我懂！可是聪明的你，真的弄懂了吗？本文就来捋一捋，这是何方神圣，究竟有什么用处？？

通过阅读本文，你将会掌握以下主要知识点：

• 隐马模型的初步理解
• 隐马模型的三个重要问题
• 什么是维特比算法
• Baum welch 算法
• EM算法
• HMM的实战演练

1. 隐马模型的初步理解

1.1 马尔科夫模型

首先，我们要知道，隐马尔可夫模型是一种概率图模型，可以看成是概率图模型里面的动态贝叶斯网络，为什么说是动态的贝叶斯网络呢，就是由于它是随着时间节点，网络结构是不断变化的。

先说马尔科夫模型，马尔科夫模型的意思是说：一个系统中，某个时刻t的状态，取决于它在时间1,2,3，…，t-1时刻的状态，就是说是从前面的状态转移过来的。好了，既然有这个定义，那么我们就可以用数学表达式写一下，时刻t状态为Si的概率是：
​​​​
其中，$q_t$是t时刻的状态变量，$S_j$表示第j个实际状态（tips：一定要注意本文定义的符号的含义，不然到后面就整懵了）

好的，这样一个式子，想一想，第100个时刻的状态，就和前面99个时刻的状态有关，那1000呢，10000呢，是不是已经复杂到不想计算了，因此，为了简化这样一个问题，我们引入一阶马尔可夫链，官方定义：

如果在特定情况下，系统在时间 t 的状态只与其在时间 t-1 的状态相关，则该系统构成一个离散的一阶马尔柯夫链：

看到没，已经将除前一个时刻以外的所有时刻的状态都忽视了，是不是看起来清爽了许多，其实这样一个假设在处理这类问题的时候是非常常用的，比如说语言模型的n-gram，不也是用的这样一个马尔可夫性的假设，来简化n元语法吗。

但是还有一个问题，就是在t时刻状态由$S_i$转移到$S_j$时，和在t’时刻，状态由$S_i$转移到$S_j$，概率是一样的吗？？如果不一样，这还是一个依赖于时间的模型，会更加复杂，所以我们还要做一个假设，就是状态转移概率独立于时间。官方定义：

如果只考虑公式(5.15)独立于时间 t 的随机过程，即所谓的不动性假设，状态与时间无关，那么：

以上这个随机过程，就被称为是马尔柯夫模型 (Markov Model)，当然不用说，要满足转移概率大于等于0，从某时刻转移出去的概率之和要等于1,这些应该都很好理解。所以：

（马尔柯夫模型又可视为随机有限状态自动机，该有限状态自动机的每一个状态转换过程都有一个相应的概率，该概率表示自动机采用这一状态转换的可能性。）

如果还是不好理解的话，画一个图来表示马尔科夫模型的转移：

现在问一下，如果是一阶的马尔科夫模型，状态序列$S_1,S_2,...,S_T$的概率应该怎么写呢？很简单对吧：

这里的$p_i$是初始状态的概率。

所有刚才那个转移图，t,i,p序列出现的概率是：

1.2 隐马尔可夫模型：

刚才是从随机过程到马尔科夫模型，好像都很容易理解的样子，现在我们讨论一下隐马尔可夫模型，隐，隐在什么地方？

要弄清楚这个，就先得分清楚：状态序列和观测序列的概念。状态序列，可以理解为刚才的马尔科夫模型的状态；观测序列是指：在每个时刻，当前时刻的状态值，以一定概率产生这样一个观测值，那是不是在每个时刻都有一个观测值，所有观测值连起来，就是一个观测序列。那定义这样一个东西有什么用呢？是由于，在实际的情况里，状态值也许是不可见的，我们能看到的，只是观测值，要挖掘这样一个隐藏的关系，则就用到了我们说的隐马尔可夫模型。

好了，你说还是不明白，太抽象，那么我们举个业界经典的例子。

• 某赌场在投骰子，根据点数决定胜负。在多次投掷骰子的时候采取了如下手段进行作弊：准备了两个骰子A和B，其中A为正常骰子，B为灌铅骰子，由于怕被发现，所有连续投掷的时候偶尔使用一下B，A和B之间转换的概率如下：

所以，我们要解决的问题，就是

• 1.通过抛出的观测序列和给定模型，评估产生这个观测序列的概率是多大？
• 2.通过抛出的点数序列，解码出最可能是什么时候作的弊？
• 3.给定观测序列，求出模型参数，即，正常骰子到灌铅骰子的转移概率，以及各自产生点数的概率？

好了，假设你已经懂了隐马尔可夫模型的基本定义，那在隐马尔可夫模型里面，我们根据上面提出，抽象出来三类基本问题：

2. 隐马模型的三个重要问题

• 1.评估问题

• 2.解码问题

• 3.学习问题

2.1 评估问题

各个击破，先是评估问题，评估问题很容易理解，就是给出模型参数和一段观测序列，然后评估产生这个序列的概率有多大。规范描述：

对于这样一个问题，首先要搞清楚的是，对于某个时刻的任意状态，都有可能产生O_t。那么我们把概率改写成如下：

上式应该好理解，添加了一个中间变量Q（状态），这是由于O必须由Q产生，所以把过程分解一下，得到这样一个式子。
然后分解来看，

这个没有问题，但是问题是，5.22和5.23，只计算了从头到尾的某一个状态序列。你想想，如果状态有N个，时间有T个，产生可能的序列是$N^T$个！！因为每个时刻都可能产生N个状态。这还算个毛，没法进行下去了。对于这样一个问题，我们隐隐约约感觉，问题是可以简化的，就是说，计算是有冗余的，而且我们还隐隐约约记得，老师曾经讲过一个叫动态规划的方法。

好的，不哔哔，拿出动态规划大法，我们知道，动态规划的问题关键是找到状态转移方程，就是说问题可以化成子问题解决。
这里就要想想了，如何定义一个变量，从这个变量推导下一个变量。这里由于是N个状态，T个时间

那么可以这样定义一个变量

这个alpha_t(i)的含义是，在前面t个时刻，观测分别是O1-Ot的时候，在t时刻状态取到S_i时的概率是多少，注意哦，只考虑t时刻状态S_i。那么我们想，在t时刻状态有N个，当然这个alpha_t(i)和我们最终的P有如下关系：

（就是把t时刻所有的状态加起来）
那么定义这个变量之后，如何列转移方程呢？即如何从alpha_t(i)到alpha_t+1(j)？当然要先成转移概率，乘以产生观测值概率，再把t时刻所有的状态加起来：

算到最后，再带入5.25,大功告成。这样，是不是就可以把所有时刻，所有状态遍历一次，就得到了最终的结果？是不是高效又牛逼？
算法的时间复杂性：
每计算一个 $a_t(i)$ 必须考虑从 t-1 时的所有 N 个状态转移到t时刻的 $S_i$状态的可能性，时间复杂性为 O(N)，对应时刻 t，要计算 N 个向前变量： $a_t(1),a_t(2),\dots ,a_t(N)$，所以时间复杂性为： $O(N)\times N=O(N^2)$。 又因 t = 1, 2, … , T，所以向前算法总的复杂性为： $O(N^{2}T)$。

以上叫前向算法，是由于它从前往后推出来的，下面后向算法是异曲同工的，大致过一下：
定义变量${\beta }_t(i)$：


算法描述：

第一个问题解决了，第二个问题：如何发现“最优”状态序列能够“最好地解释”观察序列。

2.2 解码问题

这个是什么意思呢，就是说产生这样一个观测序列的状态序列是很多个的，哪一个才是最好的呢？

这也是一个运算爆炸的问题，你可能会这么想：找最可能的隐状态序列嘛，那我对每个时刻都找概率最大产生观测值的，仔细想想，这肯定是不行的，因为仅考虑了生成概率，没有考虑转移概率，两个状态可能不能相互转移也说不定呢？另外每步的贪心算法肯定是得不到最优解的，为什么？一句话：因为在这一步选择次优的方案下一步可能得到更好结局！
说了这么多，就是要用到维特比算法，这是一个什么算法如此神奇？

维特比算法

思路是这样的，对于某时刻$t$，状态$q_i$，我如果保存在此状态下产生这个已知观测序列的最大概率，（同时保存这个概率对应的路径），由于有$N$个状态，即在t时刻保存N个最大概率，那么下一个时刻的状态$q_j$产生观测序列的最大概率是多少？当然是从前一个时刻所有状态过来的，而我们已经保存了前一个时刻产生观测序列的最大概率，当然就依次计算，得到t+1时刻的最大概率，一直到最后一步，就得到了那个状态路径！这就是这个算法的核心。
​

好的，第二个问题也收官。现在考虑第三个问题，即学习问题，给定一个观测序列，如何学习到转移概率等参数？

2.3 学习问题

这里得分类讨论：1.状态序列是已知的。2.状态序列是未知的

1. 状态序列已知

如果状态序列已知，那么根据极大似然估计，通俗一点的说，可以把频率数出来，用这个频率代替转移概率。

2. 状态序列未知

如果状态是不知道的，只能用EM算法。

又到了被EM支配的时候，EM算法又称期望最大化，通过不断地迭代，最终能够收敛的一个算法。

EM算法思路是这样的，先随机初始化模型参数，那么由这一组参数，当然就可以得到求出某一状态转移到另一个状态的期望次数，当然也能得到一个状态输出一个观测值的期望次数。现在用这个期望次数替换5.36的实际次数（由于现在状态序列未知，实际次数无法统计了），然后你看5.36，我们是在算什么？我们是在计算模型参数，那么由这样一个参数更新我们初始化的随机参数。再由参数求期望次数，不断循环迭代，参数收敛于最大似然估计值，就能够估计出最终的模型参数。

看到上式的分子，知道为什么叫前向后向算法了吧，$\alpha$是前向，$\beta$是后向。注意，这里的$\alpha$和$\beta$和前面的$\alpha 、\beta$的含义是一致的。用图来表示这样一个过程：

最后，我们发现Baum Welch算法中${\xi }_t(i,j)和{\gamma }_t(i)$这两个变量的设置，为这个EM算法注入了灵魂。为什么这么说呢，我们知道，现在是只有观测序列，而没有状态序列，如果有状态序列的话，这个事情就好办多了，可以通过最简单的统计方法用极大似然估计来做；现在不行，所以只能自己凭空造初始参数，通过初始参数搞到一组统计数据（求期望的过程），比如从状态$q_i$到状态$q_j$的次数和总共从$q_i$转移出去的次数，然后来估计参数$a_{ij}$。那么这组造的数据如何，就决定了估计出来的参数质量怎么样，所以要想办法不断调优这个期望值，那么怎么调优呢？构造的变量要包含$a_{ij}和b_j(k)$，这样估计出来的新参数，才能用来调整这个期望，而这个期望求出来了，又能来调整估计的参数。
那么如何构造${\xi }_t(i,j)和{\gamma }_t(i)$，就很重要了：首先里面要含有$a_{ij}和b_j(k)$，然后我们的观测序列$O$肯定要在里面用上，所以最后我们的${\xi }_t(i,j)$表示一个变量，是在给定参数和观测序列下，在$t$时刻是状态$s_i$，下一个时刻是$s_j$ 的概率，有人说这不是$a_{ij}$的定义吗，注意这不一样！${\xi }_t(i,j)$还要除以一个在$t$时刻为状态$s_i$的概率，也就是除以一个${\gamma }_t(i)$才是$a_{ij}$。所以最后构造出来如式子5.37的形式。

至于EM算法，我想单独写一篇来总结一下，欢迎关注后文。

至此，隐马尔可夫模型大致说完了，这个模型可以用来做句子的词性标注，比如把词性看做是状态，把句子看成观测序列，相当于是一个解码的问题。或者可以用来做分词，但是做分词，是从多个候选分词结果中选一个结果最好的，也是把词性作为状态，对多个结果进行评估，选最佳结果。

3. HMM实战案例

没有实战总感觉是纸上谈兵，缺少仪式感。

• 导入必要的包：

# import python modules -- this cell needs to be run again if you make changes to any of the files
import matplotlib.pyplot as plt
import numpy as np

from helpers import show_model
from pomegranate import State, HiddenMarkovModel, DiscreteDistribution

• 建立模型

# create the HMM model
model = HiddenMarkovModel(name="Example Model")

# emission probability distributions, P(umbrella | weather)
sunny_emissions = DiscreteDistribution({"yes": 0.1, "no": 0.9})
sunny_state = State(sunny_emissions, name="Sunny")

# TODO: create a discrete distribution for the rainy emissions from the probability table
# above & use that distribution to create a state named Rainy
rainy_emissions = DiscreteDistribution({"yes": 0.8, "no": 0.2})
rainy_state = State(rainy_emissions, name="Rainy")

# add the states to the model
model.add_states(sunny_state, rainy_state)

assert rainy_emissions.probability("yes") == 0.8, "The director brings his umbrella with probability 0.8 on rainy days"

• 添加转移概率

# create edges for each possible state transition in the model
# equal probability of a sequence starting on either a rainy or sunny day
model.add_transition(model.start, sunny_state, 0.5)
model.add_transition(model.start, rainy_state, 0.5)

# add sunny day transitions (we already know estimates of these probabilities
# from the problem statement)
model.add_transition(sunny_state, sunny_state, 0.8)  # 80% sunny->sunny
model.add_transition(sunny_state, rainy_state, 0.2)  # 20% sunny->rainy

# TODO: add rainy day transitions using the probabilities specified in the transition table
model.add_transition(rainy_state, sunny_state, 0.4)  # 40% rainy->sunny
model.add_transition(rainy_state, rainy_state, 0.6)  # 60% rainy->rainy

# finally, call the .bake() method to finalize the model
model.bake()

assert model.edge_count() == 6, "There should be two edges from model.start, two from Rainy, and two from Sunny"
assert model.node_count() == 4, "The states should include model.start, model.end, Rainy, and Sunny"
print("Great! You've finished the model.")

• 看一眼模型

show_model(model, figsize=(5, 5), filename="example.png", overwrite=True, show_ends=False)

• 重新调整顺序

column_order = ["Example Model-start", "Sunny", "Rainy", "Example Model-end"]  # Override the Pomegranate default order
column_names = [s.name for s in model.states]
order_index = [column_names.index(c) for c in column_order]

# re-order the rows/columns to match the specified column order
transitions = model.dense_transition_matrix()[:, order_index][order_index, :]
print("The state transition matrix, P(Xt|Xt-1):\n")
print(transitions)
print("\nThe transition probability from Rainy to Sunny is {:.0f}%".format(100 * transitions[2, 1]))

• 计算某观测序列产生的概率

# TODO: input a sequence of 'yes'/'no' values in the list below for testing
observations = ['yes', 'no', 'yes']

assert len(observations) > 0, "You need to choose a sequence of 'yes'/'no' observations to test"

# TODO: use model.forward() to calculate the forward matrix of the observed sequence,
# and then use np.exp() to convert from log-likelihood to likelihood
forward_matrix = np.exp(model.forward(observations))

# TODO: use model.log_probability() to calculate the all-paths likelihood of the
# observed sequence and then use np.exp() to convert log-likelihood to likelihood
probability_percentage = np.exp(model.log_probability(observations))

# Display the forward probabilities
print("         " + "".join(s.name.center(len(s.name)+6) for s in model.states))
for i in range(len(observations) + 1):
    print(" <start> " if i==0 else observations[i - 1].center(9), end="")
    print("".join("{:.0f}%".format(100 * forward_matrix[i, j]).center(len(s.name) + 6)
                  for j, s in enumerate(model.states)))

print("\nThe likelihood over all possible paths " + \
      "of this model producing the sequence {} is {:.2f}%\n\n"
      .format(observations, 100 * probability_percentage))

• 维特比算法

# TODO: input a sequence of 'yes'/'no' values in the list below for testing
observations = ['yes', 'no', 'yes']

# TODO: use model.viterbi to find the sequence likelihood & the most likely path
viterbi_likelihood, viterbi_path = model.viterbi(observations)

print("The most likely weather sequence to have generated " + \
      "these observations is {} at {:.2f}%."
      .format([s[1].name for s in viterbi_path[1:]], np.exp(viterbi_likelihood)*100)
)

• 对比两种算法

from itertools import product

observations = ['no', 'no', 'yes']

p = {'Sunny': {'Sunny': np.log(.8), 'Rainy': np.log(.2)}, 'Rainy': {'Sunny': np.log(.4), 'Rainy': np.log(.6)}}
e = {'Sunny': {'yes': np.log(.1), 'no': np.log(.9)}, 'Rainy':{'yes':np.log(.8), 'no':np.log(.2)}}
o = observations
k = []
vprob = np.exp(model.viterbi(o)[0])
print("The likelihood of observing {} if the weather sequence is...".format(o))
for s in product(*[['Sunny', 'Rainy']]*3):
    k.append(np.exp(np.log(.5)+e[s[0]][o[0]] + p[s[0]][s[1]] + e[s[1]][o[1]] + p[s[1]][s[2]] + e[s[2]][o[2]]))
    print("\t{} is {:.2f}% {}".format(s, 100 * k[-1], " <-- Viterbi path" if k[-1] == vprob else ""))
print("\nThe total likelihood of observing {} over all possible paths is {:.2f}%".format(o, 100*sum(k)))