一文读懂NLP之隐马尔科夫模型（HMM）详解加python实现

1 隐马尔科夫模型
- 1.1 HMM解决的问题
- 1.2 HMM模型的定义
- - 1.2.1HMM的两个假设
  - 1.2.2 HMM模型
- 1.3 HMM模型的三个基本问题
2 概率计算问题及算法
- 2.1 直接计算法
- 2.2 前向算法
- 2.3 后向算法
- 2.4 一些概率与期望值的计算
3 模型训练问题及算法
- 3.1 监督学习——最大似然估计
- 3.2 非监督学习——EM算法
- 3.3 Baum-Welch算法
4 序列预测问题及算法
- 4.1 近似解法
- 4.2 维特比算法（Viterbi algorithm）
5 hmmlearn使用
- 5.1 pythonceshi工具unnitest
- 5.2 离散HMM测试
- 5.3 HMM测试

1 隐马尔科夫模型

本文主要介绍NLP领域中很重要的一个模型——隐马尔科夫模型（Hidden Markov Model，HMM）。希望读完本文后，大家能够清楚地理解HMM并能够应用到实际中。

隐马尔科夫模型是结构最简单的动态贝叶斯网（dynamic Bayesian network，也被称作有向图模型），HMM是可以用于标注问题的统计数学模型，描述由隐藏的马尔科夫链随机生成观测序列的过程，属于生成模型。HMM模型在语音识别、自然语言处理、生物信息、模式识别等领域有广泛的应用。

1.1 HMM解决的问题

首先看看什么样的问题可以使用HMM模型解决。

使用HMM模型来解决的问题一般有两个特征：

1）问题是基于序列的，比如时间序列、状态序列。
2 ）问题中有两类数据，一类序列数据是可以观测到的，即观测序列；而另一类数据是不能观察到的，即隐藏状态序列，简称状态序列。

如果问题有了这两个特征，那么这个问题一般可以使用HMM模型尝试解决，这样的问题在生活中是很多的。例如，说一句话，发出的声音是观测序列，想表达的意思是隐藏状态序列；写文章的过程可以想象为先在脑海中构思好一个满足语法的词性序列，然后再将每个词性填充为具体的词语。

1.2 HMM模型的定义

隐马尔科夫模型是关于时序的概率模型，描述由一个隐藏的马尔科夫链随机生成不可观测的状态随机序列，再由各个状态生成一个观测而产生观测随机序列的过程（摘自《统计学习方法》）。隐藏的马尔科夫链随机生成的状态的序列，称为状态序列（state sequence），记作y\boldsymbol{y}y；每个状态生成一个观测，而由此产生的观测的随机序列，称为观测序列（observation sequence），记作x\boldsymbol{x}x。序列的每一个位置又可以看作一个时刻。

HMM模型示意图

1.2.1HMM的两个假设

首先介绍一下马尔科夫假设：每个事件的发生概率只取决于前一个事件。将满足该假设的连续多个事件串联在一起，就构成了马尔科夫链。在NLP语境下，可以将事件具象为单词，于是马尔科夫模型就是二元语法。

第一个假设：将马尔科夫假设作用于状态序列，假设当前状态yty_tyt仅仅依赖于前一个状态yt−1y_{t-1}yt−1，连续多个状态构成隐马尔科夫链y\boldsymbol{y}y。数学表达式为：
p(yt∣yt−1,xt−1,yt−2,xt−2,yt−3,xt−3,...,y1,x1)=p(yt∣yt−1),t=1,2,3,...,Tp(y_t|y_{t-1},x_{t-1},y_{t-2},x_{t-2},y_{t-3},x_{t-3},...,y_{1},x_1)=p(y_t|y_{t-1}), t=1,2,3,...,Tp(yt∣yt−1,xt−1,yt−2,xt−2,yt−3,xt−3,...,y1,x1)=p(yt∣yt−1),t=1,2,3,...,T

有了隐马尔科夫链，如何建立与观测序列x\boldsymbol{x}x的联系呢？HMM做了第二个假设:

第二个假设：任意时刻的观测xtx_txt只依赖于该时刻的状态yty_tyt，与其他时刻的状态或观测独立无关。数学表达式为：
p(xt∣yT,xT,yt−1,xt−1,...,yt+1,xt+1,yt,yt−1,xt−1,...y1,x1)=p(xt∣yt),t=1,2,3,...,Tp(x_t|y_T,x_T,y_{t-1},x_{t-1},...,y_{t+1},x_{t+1},y_t,y_{t-1},x_{t-1},...y_1,x_1)=p(x_t|y_t),t=1,2,3,...,Tp(xt∣yT,xT,yt−1,xt−1,...,yt+1,xt+1,yt,yt−1,xt−1,...y1,x1)=p(xt∣yt),t=1,2,3,...,T

1.2.2 HMM模型

设QQQ是所有可能的状态的集合，VVV是所有可能的观测的集合。
Q={q1,q2,...,qN}，V={v1,v2,v3,...,vN}Q=\{q_1,q_2,...,q_N\}，V=\{v_1,v_2,v_3,...,v_N\}Q={q1,q2,...,qN}，V={v1,v2,v3,...,vN}
其中，N是可能的状态数，M是可能的观测数。
HMM模型用三元组来表示λ=(π,A,B)\lambda=(\boldsymbol{\pi},A,B)λ=(π,A,B)：

π\boldsymbol{\pi}π : 初始状态概率向量
A：状态转移概率矩阵
B：发射概率矩阵

系统怎么从零开始呢？ 观测值是由隐藏状态产生的，所以系统最初应该是生成隐藏状态。

初始概率向量指的是系统启动时进入的第一个状态y1y_1y1成为初始状态，y\boldsymbol{y}y有NNN种取值，从QQQ集合中选取一个，即y∈{q1,q2,...,qN}\boldsymbol{y} \in \{q_1,q_2,...,q_N\}y∈{q1,q2,...,qN}。由于y1y_1y1是第一个状态，是一个独立的离散随机变量，可以用p(y1∣π)p(y_1|\boldsymbol{\pi})p(y1∣π)来描述，y1y_1y1只由π\boldsymbol{\pi}π来控制，其中π=(π1,π2,π3,...,πN)T,0≤πi≤1,∑i=1Nπi\boldsymbol{\pi}=(\pi_1,\pi_2,\pi_3,...,\pi_N)^T,0\leq\pi_i\leq1,\sum\limits^{N}_{i=1}{\pi_i}π=(π1,π2,π3,...,πN)T,0≤πi≤1,i=1∑Nπi=1。π\boldsymbol{\pi}π是概率分布的参数向量，称为初始状态概率向量。给定π\boldsymbol{\pi}π，初始状态y1y_1y1的取值分布就确定了。以中文分词问题为例，采用{B,M,E,S}标记时，其中B代表一个词的第一个字，M代表词的中间字，E代表词的末尾字，S代表单字成词。y1y_1y1所有可能的取值及对应的概率如下：

p(y1=B)=0.7p(y_1=B)=0.7p(y1=B)=0.7

p(y1=M)=0p(y_1=M)=0p(y1=M)=0

p(y1=E)=0p(y_1=E)=0p(y1=E)=0

p(y1=s)=0.3p(y_1=s)=0.3p(y1=s)=0.3

则π=[0.7,0,0,0.3]\pi=[0.7,0,0,0.3]π=[0.7,0,0,0.3]，也就是说句子第一个字可能是一个单字或者一个词的首字，不可能是一个词的中间或者尾字。

y1y_1y1确定之后，怎么产生x1x_1x1呢？如何确定x1x_1x1的概率分布呢？

根据第二个假设：当前观测值x1x_1x1仅取决于当前的状态y1y_1y1，对于从QQQ中取出的每一种状态y1y_1y1，x1x_1x1都可以从VVV集合中的MMM个值选一个，所以对于每一个y,xy,xy,x都是一个独立的离散随机变量，其概率参数对应一个向量，维度为MMM，即x∈{v1,v2,...,vN}\boldsymbol{x} \in \{v_1,v_2,...,v_N\}x∈{v1,v2,...,vN}。由于一共有NNN种yyy，所以这些向量构成了一个N×MN\times MN×M矩阵，称为发射概率矩阵B\boldsymbol{B}B。
B=[bij]N×M=[p(xt=vi∣yt=qj)]N×M\boldsymbol{B}=[b_{ij}]_{N\times M}=[p(x_t=v_i|y_t=q_j)]_{N\times M}B=[bij]N×M=[p(xt=vi∣yt=qj)]N×M

其中，p(xt=vi∣yt=qj)p(x_t=v_i|y_t=q_j)p(xt=vi∣yt=qj)代表t时刻，隐藏状态yty_tyt是qjq_jqj，由这个状态产生的观测值xtx_txt等于viv_ivi的概率。

状态	观测值1	观测值2	…	观测值M
状态1
状态2
…
状态N

发射概率矩阵是一个非常形象的术语：可以将y\boldsymbol{y}y想象成为不同的彩弹枪，x\boldsymbol{x}x为不同颜色的子弹，每把彩弹枪内的颜色子弹比例不一样，导致有的彩弹枪红色子弹较多比较容易发射红色彩弹，一些彩弹枪绿色子弹较多更容易发射绿色彩弹。
发射概率在中文分词中也具有实际意义，有些字符构词的位置比较固定，比如一把作为词首的枪，不容易发射出“餮”，因为“餮”一般作为“饕餮”的词尾出现。通过给p(x1=餮∣y1=B)p(x_1=餮|y_1=B)p(x1=餮∣y1=B)较低的概率，HMM模型可以有效的防止“饕餮”被错误的切分。

y1y_1y1确定之后，如何转移状态到y2y_2y2?乃至yny_nyn？

根据HMM模型的第一个假设：t+1t+1t+1时刻的状态仅仅取决于ttt时刻的状态。类似发射概率矩阵，对于ttt时刻的每一种状态，yt+1y_{t+1}yt+1是一个离散的随机变量，取值有NNN种。ttt时刻一共可能有NNN种状态，所以从ttt时刻到t+1t+1t+1时刻的状态转移矩阵为N×NN\times NN×N的方阵，称为状态转移概率矩阵A\boldsymbol{A}A:
A=[aij]N×N=[p(yt+1=vj∣yt=vi)]N×N\boldsymbol{A}=[a_{ij}]_{N\times N}=[p(y_{t+1}=v_j|y_t=v_i)]_{N\times N}A=[aij]N×N=[p(yt+1=vj∣yt=vi)]N×N

其中，p(yt+1=sj∣yt=si)p(y_{t+1}=s_j|y_t=s_i)p(yt+1=sj∣yt=si)表示ttt时刻的状态为viv_ivi，t+1t+1t+1时刻的状态为vjv_jvj的概率。

当前状态	下一状态是状态1	状态2	…	状态N
状态1
状态2
…
状态N

例如，在中文分词中，标签B后面不可能是S，于是给p(yt+1=S∣yt=B)=0p(y_{t+1}=S|y_t=B)=0p(yt+1=S∣yt=B)=0就可以防止B后面接S的情况出现。

初始状态概率向量、状态转移概率矩阵与发射概率矩阵称为HMM模型的三元组λ=(π,A,B)\lambda=(\boldsymbol{\pi},A,B)λ=(π,A,B)，只要三元组确定了，HMM模型就确定了。

举一个经典的例子：
假设有四个盒子，每个盒子里都装有红白两种颜色的球，如下表：

盒子	1	2	3	4
红球数	5	3	6	8
白球数	5	7	4	2

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

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

一共抽取出来5个球，得到一个球的颜色观察序列：
x={红，红，白，白，红}\boldsymbol{x}=\{\color{red}红，红，\color{black}白，白，\color{red}红\}x={红，红，白，白，红}

在这个过程中，观察者只能观测到球的颜色序列，观测不到球是从哪个盒子取出的，即观察不到盒子的序列。
这个例子中有两个随机序列，一个是盒子的序列（状态序列），一个是球的颜色的观测序列（观测序列），前者是隐藏的，后者是可以观测的。根据所给的条件可以明确状态集合、观测集合、序列长度以及模型的三要素。

状态集合是Q={盒子1，盒子2，盒子3，盒子4}，N=4Q=\{盒子1，盒子2，盒子3，盒子4\}，N=4Q={盒子1，盒子2，盒子3，盒子4}，N=4。
观测集合是V={红，白}V=\{\color{red}红，\color{black}白\}V={红，白}。
状态序列和观测序列的长度T=5T=5T=5，原因是一共观测了5次。
初始状态概率向量为π=(0.25,0.25,0.25,0.25)T\boldsymbol{\pi}=(0.25,0.25,0.25,0.25)^Tπ=(0.25,0.25,0.25,0.25)T，原因是第一次是等概率随机抽取一个盒子。
状态转移概率矩阵A=[01000.400.6000.400.6000.50.5]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}A=⎣⎢⎢⎡00.400100.4000.600.5000.60.5⎦⎥⎥⎤
发射概率矩阵B=[0.50.50.30.70.60.40.80.2]B=\begin{bmatrix} 0.5 & 0.5 \\ 0.3 & 0.7 \\ 0.6 & 0.4 \\ 0.8 & 0.2 \\ \end{bmatrix}B=⎣⎢⎢⎡0.50.30.60.80.50.70.40.2⎦⎥⎥⎤

1.3 HMM模型的三个基本问题

有了HMM模型之后，如何使用模型呢？HMM模型一个可以解决三个问题：

概率计算问题：给定模型参数λ=(π,A,B)\lambda=(\boldsymbol{\pi},A,B)λ=(π,A,B)，和一个观测序列x\boldsymbol{x}x,，计算在这个模型参数λ\lambdaλ下，观测序列出现的最大概率，即p(x∣λ)p(\boldsymbol{x}|\lambda)p(x∣λ)的最大值。
模型训练问题：给定训练集(x(i),y(i))(\boldsymbol{x}^{(i)},\boldsymbol{y}^{(i)})(x(i),y(i))，估计模型参数λ=(π,A,B)\lambda=(\boldsymbol{\pi},A,B)λ=(π,A,B)，使得在该模型下观测序列概率p(x∣λ)p(\boldsymbol{x}|\lambda)p(x∣λ)最大，即使用极大似然估计得方法估计参数。
序列预测问题：也称为解码问题，已知模型参数λ=(π,A,B)\lambda=(\boldsymbol{\pi},A,B)λ=(π,A,B)，给定观测序列x\boldsymbol{x}x，求最有可能的状态序列y\boldsymbol{y}y，即求p(y∣x)p(\boldsymbol{y}|\boldsymbol{x})p(y∣x)的最大值。

2 概率计算问题及算法

概率计算问题，也就是在给定的模型参数三元组的条件生成观测序列的过程。给定模型参数λ=(π,A,B)\lambda=(\boldsymbol{\pi},A,B)λ=(π,A,B)和一个观测序列x\boldsymbol{x}x,，计算在这个模型参数λ\lambdaλ下，观测序列出现的最大概率，即p(x∣λ)p(\boldsymbol{x}|\lambda)p(x∣λ)的最大值。先介绍概念上可行但计算上不行的直接计算法（暴力解法），然后介绍前向算法与后向算法。

2.1 直接计算法

给定模型参数λ=(π,A,B)\boldsymbol{\lambda}=(\boldsymbol{\pi},A,B)λ=(π,A,B)和一个观测序列x={x1,x2,x3,...xT}\boldsymbol{x}=\{x_1,x_2,x_3,...x_T\}x={x1,x2,x3,...xT}，计算观测序列出现的概率p(x∣λ)p(\boldsymbol{x}|\lambda)p(x∣λ)，最直接的方法就是按照概率公式直接计算，通过列举所有可能的长度为TTT的状态序列y={y1,y2,y3,...,yT}\boldsymbol{y}=\{y_1,y_2,y_3,...,y_T\}y={y1,y2,y3,...,yT}，求各个状态序列y\boldsymbol{y}y与观测序列x=(x1,x2,x3,...,xT)\boldsymbol{x}=(x_1,x_2,x_3,...,x_T)x=(x1,x2,x3,...,xT)的联合概率p(x,y∣λ)p(\boldsymbol{x},\boldsymbol{y}|\boldsymbol{\lambda})p(x,y∣λ),，然后对所有可能的状态序列求和，得到p(x∣λ)p(\boldsymbol{x}|\boldsymbol{\lambda})p(x∣λ)。
状态序列y=(y1,y2,y3,...,yT)\boldsymbol{y}=(y_1,y_2,y_3,...,y_T)y=(y1,y2,y3,...,yT)发生的概率是：
p(y∣λ)=πy1ay1y2ay2y3...ayT−1yTp(\boldsymbol{y}|\boldsymbol{\lambda})=\pi_{y_1}a_{y_1y_2}a_{y_2y_3}...a_{y_{T-1}y_T}p(y∣λ)=πy1ay1y2ay2y3...ayT−1yT

对固定的状态序列y=(y1,y2,y3,...,yT)\boldsymbol{y}=(y_1,y_2,y_3,...,y_T)y=(y1,y2,y3,...,yT)并且观测序列为x=(x1,x2,x3,...,xT)\boldsymbol{x}=(x_1,x_2,x_3,...,x_T)x=(x1,x2,x3,...,xT)的概率是：
p(x∣y,λ)=by1x1by2x2by3x3...byTxTp(\boldsymbol{x}|\boldsymbol{y},\boldsymbol{\lambda})=b_{y_1x_1}b_{y_2x_2}b_{y_3x_3}...b_{y_Tx_T}p(x∣y,λ)=by1x1by2x2by3x3...byTxT

在给定HMM模型参数λ\boldsymbol{\lambda}λ的条件下，x,y\boldsymbol{x},\boldsymbol{y}x,y同时出现的联合概率为：
p(x,y∣λ)=p(x∣y,λ)p(y∣λ)=πy1by1x1ay1y2by2x2ay2y3by3x3...ayT−1yTbyTxT\begin{aligned}p(\boldsymbol{x},\boldsymbol{y}|\boldsymbol{\lambda})&=p(\boldsymbol{x}|\boldsymbol{y},\boldsymbol{\lambda})p(\boldsymbol{y}|\boldsymbol{\lambda})\\ &=\pi_{y_1}b_{y_1x_1}a_{y_1y_2}b_{y_2x_2}a_{y_2y_3}b_{y_3x_3}...a_{y_{T-1}y_T}b_{y_Tx_T} \end{aligned}p(x,y∣λ)=p(x∣y,λ)p(y∣λ)=πy1by1x1ay1y2by2x2ay2y3by3x3...ayT−1yTbyTxT

然后，对所有可能的状态序列y\boldsymbol{y}y求和，得到观测序列x\boldsymbol{x}x的概率p(x∣λ)p(\boldsymbol{x}|\lambda)p(x∣λ)，即：
p(x∣λ)=∑yp(x,y∣λ)=∑yp(x∣y,λ)p(y∣λ)=∑y1,y2,y3,...,yTπy1by1x1ay1y2by2x2ay2y3by3x3...ayT−1yTbyTxT\begin{aligned} p(\boldsymbol{x}|\lambda)&=\sum\limits_\boldsymbol{y}p(\boldsymbol{x},\boldsymbol{y}|\boldsymbol{\lambda})\\ &=\sum\limits_\boldsymbol{y}p(\boldsymbol{x}|\boldsymbol{y},\boldsymbol{\lambda})p(\boldsymbol{y}|\boldsymbol{\lambda})\\ &=\sum\limits_{y_1,y_2,y_3,...,y_T}\pi_{y_1}b_{y_1x_1}a_{y_1y_2}b_{y_2x_2}a_{y_2y_3}b_{y_3x_3}...a_{y_{T-1}y_T}b_{y_Tx_T} \end{aligned}p(x∣λ)=y∑p(x,y∣λ)=y∑p(x∣y,λ)p(y∣λ)=y1,y2,y3,...,yT∑πy1by1x1ay1y2by2x2ay2y3by3x3...ayT−1yTbyTxT

利用这个公式计算量很大，复杂度是O(TNT)O(TN^T)O(TNT)，在实际中是不可行的。

2.2 前向算法

首先需要定义前向概率：给定HMM模型λ\boldsymbol{\lambda}λ，定义时刻ttt部分观测序列为x=(x1,x2,x3,...,xt)\boldsymbol{x}=(x_1,x_2,x_3,...,x_t)x=(x1,x2,x3,...,xt)且状态为qiq_iqi的概率为前向概率，记作：
αt(i)=p(x1,x2,x3,...,xt,yt=qi∣λ)\alpha_t(i)=p(x_1,x_2,x_3,...,x_t,y_t=q_i|\boldsymbol{\lambda})αt(i)=p(x1,x2,x3,...,xt,yt=qi∣λ)

前向算法：

输入：HMM模型参数λ\boldsymbol{\lambda}λ，观测序列x=(x1,x2,x3,...,xt)\boldsymbol{x}=(x_1,x_2,x_3,...,x_t)x=(x1,x2,x3,...,xt)；
输出：观测序列概率p(x∣λ)p(\boldsymbol{x}|\boldsymbol{\lambda})p(x∣λ)。

算法步骤：

初始：
根据π\boldsymbol{\pi}π生成t=1t=1t=1时刻的状态y1=qiy_1=q_iy1=qi，概率为πi\pi_iπi，并且根据发射概率矩阵B\boldsymbol{B}B由y1=qiy_1=q_iy1=qi生成x1x_1x1，概率为bix1b_{ix_1}bix1，则：
α1(i)=πibix1\alpha_1(i)=\pi_ib_{ix_1}α1(i)=πibix1
递推：
当t=2t=2t=2时，根据状态转移概率矩阵A\boldsymbol{A}A，系统的状态由y1=qjy_1=q_jy1=qj变为y2=qiy_2=q_iy2=qi，概率为ajia_{ji}aji，p(y2=qi,y1=qj,x1∣λ)=α1(j)ajip(y_2=q_i,y_1=q_j,x_1|\boldsymbol{\lambda})=\alpha_1(j)a_{ji}p(y2=qi,y1=qj,x1∣λ)=α1(j)aji。不管y1y_1y1为什么状态，都可能转移到状态y2=qiy_2=q_iy2=qi，所以需要对y1y_1y1求和:∑j=1Nα1(j)aji\sum\limits_{j=1}^N\alpha_1(j)a_{ji}j=1∑Nα1(j)aji。根据αt(i)\alpha_t(i)αt(i)的定义，t=2t=2t=2时刻，由状态y2=qiy_2=q_iy2=qi产生观测值x2x_2x2的概率为bix2b_{ix_2}bix2，所以：
α2(i)=[∑j=1Nα1(j)aji]bix2\alpha_2(i)=[\sum\limits_{j=1}^N\alpha_1(j)a_{ji}]b_{ix_2}α2(i)=[j=1∑Nα1(j)aji]bix2

对t=1,2,3,...T−1t=1,2,3,...T-1t=1,2,3,...T−1，有：
αt+1(i)=[∑j=1Nαt(j)aji]bixt+1,i=1,2,3,...N\alpha_{t+1}(i)=[\sum\limits^N_{j=1}\alpha_t(j)a_{ji}]b_{ix_{t+1}},\quad i=1,2,3,...Nαt+1(i)=[j=1∑Nαt(j)aji]bixt+1,i=1,2,3,...N

3. 终止：
时刻t=Tt=Tt=T：
根据定义可知：αT(i)\alpha_T(i)αT(i)表示，t=Tt=Tt=T时刻，处于状态yT=qiy_T=q_iyT=qi，并且观察到的序列为x=(x1,x2,x3,...,xT)\boldsymbol{x}=(x_1,x_2,x_3,...,x_T)x=(x1,x2,x3,...,xT)的概率：
αT(i)=p(x1,x2,x3,...,xT,yT=qi∣λ)=p(x,yT=qi∣λ)\begin{aligned} \alpha_T(i)&=p(x_1,x_2,x_3,...,x_T,y_T=q_i|\boldsymbol{\lambda})\\ &=p(\boldsymbol{x},y_T=q_i|\boldsymbol{\lambda}) \end{aligned}αT(i)=p(x1,x2,x3,...,xT,yT=qi∣λ)=p(x,yT=qi∣λ)

所以：

p(x∣λ)=∑i=1NαT(i)p(\boldsymbol{x}|\boldsymbol{\lambda})=\sum\limits^N_{i=1}\alpha_T(i)p(x∣λ)=i=1∑NαT(i)

前向算法实际是基于“状态序列的路径结构”递推计算p(x∣λ)p(\boldsymbol{x}|\boldsymbol{\lambda})p(x∣λ)的算法。前向算法高效的关键是其局部计算前向概率，然后利用路径结构将前向概率递推到全局，得到p(x∣λ)p(\boldsymbol{x}|\boldsymbol{\lambda})p(x∣λ)。在t=1t=1t=1时刻计算α1(i)\alpha_1(i)α1(i)的NNN个值，而当t≥2t\geq2t≥2时，计算每一个\alpha_t(i)的值都会利用前NNN个αt−1(i)\alpha_{t-1}(i)αt−1(i)的值，减少计算量的原因在于每一次计算都直接引用前一个时刻的计算结果，避免了概率的重复计算，复杂度为O(N∗N∗T)O(N*N*T)O(N∗N∗T)，而不是直接算法的O(TNT)O(TN^T)O(TNT)。

例：上文中的盒子和球模型，状态集合为Q={1,2,3}Q=\{1,2,3\}Q={1,2,3}，观测集合为V={红，白}V=\{红，白\}V={红，白}，
A=[0.50.20.30.30.50.20.20.30.5]A=\begin{bmatrix} 0.5 & 0.2 & 0.3\\ 0.3 & 0.5 & 0.2\\ 0.2 & 0.3 & 0.5 \end{bmatrix}A=⎣⎡0.50.30.20.20.50.30.30.20.5⎦⎤ ， B=[0.50.50.40.60.70.3]B=\begin{bmatrix} 0.5 & 0.5 \\ 0.4 & 0.6\\ 0.7 & 0.3 \end{bmatrix}B=⎣⎡0.50.40.70.50.60.3⎦⎤ ， π=(0.2,0.4,0.4)T\pi=(0.2,0.4,0.4)^Tπ=(0.2,0.4,0.4)T
设T=3T=3T=3，x=(红，白，红)\boldsymbol{x}=(红，白，红)x=(红，白，红)，用前向算法来计算p(x∣λ)p(\boldsymbol{x}|\boldsymbol{\lambda})p(x∣λ)。
解：

计算初值：
α1(1)=π1b1x1=0.2×0.5=0.10\alpha_1(1)=\pi_1b_{1x_1}=0.2\times0.5=0.10α1(1)=π1b1x1=0.2×0.5=0.10

α1(2)=π2b2x1=0.4×0.4=0.16\alpha_1(2)=\pi_2b_{2x_1}=0.4\times0.4=0.16α1(2)=π2b2x1=0.4×0.4=0.16

α1(3)=π3b3x1=0.4×0.7=0.28\alpha_1(3)=\pi_3b_{3x_1}=0.4\times0.7=0.28α1(3)=π3b3x1=0.4×0.7=0.28
2. 递推计算：
α2(1)=[∑i=13α1(i)ai1]b1x2=[0.10×0.5+0.16×0.3+0.28×0.2)]×0.5=0.077\alpha_2(1)=[\sum\limits_{i=1}^3\alpha_1(i)a_{i1}]b_{1x_2}=[0.10\times0.5+0.16\times0.3+0.28\times0.2)]\times0.5=0.077α2(1)=[i=1∑3α1(i)ai1]b1x2=[0.10×0.5+0.16×0.3+0.28×0.2)]×0.5=0.077

α2(2)=[∑i=13α1(i)ai2]b2x2=[0.10×0.2+0.16×0.5+0.28×0.3)]×0.5=0.1104\alpha_2(2)=[\sum\limits_{i=1}^3\alpha_1(i)a_{i2}]b_{2x_2}=[0.10\times0.2+0.16\times0.5+0.28\times0.3)]\times0.5=0.1104α2(2)=[i=1∑3α1(i)ai2]b2x2=[0.10×0.2+0.16×0.5+0.28×0.3)]×0.5=0.1104

α2(3)=[∑i=13α1(i)ai3]b3x2=[0.10×0.3+0.16×0.2+0.28×0.5)]×0.5=0.0606\alpha_2(3)=[\sum\limits_{i=1}^3\alpha_1(i)a_{i3}]b_{3x_2}=[0.10\times0.3+0.16\times0.2+0.28\times0.5)]\times0.5=0.0606α2(3)=[i=1∑3α1(i)ai3]b3x2=[0.10×0.3+0.16×0.2+0.28×0.5)]×0.5=0.0606

α3(1)=[∑i=13α2(i)ai1]b1x3=0.04187\alpha_3(1)=[\sum\limits_{i=1}^3\alpha_2(i)a_{i1}]b_{1x_3}=0.04187α3(1)=[i=1∑3α2(i)ai1]b1x3=0.04187

α3(2)=[∑i=13α2(i)ai2]b2x3=0.03551\alpha_3(2)=[\sum\limits_{i=1}^3\alpha_2(i)a_{i2}]b_{2x_3}=0.03551α3(2)=[i=1∑3α2(i)ai2]b2x3=0.03551

α3(3)=[∑i=13α2(i)ai3]b3x3=0.05284\alpha_3(3)=[\sum\limits_{i=1}^3\alpha_2(i)a_{i3}]b_{3x_3}=0.05284α3(3)=[i=1∑3α2(i)ai3]b3x3=0.05284

终止：

p(x∣λ)=∑i=13α3(i)=0.13022p(\boldsymbol{x}|\boldsymbol{\lambda})=\sum\limits_{i=1}^3\alpha_3(i)=0.13022p(x∣λ)=i=1∑3α3(i)=0.13022

2.3 后向算法

类似前向算法，定义：给定HMM模型λ\boldsymbol{\lambda}λ，定义时刻ttt状态为qiq_iqi的条件下，从t+1到Tt+1到Tt+1到T的部分观测序列为xt+1,xt+2,xt+3,...,xTx_{t+1},x_{t+2},x_{t+3},...,x_Txt+1,xt+2,xt+3,...,xT的概率为后向概率，记作：
βt(i)=p(xt+1,xt+2,xt+3,...,xT∣yt=qi,λ)\beta_t(i)=p(x_{t+1},x_{t+2},x_{t+3},...,x_T|y_t=q_i,\boldsymbol{\lambda})βt(i)=p(xt+1,xt+2,xt+3,...,xT∣yt=qi,λ)
后向算法：

输入：HMM模型参数λ\boldsymbol{\lambda}λ，观测序列x=(x1,x2,x3,...,xt)\boldsymbol{x}=(x_1,x_2,x_3,...,x_t)x=(x1,x2,x3,...,xt)；
输出：观测序列概率p(x∣λ)p(\boldsymbol{x}|\boldsymbol{\lambda})p(x∣λ)。

当t=Tt=Tt=T时:
βT(i)=p(这里没东西∣yT=qi,λ)=1\beta_T(i)=p(这里没东西|y_T=q_i,\boldsymbol{\lambda})=1βT(i)=p(这里没东西∣yT=qi,λ)=1
这里可以这么理解：按理说好需要看看T时刻过后是什么观测值，但是T时刻之后没有观测值了，不管T时刻状态是什么，之后就是没有观测值，所以没有观测值的概率为1，且与状态无关。
对t=T−1t=T-1t=T−1:
已知βT(j)\beta_T(j)βT(j)，根据HMM模型可知，yT=qjy_T=q_jyT=qj是由yT−1y_{T-1}yT−1转移而来的，假设yT−1=qiy_{T-1}=q_iyT−1=qi，转移的概率为aija_{ij}aij:
根据βT−1(i)\beta_{T-1}(i)βT−1(i)的定义，时刻T−1T-1T−1的状态yT−1=qiy_{T-1}=q_iyT−1=qi，从T到T的观测序列为xTx_TxT，T时刻状态yT=qjy_T=q_jyT=qj生成xTx_TxT的概率为bjxTb_{jx_T}bjxT，则：
βT−1(i)=∑j=1NaijbjxTβT(j)\beta_{T-1}(i)=\sum\limits^N_{j=1}a_{ij}b_{jx_T}\beta_T(j)βT−1(i)=j=1∑NaijbjxTβT(j)
对t=T−1,T−2,...,1t=T-1,T-2,...,1t=T−1,T−2,...,1:
bjxt=p(xt∣yt=qj,λ))b_{jx_t}=p(x_t|y_t=q_j,\boldsymbol{\lambda}))bjxt=p(xt∣yt=qj,λ))

aij=p(yt+1=qj∣yt=qi,λ))a_{ij}=p(y_{t+1}=q_j|y_t=q_i,\boldsymbol{\lambda}))aij=p(yt+1=qj∣yt=qi,λ))

βt+1(j)=p(xt+2,xt+3,...,xT∣yt+1=qj,λ)\beta_{t+1}(j)=p(x_{t+2},x_{t+3},...,x_T|y_{t+1}=q_j,\boldsymbol{\lambda})βt+1(j)=p(xt+2,xt+3,...,xT∣yt+1=qj,λ)

由状态yt+1=qjy_{t+1}=q_jyt+1=qj生成t+1t+1t+1时刻的观测值xt+1x_{t+1}xt+1:
bjxt+1βt+1(j)=p(xt+1∣yt+1=qj,λ)p(xt+2,xt+3,...,xT∣yt+1=qj,λ)=p(xt+1,xt+2,...,xT∣yt+1=qj,λ)\begin{aligned} b_{jx_{t+1}}\beta_{t+1}(j)&=p(x_{t+1}|y_{t+1}=q_j,\boldsymbol{\lambda})p(x_{t+2},x_{t+3},...,x_T|y_{t+1}=q_j,\boldsymbol{\lambda})\\ &=p(x_{t+1},x_{t+2},...,x_T|y_{t+1}=q_j,\boldsymbol{\lambda})\end{aligned}bjxt+1βt+1(j)=p(xt+1∣yt+1=qj,λ)p(xt+2,xt+3,...,xT∣yt+1=qj,λ)=p(xt+1,xt+2,...,xT∣yt+1=qj,λ)

按照条件概率来理解：在t+1t+1t+1时刻状态为yt+1=qjy_{t+1}=q_jyt+1=qj的条件下，观察到xt+1,xt+2,...,xTx_{t+1},x_{t+2},...,x_Txt+1,xt+2,...,xT的概率为bjxt+1βt+1(j)b_{jx_{t+1}}\beta_{t+1}(j)bjxt+1βt+1(j)，由于βt(i)\beta_t(i)βt(i)与ttt时刻的状态yty_tyt有关，根据HMM模型的第一个假设，yt+1y_{t+1}yt+1仅仅与yty_tyt有关，且概率aija_{ij}aij由状态转移矩阵矩阵A提供。βt(i)\beta_t(i)βt(i)表示在ttt时刻状态为yt=qiy_t=q_iyt=qi的条件下，观察到xt+1,xt+2,...,xTx_{t+1},x_{t+2},...,x_Txt+1,xt+2,...,xT的概率，这个概率若要用βt+1(j)\beta_{t+1}(j)βt+1(j)来表示，针对每一个yt+1=qjy_{t+1}=q_jyt+1=qj，都要乘以从yt=qiy_t=q_iyt=qi到yt+1=qjy_{t+1}=q_jyt+1=qj的状态转移概率aija_{ij}aij，yt+1y_{t+1}yt+1一共有N种状态，所以：
aijbjxt+1βt+1(j)=p(xt+1,xt+2,...,xT∣yt+1=qj,λ)aij=p(xt+1,xt+2,...,xT∣yt+1=qj,λ)p(yt+1=qj∣yt=qi,λ))=p(xt+1,xt+2,...,xT,yt+1=qj,∣yt=qi,λ)\begin{aligned} a_{ij}b_{jx_{t+1}}\beta_{t+1}(j)&=p(x_{t+1},x_{t+2},...,x_T|y_{t+1}=q_j,\boldsymbol{\lambda})a_{ij}\\ &=p(x_{t+1},x_{t+2},...,x_T|y_{t+1}=q_j,\boldsymbol{\lambda})p(y_{t+1}=q_j|y_t=q_i,\boldsymbol{\lambda}))\\ &=p(x_{t+1},x_{t+2},...,x_T,y_{t+1}=q_j,|y_t=q_i,\boldsymbol{\lambda}) \end{aligned}aijbjxt+1βt+1(j)=p(xt+1,xt+2,...,xT∣yt+1=qj,λ)aij=p(xt+1,xt+2,...,xT∣yt+1=qj,λ)p(yt+1=qj∣yt=qi,λ))=p(xt+1,xt+2,...,xT,yt+1=qj,∣yt=qi,λ)

βt(i)=p(xt+1,xt+2,...,xT∣yt=qi,λ)=∑yt+1p(xt+1,xt+2,...,xT,yt+1=qj,∣yt=qi,λ)=∑j=1Naijbjxt+1βt+1(j)\begin{aligned}\beta_t(i)&=p(x_{t+1},x_{t+2},...,x_T|y_t=q_i,\boldsymbol{\lambda})\\ &=\sum\limits_{y_{t+1}}p(x_{t+1},x_{t+2},...,x_T,y_{t+1}=q_j,|y_t=q_i,\boldsymbol{\lambda})\\ &=\sum\limits_{j=1}^Na_{ij}b_{jx_{t+1}}\beta_{t+1}(j) \end{aligned}βt(i)=p(xt+1,xt+2,...,xT∣yt=qi,λ)=yt+1∑p(xt+1,xt+2,...,xT,yt+1=qj,∣yt=qi,λ)=j=1∑Naijbjxt+1βt+1(j)
所以：
βt(i)=∑j=1Naijbjxt+1βt+1(j),i=1,2,3,...,N\beta_t(i)=\sum\limits_{j=1}^Na_{ij}b_{jx_{t+1}}\beta_{t+1}(j),\quad i=1,2,3,...,Nβt(i)=j=1∑Naijbjxt+1βt+1(j),i=1,2,3,...,N

3. 当t=1t=1t=1时，β1(i)=p(x2,x3,x4,...,xT∣y1=qi,λ)\beta_1(i)=p(x_2,x_3,x_4,...,x_T|y_1=q_i,\boldsymbol{\lambda})β1(i)=p(x2,x3,x4,...,xT∣y1=qi,λ)；按照步骤2的思想，针对t=1t=1t=1时刻的每一种状态y1=qiy_1=q_iy1=qi，都需要乘上初始状态概率向量和相应的发射概率矩阵。
p(x∣λ)=∑i=1Nπibix1β1(i)p(\boldsymbol{x}|\boldsymbol{\lambda})=\sum\limits_{i=1}^N\pi_ib_{ix_1}\beta_{1}(i)p(x∣λ)=i=1∑Nπibix1β1(i)

2.4 一些概率与期望值的计算

利用前向概率和后向概率，可以得到关于单个状态和两个状态概率的计算公式。

给定模型λ\boldsymbol{\lambda}λ和观测x\boldsymbol{x}x，在时刻ttt处于状态qiq_iqi的概率，记作γt(i)\gamma_t(i)γt(i)：
γt(i)=p(yt=qi∣x,λ)\gamma_t(i)=p(y_t=q_i|\boldsymbol{x},\boldsymbol{\lambda})γt(i)=p(yt=qi∣x,λ)

可以用过前向和后向算法来算：

γt(i)=p(yt=qi∣x,λ)=p(yt=qi,x∣λ)p(x∣λ)\begin{aligned} \gamma_t(i)=p(y_t=q_i|\boldsymbol{x},\boldsymbol{\lambda})=\frac{p(y_t=q_i,\boldsymbol{x}|\boldsymbol{\lambda})}{p(\boldsymbol{x}|\boldsymbol{\lambda})} \end{aligned}γt(i)=p(yt=qi∣x,λ)=p(x∣λ)p(yt=qi,x∣λ)

αt(i)\alpha_t(i)αt(i)定义为时刻ttt部分观测序列为x=(x1,x2,x3,...,xt)\boldsymbol{x}=(x_1,x_2,x_3,...,x_t)x=(x1,x2,x3,...,xt)且状态为qiq_iqi的概率，βt(i)\beta_t(i)βt(i)定义为时刻ttt状态为qiq_iqi的条件下，从t+1到Tt+1到Tt+1到T的部分观测序列为xt+1,xt+2,xt+3,...,xTx_{t+1},x_{t+2},x_{t+3},...,x_Txt+1,xt+2,xt+3,...,xT的概率，则两者相乘表示，观察序列为x=(x1,x2,x3,...,xT)\boldsymbol{x}=(x_1,x_2,x_3,...,x_T)x=(x1,x2,x3,...,xT)且状态为qiq_iqi的概率，数学表示为：
αt(i)βt(i)=p(yt=qi,x∣λ)\alpha_t(i)\beta_t(i)=p(y_t=q_i,\boldsymbol{x}|\boldsymbol{\lambda})αt(i)βt(i)=p(yt=qi,x∣λ)

p(x∣λ)=∑ytp(yt=qi,x∣λ)=∑j=1Nαt(i)βt(i)p(\boldsymbol{x}|\boldsymbol{\lambda})=\sum\limits_{y_t}p(y_t=q_i,\boldsymbol{x}|\boldsymbol{\lambda})=\sum\limits_{j=1}^N\alpha_t(i)\beta_t(i)p(x∣λ)=yt∑p(yt=qi,x∣λ)=j=1∑Nαt(i)βt(i)

所以：

γt(i)=αt(i)βt(i)p(x∣λ)=αt(i)βt(i)∑j=1Nαt(i)βt(i)\gamma_t(i)=\frac{\alpha_t(i)\beta_t(i)}{p(\boldsymbol{x}|\boldsymbol{\lambda})}=\frac{\alpha_t(i)\beta_t(i)}{\sum\limits_{j=1}^N\alpha_t(i)\beta_t(i)}γt(i)=p(x∣λ)αt(i)βt(i)=j=1∑Nαt(i)βt(i)αt(i)βt(i)
2. 给定HMM模型参数λ\boldsymbol{\lambda}λ和观测序列x\boldsymbol{x}x，在时刻ttt处于状态qiq_iqi且在时刻t+1t+1t+1处于状态qjq_jqj的概率，记作ξt(i,j)\xi_t(i,j)ξt(i,j)：
ξt(i,j)=p(yt=qi,yt+1=qj∣x,λ)\xi_t(i,j)=p(y_t=q_i,y_{t+1}=q_j|\boldsymbol{x},\boldsymbol{\lambda})ξt(i,j)=p(yt=qi,yt+1=qj∣x,λ)

使用前向和后向算法来计算：
ξt(i,j)=p(yt=qi,yt+1=qj∣x,λ)=p(yt=qi,yt+1=qj,x∣λ)p(x,λ)=p(yt=qi,yt+1=qj,x∣λ)∑i=1N∑j=1Np(yt=qi,yt+1=qj,x∣λ)\begin{aligned} \xi_t(i,j)&=p(y_t=q_i,y_{t+1}=q_j|\boldsymbol{x},\boldsymbol{\lambda})\\ &=\frac{p(y_t=q_i,y_{t+1}=q_j,\boldsymbol{x}|\boldsymbol{\lambda})}{p(\boldsymbol{x},\boldsymbol{\lambda)}}\\ &=\frac{p(y_t=q_i,y_{t+1}=q_j,\boldsymbol{x}|\boldsymbol{\lambda})}{\sum\limits_{i=1}^N\sum\limits_{j=1}^Np(y_t=q_i,y_{t+1}=q_j,\boldsymbol{x}|\boldsymbol{\lambda})} \end{aligned}ξt(i,j)=p(yt=qi,yt+1=qj∣x,λ)=p(x,λ)p(yt=qi,yt+1=qj,x∣λ)=i=1∑Nj=1∑Np(yt=qi,yt+1=qj,x∣λ)p(yt=qi,yt+1=qj,x∣λ)

而p(yt=qi,yt+1=qj,x∣λ)p(y_t=q_i,y_{t+1}=q_j,\boldsymbol{x}|\boldsymbol{\lambda})p(yt=qi,yt+1=qj,x∣λ)表示在模型参数下，观测序列x\boldsymbol{x}x以及ttt时刻状态为qiq_iqi且时刻t+1t+1t+1处于状态qjq_jqj的概率，αt(i)=p(yt=qi,x1,x2,x3,...,xt∣λ)\alpha_t(i)=p(y_t=q_i,x_1,x_2,x_3,...,x_t|\boldsymbol{\lambda})αt(i)=p(yt=qi,x1,x2,x3,...,xt∣λ)表示在模型参数下，观测序列x1,x2,x3,...,xtx_1,x_2,x_3,...,x_tx1,x2,x3,...,xt以及ttt时刻状态为qiq_iqi的概率，βt+1(j)=p(xt+2,xt+3,...,xT∣yt+1=qj,λ)\beta_{t+1}(j)=p(x_{t+2},x_{t+3},...,x_T|y_{t+1}=q_j,\boldsymbol{\lambda})βt+1(j)=p(xt+2,xt+3,...,xT∣yt+1=qj,λ)表示在模型参数和t+1t+1t+1时刻状态为qjq_jqj的条件下，观察序列为xt+2,xt+3,...,xTx_{t+2},x_{t+3},...,x_Txt+2,xt+3,...,xT，两者之间差一个从时刻t到t+1t到t+1t到t+1的状态转移概率以及时刻t+1t+1t+1状态qjq_jqj产生序列xt+1x_{t+1}xt+1的概率：
p(yt=qi,yt+1=qj,x∣λ)=p(yt=qi,x1,x2,...,xt∣λ)p(yt+1∣yt,λ)p(xt+1∣yt+1,λ)p(xt+2,xt+3,...,xT∣yt+1=qj,λ)=αt(i)aijbjxt+1βt+1(j)\begin{aligned} p(y_t=q_i,y_{t+1}=q_j,\boldsymbol{x}|\boldsymbol{\lambda})&=p(y_t=q_i,x_1,x_2,...,x_t|\boldsymbol{\lambda})p(y_{t+1}|y_t,\boldsymbol{\lambda})\\&\qquad p(x_{t+1}|y_{t+1},\boldsymbol{\lambda})p(x_{t+2},x_{t+3},...,x_T|y_{t+1}=q_j,\boldsymbol{\lambda})\\ &=\alpha_t(i)a_{ij}b_{jx_{t+1}}\beta_{t+1}(j) \end{aligned}p(yt=qi,yt+1=qj,x∣λ)=p(yt=qi,x1,x2,...,xt∣λ)p(yt+1∣yt,λ)p(xt+1∣yt+1,λ)p(xt+2,xt+3,...,xT∣yt+1=qj,λ)=αt(i)aijbjxt+1βt+1(j)

所以：
ξt(i,j)=αt(i)aijbjxt+1βt+1(j)∑i=1N∑j=1Nαt(i)aijbjxt+1βt+1(j)\xi_t(i,j)=\frac{\alpha_t(i)a_{ij}b_{jx_{t+1}}\beta_{t+1}(j)}{\sum\limits_{i=1}^N\sum\limits_{j=1}^N\alpha_t(i)a_{ij}b_{jx_{t+1}}\beta_{t+1}(j)}ξt(i,j)=i=1∑Nj=1∑Nαt(i)aijbjxt+1βt+1(j)αt(i)aijbjxt+1βt+1(j)

3.将ξt(i,j)\xi_t(i,j)ξt(i,j)和γt(i)\gamma_t(i)γt(i)对各个时刻求和，可以得到一些有用的期望值：

在观测x\boldsymbol{x}x下，状态qiq_iqi出现的期望值：∑t=1Tγt(i)\sum\limits_{t=1}^T\gamma_t(i)t=1∑Tγt(i)
在观测x\boldsymbol{x}x下,由状态qiq_iqi转移的期望值：∑t=1T−1γt(i)\sum\limits_{t=1}^{T-1}\gamma_t(i)t=1∑T−1γt(i)
在观测x\boldsymbol{x}x下,由状态qiq_iqi转移到状态qjq_jqj的期望值：∑t=1T−1ξt(i,j)\sum\limits_{t=1}^{T-1}\xi_t(i,j)t=1∑T−1ξt(i,j)

3 模型训练问题及算法

给定训练集(x(i),y(i))(\boldsymbol{x}^{(i)},\boldsymbol{y}^{(i)})(x(i),y(i))，估计模型参数λ=(π,A,B)\lambda=(\boldsymbol{\pi},A,B)λ=(π,A,B)，使得在该模型下观测序列概率p(x∣λ)p(\boldsymbol{x}|\lambda)p(x∣λ)最大。根据训练数据是否有状态序列数据分为：完全数据和非完全数据，分别使用监督学习和非监督学习实现。

3.1 监督学习——最大似然估计

在监督学习中，我们使用极大似然法来估计HMM模型参数。
假设给定训练数据包含S个长度相同的观测序列{(x1,y1),(x2,y2),...,(xS,yS)}\{(\boldsymbol{x}^1,\boldsymbol{y}^1),(\boldsymbol{x}^2,\boldsymbol{y}^2),...,(\boldsymbol{x}^S,\boldsymbol{y}^S)\}{(x1,y1),(x2,y2),...,(xS,yS)}，使用极大似然法来估计HMM模型参数。

初始状态概率向量的估计：
统计S个样本中，初始状态为qiq_iqi的频率。
π^i=NqiS\hat{\pi}_i=\frac{N_{q_i}}{S}π^i=SNqi
其中，NqiN_{q_i}Nqi是初始状态为qiq_iqi的样本数量，S是样本的数量。

状态转移概率矩阵的估计：
设样本中时刻t处于状态qiq_iqi，时刻t+1处于状态qjq_jqj的频数为AijA_{ij}Aij，那么状态转移概率矩阵的估计为：
a^ij=Aij∑j=1NAij,j=1,2,3,...,N;i=1,2,3,...,N\hat{a}_{ij}=\frac{A_{ij}}{\sum\limits_{j=1}^NA_{ij}},\quad j=1,2,3,...,N;\quad i=1,2,3,...,Na^ij=j=1∑NAijAij,j=1,2,3,...,N;i=1,2,3,...,N

发射概率矩阵的估计：
设样本中状态为iii并观测值为jjj的频数BijB_{ij}Bij，那么状态为iii观测为jjj的概率bijb_{ij}bij的估计为：
b^ij=Bij∑j=1MBij,j=1,2,3,...,M;i=1,2,3,...,N\hat{b}_{ij}=\frac{B_{ij}}{\sum\limits_{j=1}^MB_{ij}},\quad j=1,2,3,...,M;\quad i=1,2,3,...,Nb^ij=j=1∑MBijBij,j=1,2,3,...,M;i=1,2,3,...,N

监督学习的方法就是拿频率来估计概率。

3.2 非监督学习——EM算法

假设给定训练数据只包含S个长度为T的观测序列x={x1,x2,,...,xS}\boldsymbol{x}=\{\boldsymbol{x}^1,\boldsymbol{x}^2,,...,\boldsymbol{x}^S\}x={x1,x2,,...,xS}而没有对应的状态序列，目标是学习HMM模型的参数λ=(π,A,B)\boldsymbol{\lambda}=(\boldsymbol{\pi},A,B)λ=(π,A,B)。将状态序列看作不可观测的隐数据Y\boldsymbol{Y}Y，HMM模型事实上是一个含有隐变量的概率模型：
p(X∣λ)=∑Yp(X∣Y,λ)p(Y∣λ)p(\boldsymbol{X}|\boldsymbol{\lambda})=\sum\limits_\boldsymbol{Y}p(\boldsymbol{X}|\boldsymbol{Y},\boldsymbol{\lambda})p(\boldsymbol{Y}|\boldsymbol{\lambda})p(X∣λ)=Y∑p(X∣Y,λ)p(Y∣λ)
这个参数可以由EM算法实现。

确定数据的对数似然函数
所有观测数据写成x={x1,x2,x3,...,xT}\boldsymbol{x}=\{x_1,x_2,x_3,...,x_T\}x={x1,x2,x3,...,xT}，所有的隐藏状态数据写成y={y1,y2,,...,yT}\boldsymbol{y}=\{y_1,y_2,,...,y_T\}y={y1,y2,,...,yT}，则完全数据是(x,y)=(x1,x2,x3,...,xT,y1,y2,,...,yT)(\boldsymbol{x},\boldsymbol{y})=(x_1,x_2,x_3,...,x_T,y_1,y_2,,...,y_T)(x,y)=(x1,x2,x3,...,xT,y1,y2,,...,yT)，完全数据的对数似然函数是log⁡p(x,y∣λ)\log{p(\boldsymbol{x},\boldsymbol{y}|\boldsymbol{\lambda})}logp(x,y∣λ)。
Em算法的E步：求Q函数Q(λ,λ‾)=Ey[log⁡p(y,y∣λ)∣x,λ‾]=∑ylog⁡p(x,y∣λ)p(x,y∣λ‾)\begin{aligned} Q(\boldsymbol{\lambda},\overline{\boldsymbol{\lambda}})&=E_\boldsymbol{y}[\log{p(\boldsymbol{y},\boldsymbol{y}|\boldsymbol{\lambda})}|\boldsymbol{x},\overline{\boldsymbol{\lambda}}]\\ &=\sum\limits_\boldsymbol{y}\log{p(\boldsymbol{x},\boldsymbol{y}|\boldsymbol{\lambda})}p(\boldsymbol{x},\boldsymbol{y}|\overline{\boldsymbol{\lambda}}) \end{aligned}Q(λ,λ)=Ey[logp(y,y∣λ)∣x,λ]=y∑logp(x,y∣λ)p(x,y∣λ)

其中λ‾\overline{\boldsymbol{\lambda}}λ是HMM模型参数的当前估计值，λ\boldsymbol{\lambda}λ是要极大化的HMM模型参数。
p(x,y∣λ)=πy1by1x1ay1y2by2x2⋅...⋅ayT−1yTbyTxTp(\boldsymbol{x},\boldsymbol{y}|\boldsymbol{\lambda})=\pi_{y_1}b_{y_1x_1}a_{y_1y_2}b_{y_2x_2}\cdot ...\cdot a_{y_{T-1}y_T}b_{y_Tx_T}p(x,y∣λ)=πy1by1x1ay1y2by2x2⋅...⋅ayT−1yTbyTxT
所以:
Q(λ,λ‾)=∑ylog⁡πy1p(x,y∣λ‾)+∑y[∑t=1T−1log⁡aytyt+1]p(x,y∣λ‾)+∑y[∑t=1Tlog⁡bytxt]p(x,y∣λ‾)Q(\boldsymbol{\lambda},\overline{\boldsymbol{\lambda}})=\sum\limits_\boldsymbol{y}\log{\pi_{y_1}p(\boldsymbol{x},\boldsymbol{y}|\overline{\boldsymbol{\lambda}}})+\sum\limits_\boldsymbol{y}[\sum\limits^{T-1}_{t=1}\log{a_{y_ty_{t+1}}]p(\boldsymbol{x},\boldsymbol{y}|\overline{\boldsymbol{\lambda}}})+\sum\limits_\boldsymbol{y}[\sum\limits_{t=1}^T\log{b_{y_tx_t}]p(\boldsymbol{x},\boldsymbol{y}|\overline{\boldsymbol{\lambda}}})Q(λ,λ)=y∑logπy1p(x,y∣λ)+y∑[t=1∑T−1logaytyt+1]p(x,y∣λ)+y∑[t=1∑Tlogbytxt]p(x,y∣λ)

EM算法的M步：极大化Q函数求模型的参数
第一项：
∑ylog⁡πy1p(x,y∣λ‾)=∑i=1Nlog⁡πip(x,y1=qi∣λ‾)\sum\limits_\boldsymbol{y}\log{\pi_{y_1}p(\boldsymbol{x},\boldsymbol{y}|\overline{\boldsymbol{\lambda}}})=\sum\limits_{i=1}^N\log{\pi_ip(\boldsymbol{x},y_1=q_i|\overline{\boldsymbol{\lambda}})}y∑logπy1p(x,y∣λ)=i=1∑Nlogπip(x,y1=qi∣λ)
注意到∑i=1Nπi=1\sum\limits^{N}_{i=1}{\pi_i}=1i=1∑Nπi=1，利用拉格朗日乘子法，写出拉格朗日函数：

∑i=1Nlog⁡πip(x,y1=qi∣λ‾)+γ(∑i=1N−1)\sum\limits_{i=1}^N\log{\pi_ip(\boldsymbol{x},y_1=q_i|\overline{\boldsymbol{\lambda}})}+\gamma(\sum\limits^{N}_{i=1}-1)i=1∑Nlogπip(x,y1=qi∣λ)+γ(i=1∑N−1)

求偏导并令其等于0：

∂∂πi[∑i=1Nlog⁡πip(x,y1=qi∣λ‾)+γ(∑i=1N−1)]=0\frac{\partial}{\partial \pi_i}[\sum\limits_{i=1}^N\log{\pi_ip(\boldsymbol{x},y_1=q_i|\overline{\boldsymbol{\lambda}})}+\gamma(\sum\limits^{N}_{i=1}-1)]=0∂πi∂[i=1∑Nlogπip(x,y1=qi∣λ)+γ(i=1∑N−1)]=0

得到：

p(x,y1=qi∣λ‾)+γπi=0p(\boldsymbol{x},y_1=q_i|\overline{\boldsymbol{\lambda}})+\gamma\pi_i=0p(x,y1=qi∣λ)+γπi=0

对i求和：
γ=−p(x∣λ‾)\gamma=-p(\boldsymbol{x}|\overline{\boldsymbol{\lambda}})γ=−p(x∣λ)

所以，得到πi\pi_iπi：

πi=p(x,y1=qi∣λ‾)p(x∣λ‾)\pi_i=\frac{p(\boldsymbol{x},y_1=q_i|\overline{\boldsymbol{\lambda}})}{p(\boldsymbol{x}|\overline{\boldsymbol{\lambda}})}πi=p(x∣λ)p(x,y1=qi∣λ)

第二项：
∑y[∑t=1T−1log⁡aytyt+1]p(x,y∣λ‾)=∑i=1N∑j=1N∑t=1T−1log⁡aijp(x,yt=qi,yt+1=qj∣λ‾)\sum\limits_\boldsymbol{y}[\sum\limits^{T-1}_{t=1}\log{a_{y_ty_{t+1}}]p(\boldsymbol{x},\boldsymbol{y}|\overline{\boldsymbol{\lambda}}})=\sum\limits_{i=1}^N\sum\limits_{j=1}^N\sum\limits_{t=1}^{T-1}\log{a_{ij}p(\boldsymbol{x},y_t=q_i,y_{t+1}=q_j|\overline{\boldsymbol{\lambda}})}y∑[t=1∑T−1logaytyt+1]p(x,y∣λ)=i=1∑Nj=1∑Nt=1∑T−1logaijp(x,yt=qi,yt+1=qj∣λ)

类似第一项，应用∑j=1N=1\sum\limits_{j=1}^N=1j=1∑N=1的拉格朗日乘子法可以求出：
aij=∑t=1T−1p(x,yt=qi,yt+1=qj∣λ‾)∑t=1T−1p(x,yt=qi∣λ‾)a_{ij}=\frac{\sum\limits_{t=1}^{T-1}p(\boldsymbol{x},y_t=q_i,y_{t+1}=q_j|\overline{\boldsymbol{\lambda}})}{\sum\limits_{t=1}^{T-1}p(\boldsymbol{x},y_t=q_i|\overline{\boldsymbol{\lambda}})}aij=t=1∑T−1p(x,yt=qi∣λ)t=1∑T−1p(x,yt=qi,yt+1=qj∣λ)

第三项：
∑y[∑t=1Tlog⁡bytxt]p(x,y∣λ‾)=∑j=1N∑t=1Tlog⁡bjxtp(x,yt=qj∣λ‾)\sum\limits_\boldsymbol{y}[\sum\limits_{t=1}^T\log{b_{y_tx_t}]p(\boldsymbol{x},\boldsymbol{y}|\overline{\boldsymbol{\lambda}}})=\sum\limits_{j=1}^N\sum\limits_{t=1}^{T}\log{b_{jx_t}p(\boldsymbol{x},y_t=q_j|\overline{\boldsymbol{\lambda}})}y∑[t=1∑Tlogbytxt]p(x,y∣λ)=j=1∑Nt=1∑Tlogbjxtp(x,yt=qj∣λ)

同样使用拉格朗日乘子法，求得：
bjk=∑t=1Tp(x,yt=qj∣λ‾)I(xt=vk)∑t=1Tp(x,yt=qj∣λ‾)b_{jk}=\frac{\sum\limits_{t=1}^Tp(\boldsymbol{x},y_t=q_j|\overline{\boldsymbol{\lambda}})I(x_t=v_k)}{\sum\limits_{t=1}^Tp(\boldsymbol{x},y_t=q_j|\overline{\boldsymbol{\lambda}})}bjk=t=1∑Tp(x,yt=qj∣λ)t=1∑Tp(x,yt=qj∣λ)I(xt=vk)

3.3 Baum-Welch算法

将EM算法的参数式子分别用前向和后向概率算出的γt(i),ξt(i,j)\gamma_t(i),\xi_t(i,j)γt(i),ξt(i,j)表示则:
aij=∑t=1T−1ξt(i,j)∑t=1T−1γt(i)a_{ij}=\frac{\sum\limits_{t=1}^{T-1}\xi_t(i,j)}{\sum\limits_{t=1}^{T-1}\gamma_t(i)}aij=t=1∑T−1γt(i)t=1∑T−1ξt(i,j)

bij=∑t=1,xt=vjTγt(i)∑t=1Tγt(i)b_{ij}=\frac{\sum\limits_{t=1,x_t=v_j}^T\gamma_t(i)}{\sum\limits_{t=1}^T\gamma_t(i)}bij=t=1∑Tγt(i)t=1,xt=vj∑Tγt(i)

πi=γ1(i)\pi_i=\gamma_1(i)πi=γ1(i)

4 序列预测问题及算法

序列预测问题就是已知模型参数λ=(π,A,B)\lambda=(\boldsymbol{\pi},A,B)λ=(π,A,B)，给定观测序列x\boldsymbol{x}x，求最有可能的状态序列y\boldsymbol{y}y，即求p(y∣x)p(\boldsymbol{y}|\boldsymbol{x})p(y∣x)的最大值。一般有两种解法：近似解法与维特比解法。

4.1 近似解法

近似算法的思想是，在每个时刻ttt选择该时刻最有可能出现的状态yt∗y_t^*yt∗，从而得到一个状态序列y=(yi∗,y2∗,...,yt∗)\boldsymbol{y}=(y_i^*,y_2^*,...,y_t^*)y=(yi∗,y2∗,...,yt∗)，将它作为预测的结果。
给定模型λ\boldsymbol{\lambda}λ和观测x\boldsymbol{x}x，在时刻ttt处于状态qiq_iqi的概率:
γt(i)==αt(i)βt(i)p(x∣λ)=αt(i)βt(i)∑j=1Nαt(i)βt(i)\gamma_t(i)==\frac{\alpha_t(i)\beta_t(i)}{p(\boldsymbol{x}|\boldsymbol{\lambda})}=\frac{\alpha_t(i)\beta_t(i)}{\sum\limits_{j=1}^N\alpha_t(i)\beta_t(i)}γt(i)==p(x∣λ)αt(i)βt(i)=j=1∑Nαt(i)βt(i)αt(i)βt(i)

在每一个时刻最有可能的状态yt∗y_t^*yt∗是：
yt∗=arg max⁡1≤i≤N[γt(i)],t=1,2,3...,Ty_t^*=\argmax_{1\leq i\leq N}[\gamma_t(i)],\quad t=1,2,3...,Tyt∗=1≤i≤Nargmax[γt(i)],t=1,2,3...,T
从而得到整个序列。
近似算法的优点是计算简单，其缺点是不能保证预测的状态序列整体是最有可能的状态序列，因为预测的状态序列可能有实际不发生的部分。事实上，上述方法得到的状态序列中有可能存在转移概率为0的相邻状态。尽管如此，近似算法仍然是有用的。

4.2 维特比算法（Viterbi algorithm）

维特比算法实际是用动态规划解HMM模型序列预测问题，用动态规划求解概率最大路径（最优路径），一条路径对应一个状态序列。
根据图论，假设最优路径为y∗\boldsymbol{y}^*y∗，其中从起点到时刻ttt的一段最优路径是(y1∗,y2∗,...,yt∗)(y_1^*,y_2^*,...,y_t^*)(y1∗,y2∗,...,yt∗)，则这部分路径对于后序最优路径（yt∗,yt+1∗,yt+2∗,...,yT∗y_t^*,y_{t+1}^*,y_{t+2}^*,...,y_T^*yt∗,yt+1∗,yt+2∗,...,yT∗）的选取来说一定也是最优的。可以使用反证法来证明：假设存在另一条局部路径(y1,,y2,,...,yt,)(y_1^,,y_2^,,...,y_t^,)(y1,,y2,,...,yt,)要优于(y1∗,y2∗,...,yt∗)(y_1^*,y_2^*,...,y_t^*)(y1∗,y2∗,...,yt∗)，那么它与（yt∗,yt+1∗,yt+2∗,...,yT∗y_t^*,y_{t+1}^*,y_{t+2}^*,...,y_T^*yt∗,yt+1∗,yt+2∗,...,yT∗）拼接起来蝴蝶刀另一条更优的全局最优路径，与定义矛盾。
根据HMM模型的第一个假设，yt+1y_{t+1}yt+1仅仅只与yty_tyt相关，所以网状图可以动态规划地搜索。定义：二维数组δt(i)\delta_t(i)δt(i)表示在时刻ttt以qjq_jqj结尾的所有局部路径的最大概率。从第一步推到第T步，每次递推都在上一次的N条局部路径中挑选，所以复杂度为O(TN)O(TN)O(TN)。为了得到路径，还需要定义一个二维数组ψt(i)\psi_t(i)ψt(i)，记录每个状态的前驱。

δt(i)=max⁡y1,y2,...,yt−1p(yt=qi,yt−1,...y1,xt,xt−1,...,x1∣λ)\delta_t(i)=\max_{y_1,y_2,...,y_{t-1}}p(y_t=q_i,y_{t-1},...y_1,x_t,x_{t-1},...,x_1|\boldsymbol{\lambda})δt(i)=y1,y2,...,yt−1maxp(yt=qi,yt−1,...y1,xt,xt−1,...,x1∣λ)

ψt(i)=arg⁡max⁡1≤j≤N[δt−1(j)aji]\psi_t(i)=\arg\max_{1\leq j \leq N}[\delta_{t-1}(j)a_{ji}]ψt(i)=arg1≤j≤Nmax[δt−1(j)aji]

维特比算法：

初始化，当t=1t=1t=1时，最优路径的备选由N个状态组成，它的前驱为空：
δ1(i)=πibix1,i=1,2,3,..,N\delta_1(i)=\pi_ib_{ix_1},\quad i=1,2,3,..,Nδ1(i)=πibix1,i=1,2,3,..,N

ψ1(i)=0,i=1,2,3,..,N\psi_1(i)=0,\quad i=1,2,3,..,Nψ1(i)=0,i=1,2,3,..,N

递推，当t≥2t\geq 2t≥2时，每条备选的路径像贪吃蛇一样吃入一个状态，长度增加一个单位，根据状态转移概率和发射概率计算花费。找出新的局部最优路径，更新两个数组。
δt(i)=max⁡1≤j≤N[δt−1(j)aji]bixt,i=1,2,3,..,N\delta_t(i)=\max_{1\leq j \leq N}[\delta_{t-1}(j)a_{ji}]b_{ix_t},\quad i=1,2,3,..,Nδt(i)=1≤j≤Nmax[δt−1(j)aji]bixt,i=1,2,3,..,N

ψt(i)=arg⁡max⁡1≤j≤N(δt(j)aji),i=1,2,3,..,N\psi_t(i)=\arg\max_{1\leq j \leq N}(\delta_t(j)a_{ji}),\quad i=1,2,3,..,Nψt(i)=arg1≤j≤Nmax(δt(j)aji),i=1,2,3,..,N

终止，找出最终时刻(δt(i)(\delta_t(i)(δt(i)数组中最大概率p∗p^*p∗，以及相应的结尾状态下表yT∗y_T^*yT∗
P∗=max⁡1≤j≤NδT(i)P^*=\max_{1\leq j \leq N}\delta_T(i)P∗=1≤j≤NmaxδT(i)

yT∗=arg⁡max⁡1≤j≤NδT(i)y_T^*=\arg\max_{1\leq j \leq N}\delta_T(i)yT∗=arg1≤j≤NmaxδT(i)

回溯，根据前驱数组ψt\psi_tψt回溯前驱状态，取得最优路径状态下标y∗=(y1∗,y2∗,...,yT∗)\boldsymbol{y}^*=(y_1^*,y_2^*,...,y_T^*)y∗=(y1∗,y2∗,...,yT∗)。
yt∗=ψt+1(yt+1∗),t=T−1,T−2,...,1y_t^*=\psi_{t+1}(y_{t+1}^*),\quad t=T-1,T-2,...,1yt∗=ψt+1(yt+1∗),t=T−1,T−2,...,1

举个例子：
上文中的盒子和球模型，状态集合为Q={1,2,3}Q=\{1,2,3\}Q={1,2,3}，观测集合为V={红，白}V=\{红，白\}V={红，白}，
A=[0.50.20.30.30.50.20.20.30.5]A=\begin{bmatrix} 0.5 & 0.2 & 0.3\\ 0.3 & 0.5 & 0.2\\ 0.2 & 0.3 & 0.5 \end{bmatrix}A=⎣⎡0.50.30.20.20.50.30.30.20.5⎦⎤ ， B=[0.50.50.40.60.70.3]B=\begin{bmatrix} 0.5 & 0.5 \\ 0.4 & 0.6\\ 0.7 & 0.3 \end{bmatrix}B=⎣⎡0.50.40.70.50.60.3⎦⎤ ， π=(0.2,0.4,0.4)T\pi=(0.2,0.4,0.4)^Tπ=(0.2,0.4,0.4)T
设T=3T=3T=3，x=(红，白，红)\boldsymbol{x}=(红，白，红)x=(红，白，红)，求最优状态序列，即最优路径y∗=(y1∗,y2∗,...,yT∗)\boldsymbol{y}^*=(y_1^*,y_2^*,...,y_T^*)y∗=(y1∗,y2∗,...,yT∗)。
解：

初始化，在t=1t=1t=1时，对每一个状态qi，i=1,2,3q_i，i=1,2,3qi，i=1,2,3，求状态为qiq_iqi观测为x1x_1x1为红的概率，记此概率为δ1(i)\delta_1(i)δ1(i)。
δ1(i)=πibix1=πibi红，i=1,2,3\delta_1(i)=\pi_ib_{ix_1}=\pi_ib_{i红}，\quad i=1,2,3δ1(i)=πibix1=πibi红，i=1,2,3

δ1(1)=0.2×0.5=0.1,δ1(2)=0.4×0.4=0.16,δ1(3)=0.4×0.7=0.28\delta_1(1)=0.2\times0.5=0.1,\delta_1(2)=0.4\times0.4=0.16,\delta_1(3)=0.4\times0.7=0.28δ1(1)=0.2×0.5=0.1,δ1(2)=0.4×0.4=0.16,δ1(3)=0.4×0.7=0.28

ψ1(i)=0,i=1,2,3\psi_1(i)=0,\quad i=1,2,3ψ1(i)=0,i=1,2,3
2. 在t=2t=2t=2时，对每个状态i，i=1,2,3i，i=1,2,3i，i=1,2,3:
δ2(i)=max⁡1≤j≤3[δ1(j)aji]bi白,i=1,2,3,..,N\delta_2(i)=\max_{1\leq j \leq 3}[\delta_{1}(j)a_{ji}]b_{i白},\quad i=1,2,3,..,Nδ2(i)=1≤j≤3max[δ1(j)aji]bi白,i=1,2,3,..,N

δ2(1)=max⁡1≤j≤3[δ1(j)aji]bi白=max⁡j{o.1×0.5,0.16×0.3,0.28×0.2}×0.5=0.028\begin{aligned} \delta_2(1)&=\max_{1\leq j \leq 3}[\delta_{1}(j)a_{ji}]b_{i白}\\ &=\max_j\{o.1\times0.5,0.16\times0.3,0.28\times0.2\}\times0.5\\ &=0.028 \end{aligned}δ2(1)=1≤j≤3max[δ1(j)aji]bi白=jmax{o.1×0.5,0.16×0.3,0.28×0.2}×0.5=0.028

ψ2(1)=3\psi_2(1)=3ψ2(1)=3

同理：

δ2(2)=0.0504,ψ2(3)=3\delta_2(2)=0.0504,\quad\psi_2(3)=3δ2(2)=0.0504,ψ2(3)=3

δ2(3)=0.042,ψ2(3)=3\delta_2(3)=0.042,\quad\psi_2(3)=3δ2(3)=0.042,ψ2(3)=3

同样，t=3t=3t=3：

δ3(1)=0.00756,ψ3(1)=2\delta_3(1)=0.00756,\quad\psi_3(1)=2δ3(1)=0.00756,ψ3(1)=2

δ3(2)=0.01008,ψ3(2)=2\delta_3(2)=0.01008,\quad\psi_3(2)=2δ3(2)=0.01008,ψ3(2)=2

δ3(3)=0.0417,ψ3(3)=3\delta_3(3)=0.0417,\quad\psi_3(3)=3δ3(3)=0.0417,ψ3(3)=3

最优路径概率:
P∗=max⁡1≤i≤3δ3(i)=0.0147P^*=\max_{1\leq i\leq 3}\delta_3(i)=0.0147P∗=1≤i≤3maxδ3(i)=0.0147

y3∗=arg max⁡1≤j≤3δ3(i)=3y_3^*=\argmax_{1\leq j \leq 3}\delta_3(i)=3y3∗=1≤j≤3argmaxδ3(i)=3

回溯：
t=2,y2∗=ψ3(y3∗)=ψ3(3)=3t=2,\quad y_2^*=\psi_3(y_3^*)=\psi_3(3)=3t=2,y2∗=ψ3(y3∗)=ψ3(3)=3

t=1,y1∗=ψ2(y2∗)=ψ2(3)=3t=1,\quad y_1^*=\psi_2(y_2^*)=\psi_2(3)=3t=1,y1∗=ψ2(y2∗)=ψ2(3)=3

所以，最优路径为y∗=(y1∗,y2∗,y3∗)=(3,3,3)\boldsymbol{y}^*=(y_1^*,y_2^*,y_3^*)=(3,3,3)y∗=(y1∗,y2∗,y3∗)=(3,3,3)。

5 hmmlearn使用

hmmlearn是一个实现了hmm的python库，安装很简单，使用pip install hmmlearn就行。
hmmlearn实现了三种HMM模型，分成两类：

针对观测状态是连续的：GaussianHMM和GMMHMM（广泛用于语音识别）
针对观测状态是离散的：MultinomialHMM，也就是上文中提到的。

对于MultinomialHMM的模型，使用比较简单，"startprob_"参数对应我们的初始状态概率向量π\piπ， "transmat_"对应我们的状态转移矩阵AAA, "emissionprob_"对应我们的发射概率矩阵BBB。

对于连续观测状态的HMM模型，GaussianHMM类假设观测状态符合高斯分布，而GMMHMM类则假设观测状态符合混合高斯分布。一般情况下我们使用GaussianHMM即高斯分布的观测状态即可。以下对于连续观测状态的HMM模型，我们只讨论GaussianHMM类。

在GaussianHMM类中，"startprob_"参数对应我们的隐藏状态初始分布π\piπ , "transmat_"对应我们的状态转移矩阵A, 比较特殊的是发射概率的表示方法，此时由于观测状态是连续值，我们无法像MultinomialHMM一样直接给出矩阵B。而是采用给出各个隐藏状态对应的观测状态高斯分布的概率密度函数的参数。

如果观测序列是一维的，则观测状态的概率密度函数是一维的普通高斯分布。如果观测序列是
N维的，则隐藏状态对应的观测状态的概率密度函数是N维高斯分布。高斯分布的概率密度函数参数可以用μ表示高斯分布的期望向量，Σ表示高斯分布的协方差矩阵。在GaussianHMM类中，“means”用来表示各个隐藏状态对应的高斯分布期望向量μ形成的矩阵，而“covars”用来表示各个隐藏状态对应的高斯分布协方差矩阵Σ形成的三维张量。

参考：深度剖析HMM

5.1 pythonceshi工具unnitest

测试工具unnitest非常容易使用，首先是建立一个继承自TestCase的测试类，然后通过覆盖setUp()完成相关初始化，最后通过覆盖tearDown()方法清除测试中产生的数据，为以后的TestCase留下一个干净的环境。我们需要在测试类中编写以test_开头的测试函数，unnitest会在测试中自动执行以test_开头的测试函数，unnitest的使用框架：

import unittest  class XXXXX(unittest.TestCase):  def setUp(self):  # 自行编写passif __name__ == '__main__':  unittest.main()

5.2 离散HMM测试

对离散HMM模型的测试代码：

# 计算平方误差
def s_error(A, B):return sqrt(np.sum((A-B)*(A-B)))/np.sum(B)class DiscreteHMM_Test(unittest.TestCase):def setUp(self):# 建立两个HMM，隐藏状态个数为4，X可能分布为10类n_state =4n_feature = 10X_length = 1000n_batch = 100 # 批量数目self.n_batch = n_batchself.X_length = X_lengthself.test_hmm = hmm.DiscreteHMM(n_state, n_feature)self.comp_hmm = ContrastHMM(n_state, n_feature)self.X, self.Z = self.comp_hmm.module.sample(self.X_length*10)self.test_hmm.train(self.X, self.Z)def test_train_batch(self):X = []Z = []for b in range(self.n_batch):b_X, b_Z = self.comp_hmm.module.sample(self.X_length)X.append(b_X)Z.append(b_Z)batch_hmm = hmm.DiscreteHMM(self.test_hmm.n_state, self.test_hmm.x_num)batch_hmm.train_batch(X, Z)# 判断概率参数是否接近# 初始概率判定没有通过！！！self.assertAlmostEqual(s_error(batch_hmm.start_prob, self.comp_hmm.module.startprob_), 0, 1)self.assertAlmostEqual(s_error(batch_hmm.transmat_prob, self.comp_hmm.module.transmat_), 0, 1)self.assertAlmostEqual(s_error(batch_hmm.emission_prob, self.comp_hmm.module.emissionprob_), 0, 1)def test_train(self):# 判断概率参数是否接近# 单批量的初始概率一定是不准的# self.assertAlmostEqual(s_error(self.test_hmm.start_prob, self.comp_hmm.module.startprob_), 0, 1)self.assertAlmostEqual(s_error(self.test_hmm.transmat_prob, self.comp_hmm.module.transmat_), 0, 1)self.assertAlmostEqual(s_error(self.test_hmm.emission_prob, self.comp_hmm.module.emissionprob_), 0, 1)def test_X_prob(self):X,_ = self.comp_hmm.module.sample(self.X_length)prob_test = self.test_hmm.X_prob(X)prob_comp = self.comp_hmm.module.score(X)self.assertAlmostEqual(s_error(prob_test, prob_comp), 0, 1)def test_predict(self):X, _ = self.comp_hmm.module.sample(self.X_length)prob_next = self.test_hmm.predict(X,np.random.randint(0,self.test_hmm.x_num-1))self.assertEqual(prob_next.shape,(self.test_hmm.n_state,))def test_decode(self):X,_ = self.comp_hmm.module.sample(self.X_length)test_decode = self.test_hmm.decode(X)_, comp_decode = self.comp_hmm.module.decode(X)self.assertAlmostEqual(s_error(test_decode, comp_decode), 0, 1)if __name__ == '__main__':unittest.main()

5.3 HMM测试

利用hmmlearn初始化一个高斯模型

class ContrastHMM():def __init__(self, n_state, n_feature):self.module = hmmlearn.hmm.GaussianHMM(n_components=n_state,covariance_type="full")# 初始概率self.module.startprob_ = np.random.random(n_state)self.module.startprob_ = self.module.startprob_ / np.sum(self.module.startprob_)# 转换概率self.module.transmat_ = np.random.random((n_state,n_state))self.module.transmat_ = self.module.transmat_ / np.repeat(np.sum(self.module.transmat_, 1),n_state).reshape((n_state,n_state))# 高斯发射概率self.module.means_ = np.random.random(size=(n_state,n_feature))*10self.module.covars_ = .5 * np.tile(np.identity(n_feature), (n_state, 1, 1))

利用hmmlearn初始化一个离散HMM模型:

class ContrastHMM():def __init__(self, n_state, n_feature):self.module = hmmlearn.hmm.MultinomialHMM(n_components=n_state)# 初始概率self.module.startprob_ = np.random.random(n_state)self.module.startprob_ = self.module.startprob_ / np.sum(self.module.startprob_)# print self.module.startprob_# 转换概率self.module.transmat_ = np.random.random((n_state,n_state))self.module.transmat_ = self.module.transmat_ / np.repeat(np.sum(self.module.transmat_, 1),n_state).reshape((n_state,n_state))# print self.module.transmat_# 发射概率self.module.emissionprob_ = np.random.random(size=(n_state,n_feature))self.module.emissionprob_ = self.module.emissionprob_ / np.repeat(np.sum(self.module.emissionprob_, 1),n_feature).reshape((n_state,n_feature))# print self.module.emissionprob_

代码详情请见我的github，请大家帮忙点个star。

一文读懂NLP之隐马尔科夫模型（HMM）详解加python实现相关推荐

隐马尔科夫模型HMM详解(1)
目录隐马尔科夫模型基本概念隐马尔科夫模型的三个基本问题概率计算预测算法-Viterbi算法 HMM学习算法参考下篇文章代码地址:https://gitee.com/liangcd/speec ...
隐马尔科夫模型HMM之Baum-Welch算法Python代码实现
☕️ 本文系列文章汇总: (1)HMM开篇:基本概念和几个要素 (2)HMM计算问题:前后向算法代码实现 (3)HMM学习问题:Baum-Welch算法 (4) HMM预测问题:维特比算法本篇算 ...
【NLP】用于语音识别、分词的隐马尔科夫模型HMM
大家好,今天介绍自然语言处理中经典的隐马尔科夫模型(HMM).HMM早期在语音识别.分词等序列标注问题中有着广泛的应用. 了解HMM的基础原理以及应用,对于了解NLP处理问题的基本思想和技术发展脉络有 ...
隐马尔科夫模型HMM自学（3）
Viterbi Algorithm 本来想明天再把后面的部分写好,可是睡觉今天是节日呢?一时情不自禁就有打开电脑.......... 找到可能性最大的隐含状态序列崔晓源翻译多数情况下,我们都希望 ...
隐马尔科夫模型 HMM 与语音识别 speech recognition (1)：名词解释
0.引言想在 CSDN 上看一下隐马尔科夫模型,简称HMM(Hidden Markov Model)的例子,找了几篇博文,却发现大部分都是转载的,转载的还没有出处,文中的表述与逻辑也看的人晕头转向, ...
隐马尔科夫模型(HMM)择时应用的量化策略
HMM模型隐马尔科夫模型(HMM)择时应用的量化策略. 仅为研究学习使用, 不作为任何投资策略建议. 文章内容从各处整理汇总而成, 感谢各位大神分享. 具体策略代码均调试通过. 一.从大奖章讲起 ...
隐马尔科夫模型HMM自学（2）
HMM 定义崔晓源翻译 HMM是一个三元组 (,A,B). the vector of the initial state probabilities; the state transitio ...
隐马尔科夫模型HMM自学（1）
介绍崔晓源翻译我们通常都习惯寻找一个事物在一段时间里的变化规律.在很多领域我们都希望找到这个规律,比如计算机中的指令顺序,句子中的词顺序和语音中的词顺序等等.一个最适用的例子就是天气的预测. 首 ...
python地图匹配_基于隐马尔科夫模型(HMM)的地图匹配(Map-Matching)算法
1. 摘要本篇博客简单介绍下用隐马尔科夫模型(Hidden Markov Model, HMM)来解决地图匹配(Map-Matching)问题.转载请注明网址. 2. Map-Matching(MM ...

一文读懂NLP之隐马尔科夫模型（HMM）详解加python实现

一文读懂NLP之隐马尔科夫模型（HMM）详解加python实现

1 隐马尔科夫模型

1.1 HMM解决的问题

1.2 HMM模型的定义

1.2.1HMM的两个假设

1.2.2 HMM模型

1.3 HMM模型的三个基本问题

2 概率计算问题及算法

2.1 直接计算法

2.2 前向算法

2.3 后向算法

2.4 一些概率与期望值的计算

3 模型训练问题及算法

3.1 监督学习——最大似然估计

3.2 非监督学习——EM算法

3.3 Baum-Welch算法

4 序列预测问题及算法

4.1 近似解法

4.2 维特比算法（Viterbi algorithm）

5 hmmlearn使用

5.1 pythonceshi工具unnitest

5.2 离散HMM测试

5.3 HMM测试

一文读懂NLP之隐马尔科夫模型（HMM）详解加python实现相关推荐

最新文章

热门文章