隐马尔科夫模型（HMM）笔记（公式+代码）

Michael阿明

发布于 2020-07-13 16:28:57

5.9K0

文章被收录于专栏：Michael阿明学习之路Michael阿明学习之路

隐马尔科夫模型（hidden Markov model，HMM）是可用于标注问题的统计学习模型，描述由隐藏的马尔可夫链随机生成观测序列的过程，属于生成模型。隐马尔可夫模型在语音识别、自然语言处理、生物信息、模式识别等领域有着广泛的应用。

本文内容部分引用于李航《统计学习方法》

1. 基本概念

1.1 HMM模型定义

隐马尔可夫模型是关于时序的概率模型，描述由一个隐藏的马尔可夫链随机生成不可观测的状态随机序列，再由各个状态生成一个观测从而产生观测随机序列的过程。
隐藏的马尔可夫链随机生成的状态的序列，称为状态序列（state sequence）；
每个状态生成一个观测，而由此产生的观测的随机序列，称为观测序列（observation sequence）。
序列的每一个位置又可以看作是一个时刻。

一个不一定恰当的例子：你只知道一个人今天是散步了，逛街了，还是在家打扫卫生，推测今天的天气；这个人在下雨天可能在家打扫卫生，也有可能雨中散步，晴天可能去逛街，也有可能散步，或者打扫卫生隐藏状态 Y（不可见）：下雨，晴天观察状态 X（可见）：散步，逛街，打扫房间

假设我们观测到的状态是

序列，

X=(x_1,x_2,x_3,...x_n)

, 隐藏状态是

序列，

Y=(y_1,y_2,y_3,...y_n)

我们在知道观测序列

的情况下，求有多大概率是该隐藏状态序列

贝叶斯公式

P(B|A)=\frac {P(AB)}{P(A)}=\frac{P(A|B)P(B)}{P(A)}

现在我们要求的就是

P(Y|X)

\begin{aligned} P(Y|X) & = P(y_1,y_2,y_3,\cdots,y_n|x_1,x_2,x_3,\cdots,x_n)\\ & =\frac{P(x_1,x_2,x_3,\cdots,x_n|y_1,y_2,y_3,\cdots,y_n)P(y_1,y_2,y_3,\cdots,y_n)}{P(x_1,x_2,x_3,\cdots,x_n)}\quad\quad(0)\\ P(y_1,y_2,y_3,\cdots,y_n) & =P(y_1){\color{red}P(y_2,y_3,y_4,\cdots,y_n|y_1)}\\ & =P(y_1)P(y_2|y_1){\color{red}P(y_3,y_4,\cdots,y_n|y_1,y_2)}\\ & =P(y_1)P(y_2|y_1)P(y_3|y_1,y_2){\color{red}P(y_4,\cdots,y_n|y_1,y_2,y_3)}\\ \vdots\\ & =P(y_1)P(y_2|y_1)P(y_3|y_1,y_2)\cdots P(y_{n-1}|y_1,y_2,\cdots,y_{n-2})P(y_n|y_1,y_2,\cdots,y_{n-1})\\ & =P(y_1)\prod_{i=2}^{n}P(y_i|y_1,y_2,y_3,\cdots,y_{i-1})\\ & =\color{blue}{P(y_1)\prod_{i=2}^{n}P(y_i|y_{i-1})}\quad\quad\quad(1)\\ \end{aligned}

隐马尔可夫模型 两个假设

上式最后一步是隐马尔可夫的 齐次性 假设：当前状态

y_i

仅依赖于前一个状态

y_{i-1}

下面式子最后一步是隐马尔可夫的 观测独立性 假设：观测值

x_i

只跟他的隐含状态

y_i

1.2 盒子和球模型

假设有4个盒子，每个盒子里都装有红、白两种颜色的球，盒子里的红、白球数由表列出。

按下面的方法抽球，产生一个球的颜色的观测序列

：

等概率随机选1个盒子，从这个盒子里随机抽出1个球，记录颜色，放回；
从当前盒子随机转移到下一个盒子，规则是：如果当前盒子是盒子1，那么下一盒子一定是盒子2；如果当前是盒子2或3，那么分别以概率0.4和0.6转移到左边或右边的盒子；如果当前是盒子4，那么各以0.5的概率停留在盒子4或转移到盒子3；
确定转移的盒子后，再从这个盒子里随机抽出1个球，记录颜色，放回；
重复5次，得到一个球的颜色的观测序列

各个隐含状态的初始概率

P(y_i)

：

\pi = (0.25,0.25,0.25,0.25)^T

各个隐含状态间的转移概率

P(y_i|y_{i-1})

：

A = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0.4 & 0 & 0.6 & 0 \\ 0 & 0.4 & 0 & 0.6 \\ 0 & 0 & 0.5 & 0.5 \\ \end{bmatrix}

观测(发射)概率

P(x_i|y_i)

：

B = \begin{bmatrix} 0.5 & 0.5 \\ 0.3 & 0.7 \\ 0.6 & 0.4 \\ 0.8 & 0.2 \\ \end{bmatrix}

1.3 观测序列生成过程

输入：隐马模型

\lambda=（A，B，\pi）

，观测序列长度

(该例子为 5)

输出：观测序列

X = (x_1,x_2,x_3,x_4,x_5)

, (红、白球)

跟据

\pi

随机产生一个

y_i

状态(盒子), 序列长度 t = 0

y_i

中，按照其观测概率

对应的行，产生一个观测序列

的元素

x_i

（红球 or 白球），计数 t++

根据

y_i

的状态转移概率

对应的行，产生下一个状态

y_{i+1}

, while(t < T), 重复2,3步骤

1.4 HMM模型3个基本问题

概率计算问题：给定模型

\lambda=（A，B，\pi）

和观测序列

X = (x_1,x_2.....,x_n)

, 计算在模型

\lambda

下，观测序列

出现的概率

P(X|\lambda)

学习问题：已知观测序列

X = (x_1,x_2.....,x_n)

，估计模型

\lambda=（A，B，\pi）

的参数，使得在该模型下，观测序列概率

P(X|\lambda)

最大。极大似然估计的方法估计参数

预测问题：也称解码(decoding)问题。已知模型

\lambda=（A，B，\pi）

和观测序列

X = (x_1,x_2.....,x_n)

，求对给定观测序列条件概率

P(Y|X)

最大的状态序列

Y = (y_1,y_2......,y_n)

。即给定观测序列

，求最有可能的对应隐含状态序列

2. 概率计算问题

给定模型

\lambda=（A，B，\pi）

和观测序列

X = (x_1,x_2.....,x_n)

, 计算在模型

\lambda

下，观测序列

出现的概率

P(X|\lambda)

2.1 直接计算法

列举所有的长度为

的状态序列

Y = (y_1,y_2......,y_t)

求各个状态序列

与观测序列

X = (x_1,x_2.....,x_t)

的联合概率

P(XY|\lambda)

然后对上面所有的状态序列求和

\Sigma

，得到

P(X|\lambda)

状态序列

Y = (y_1,y_2......,y_t)

的概率是：

P(Y|\lambda)=\pi_{y_1}A_{y_1,y_2}A_{y_2,y_3}......A_{y_{t-1},y_t}

对固定的状态序列

Y = (y_1,y_2......,y_t)

，观测序列

X = (x_1,x_2.....,x_t)

的概率是：

P(X|Y,\lambda)=B_{y_1,x_1}B_{y_2,x_2}......B_{y_t,x_t}

X,Y

同时出现的联合概率为：

P(XY|\lambda)=P(X|Y,\lambda)P(Y|\lambda)=\pi_{y_1}B_{y_1,x_1}A_{y_1,y_2}B_{y_2,x_2}......A_{y_{t-1},y_t}B_{y_t,x_t}

对上式在所有可能的状态序列

Y_i

的情况下求和

\Sigma

，即可得到

P(X|\lambda)=\sum_{Y_1}^{Y_n} P(X|Y,\lambda)P(Y|\lambda)

但是上面计算量很大，复杂度

O(TN^T)

，不可行

2.2 前向算法

前向概率 概念：给定HMM模型

\lambda

，定义到时刻 t 部分观测序列

X_{part} = (x_1,x_2.....,x_t)

且 t 时刻的状态

y_t

为

q_i

的概率为前向概率，记为：

a_t(i)=P(x_1,x_2,...,x_t,y_t=q_i | \lambda)

递推求解前向概率

a_t(i)

及观测序列概率

P(X|\lambda)

输入：给定HMM模型

\lambda

, 观测序列

输出：观测序列概率

P(X|\lambda)

初值：

a_1(i) = \pi_iB_{i,x_1}\quad\quad i = 1,2,...,N

，

为盒子个数

递推：对

t=1,2,...,T-1

，

a_{t+1}(i)= \begin{bmatrix} \sum_{j=1}^N a_t(j)A_{j,i} \end{bmatrix}B_{i,x_{t+1}}\quad i = 1,2,...,N

终止：

P(X|\lambda) = \sum_{i=1}^N a_T(i)

（看前向概率定义，全概率公式）

算法解释：

初始时刻

t=1

，观测到的是

x_1

, 可能的状态是盒子的编号

i = 1,2,...,N

，概率为

\pi_i

，在

盒子发射出

x_1

颜色球的概率为

B_{i,x_1}

，所以

a_1(i) = \pi_iB_{i,x_1}

递推：上一时刻

的 N 种情况下都可能转移到

t+1

时的

q_i

状态，对应的前向概率乘以转移概率，并求和，得到状态

q_i

的概率，再乘以发射概率

B_{i,x_{t+1}}

，就是

t+1

时的前向概率

最后一个时刻

时的所有前向概率求和就是

P(X|\lambda) = \sum_{i=1}^N a_T(i)

前向算法是基于路径的，

t+1

时刻，直接用

时刻的结果，时间复杂度

O(TN^2)

2.2.1 前向公式证明

首先有公式联合概率

P(ABC) = P(A)P(B|A)P(C|AB)

,对任意多个项都成立递推公式证明：

\begin{aligned} a_t(j)A_{j,i} &= P(x_1,x_2,...,x_t,y_t=q_j|\lambda)P(y_{t+1}=q_i|y_t=q_j,\lambda)\\ &={\color{red}P(x_1,x_2,...,x_t|y_t=q_j,\lambda)}P(y_t=q_j|\lambda)P(y_{t+1}=q_i|y_t=q_j,\lambda)\\ &= {\color{red}P(x_1,x_2,...,x_t|y_t=q_j,{\color{blue}y_{t+1}=q_i},\lambda)}P(y_t=q_j,y_{t+1}=q_i|\lambda)\\ &=P(x_1,x_2,...,x_t,y_t=q_j,y_{t+1}=q_i|\lambda)\\ \{\sum_{j=1}^N a_t(j)A_{j,i}\} B_{i,x_{t+1}} &= P(x_1,x_2,...,x_t,y_{t+1}=q_i|\lambda)P(x_{t+1}|y_{t+1}=q_i,\lambda)\\ &={\color{red}P(x_1,x_2,...,x_t|y_{t+1}=q_i,\lambda)}P(y_{t+1}=q_i|\lambda)P(x_{t+1}|y_{t+1}=q_i,\lambda)\\ &={\color{red}P(x_1,x_2,...,x_t|y_{t+1}=q_i,{\color{blue}x_{t+1}},\lambda)}P(x_{t+1},y_{t+1}=q_i|\lambda)\\ &=P(x_1,x_2,...,x_t,x_{t+1},y_{t+1}=q_i|\lambda)=a_{t+1}(i) \end{aligned}

第一个蓝色处：前

个观测序列，显然跟

t+1

时刻的状态

y_{t+1}

无关，第二个蓝色处：观测独立性

终止公式证明(全概率公式)：

\begin{aligned} \sum_{i=1}^N a_T(i)&=P(x_1,x_2,...,x_T,y_T=q_1|\lambda)+P(x_1,x_2,...,x_T,y_T=q_2|\lambda)+P(x_1,x_2,...,x_T,y_T=q_N|\lambda)\\ &= P(x_1,x_2,...,x_T|\lambda)=P(X|\lambda) \end{aligned}

2.2.2 盒子和球例子

考虑盒子和球模型

\lambda = (A,B,\pi)

，状态集合（盒子的编号）

Q=\{1,2,3\}

，观测集合

V=\{红，白\}

\pi=\begin{bmatrix} 0.2 \\ 0.4 \\ 0.4 \\ \end{bmatrix}, A = \begin{bmatrix} 0.5 & 0.2 & 0.3 \\ 0.3 & 0.5 & 0.2 \\ 0.2 & 0.3 & 0.5 \\ \end{bmatrix}, B = \begin{bmatrix} 0.5 & 0.5 \\ 0.4 & 0.6 \\ 0.7 & 0.3 \\ \end{bmatrix}

设总的时间长度

T=3

, 观测序列

X=(红，白，红)

，用前向算法计算

P(X|\lambda)

解：

计算前向概率初值：

a_1(i) = \pi_iB_{i,x_1}

a_1(1) = \pi_1B_{1,1}=0.2*0.5=0.10

，（1号盒子，时刻1摸出红色(第一列)的概率）

a_1(2) = \pi_2B_{2,1}=0.4*0.4=0.16

，（2号盒子，时刻1摸出红色的概率）

a_1(3) = \pi_3B_{3,1}=0.4*0.7=0.28

，（3号盒子，时刻1摸出红色的概率）

递推计算：对

t=1,2,...,T-1

，

a_{t+1}(i)= \begin{bmatrix} \sum_{j=1}^N a_t(j)A_{j,i} \end{bmatrix}B_{i,x_{t+1}}\quad i = 1,2,...,N

a_2(1)=\{ \sum_{j=1}^3 a_1(j)A_{j,1} \} B_{1,2}=(0.10*0.5+0.16*0.3+0.28*0.2)*0.5=0.077

（时刻2时，从1号盒子摸出白色的概率）

a_2(2)=\{ \sum_{j=1}^3 a_1(j)A_{j,2} \} B_{2,2}=(0.10*0.2+0.16*0.5+0.28*0.3)*0.6=0.1104

（时刻2时，从2号盒子摸出白色的概率）

a_2(3)=\{ \sum_{j=1}^3 a_1(j)A_{j,3} \} B_{3,2}=(0.10*0.3+0.16*0.2+0.28*0.5)*0.3=0.0606

（时刻2时，从3号盒子摸出白色的概率）

a_3(1)=\{ \sum_{j=1}^3 a_2(j)A_{j,1} \} B_{1,1}=(0.077*0.5+0.1104*0.3+0.0606*0.2)*0.5=0.04187

（时刻3时，从1号盒子摸出红色的概率）

a_3(2)=\{ \sum_{j=1}^3 a_2(j)A_{j,2} \} B_{2,1}=(0.077*0.2+0.1104*0.5+0.0606*0.3)*0.4=0.035512

（时刻3时，从2号盒子摸出红色的概率）

a_3(3)=\{ \sum_{j=1}^3 a_2(j)A_{j,3} \} B_{3,1}=(0.077*0.3+0.1104*0.2+0.0606*0.5)*0.7=0.052836

（时刻3时，从3号盒子摸出红色的概率）

终止：

P(X|\lambda) = \sum_{i=1}^N a_T(i) = \sum_{i=1}^3 a_3(i)=0.04187+0.035512+0.052836=0.130218‬

2.2.3 前向算法Python代码

借鉴他人代码，手敲一遍并理解

# -*- coding:utf-8 -*-
# python3.7
# @Time: 2019/12/17 22:33
# @Author: Michael Ming
# @Website: https://michael.blog.csdn.net/
# @File: hiddenMarkov.py
import numpy as np
class HiddenMarkov:
    def forward(self, Q, V, A, B, X, PI):   # 前向算法
        N = len(Q)  # 隐藏状态数量
        M = len(X)  # 观测序列大小
        alphas = np.zeros((N, M))   # 前向概率 alphas 矩阵
        T = M   # 时刻长度，即观测序列长度
        for t in range(T):  # 遍历每个时刻，计算前向概率 alphas
            indexOfXi = V.index(X[t])   # 观测值Xi对应的V中的索引
            for i in range(N):  # 对每个状态进行遍历
                if t == 0: # 计算初值
                    alphas[i][t] = PI[t][i] * B[i][indexOfXi]   # a1(i)=πi B(i,x1)
                    print("alphas1(%d)=p%dB%d(x1)=%f" %(i,i,i,alphas[i][t]))
                else:
                    alphas[i][t] = np.dot(
                        [alpha[t-1] for alpha in alphas],   #  取的alphas的t-1列
                        [a[i] for a in A]) *B[i][indexOfXi]    # 递推公式
                    print("alpha%d(%d)=[sigma alpha%d(i)ai%d]B%d(x%d)=%f"
                        %(t, i, t-1, i, i, t, alphas[i][t]))
                    print(alphas)
        P = sum(alpha[M-1] for alpha in alphas)  # 最后一列概率的和
        print("P=%f" % P)

Q = [1, 2, 3]
V = ['红', '白']
A = [[0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]]
B = [[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]]
# X = ['红', '白', '红', '红', '白', '红', '白', '白']
X = ['红', '白', '红']    # 书上的例子
PI = [[0.2, 0.4, 0.4]]

hmm = HiddenMarkov()
hmm.forward(Q, V, A, B, X, PI)

关于numpy的一些操作：

>>> a
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])
>>> b
array([[ 7,  8],
       [ 9, 10],
       [11, 12]])
>>> np.dot(a[0],[B[1] for B in b])
64
>>> np.dot([A[0] for A in a],[B[1] for B in b])
132
>>> np.multiply(a[0],[B[1] for B in b])
array([ 8, 20, 36])
>>> np.multiply([A[0] for A in a],[B[1] for B in b])
array([ 8, 40, 84])

运行结果：（跟上面计算一致）

alphas1(0)=p0B0(x1)=0.100000
alphas1(1)=p1B1(x1)=0.160000
alphas1(2)=p2B2(x1)=0.280000
alpha1(0)=[sigma alpha0(i)ai0]B0(x1)=0.077000
[[0.1   0.077 0.   ]
 [0.16  0.    0.   ]
 [0.28  0.    0.   ]]
alpha1(1)=[sigma alpha0(i)ai1]B1(x1)=0.110400
[[0.1    0.077  0.    ]
 [0.16   0.1104 0.    ]
 [0.28   0.     0.    ]]
alpha1(2)=[sigma alpha0(i)ai2]B2(x1)=0.060600
[[0.1    0.077  0.    ]
 [0.16   0.1104 0.    ]
 [0.28   0.0606 0.    ]]
alpha2(0)=[sigma alpha1(i)ai0]B0(x2)=0.041870
[[0.1     0.077   0.04187]
 [0.16    0.1104  0.     ]
 [0.28    0.0606  0.     ]]
alpha2(1)=[sigma alpha1(i)ai1]B1(x2)=0.035512
[[0.1      0.077    0.04187 ]
 [0.16     0.1104   0.035512]
 [0.28     0.0606   0.      ]]
alpha2(2)=[sigma alpha1(i)ai2]B2(x2)=0.052836
[[0.1      0.077    0.04187 ]
 [0.16     0.1104   0.035512]
 [0.28     0.0606   0.052836]]
P=0.130218

2.3 后向算法

后向概率 概念：给定HMM模型

\lambda

，定义到时刻

状态

y_t

为

q_i

的条件下，从

t+1

到

的部分观测序列

X_{part} = (x_{t+1},x_{t+2}.....,x_T)

的概率为后向概率，记为：

\beta_t(i)=P(x_{t+1},x_{t+2}.....,x_T|y_t=q_i , \lambda)

递推求解后向概率

\beta_t(i)

及观测序列概率

P(X|\lambda)

输入：给定HMM模型

\lambda

, 观测序列

输出：观测序列概率

P(X|\lambda)

初值：

\beta_T(i) = 1\quad\quad i = 1,2,...,N

，

为盒子个数

递推：对

t=T-1,T-2,...,1

，

\beta_{t}(i)=\sum_{j=1}^N A_{i,j}B_{j,x_{t+1}}\beta_{t+1}(j)\quad i = 1,2,...,N

终止：

P(X|\lambda) = \sum_{i=1}^N \pi_iB_{i,x_1}\beta_1(i)

算法解释：

2.3.1 后向公式证明

递推公式证明：

\begin{aligned} \sum_{j=1}^N A_{i,j}B_{j,x_{t+1}}\beta_{t+1}(j) &= \sum_{j=1}^N A_{i,j}{\color{red}P(x_{t+1}|y_{t+1}=q_j,\lambda)}P(x_{t+2},...,x_T|y_{t+1}=q_j,\lambda)\\ &= \sum_{j=1}^N A_{i,j}{\color{red}P(x_{t+1}|{\color{blue}x_{t+2},...,x_T},y_{t+1}=q_j,\lambda)}P(x_{t+2},...,x_T|y_{t+1}=q_j,\lambda)\\ &= \sum_{j=1}^N A_{i,j}P(x_{t+1},x_{t+2},...,x_T|{\color{red}y_{t+1}=q_j},\lambda)\\ &= \sum_{j=1}^N A_{i,j}P(x_{t+1},x_{t+2},...,x_T|{\color{red}y_{t+1}=q_j},{\color{blue}y_t=q_i},\lambda)\\ &= \sum_{j=1}^N P(y_{t+1}=q_j|y_t=q_i,\lambda)P(x_{t+1},x_{t+2},...,x_T|{\color{red}y_{t+1}=q_j},{\color{blue}y_t=q_i},\lambda)\\ &= \sum_{j=1}^N P(x_{t+1},x_{t+2},...,x_T,{\color{red}y_{t+1}=q_j}|{\color{blue}y_t=q_i},\lambda)\\ &= P(x_{t+1},x_{t+2},...,x_T|{\color{blue}y_t=q_i},\lambda)=\beta_{t}(i)\\ \end{aligned}\\

第一个蓝色处：观测独立性；第二个蓝色处：观测独立性（

x_{t+1},...x_T

都与

y_t

无关）

终止公式证明：

\begin{aligned} \pi_iB_{i,x_1}\beta_1(i)&=P(y_1=q_i|\lambda)P(x_1|y_1=q_i,\lambda)\color{red}P(x_2,x_3,...,x_T|y_1=q_i,\lambda)\\ &=P(x_1,y_1=q_i|\lambda)\color{red}P(x_2,x_3,...,x_T|{\color{blue}{x_1}},y_1=q_i,\lambda)\\ &=P(x_1,x_2,...,x_T,y_1=q_i|\lambda)\\ \sum_{i=1}^N \pi_iB_{i,x_1}\beta_1(i) &= \sum_{i=1}^N P(x_1,x_2,...,x_T,y_1=q_i|\lambda)=P(x_1,x_2,...,x_T|\lambda)=P(X|\lambda)\\ \end{aligned}

利用前向和后向概率的定义，可以将观测序列概率

P(X|\lambda)

写成：

P(X|\lambda)=\sum_{i=1}^N\sum_{j=1}^N \alpha_t(i)A_{i,j}B_{j,x_{t+1}}\beta_{t+1}(j),\quad t= 1,2,...,T-1

证明如下：

\begin{aligned} \sum_{i=1}^N\sum_{j=1}^N \alpha_t(i)A_{i,j}B_{j,x_{t+1}}\beta_{t+1}(j) &= \sum_{i=1}^N\sum_{j=1}^N P(x_1,x_2,...,x_t,y_t=q_i|\lambda)P(y_{t+1}=q_j|y_t=q_i,\lambda)\\&\quad\quad\quad\quad\quad*P(x_{t+1}|y_{t+1}=q_j,\lambda)P(x_{t+2},x_{t+3},...,x_T|y_{t+1}=q_j,\lambda)\\ &= \sum_{i=1}^N\sum_{j=1}^N P(x_1,x_2,...,x_t|y_t=q_i,\lambda)P(y_t=q_i|\lambda)P(y_{t+1}=q_j|y_t=q_i,\lambda)\\&\quad\quad\quad\quad\quad*P(x_{t+1}|{\color{red}x_{t+2},x_{t+3},...,x_{T}},y_{t+1}=q_j,\lambda)P(x_{t+2},x_{t+3},...,x_T|y_{t+1}=q_j,\lambda)\\ &= \sum_{i=1}^N\sum_{j=1}^N P(x_1,x_2,...,x_t|y_t=q_i,{\color{blue}y_{t+1}=q_j},\lambda)P(y_t=q_i|\lambda)P(y_{t+1}=q_j|y_t=q_i,\lambda)\\&\quad\quad\quad\quad\quad*P(x_{t+1}|{x_{t+2},x_{t+3},...,x_{T}},y_{t+1}=q_j,\lambda)P(x_{t+2},x_{t+3},...,x_T|y_{t+1}=q_j,\lambda)\\ &= \sum_{i=1}^N\sum_{j=1}^N P(x_1,x_2,...,x_t|y_t=q_i,{\color{blue}y_{t+1}=q_j},\lambda)P(y_t=q_i|\lambda)P(y_{t+1}=q_j|y_t=q_i,\lambda)\\&\quad\quad\quad\quad\quad*P(x_{t+1},{x_{t+2},x_{t+3},...,x_{T}}|y_{t+1}=q_j,\lambda)\\ &= \sum_{i=1}^N\sum_{j=1}^N P(x_1,x_2,...,x_t,{\color{blue}y_{t+1}=q_j}|y_t=q_i,\lambda)P(y_t=q_i|\lambda)\\&\quad\quad\quad\quad\quad*P(x_{t+1},{x_{t+2},x_{t+3},...,x_{T}}|{\color{red}x_1,x_2,...,x_t},y_{t+1}=q_j,{\color{blue}y_t=q_i},\lambda)\\ &= \sum_{i=1}^N\sum_{j=1}^N P(y_t=q_i|\lambda)P({\color{orangered}x_1,x_2,...,x_t,y_{t+1}=q_j}|y_t=q_i,\lambda)\\&\quad\quad\quad\quad\quad*P(x_{t+1},{x_{t+2},x_{t+3},...,x_{T}}|{\color{orangered}x_1,x_2,...,x_t,y_{t+1}=q_j},{\color{blue}y_t=q_i},\lambda)\\ &= \sum_{i=1}^N\sum_{j=1}^N P(y_t=q_i|\lambda)P({\color{orangered}x_1,x_2,...,x_t,y_{t+1}=q_j},x_{t+1},{x_{t+2},x_{t+3},...,x_{T}}|{\color{blue}y_t=q_i},\lambda)\\ &= \sum_{i=1}^N P(y_t=q_i|\lambda)P({\color{orangered}x_1,x_2,...,x_t},x_{t+1},{x_{t+2},x_{t+3},...,x_{T}}|{\color{blue}y_t=q_i},\lambda)\\ &= \sum_{i=1}^N P({\color{orangered}x_1,x_2,...,x_t},x_{t+1},{x_{t+2},x_{t+3},...,x_{T}},{\color{blue}y_t=q_i}|\lambda)\\ &= P(x_1,x_2,...,x_T|\lambda)=P(X|\lambda)\\ \end{aligned}

2.3.2 后向算法Python代码

import numpy as np
class HiddenMarkov:
    def backward(self, Q, V, A, B, X, PI):  # 后向算法
        N = len(Q)  # 隐藏状态数量
        M = len(X)  # 观测序列大小
        betas = np.ones((N, M))    # 后向概率矩阵 betas
        for i in range(N):
            print("beta%d(%d)=1" %(M, i))   # 后向概率初始为1
        for t in range(M-2, -1, -1):
            indexOfXi = V.index(X[t+1])     # 观测值Xi对应的V中的索引
            for i in range(N):
                betas[i][t] = np.dot(
                    np.multiply(A[i], [b[indexOfXi] for b in B]),   # A的i行 B的Xi列 = 行
                    [beta[t+1] for beta in betas])  # 递推公式
                realT = t+1 # 真实的时间，t从1开始
                realI = i+1 # 状态标号，从1开始
                print("beta%d(%d)=[sigma A%djBj(x%d)beta%d(j)]=("
                      %(realT, realI, realI, realT+1, realT+1), end='') # end,表示以该符号空开，不换行
                for j in range(N):
                    print("%.2f*%.2f*%.2f+" %(A[i][j], B[j][indexOfXi],betas[j][t+1]), end=' ')
                print("0)=%.3f" % betas[i][t])
        print(betas)
        indexOfXi = V.index(X[0])
        P = np.dot(
            np.multiply(PI, [b[indexOfXi] for b in B]),
            [beta[0] for beta in betas])
        print("P(X|lambda)=", end=" ")
        for i in range(N):
            print("%.1f*%.1f*%.5f+" % (PI[0][i], B[i][indexOfXi], betas[i][0]), end=" ")
        print("0=%f" % P)

Q = [1, 2, 3]
V = ['红', '白']
A = [[0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]]
B = [[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]]
# X = ['红', '白', '红', '红', '白', '红', '白', '白']
X = ['红', '白', '红']    # 书上的例子
PI = [[0.2, 0.4, 0.4]]

hmm = HiddenMarkov()
hmm.backward(Q, V, A, B, X, PI)

运行结果：

beta3(0)=1
beta3(1)=1
beta3(2)=1
beta2(1)=[sigma A1jBj(x3)beta3(j)]=(0.50*0.50*1.00+ 0.20*0.40*1.00+ 0.30*0.70*1.00+ 0)=0.540
beta2(2)=[sigma A2jBj(x3)beta3(j)]=(0.30*0.50*1.00+ 0.50*0.40*1.00+ 0.20*0.70*1.00+ 0)=0.490
beta2(3)=[sigma A3jBj(x3)beta3(j)]=(0.20*0.50*1.00+ 0.30*0.40*1.00+ 0.50*0.70*1.00+ 0)=0.570
beta1(1)=[sigma A1jBj(x2)beta2(j)]=(0.50*0.50*0.54+ 0.20*0.60*0.49+ 0.30*0.30*0.57+ 0)=0.245
beta1(2)=[sigma A2jBj(x2)beta2(j)]=(0.30*0.50*0.54+ 0.50*0.60*0.49+ 0.20*0.30*0.57+ 0)=0.262
beta1(3)=[sigma A3jBj(x2)beta2(j)]=(0.20*0.50*0.54+ 0.30*0.60*0.49+ 0.50*0.30*0.57+ 0)=0.228
[[0.2451 0.54   1.    ]
 [0.2622 0.49   1.    ]
 [0.2277 0.57   1.    ]]
P(X|lambda)= 0.2*0.5*0.24510+ 0.4*0.4*0.26220+ 0.4*0.7*0.22770+ 0=0.130218
# 概率跟前向算法计算的是一致的

2.4 一些概率与期望值

结论1：

P(X,y_t=q_i|\lambda)=\alpha_t(i)\beta_t(i)

证明：

\begin{aligned} P(X,y_t=q_i|\lambda)&=P(x_1,x_2,...,x_t,x_{t+1},...,x_T|y_t=q_i,\lambda)P(y_t=q_i|\lambda)\\ &=P(x_1,x_2,...,x_t|y_t=q_i,\lambda)P(x_{t+1},...,x_T|y_t=q_i,\lambda)P(y_t=q_i|\lambda)\\ &=P(x_1,x_2,...,x_t,y_t=q_i | \lambda)P(x_{t+1},...,x_T|y_t=q_i,\lambda)\\ &= \alpha_t(i)\beta_t(i) \end{aligned}

结论2：给定模型

\lambda

和观测序列

，在时刻

处于状态

q_i

的概率：

\gamma_t(i)=P(y_t=q_i|X,\lambda)=\frac{ \alpha_t(i)\beta_t(i)}{\sum_{i=1}^N \alpha_t(i)\beta_t(i)}

证明：

\begin{aligned} \gamma_t(i)&=P(y_t=q_i|X,\lambda)=\frac{P(y_t=q_i,X|\lambda)}{P(X|\lambda)}\\ &=\frac{ \alpha_t(i)\beta_t(i)}{\sum_{i=1}^N P(y_t=q_i,X|\lambda)}\\ &=\frac{ \alpha_t(i)\beta_t(i)}{\sum_{i=1}^N \alpha_t(i)\beta_t(i)} \end{aligned}

结论3：给定模型

\lambda

和观测序列

，在时刻

处于状态

q_i

且在时刻

t+1

处于

q_j

的概率：

\xi_t(i,j)=P(y_t=q_i,y_{t+1}=q_j |X, \lambda)=\frac{\alpha_t(i)A_{i,j}B_{j,x_{t+1}}\beta_{t+1}(j)}{\sum_{i=1}^N\sum_{j=1}^N \alpha_t(i)A_{i,j}B_{j,x_{t+1}}\beta_{t+1}(j)}

证明：由上面可知：

P(X|\lambda)=\sum_{i=1}^N\sum_{j=1}^N \alpha_t(i)A_{i,j}B_{j,x_{t+1}}\beta_{t+1}(j)

\alpha_t(i)A_{i,j}B_{j,x_{t+1}}\beta_{t+1}(j)=P(X,y_t=q_i,y_{t+1}=q_j | \lambda)

\begin{aligned} \xi_t(i,j)&=P(y_t=q_i,y_{t+1}=q_j |X, \lambda)\\ &=\frac{P(X,y_t=q_i,y_{t+1}=q_j | \lambda)}{P(X|\lambda)}\\ &=\frac{\alpha_t(i)A_{i,j}B_{j,x_{t+1}}\beta_{t+1}(j)}{\sum_{i=1}^N\sum_{j=1}^N \alpha_t(i)A_{i,j}B_{j,x_{t+1}}\beta_{t+1}(j)} \end{aligned}

结论4：

在观测序列

下状态

q_i

出现的期望值：

\sum_{t=1}^T \gamma_t(i)

在观测序列

下由状态

q_i

转移的期望值：

\sum_{t=1}^{T-1} \gamma_t(i)

在观测序列

下由状态

q_i

转移到状态

q_j

的期望值：

\sum_{t=1}^{T-1} \xi_t(i,j)

3. 学习算法

隐马尔可夫模型的学习，根据训练数据进行分类：

数据	学习方法
观测序列 XXX + 对应的状态序列 YYY	监督学习
观测序列 XXX	无监督学习（BaumBaumBaum-WelchWelchWelch / EMEMEM 算法）

+ 对应的状态序列

监督学习观测序列

无监督学习（

Baum

Welch

算法）

3.1 监督学习方法

转移概率

\quad i->j

：

到

状态的频数

转移到所有状态的频数和

观察概率

\quad i

观测

：

状态观测到

的频数

状态观测到所有

x_i

的频数和

初始状态概率

\pi

\quad y_i

在样本中的初始概率

由于监督学习需要使用标注的训练数据，而人工标注训练数据往往代价很高，需要用无监督学习的方法。

3.2 无监督

Baum

Welch

算法

(先跳过)

4. 预测算法

4.1 近似算法

根据 2.4 节结论2：每个时刻 t 处于状态

q_i

的概率可以计算

那么每一时刻 t 的最有可能的状态就是概率最大的那个状态。

优点是计算简单
缺点是不能保证预测的状态序列整体是最有可能的状态序列，因为预测的状态序列可能有实际不发生的部分。事实上，上述方法得到的状态序列中有可能存在转移概率为0的相邻状态，尽管如此，近似算法仍然是有用的。

4.2 维特比

Viterbi

算法

维特比算法实际是用动态规划（dynamic programming）解隐马尔可夫模型预测问题，即用动态规划求概率最大路径（最优路径）。这时一条路径对应着一个状态序列。

2个定义：

在时刻 t 状态为 i 的所有单个路径中概率最大值为

\delta_t(i) = \max \limits_{y_1,y_2,...,y_{t-1}}P(y_t=i,y_{t-1},...,y_1,x_t,...,x_1|\lambda)

由定义可得变量

\delta

的递推公式：

\begin{aligned} \delta_{t+1}(i) &= \max \limits_{y_1,y_2,...,y_{t}}P(y_{t+1}=i,y_{t},...,y_1,x_{t+1},...,x_1|\lambda)\\ &= \max \limits_{1 \leq j \leq N }[\delta_t(j)A_{j,i}]B_i(x_{t+1}) \end{aligned}

在时刻 t 状态为 i 的所有单个路径中概率最大的路径的第 t-1 个结点为

\Psi_t(i)=arg \quad \max \limits_{1 \leq j \leq N }[\delta_{t-1}(j)A_{j,i}]

viterbi 算法：

输入：模型

\lambda=(A,B,\pi)

，观测序列

X=(x_1,x_2,...,x_T)

输出：最优路径

Y^*=(y_1^*,y_2^*,...,y_T^*)

初始化

\delta_1(i)=\pi_iB_i(x_1)

；

\Psi_1(i)=0

递推，对

t=2,3,...,T

\begin{aligned} \delta_{t}(i) &= \max \limits_{1 \leq j \leq N }[\delta_{t-1}(j)A_{j,i}]B_i(x_{t})\\ \Psi_t(i)&=arg \quad \max \limits_{1 \leq j \leq N }[\delta_{t-1}(j)A_{j,i}] \end{aligned}

终止
最优路径回溯

y_t^*=\Psi_{t+1}(y_{t+1}^*)

例题：考虑盒子和球模型

\lambda = (A,B,\pi)

，状态集合（盒子的编号）

Q=\{1,2,3\}

，观测集合

V=\{红，白\}

已知观测序列

X=(红，白，红)

，求最优状态序列

Y^*

解：

（1）初始化在

t=1

时刻，对每一个状态 i ，i = 1,2,3, 求状态为 i 观测

x_1

为红的概率，记为

\delta_1(i)

\delta_1(i)=\pi_iB_i(x_1)=\pi_iB_i(红)，i=1,2,3

\delta_1(1)=0.2*0.5=0.10

\delta_1(2)=0.4*0.4=0.16

\delta_1(3)=0.4*0.7=0.28

\Psi_1(i)=0，i=1,2,3

（2）在

t=2

时刻，对每个状态 i ，i = 1,2,3, 求在

t=1

时状态为

观测为红并在

t=2

时刻状态为 i 观测为

x_2

为白的路径的最大概率，记此最大概率为

\delta_2(i)

\delta_2(i)=\max \limits_{1 \leq j \leq 3 }[\delta_{1}(j)A_{j,i}]B_i(x_{2})

同时，对每个状态

i, i=1,2,3,

记录概率最大路径的前一个状态

\Psi_2(i) = arg \max \limits_{1 \leq j \leq 3}[\delta_1(j)A_{j,i}],\quad i=1,2,3

计算：

\delta_2(1)=\max \limits_{1 \leq j \leq 3 }[\delta_{1}(j)A_{j,1}]B_1(x_{2}) = \max \limits_j \{0.10*0.5, 0.16*0.3, 0.28*0.2\}*0.5=0.056*0.5=0.028

\Psi_2(1) = 3

\delta_2(2)=\max \limits_{1 \leq j \leq 3 }[\delta_{1}(j)A_{j,2}]B_2(x_{2}) = \max \limits_j \{0.10*0.2, 0.16*0.5, 0.28*0.3\}*0.6=0.084*0.6=0.0504

\Psi_2(2) = 3

\delta_2(3)=\max \limits_{1 \leq j \leq 3 }[\delta_{1}(j)A_{j,3}]B_3(x_{2}) = \max \limits_j \{0.10*0.3, 0.16*0.2, 0.28*0.5\}*0.3=0.14*0.3=0.042

\Psi_2(3) = 3

在

t=3

时刻,

\delta_3(i)=\max \limits_{1 \leq j \leq 3 }[\delta_{2}(j)A_{j,i}]B_i(x_{3})

\Psi_3(i) = arg \max \limits_{1 \leq j \leq 3 } [\delta_2(j)A_{j,i}]

\delta_3(1)=\max \limits_{1 \leq j \leq 3 }[\delta_{2}(j)A_{j,1}]B_1(x_{3}) = \max \limits_j \{0.028*0.5, 0.0504*0.3, 0.042*0.2\}*0.5=0.01512*0.5=0.00756

\Psi_3(1) = 2

\delta_3(2)=\max \limits_{1 \leq j \leq 3 }[\delta_{2}(j)A_{j,2}]B_2(x_{3}) = \max \limits_j \{0.028*0.2, 0.0504*0.5, 0.042*0.3\}*0.4=0.0252*0.4=0.01008

\Psi_3(2) = 2

\delta_3(3)=\max \limits_{1 \leq j \leq 3 }[\delta_{2}(j)A_{j,3}]B_3(x_{3}) = \max \limits_j \{0.028*0.3, 0.0504*0.2, 0.042*0.5\}*0.7=0.021*0.7=0.0147

\Psi_3(3) = 3

（3）以

Y^*

表示最优路径的概率，则

Y^* = \max \limits_{1 \leq i \leq 3}\delta_3(i)=0.0147

最优路径的终点是

i_3^*=arg \max \limits_i [\delta_3(i)]=3

（4）由最优路径的终点

i_3^*

，逆向找到前面的路径

t=2, \quad i_2^* = \Psi_3(i_3^*)=\Psi_3(3)=3

t=1, \quad i_1^* = \Psi_2(i_2^*)=\Psi_2(3)=3

于是求得最优路径，最优状态序列

Y^*=(i_1^*,i_2^*,i_3^*)=(3,3,3)

4.3 维特比算法Python代码

import numpy as np
class HiddenMarkov:
    def viterbi(self, Q, V, A, B, X, PI):   # 维特比算法，求最优状态序列
        N = len(Q)  # 隐藏状态数量
        M = len(X)  # 观测序列大小
        deltas = np.zeros((N, M))
        psis = np.zeros((N, M))
        Y = np.zeros((1, M))    # 最优（概率最大）状态序列
        for t in range(M):
            realT = t+1 # 真实时刻，从1开始
            indexOfXi = V.index(X[t])   # 观测值Xi对应的V中的索引
            for i in range(N):
                realI = i+1 # 状态序号，从1开始
                if t == 0:
                    deltas[i][t] = PI[0][i]* B[i][indexOfXi]    # 初始值
                    psis[i][t] = 0
                    print("delta1(%d)=pi_%d * B%d(x1)=%.2f * %.2f=%.2f" %
                          (realI, realI, realI, PI[0][i], B[i][indexOfXi], deltas[i][t]))
                    print('psis1(%d)=0' % realI)
                else:
                    deltas[i][t] = np.max(
                        np.multiply([delta[t-1] for delta in deltas],[a[i] for a in A])) \
                                   * B[i][indexOfXi]    # 递推公式
                    print("delta%d(%d)=max[delta%d(j)Aj%d]B%d(x%d)=%.2f*%.2f=%.5f"
                          % (realT, realI, realT-1, realI, realI, realT,
                             np.max(np.multiply([delta[t-1] for delta in deltas], [a[i] for a in A])),
                             B[i][indexOfXi], deltas[i][t]))
                    psis[i][t] = np.argmax(
                        np.multiply([delta[t-1] for delta in deltas],[a[i] for a in A])) + 1    # 序号从1开始
                    print("psis%d(%d)=argmax[delta%d(j)Aj%d]=%d" %
                          (realT, realI, realT-1, realI, psis[i][t]))
        print(deltas)
        print(psis)
        Y[0][M-1] = np.argmax([delta[M-1] for delta in deltas]) +1  # 序号从1开始
        print("Y%d=argmax[deltaT(i)]=%d" %(M, Y[0][M-1]))   # 最优路径的终点
        for t in range(M-2, -1, -1):    # 逆序推导前面的路径
            Y[0][t] = psis[int(Y[0][t+1])-1][t+1] # 寻找前面路径
            print("Y%d=psis%d(Y%d)=%d" %(t+1, t+2, t+2, Y[0][t]))
        print("最大概率的状态序列Y是: ", Y)

Q = [1, 2, 3]
V = ['红', '白']
A = [[0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]]
B = [[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]]
# X = ['红', '白', '红', '红', '白', '红', '白', '白']
X = ['红', '白', '红']    # 书上的例子
PI = [[0.2, 0.4, 0.4]]

hmm = HiddenMarkov()
hmm.viterbi(Q, V, A, B, X, PI)

运行结果：（与上面计算一致）

delta1(1)=pi_1 * B1(x1)=0.20 * 0.50=0.10
psis1(1)=0
delta1(2)=pi_2 * B2(x1)=0.40 * 0.40=0.16
psis1(2)=0
delta1(3)=pi_3 * B3(x1)=0.40 * 0.70=0.28
psis1(3)=0
delta2(1)=max[delta1(j)Aj1]B1(x2)=0.06*0.50=0.02800
psis2(1)=argmax[delta1(j)Aj1]=3
delta2(2)=max[delta1(j)Aj2]B2(x2)=0.08*0.60=0.05040
psis2(2)=argmax[delta1(j)Aj2]=3
delta2(3)=max[delta1(j)Aj3]B3(x2)=0.14*0.30=0.04200
psis2(3)=argmax[delta1(j)Aj3]=3
delta3(1)=max[delta2(j)Aj1]B1(x3)=0.02*0.50=0.00756
psis3(1)=argmax[delta2(j)Aj1]=2
delta3(2)=max[delta2(j)Aj2]B2(x3)=0.03*0.40=0.01008
psis3(2)=argmax[delta2(j)Aj2]=2
delta3(3)=max[delta2(j)Aj3]B3(x3)=0.02*0.70=0.01470
psis3(3)=argmax[delta2(j)Aj3]=3
[[0.1     0.028   0.00756]
 [0.16    0.0504  0.01008]
 [0.28    0.042   0.0147 ]]
[[0. 3. 2.]
 [0. 3. 2.]
 [0. 3. 3.]]
Y3=argmax[deltaT(i)]=3
Y2=psis3(Y3)=3
Y1=psis2(Y2)=3
最大概率的状态序列Y是:  [[3. 3. 3.]]

完整代码

5. PosTagging词性标注实践

基于HMM的中文词性标注 POSTagging

本文参与腾讯云自媒体同步曝光计划，分享自作者个人站点/博客。

原始发表：2019/12/07 ，如有侵权请联系 cloudcommunity@tencent.com 删除

serverless

python

本文分享自作者个人站点/博客前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

serverless

python