机器学习之隐马尔科夫（HMM）模型
机器学习之隐马尔科夫模型（HMM）
- 1、对隐马尔可夫模型进行概述
- 2、阐述其数学理论基础
- 3、详细说明Python代码实现过程
- 4、归纳总结其核心内容和应用价值
隐马尔可夫模型介绍
隐式马尔可夫模型（Hidden Markov Model, HMM）是一种基于时间序列的概率建模方法。该方法主要研究隐藏的马尔可夫过程如何通过产生不可观察的状态序列来影响可观测数据序列的统计特性，并且其核心在于通过每个状态产生观测数据进而形成观测数据序列的过程。这种建模方法主要采用一种生成性建模策略
下面我们来从概率学角度定义马尔科夫模型，从一个典型例子开始：
假设有四个不同的容器（如四个盒子），每个容器内分别装有两种类型的乒乓球——红色与白色的小球。这些容器内的乒乓球数目各不相同。具体数据列于下表中：
| 盒子编号 | 1 | 2 | 3 | 4 |
|---|---|---|---|---|
| 红球数 | 5 | 3 | 6 | 8 |
| 白球数 | 5 | 7 | 4 | 2 |
从这些盒子中依次取出T个球，在这一过程中遵循以下取样规则：每次随机选择一个盒子并从中取出一球，并记录其颜色后将其放回原处。在这一过程中我们无法得知具体来自哪个盒子的信息即无法得知a box sequence, also known as a state sequence; the color sequence observed is referred to as an observation or measurement sequence here.这两个随机序列共同构成了一个马尔科夫过程的例子。
称\mathbf{Q}为系统中所有可能状态构成的集合，则V则代表所有可能观测构成的集合：
\mathbf{Q} = \{\mathbf{q}_1, \mathbf{q}_2, \dots, \mathbf{q}_N\}, \quad V = \{v_1, v_2, \dots, v_M\}
其中,N代表系统中的状态总数,M则代表系统中观测的数量,比如在上述例子中可得N=4,M=2。
\mathbf{I}被定义为长度为T的状态序列；而\mathbf{O}则表示相应的观测序列：
B被定义为观测概率矩阵：
B = [b_j(k)]_{N \times M}
其中，
b_j(k)表示在时刻t处于状态q_j的条件下生成观测v_k的概率，
k的取值范围为1到M，
j的取值范围为1到N。
\pi表示初始状态的概率向量：
\pi = (\pi_i)
其中，\pi_i = P(i_1 = q_i), i=1,2,\cdots,N 表示在时刻t=1处于状态q_i的概率。
由此可知，隐马尔科夫模型\lambda采用三元符号体系其三元符号体系是什么即
\lambda = (A,B,\pi)
其中A、B和\pi被定义为隐马尔科夫模型的基本要素
由定义可知隐马尔可夫模型做了两个基本假设：
(1) 基于齐次马尔科夫性质假说，在任意时刻t的状态仅依赖于其前一个状态；
P(i_t \mid i_{t-1}, o_{t-1}, \dots, i_1, o_1) = P(i_t \mid i_{t-1}), \quad t=1,2,\dots,T
(2) 根据观测独立性假说，在各个时间点上观测结果仅与当前时刻的状态相关；
如前所述，在盒子取球的例子中
状态集合：Q = {盒子1，盒子2，盒子3，盒子4}， N=4
观测集合：V = {红球，白球} M=2
初始化概率分布：
\pi = (0.25,0.25,0.25,0.25)^T
状态转移矩阵：
A = \left[ \begin{matrix} 0 & 1 & 0 & 0 \\ 0.4 & 0 & 0.6 & 0 \\ 0 & 0.4 & 0 & 0.6 \\ 0 & 0 & 0.5 & 0.5 \end{matrix} \right]
观测矩阵:
B = \left[ \begin{matrix} 0.5 & 0.5 \\ 0.3 & 0.7 \\ 0.6 & 0.4 \\ 0.8 & 0.2 \end{matrix} \right]
隐马尔可夫模型的三个基本问题
1、概率计算问题
给定：\lambda = (A,B,\pi) O=(o_1,o_2,\cdots,o_T)
计算：P(O|\lambda)
2、学习问题
已知：O=(o_1,o_2,\cdots,o_T)
估计：\lambda = (A,B,\pi),使P(O|\lambda)最大
3、预测问题（解码）
已知：\lambda = (A,B,\pi) O=(o_1,o_2,\cdots,o_T)
求：使P(I|O)最大的状态序列I=(i_1,i_2,\cdots,i_T)
下面我们使用python代码写一个HMM模型生成序列O的示例代码
import numpy as np
class HMM(object):
def __init__(self, N, M, pi=None, A=None, B=None):
self.N = N
self.M = M
self.pi = pi
self.A = A
self.B = B
def get_data_with_distribute(self, dist): # 根据给定的概率分布随机返回数据（索引）
r = np.random.rand()
for i, p in enumerate(dist):
if r < p: return i
r -= p
def generate(self, T: int):
'''
根据给定的参数生成观测序列
T: 指定要生成观测序列的长度
'''
result = []
for ind in range(T): # 依次生成状态和观测数据
if ind==0:
i = self.get_data_with_distribute(self.pi)
else:
i = self.get_data_with_distribute(self.A[i])
o = self.get_data_with_distribute(self.B[i])
result.append(o)
return result
if __name__ == "__main__":
pi = np.array([0.25, 0.25, 0.25, 0.25])
A = np.array([
[0, 1, 0, 0],
[0.4, 0, 0.6, 0],
[0, 0.4, 0, 0.6],
[0, 0, 0.5, 0.5]])
B = np.array([
[0.5, 0.5],
[0.3, 0.7],
[0.6, 0.4],
[0.8, 0.2]])
hmm = HMM(4, 2, pi, A, B)
print(hmm.generate(10)) # 生成10个数据
# 生成结果如下
[0, 0, 0, 0, 1, 1, 1, 1, 0, 0] # 0代表红球，1代表白球隐马尔可夫模型数学原理
下面我们从隐马尔可夫模型的三个基本问题出发，逐个解释问题解法：
1、概率计算问题
* 直接计算法
* 前向算法
* 后向算法直接计算法：
状态序列I=(i_1,i_2,\cdots,i_T)概率：P(I|\lambda)=\pi_{i_1} a_{i_1 i_2} a_{i_2 i_3}\cdots a_{i_{T-1}i_T}
给定固定的观测序列I的概率度：P(O|I,\lambda)
给定固定的观测序列I的概率度为各时刻条件概率的连乘积。
表示O与I共同出现的概率为：
P(O,I|\lambda) = P(O|I,\lambda)P(I|\lambda) = \pi_{i_1}b_{i_1}(o_1)a_{i_1 i_2} b_{i_2}(o_2) \cdots a_{i_{T-1}i_T}b_{i_T}(o_T)
其中每个因子对应于时间序列中的一个状态转移过程
所有可能的状态序列I被求和以计算观测序列O的概率：
公式保持不变：
P(O|\pi) = \sum_I{P(O|i,\lambda)P(I|\lambda)} = \sum_{i_1,i_2,\cdots,i_T}\pi_{i_1}b_{i_1}(o_1)a_{i_1 i_2} b_{i_2}(o_2) \cdots a_{i_{T-1}i_T}b_{i_T}(o_T)
此算法的复杂度为O(TN^T)，计算量太大，不可行。
前向算法：
前向概率定义：基于隐马尔可夫模型\lambda，在时刻t的部分观测序列o_1,o_2,\cdots,o_t中状态处于q_i的概率作为前向概率记作\alpha_t(i)。即：
\alpha_{t}(i) = P(o_1,o_2,\cdots,o_{t},i_{t}=q_{i}|\lambda)
其中，
输入包括：
\lambda \text{ 隐马尔可夫模型}
以及：
O \text{ 观测序列}
输出：观测序列概率P(O|\lambda)
初值：\alpha_1(i) = \pi_ib_i(o_1)，i=1,2,\cdots,N
递推：\alpha_{t+1}(i) = \left[ \sum_{j=1}^N \alpha_t(j) a_{_ji} \right]b_i(o_{t+1}),i=1,2,\cdots,N
终止：P(O|\lambda) = \sum_{i=1}^N \alpha_T(i)
算法复杂度O(N^2T)，算法复杂度大大降低
因为：\alpha_T(i) = P(o_1,o_2,\cdots,o_T,i_T=q_i|\lambda)
所以：P(O|\lambda) = \sum_{i=1}^N \alpha_T(i)
以下通过实例说明前向算法的操作流程：
给定状态转移概率矩阵A = \left[ \begin{matrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{matrix} \right]，
观测概率矩阵B = [b_j(o)]，
初始状态概率分布\pi=[\pi_1,\pi_2,\pi_3]，
其中观测序列O=(红,白,红)，模型参数包括状态转移矩阵A、观测概率矩阵B以及初始状态分布π。
模型λ=(A,B,π)，状态集合Q={1,2,3}, 观测集合V={红、白}
假设时间步数T=3
运用前向算法求得观测序列概率P(O|λ)
（1）计算初值
\alpha_1(1) = \pi_1 b_1(o_1) = 0.2\times0.5 = 0.1 \\ \alpha_1(2) = \pi_2 b_2(o_1) = 0.4\times0.4 = 0.16 \\ \alpha_1(3) = \pi_3 b_3(o_1) = 0.4\times0.7 = 0.28
（2）递推计算
\begin{aligned} \alpha_2(1) & = \left[ \begin{matrix} \sum_{i=1}^3\alpha_1(i)a_{i1} \end{matrix} \right]b_1(o_2) = (0.1\times0.5+0.16\times0.3+0.28\times0.2)\times0.5 = 0.154\times0.5 = 0.077 \\ \alpha_1(2) & = \left[ \begin{matrix} \sum_{i=1}^3\alpha_1(i)a_{i2} \end{matrix} \right]b_2(o_2) = 0.184\times0.6 = 0.1104 \\ \alpha_1(3) & = \left[ \begin{matrix} \sum_{i=1}^3\alpha_1(i)a_{i3} \end{matrix} \right]b_3(o_2) = 0.202\times0.3 = 0.0606 \\ \alpha_3(1) & = \left[ \begin{matrix} \sum_{i=1}^3\alpha_2(i)a_{i1} \end{matrix} \right]b_1(o_3) = 0.04187 \\ \alpha_3(2) & = \left[ \begin{matrix} \sum_{i=1}^3\alpha_2(i)a_{i2} \end{matrix} \right]b_2(o_3) = 0.03551 \\ \alpha_3(3) & = \left[ \begin{matrix} \sum_{i=1}^3\alpha_2(i)a_{i3} \end{matrix} \right]b_3(o_3) = 0.05284 \end{aligned}
（3）终止
P(O|\lambda) = \sum_{i=1}^3 \alpha_3(i) = 0.13022
后向算法：
后向概率定义：基于隐马尔可夫模型λ，在时刻t状态下处于状态qi的情况下（即i_t = qi），从时刻t+1到T的部分观测序列为o_{i+1}, o_{i+2}, ..., o_T的概率被定义为后向概率，并将其记作：
\beta_t(i) = P(o_{t+1}, o_{t+2}, \ldots, o_T \mid i_t = q_i, \lambda)
同样地，在这种情况下也可以通过建立递推关系式来计算相应的后向概率\beta_t(i)以及完整的观测序列的概率P(O|\lambda)。
输入：隐马尔科夫模型\lambda，观测序列O；
输出：观测序列概率P(O|\lambda).
（1） \beta_T(i) = 1, i=1,2,\cdots,N
（2）对t=T-1,T-2,\cdots,1
\beta_t(i)=\sum_{j=1}^N a_{ij}b_j(o_{t+1})\beta_{t+1}(j), i=1,2,\cdots,N
（3） P(O|\lambda)=\sum_{i=1}^N \pi_ib_i(o_1)\beta_1(i)
前后向统一写为：（t=1和t=T-1分别对应）
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), t=1,2,\cdots,T-1
一些概率和期望值计算方式：
1、给定模型\lambda和观测O，在时刻t处于状态q_i的概率，记\gamma_t(i)=P(i_t=q_t|O,\lambda)
\gamma_t(i)=P(i_t=q_t|O,\lambda)=\frac{P(i_t=q_t,O|\lambda)}{P(O|\lambda)}=\frac{\alpha_t(i)\beta_t(i)}{\sum_{j=1}^N \alpha_t(j)\beta_t(j)}
2、给定模型\lambda和观测O，在时刻t处于状态q_i且在时刻t+1处于q_j概率，记\xi_t(i,j)=P(i_t=q_i,i_{t+1}=q_j|O,\lambda)
\xi_t(i,j)=P(i_t=q_i,i_{t+1}=q_j,O|\lambda)=\frac{P(i_t=q_t,i_{t+1}=q_j,O|\lambda)}{P(O|\lambda)}=\frac{P(i_t=q_t,i_{t+1}=q_j,O|\lambda)}{\sum_{i=1}^N \sum_{j=1}^N P(i_t=q_i,i_{t+1}=q_j,O|\lambda)} = \frac{\alpha_t(i)a_{ij}b_j(o_{t+1}) \beta_{t+1}(j)}{\sum_{i=1}^N \sum_{j=1}^N \alpha_t(i)a_{ij}b_j(o_{t+1}) \beta_{t+1}(j)}
3、由1和2可得，(1)在观测O下状态i出现的期望值：\sum_{t=1}^T \gamma_t(i)，（2）在观测O下由状态i转移的期望值：\sum_{t=1}^{T-1} \gamma_t(i)，（3）在观测O下由状态i转移到状态j的期望值：\sum_{t=1}^{T-1} \xi_t(i,j)。
2、学习问题
监督学习方法：
- 称训练数据集由观测序列O与对应的状态序列I共同构成，并表示为{(O₁,I₁),(O₂,I₂),⋯,(O_S,I_S)}
- 可以通过极大似然估计法来估计隐马尔可夫模型的参数。
已知：{(O_1,I_1)(O_2,I_2),\cdots,(O_S,I_S)}
（1）转移概率a_{ij}的估计：假设样本中时刻t处于状态i，时刻t+1转移到状态j的频数为A_{ij}那么转台转移概率a_{ij}的估计是：
\hat{a}_{ij} = \frac{A_{ij}}{\sum_{j=1}^N A_{ij}}, i=1,2,\cdots,N;j=1,2,\cdots,N
（2）观测概率b_j(k)的估计：设样本中状态为j并观测为k的频数是B_j(k)那么状态j观测为k的概率
\hat{b}_j(k) = \frac{B_{jk}}{\sum_{k=1}^M B_{jk}}， j=1,2,\cdots,N;k=1,2,\cdots,M
（3）初始状态概率\pi_i的估计\hat{\pi}_i为S个样本中初始状态为q_i的频率。
非监督学习方法：
Baum-Welch算法
假设训练数据仅包含{O₁,O₂,…,Oₛ}，请估计模型参数λ=(A,B,π)，这些参数本质上代表了一个含有隐变量的概率模型；其中EM算法用于计算概率P(O|λ)，具体而言就是对所有可能的状态序列I取加权平均。
（1）确定完全数据的对数似然函数
完全数据(O,I)=(o_1,o_2,\cdots,o_T,i_1,i_2,\cdots,i_r)
完全数据的对数似然函数logP(O,I|\lambda)
（2）在EM算法中进行E步时,计算Q函数
Q(\theta,\hat\theta)=\mathbb{E}_{Z|X;\hat\theta}\left[\log P(X,Z;\theta)\right]
其中\mathbb{E}_{Z|X;\hat\theta}表示在参数\hatθ下给定观测数据X对隐变量Z的条件期望；进而得出：
J(θ)=\mathbb{E}_{Z|X;θ̂}[logP(X,Z;θ)]
其中\mathbb{E}_{Z|X;θ̂}表示在当前参数θ下给定观测数据X对隐变量Z的条件期望；针对观测数据X的具体应用场景实施相关操作
（3）EM算法的M步，极大化Q(\lambda,\bar{\lambda})求模型参数A,B,\pi
第一项\sum_I log\pi_{i_0}P(O,I|\hat{\lambda})=\sum_{i=1}^N log\pi_i P(O,i_1=i|\bar{\lambda}) ，基于约束条件：\sum_{i=1}^N \pi_i = 1借助拉格朗日乘子法构造拉格朗日函数\sum_{i=1}^N log\pi_iP(O,i_1=i|\bar{\lambda})+\gamma(\sum_{i=1}^N \pi_i -1)并对该函数关于各\pi_i求偏导数并令其等于零
对拉格朗日函数\sum_{i=1}^N log\pi_i P(O,i_1=i|\bar{\lambda})+\gamma(\sum_{i=1}^N \pi_i -1)分别关于每个概率\pi_i求偏导数得到方程组：
\frac{\partial}{\partial \pi_i}\left[\sum_{i=1}^N log\pi_i P(O,i_1=i|\bar{\lambda})+\gamma(\sum_{i=1}^N \pi_i -1)\right]=0
解此方程组得到关系式：
P(O,i_1 = i|\bar{\lambda}) + \gamma\pi_i = 0
其中常数因子满足\gamma = -P(O|\bar{\lambda}) ，进而得到概率分布：
\pi_i = \frac{P(O,i_1=i|\bar{\lambda})}{P(O|\bar{\lambda})}
第二项可写成\sum_I(\sum_{t=1}^{T-1}log a_{i_t i_{t+1}})P(O,I|\bar{\lambda})=\sum_{i=1}^N \sum_{j=1}^N \sum_{t=1}^{T-1} log a_{ij}P(O,i_t=i,i_{t+1}=j|\bar{\lambda})由约束条件\sum_{j=1}^N a_{ij}=1，拉格朗日乘子法得：
a_{ij}=\frac{\sum_{t=1}^{T-1}P(O,i_t=i,i_{t+1}=j|\bar{\lambda})}{\sum_{t=1}^{T-1}P(O,i_t=i|\bar{\lambda})}
第三项\sum_I(\sum_{t=1}^T logb_{i_t}(o_t))P(O,I|\bar{\lambda})=\sum_{j=1}^N \sum_{t=1}^T logb_j(o_t)P(O,i_t=j|\bar{\lambda})由约束条件\sum_{k=1}^M b_j(k)=1，注意，只有在o_t=v_k时b_j(o_t)对b_j(k)的偏导数才不为0，以I(o_t=v_k)表示，求得
b_j(k)=\frac{\sum_{t=1}^T P(O,i_t=j|\bar{\lambda})I(o_t=v_k)}{\sum_{t=1}^N P(O,i_t=j|\bar{\lambda})}
将以上得到的概率分别用\gamma_t(i)，\xi_t(i,j)表示
a_{ij}=\frac{\sum_{t=1}^{T-1}\xi_t(i,j)}{\sum_{t=1}^{T-1}\gamma_t(i)} \\ b_j(k)=\frac{\sum_{t=1,o_t=v_k}^T \gamma_t(j)}{\sum_{t=1}^T \gamma_t(j)} \\ \pi_i = \gamma_1(i)
Baum-Welch算法流程：
输入：观测数据O=(o_1,o_2,\cdots,o_T)；
输出：隐马尔可夫模型参数。
（1）初始化
对n=0，选取a_{ij}^{(0)}，b_j(k)^{(0)}，\pi_i^{(0)}得到模型\lambda^{(0)}=(A^{(0)},B^{(0)},\pi^{(0)})
（2）递推过程基于观测序列O=(o₁,o₂,…,o_T)和模型λ(n)=(A(n),B(n),π(n))进行计算：
对于每个状态变量i,j和时间步长n=1,2,…,
a_{ij}^{(n+1)} = \frac{\sum_{t=1}^{T-1}\xi_t(i,j)}{\sum_{t=1}^{T-1}\gamma_t(i)}
b_j(k)^{(n+1)} = \frac{\sum_{t=1,o_t=v_k}^T \gamma_t(j)}{\sum_{t=1}^T \gamma_t(j)}
\pi_i^{(n+1)} = \gamma_1(i)
（3）终止，得到模型参数\lambda^{(n+1)}=(A^{(n+1)},B^{(n+1)},\pi^{(n+1)})
3、预测问题
预测算法
（1）近似算法
想法：每个时间点t都挑选该时间点之后可能出现的状态i_t^*进行评估，并由此生成一个完整的状态序列I^*=(i_1^*,i_2^*,\cdots,i_T^*)作为预测结果；同时，在时间点t处于状态q_i的概率值\gamma_t(i)被定义为\frac{\alpha_t(i) \beta_t(i)}{P(O|\lambda)}的形式，并进一步简化为\frac{\alpha_t(i) \beta_t(i)}{\sum_{j=1}^N \alpha_t(j) \beta_t(j)}；在此基础上确定每个时间点t的最可能状态i_t^*的方法是通过最大化\gamma_t(i)值来实现的
从而得到状态序列：I_*=(i_1^*,i_2^*,\cdots,i_T^*)
得到的状态有可能实际不发生
（2）维比特算法
通过动态规划方法求解概率最大的路径问题时，每条路径都对应着一个状态序列。
最优路径具有以下特征：当最优路径在时刻t经过节点i_t^*时，则该路径从节点i_t^*延伸至终点i_T^*的所有部分中，在t到T阶段的各部分中达到最优。
只需从t=1开始逐步递推，在每个时刻t状态下计算所有可能部分的最大概率值。最终，在t=T时获得的最大概率值即为整体最优解的概率值P^*. 同时，在此阶段也确定了最终的终结点位置。
接下来，在确定了最终节点后需逆向推导出完整的最优路线。具体方法是从t=T-1开始依次向前追溯至t=1.
导入两个变量\delta和\psi，定义在时刻t状态为i的所有单个路径(i_1,i_2,\cdots,i_t)中概率最大值为：
\delta_t(i)=\underset{i_1,i_2,\cdots,i_{t-1}}{\operatorname{max}}P(i_t=i,i_{t-1},\cdots,i_1,o_t,\cdots,o_1|\lambda),i=1,2,\cdots,N
由定义可得变量\delta的递推公式：
\delta_{t+1}(i)=\underset{i_1,i_2,\cdots,i_t}{\operatorname{max}}P(i_{t+1}=i,i_t,\cdots,i_1,o_{t+1},\cdots,o_1|\lambda)=\underset{1\leq j \leq N}{\operatorname{max}}[\delta_t(j)a_{ji}]b_i(o_{t+1}), i=1,2,\cdots,N;t=1,2,\cdots,T-1
定义在时刻t状态为i的所有单个路径(i_1,i_2,\cdots,i_{t-1},i)中概率最大的路径的第t-1个结点为
\psi_t(i)=\underset{1\leq j \leq N}{\operatorname{argmax}}[\delta_{t-1}(j)a_{ji}], i=1,2,\cdots,N
输入：模型\lambda=(A,B,\pi)和观测O=(o_1,o_2,\cdots,o_T)；
输出：最优路径I^*=(i_1^*,i_2^*,\cdots,i_T^*).
①初始化
\delta_1(i)=\pi_ib_i(o_1), i=1,2,\cdots,N \\ \psi_1(i)=0, i=1,2,\cdots,N
②递推，对t=2,3,\cdots,T
\delta_t(i)=\underset{1\leq j \leq N}{\operatorname{max}}[\delta_{t-1}(j)a_{ji}]b_i(o_t), i=1,2,\cdots,N \\ \psi_t(i)=\underset{1\leq j \leq N}{\operatorname{argmax}}[\delta_{t-1}(j)a_{ji}], i=1,2,\cdots,N
③终止
P*=\underset{1\leq i \leq N}{\operatorname{max}}\delta_T(i) \\ i_T^*=\underset{1\leq i \leq N}{\operatorname{argmax}}[\delta_T(i)]
④最优路径回溯，对t=T-1,T-2,\cdots,1
i_t^*=\psi_{t+1}(i_{t+1}^*)
求得最优路径I^*=(i_1^*,i_2^*,\cdots,i_T^*)
矩阵A= \left[\begin{matrix} 0.5 & 0.2 & 0.3 \\ 0.3 & 0.5 & 0.2 \\ 0.2 & 0.3 & 0.5 \end{matrix}\right]、矩阵B = \left[\begin{matrix} 1/2 &1/2 \\4/1&6/1\\7/1&3/1\end{matrix}\right]、初始概率向量\pi = (1/5,2/5,2/5)^T已给出。请寻求最佳路径序列I^\ast=(i_1^\ast,i_2^\ast,i_3^\ast)。
① 初始化阶段，在时间步t=1时，对于每一个状态变量i（其中i=1,2,3），计算其观测值o₁为红色的概率，并将其定义为δ₁(i) = π_i b_i(o₁) = π_i b_i(红)，其中i分别取值于1、2和3。
代入实际数据：\delta_1(1)=0.10，\delta_1(2)=0.16，\delta_1(3)=0.28，记\psi_1(i)=0,i=1,2,3
当时间步t=2时，在每个状态i（i=1, 2, 3）处进行如下运算：计算从时间步t=1到t=2的状态转移概率，并考虑观察结果O_1和O_2的影响。该概率被计算为\delta_2(i)：
\delta_2(i) = \max_{1 \leq j \leq 3} [\delta_1(j)a_{ji}]b_i(o_2)
同时记录每个状态i的最佳前驱状态j：
\psi_2(i) = \arg\max_{1 \leq j \leq 3} [\delta_1(j)a_{ji}], \quad i = 1, 2, 3
计算：
\delta_2(1)=\underset{1 \leq j \leq 3}{\operatorname{max}}[\delta_1(j)a_{j1}]b_1(o_2)=\underset{j}{\operatorname{max}}\{0.10\times0.5,0.16\times0.3,0.28\times0.2\}\times0.5=0.028 \\ \psi_2(1)=3 \\ \delta_2(2)=0.0504，\psi_2(2)=3 \\ \delta_2(3)=0.042，\psi_2(3)=3
同样，在t=3时
\delta_3(i)=\underset{1 \leq j \leq 3}{\operatorname{max}}[\delta_2(j)a_{ji}]b_i(o_3) \\ \psi_3(i)=\underset{1\leq j \leq 3}{\operatorname{argmax}}[\delta_2(j)a_{ji}] \\ \delta_3(1) = 0.00756，\psi_3(1)=2 \\ \delta_3(2) = 0.01008，\psi_3(2)=2 \\ \delta_3(3) = 0.0147，\psi_3(3)=3
③以P*表示最优路径的概率：
P^*=\underset{1 \leq i \leq 3}{\operatorname{max}}\delta_3(i)=0.0147
最优路径的终点是：
i_3^* = \underset{i}{\operatorname{argmax}}[\delta_3(i)]=3
④由最优路径的终点i_3^*，逆向找到i_2^*,i_1^*
在t=2时，i_2^*=\psi_3(i_3^*)=\psi_3(3)=3 \\ 在t=1时，i_1^*=\psi_2(i_2^*)=\psi_2(3)=3 \\
于是求得最优路径，即最优状态序列：
I^*=(i_1^*,i_2^*,i_3^*)=(3,3,3)
隐马尔科夫模型Python实现
import numpy as np
class HMM(object):
def __init__(self, N, M, pi=None, A=None, B=None):
self.N = N
self.M = M
self.pi = pi
self.A = A
self.B = B
def get_data_with_distribute(self, dist): # 根据给定的概率分布随机返回数据（索引）
r = np.random.random()
for i, p in enumerate(dist):
if r < p: return i
r -= p
def generate(self, T: int):
'''
根据给定的参数生成观测序列
T: 指定要生成观测序列的长度
'''
result = []
i = 0
for ind in range(T): # 依次生成状态和观测数据
if ind==0:
i = self.get_data_with_distribute(self.pi)
else:
i = self.get_data_with_distribute(self.A[i,:])
o = self.get_data_with_distribute(self.B[i,:])
result.append(o)
return result
def forward(self,X):
'''
根据给定的参数计算条件概率(前向算法)
X：观测数据
'''
alpha = self.pi * self.B[:,X[0]]
for x in X[1:]:
#将alpha构造成(N,1)维度，使用广播机制
alpha = np.sum(self.A * alpha[...,np.newaxis],axis=0) * self.B[:,x]
return alpha.sum()
def backward(self,X):
'''
根据给定的参数计算条件概率(后向算法)
X：观测数据
'''
beta = np.ones(self.N)
for x in X[:0:-1]:
#beta默认为(1,N)维度，使用广播机制
beta = np.sum(self.A * self.B[:,x] * beta,axis=1)
return np.sum(beta * self.pi * self.B[:,X[0]],axis=0)
def get_alpha_beta(self,X):
'''
根据给定数据与参数，计算所有时刻的前向概率和后向概率
'''
T = len(X)
alpha = np.zeros((T,self.N))
alpha[0,:] = self.pi * self.B[:,X[0]]
for t in range(T-1):
alpha[t+1,:] = np.sum(self.A * alpha[t,:][...,np.newaxis],axis=0) * self.B[:,X[t+1]]
beta = np.ones((T,self.N))
for t in range(T-1,0,-1):
beta[t-1,:] = np.sum(self.A * beta[t,:] * self.B[:,X[t]],axis=1)
# print(alpha[T-1,:].sum(),np.sum(beta[0,:] * self.pi * self.B[:,X[0]],axis=0))
return alpha,beta
def fit(self,X):
'''
根据观测序列估计参数
'''
#初始化参数pi，A，B
self.pi = np.random.random(self.N)
self.pi /= np.sum(self.pi)
self.A = np.random.random((self.N,self.N))
self.A /= np.sum(self.A,axis=1)[:,np.newaxis]
self.B = np.random.random((self.N,self.M))
self.B /= np.sum(self.B,axis=1)[:,np.newaxis]
T = len(X)
for epoch in range(10):
#根据公式计算下一时刻参数值
alpha,beta = self.get_alpha_beta(X)
gamma = alpha * beta
for i in range(self.N):
for j in range(self.N):
self.A[i,j] = np.sum(alpha[:-1,i] * beta[1:,j] * self.A[i,j] * self.B[j,X[1:]],axis=0) / np.sum(gamma[:-1,i],axis=0)
for i in range(self.N):
for j in range(self.M):
self.B[i,j] = np.sum(gamma[:,i] * (X == j),axis=0) / np.sum(gamma[:,i],axis=0)
self.pi = gamma[0] / gamma[0].sum()
def decode(self,X):
'''
解码问题，参数已知，根据给定的观测序列，返回最大概率的隐藏序列
'''
T = len(X)
x = X[0]
delta = self.pi[:,np.newaxis] * self.B[:,x]
varphi = np.zeros((T,self.N),dtype=np.int)
path = [0] * T
for t in range(1,T):
# delta = delta.reshape(-1,1)
tmp = delta * self.A
varphi[t,:] = np.argmax(tmp,axis=0)
delta = np.max(tmp,axis=0) * self.B[:,X[t]]
path[-1] = np.argmax(delta)
#回溯最优路径
for t in range(T-1,0,-1):
path[t-1] = varphi[t,path[t]]
return path
if __name__ == "__main__":
pi = np.array([0.25, 0.25, 0.25, 0.25])
A = np.array([
[0, 1, 0, 0],
[0.4, 0, 0.6, 0],
[0, 0.4, 0, 0.6],
[0, 0, 0.5, 0.5]])
B = np.array([
[0.5, 0.5],
[0.3, 0.7],
[0.6, 0.4],
[0.8, 0.2]])
hmm = HMM(4, 2, pi, A, B)
print(hmm.generate(10)) # 生成10个数据
X = [0,0,1,1,0] #观测序列为红，红，白，白，红
print(hmm.forward(X))
print(hmm.backward(X))
import matplotlib.pyplot as plt
def triangle_data(T): # 生成三角波形状的序列
data = []
for x in range(T):
x = x % 2
data.append(x if x <= 1 else 2-x)
return data
data = np.array(triangle_data(15))
hmm.fit(data)
gen_obs = hmm.generate(15) # 再根据学习的参数生成数据
x = np.arange(15)
plt.scatter(x, gen_obs, marker='*', color='r')
plt.plot(x, data, color='g')
plt.show()
pi = np.array([0.2, 0.4, 0.4])
A = np.array([
[0.5, 0.2, 0.3],
[0.3, 0.5, 0.2],
[0.2, 0.3, 0.5]])
B = np.array([
[0.5, 0.5],
[0.4, 0.6],
[0.7, 0.3]])
O = [0,1,0]
hmm = HMM(3, 2, pi, A, B)
path = hmm.decode(O)
print(path) #输出(2,2,2)为索引总结
隐马尔科夫模型是一种基于序列数据的概率统计模型。其参数由初始状态概率向量π、状态转移概率矩阵A以及观测概率矩阵B共同决定。其形式可表示为λ=(A, B, π)。隐马尔科夫模型作为一种生成式语言模型，在分词任务中表现出色，并广泛应用于自然语言处理以及语音识别等多个领域均有较为广泛的运用。该方法作为机器学习领域中的经典算法之一