14. HMM
1. HMM Background
机器学习大致分为两个派别,也就是频率派和贝叶斯派。
- 频率派,称为一门独立学科叫统计机器学习,核心问题是优化问题。主要有 3 个步骤:
- model,如 f(w)=wT+b 等。
- strategy,loss function。
- algorithm:GD/SGD/ 牛顿法 / 拟牛顿法。
- 贝叶斯派,衍生出概率图模型 PGM,最终是推断,即求解后验分布 p(z|x),因为贝叶斯公式中,分母需要对所有状态空间进行积分,非常复杂,必须要解决这个数值积分问题,才能让贝叶斯网络能实际运用。这类方法代表有 MCMC,MCMC 发展出来的一系列数值积分方法把贝叶斯理论引入实际应用中。
Hidden Markov Model 隐马尔可夫模型 在概率图模型中的背景:

离散动态模型中 Kalman filter 和 Particle filter 将会在后续章节介绍。
2. HMM 算法整体介绍
1 个基本模型
隐马尔可夫模型是关于 时序 的概率模型 ,描述一个 隐藏的马尔可夫链 随机 生成不可观测的状态随机序列 ,再由各个 状态 生成一个观测从而产生 观测随机序列的过程。HMM 结构图如下:

我们引入些符号来定义模型:
- 状态变量定义为:I=i1,i2,⋯,it,如图中 I 这个序列。
- 观测变量定义为:O=o1,o2,⋯,ot, 如图中 O 这个序列。
- 状态变量值的集合用 Q=q1,q2,⋯,qN, 就是说状态变量 it 的值域,每个观测变量 it 有 N 个状态可以取。
- 观测变量值的集合用 V=v1,v2,⋯,vM, 同样 V 是观测变量 ot 的值域。
- 状态转移矩阵 A=[aij]N×N, aij=P(it+1=qj|it=qi)
这里表示的就是状态变量从前一时刻 t 到下一时刻 t + 1 的概率值,就是 I 表示的空心圈状态序列从前往后变化的概率,根据上面定义的符号这是个 N×N 的矩阵。 - 发射矩阵(观测概率矩阵)B=[bj(k)]N×M,bj(k)=P(ot=vk|it=qj), 这里表示的是在 t 时刻,从状态变量 it=qj 转移到观测变量 ot=vk 的概率值,它是 N×M 的。
- 还有个初始状态概率 π=π1,π2,⋯,πN. 在 POS 中定义为每个词性出现在句首的概率。
这样,我们可以将 HMM 表示为一个三元组 λ=(π,A,B)。
2 个基本假设
另外,我们对于 HMM 有两个基本假设 李航——统计学习方法 P194
:
- 齐次马尔可夫性假设 ,就是假设隐藏状态马尔可夫链在任意时刻只与前一时刻状态有关。
P(it|it−1,ot−1,⋯,i1,o1)=P(it|it−1)
就是满足该假设情况下,我们不用考虑隐状态和观察状态关系,也不需要考虑隐状态之前所有状态关系,只用计算前一时刻。 - 观测独立性假设:即任意时刻的观测只依赖该时刻马尔科夫链的状态,与其他观测和其他状态无关。
P(ot|it,ot−1,it−1,⋯,o2,i2,o1,i1)=P(ot|it)
从图上看,就是计算 t 时刻的观测状态概率只需要考虑该时刻的隐状态,不需要其它状态和节点。
3 个基本问题
- 概率计算问题:上面部分我们给定的 HMM λ=(π,A,B) 和观测序列 O, 计算其在模型 λ 下,序列 O 出现的概率 P(O|λ)。具体来说,就是前向后向算法(接下来会具体介绍)。
- 学习问题: 已知观测序列 O,估计模型 λ=(π,A,B) 参数,使得在该模型下观测序列概率 P(O|λ) 最大。如 POS 中就是已知:
- 在有词性的初始状态(句首某个词性的概率)π,
- 当前词性转移到下一个词性的转移矩阵 A,
- 以及从该时刻的词性发射出具体某个词的发射矩阵 B 的情况下,有具体的输出序列 O,因为已经生成了序列 O, 那么有观测序列概率 P(O|λ) 最大。
- 具体求解 λ 算法有:EM 算法和 Baum Welch 算法。
- Decoding 问题:求解后的模型 λ=(π,A,B) 和观测序列 O=o1,o2,⋯,ot, 求给定观测序列条件概率 P(I|O) 最大的状态序列 I=i1,i2,⋯,it。简单来说,就是在给定模型和观测序列后,求解隐状态序列。这分为预测和滤波:
- 预测问题:P(it+1|o1,⋯,ot), 也就是在已知当前观测变量情况下预测下一个状态,具体就是 Viterbi 算法。
- 滤波问题:P(it|o1,⋯,ot), 已知观测变量序列求 t 时刻隐变量。
接下来,我们逐个介绍 3 个问题的具体求解过程。
3. HMM 概率计算问题
由上节,我们知道概率计算问题就是在给定模型参数 λ 情况下,求解 P(O|λ),总的来说一下几种计算方法:
- 直接计算法 我们利用边缘概率和联合概率乘法公式有:
P(O|λ)=∑IP(O,I|λ)=∑IP(O|I,λ)P(I|λ)
其中 ∑I 表示所有隐状态序列求和,P(O|I,λ) 表示对应隐状态产生观测状态序列的概率,P(I|λ) 是某个隐状态出现概率。
并且利用链式法则有:
P(I|λ)=P(i1,i2⋯,iT|λ)=P(iT|i1,i2⋯,iT−1,λ)⋯P(i2|i1,λ)P(i1|λ)
而根据齐次马尔可夫假设式 1 有:
P(iT|i1,i2⋯,iT−1,λ)=P(iT|iT−1)=aiT−1iT
那么式 4 根据式 5 依次代入替换有:
P(I|λ)=P(iT|i1,i2⋯,iT−1,λ)⋯P(i2|i1,λ)P(i1|λ)=aiT−1iTaiT−2iT−1⋯ai1i2πi1=πi1T∏t=2aiT−1iT
对于公式 3 我们只要求计算出 P(O|I,λ) 就可以了,利用观测独立性假设及联合概率有:
P(O|I,λ)=P(o1,o2,⋯,oT|I,λ)=T∏i=1P(ot|I,λ)=T∏i=1bit(ot)
现在公式 3 就等于:
P(O|λ)=∑IP(O|I,λ)P(I|λ)=∑Iπi1T∏t=2aiT−1iTT∏i=1bit(ot)=∑i1∑i2⋯∑iTπi1T∏t=2aiT−1iTT∏i=1bit(ot)
公式 8 最后一步只是为了说明这将会有 O(NT) 个求和,总共计算复杂度为 O(TNT)。在实际中是不可行的,为此我们引入前向 - 后向算法(forward-backward algorithm)。
- 前向 - 后向算法
前向算法

如上图橙色阴影所示,到时刻 t,隐状态变量为 it=qi, 观测状态序列为 o1,o2,⋯,ot 的概率,为前向概率,我们引入记号 αt(i)=P(o1,o2,⋯,ot,it=qi|λ) 来表示,更普遍一点表示为 αT(i)=P(O,iT=qT|λ)。
那么 P(O|λ)=∑Ni=1P(O,it=qi|λ) , 这表示时刻 t 时隐状态为 qi, 观测状态序列为 O 的前向概率。现在我们利用上述记号来推导:
αt+1(j)=P(o1,⋯,ot,ot+1,it+1=qj|λ)=N∑i=1P(o1,⋯,ot,ot+1,it+1=qj,it=qi|λ)=N∑i=1P(ot+1|o1,⋯,ot,it+1=qj,it=qi,λ)⋅P(o1,⋯,ot,it+1=qj,it=qi,λ)=N∑i=1P(ot+1|it+1=qj)⋅P(o1,⋯,ot,it+1=qj,it=qi,λ)式 2:观测独立性假设=N∑i=1P(ot+1|it+1=qj)⋅P(it+1=qj|o1,⋯,ot,it=qi,λ)⋅P(o1,⋯,ot,it=qi|λ)=N∑i=1bj(ot+1)aijαt(i)
这表示,t 时刻隐状态为 qi, 观测状态为 O 乘以从隐状态 i 转移到 j 的概率再乘以由隐状态变量 it+1=qj 观测到 ot+1 的概率再求和。
下面借用 14 隐马尔可夫模型.pdf笔记中图来加强理解,对于前向算法来说,每个时刻 t 的输出状态都等于上一个隐状态分别计算概率再求和 。这实际上跟神经网络非常类似,就是前一状态值乘以权重。因为隐状态可取值为 N 个,乘以对应概率,在对总共的序列长度 T 求和,因此总的时间复杂度为 O(TN2),减少了大量计算。还可以看看 统计机器学习 P200 例 10.2
。

后向算法

如上图所示,我们记 βt(i)=P(ot+1,⋯,oT|it=qi,λ), 就是 t 时刻在隐状态 it=qi 条件下,从 t+1 到 T 部分观测序列为 ot+1,⋯,oT 的概率,这就是后向概率。那么 βt(1)=P(o2,⋯,oT|i1=qi,λ)。若要计算 P(O|λ),那么有:
P(O|λ)=P(o1,o2,⋯,oN|λ)=N∑i=1P(o1,o2,⋯,oN,i1=qi|λ)边缘概率=N∑i=1P(o1|o2,⋯,oN,i1=qi,λ)P(i1=qi|λ)⏟πi条件概率=∑i=1NP(o1|o2,⋯,oN,i1=qi,λ)⋅P(o2,⋯,oN,i1=qi|λ)⋅πi=N∑i=1P(o1|i1=qi,λ)⋅β1(i)⋅πi =N∑i=1bi(o1)⋅β1(i)⋅πi=N∑i=1πi⋅bi(o1)⋅β1(i)
其中,πi 为某个初始状态的概率,bi(o1) 第一个隐状态条件下生成第一个观测状态概率,β1(i) 第 1 个观测状态确定后,生成后面观测序列的概率。另外,对于上述推导过程中,可以将 λ 省去再思考推导过程。依旧借用 14 隐马尔可夫模型.pdf 中图如下,要确定输出观测序列,我们要知道初始隐状态概率,该隐状态转化为第一个观测状态概率,以及第一个观测条件下,生成后续观测序列的可能性只要将其相乘就可以了,分别对应下图蓝色,绿色,红色部分,并且可以看作一个信息传递过程。

现在,我们要推广到 βt(i) 情况,这里省去 λ,方便推导表达式。
βt(i)=P(O|λ)=P(ot+1,ot+2,⋯,oT|it=qi)=N∑j=1P(ot+1,ot+2,⋯,oT|it+1=qj,it=qi)P(ot+1,ot+2,⋯,oT,it+1=qj|it=qi)链式法则=N∑j=1P(ot+1,ot+2,⋯,oT|it+1=qj,it=qi)⋅P(it+1=qj|it=qi)齐次假设=N∑j=1P(ot+1,ot+2,⋯,oT|it+1=qj)⋅aij=N∑j=1P(ot+1|ot+2,⋯,oT,it+1=qj)⋅P(ot+2,⋯,oT|it+1=qj)⋅aij=N∑j=1bj(ot+1)⋅βt+1(j)⋅aij=N∑j=1aij⋅bj(ot+1)⋅βt+1(j)
跟上面一样的递推过程,依旧如下图的后向算法计算过程。

例:考虑盒子和球模型 λ=(π,A,B),状态集合 Q=1,2,3,观测集合 V=红,白。(李航《统计学方法》P200)
A=[0.50.20.3 0.30.50.2 0.20.30.5],B=[0.50.5 0.40.6 0.70.3],π=(0.2,0.4,0.4)T
设观测序列长度 T=3,观测序列为:O=(红,白,红),求:P(O|λ)。代码来自 HMM-Evaluation.
import numpy as np | |
def forward(lamda, O): | |
''' | |
HMM 模型:概率计算问题——前向算法,李航《统计学方法》P177 | |
step1: | |
t= 0 时刻,计算初值:alpha[t=1][j] = pai(j) * B[j][o1] | |
step2: | |
递推:t=1, 2, ..., T-1,计算 alpha[t][j]:alpha[t][j] = ∑(alpha[t-1][i] * A[i][j]) * B[j][ot] | |
step3: | |
计算:P(O|λ) = ∑alpha[T][j] | |
''' | |
A, B, pai = lamda # HMM 的参数:A, B, pai | |
N, M = B.shape # 状态集合 Q 的长度 N,观测集合 V 的长度 M | |
T = len(O) # 由观测序列得到 T 个时刻 | |
'''HMM 的前向概率矩阵:T*N''' | |
alphas = np.zeros((T, N)) | |
'''step1:t= 0 时刻,计算初值:alpha[t=0][j] = pai(j) * B[j][o1]''' | |
for j in range(N): # j 是当前状态 | |
k = V.index(O[0]) # 第 Ot 个观测对应的 V 观测集合的位置 k:O[0]=V[k] | |
alphas[0][j] = pai[j] * B[j][k] | |
'''step2:递推:t=1, 2, ..., T-1,计算:alpha[t][j] = ∑(alpha[t-1][i] * A[i][j]) * B[j][Ot]''' | |
for t in range(1, T): # j 是当前状态 | |
k = V.index(O[t]) # 第 Ot 个观测对应的 V 观测集合的位置 k:O[t]=V[k] | |
for j in range(N): | |
temp = 0 | |
for i in range(N): | |
temp += alphas[t - 1][i] * A[i][j] # 计算:∑(alpha[t-1][i] * A[i][j]) | |
alphas[t][j] = temp * B[j][k] # 计算:∑(alpha[t-1][i] * A[i][j]) * B[j][Ot] | |
'''step3:计算:P(O|λ) = ∑alpha[T][j]''' | |
prob = 0 # 求解 P(O|λ) | |
for j in range(N): # j 是当前隐状态 | |
prob += alphas[T - 1][j] | |
return prob | |
if __name__ == '__main__': | |
'''状态集合:Q = {q1, q2, ..., qN}''' | |
Q = ['盒子 1', '盒子 2', '盒子 3'] | |
'''观测集合:V = {v1, v2, ..., vM}''' | |
V = ['红', '白'] | |
''' | |
状态转移概率矩阵:A = [aij],其中 aij = P(it+1=qj | it=qi) | |
A: shape N * N (盒子之间互相转换) | |
''' | |
A = np.array([[0.5, 0.2, 0.3], | |
[0.3, 0.5, 0.2], | |
[0.2, 0.3, 0.5]]) | |
''' | |
观测概率矩阵:B=[bj(k)],其中 bj(k) = P(ot=vk | it=qj) | |
B: shape N * M (N 个盒子 M 个观测) | |
''' | |
B = np.array([[0.5, 0.5], | |
[0.4, 0.6], | |
[0.7, 0.3]]) | |
'''初始状态概率向量:pai = P(i1=qi):第一次出现每个盒子的概率分布''' | |
pai = np.array([0.2, 0.4, 0.4]) | |
'''随机序列 1——观测序列:O = {o1, o2, ..., oT}''' | |
'''随机序列 2——状态序列(隐序列,叠加态的黑盒子):I = {i1, i2, ..., iT}''' | |
O = np.array(['红', '白', '红']) | |
'''前向算法计算:P(O|λ)''' | |
P_O_or_lamda = forward(lamda=[A, B, pai], O=O) | |
print('通过 HMM 前向算法计算概率:\n 观测序列:{} 出现的概率 P(O|λ) = {:.4f}'.format(O, P_O_or_lamda)) |
4. Viterbi 算法
隐马尔可夫模型中的维特比算法,用于解码给定的观测序列和模型参数,得到最有可能的隐藏状态序列。
该代码实现了 HMM 模型的三个基本问题:
- 概率计算问题:给定模型参数和观测序列,计算这个序列出现的概率。
- 学习问题:给定观测序列,估计模型参数,使得这个观测序列出现的概率最大。
- 解码问题:给定模型参数和观测序列,找到最有可能的隐藏状态序列。
具体来说,该代码实现了维特比算法(Viterbi algorithm),用于解码问题,即在给定模型参数和观测序列的情况下,找到最有可能的隐藏状态序列。
代码中的 viterbi
函数接受两个参数:一个观测序列 obs
,一个 HMM 模型 model
。model
包括三个参数:初始状态概率向量 model['start_prob']
、状态转移概率矩阵 model['trans_prob']
、观测概率矩阵 model['emit_prob']
。
函数中的主要步骤是:
- 初始化:将起始概率向量
model['start_prob']
和第一个观测值对应的观测概率向量model['emit_prob'][:, obs[0]]
相乘,得到一个长度为状态数的向量probs
,并将其作为维特比算法的第一个状态。 - 递推:对于观测序列中的每个观测值,根据当前状态向量
probs
和状态转移概率矩阵model['trans_prob']
,计算出下一个状态向量next_probs
,并将其作为下一次迭代的当前状态向量。 - 回溯:在所有观测值都处理完之后,根据每个状态向量的最大值以及指向最大值的箭头,得到最有可能的隐藏状态序列。
import numpy as np | |
# Define the Viterbi algorithm | |
def viterbi(obs, states, obs_symbols, initial_prob, transitions, emissions): | |
""" | |
viterbi algorithm to get best_path and prob_best_path | |
:param obs: 要已经划分完毕的标注的句子, split/tokenized | |
:param states: 词性的所有状态 | |
:param obs_symbols: 所有词的列表 | |
:param initial_prob: 初始概率矩阵 | |
:param transitions: 转移矩阵 | |
:param emissions: 发射矩阵 | |
:return: | |
""" | |
# Define the T x N matrix for probabilities | |
T = len(obs) #要标注的句子长度 | |
N = len(states) #pos tag 长度 | |
probs = np.zeros((T, N)) #保存句子每个词到指定 tag 的概率 | |
# Define the T x N matrix for pointers to backtrack | |
pointers = np.zeros((T, N), dtype=int) | |
# Calculate probabilities and pointers for the first observation | |
# 初始化第一个词对应词性 prob | |
obs_index = obs_symbols.index(obs[0]) | |
probs[0] = np.multiply(initial_prob, emissions[:, obs_index]) | |
# Iterate over the remaining observations | |
for t in range(1, T): | |
#输入句子的词对应在观察序列列表的索引 | |
obs_index = obs_symbols.index(obs[t]) | |
for s in range(N): | |
#第 t 个词前一个词的概率乘以转移概率. 这里从 1 循环到 T 句子长度的过程中都是拿前面一个 | |
#来获得当前词性的概率, 已经初始化得到 0,那么 1 就是 probs[0] *transitions[:, s] | |
#迭代获取所有当前词性的概率 | |
temp_probs = np.multiply(probs[t - 1], transitions[:, s]) | |
#取最大概率的 index | |
max_prob_index = np.argmax(temp_probs) | |
#那么对应最大概率 = 发射矩阵矩阵中词 * 当前词性最大概率 | |
probs[t, s] = emissions[s, obs_index] * temp_probs[max_prob_index] | |
#为解码取路径保存索引来 这个索引表示 t 时刻对应上一时刻的最大的 prob 连接的索引 | |
pointers[t, s] = max_prob_index | |
# Find the best path | |
#最后一行的最大的索引 | |
best_path = [np.argmax(probs[-1])] | |
for t in range(T - 1, 0, -1): | |
best_path.append(pointers[t, best_path[-1]]) | |
best_path.reverse() | |
# Calculate the probability of the best path | |
prob_best_path = np.max(probs[-1]) | |
return best_path, prob_best_path | |
# Example usage | |
obs = ['walk', 'shop', 'clean'] | |
# Define states | |
states = ('Rainy', 'Sunny') | |
# Define observation symbols | |
obs_symbols = ('walk', 'shop', 'clean') | |
# Define initial probabilities | |
initial_prob = np.array([0.6, 0.4]) | |
# Define transition matrix | |
transitions = np.array([[0.7, 0.3], | |
[0.4, 0.6]]) | |
# Define emission matrix | |
emissions = np.array([[0.1, 0.4, 0.5], | |
[0.6, 0.3, 0.1]]) | |
best_path, prob_best_path = viterbi(obs, states, obs_symbols, initial_prob, transitions, emissions) | |
print("Best path:", [states[i] for i in best_path]) | |
print("Probability of the best path:", prob_best_path) |
Inference
[1] HMM 前向与后向算法
[3] HMM-Evaluation