隐马尔可夫模型（Hidden Markov Model, HMM）是可用于标注问题的模型，描述由隐藏的马尔可夫链随机生成观测序列的过程，属于生成模型。马尔可夫链不懂的可以把本科的《概率论与数理统计》找回来看一下，很简单的，就是离散状态之间的转换。下面直接定义基本概念，为后面的算法做准备。
基本概念
变量定义
HMM各个时期会处于各种状态，设$Q$是所有可能状态的集合；每个状态可以产生各种观测，设$V$是所有可能观测的集合。$Q,V$的定义如下：
$Q={q_1,q_2,dots, q_N},;;V={v_1,v_2,dots,v_M}$
其中$N$是可能的状态数，$M$是可能的观测数。注意，这里的状态和观测不一定是数字，也可以是各种具体的对象。
然后定义$S$为长度为$T$的模型状态序列，定义$O$为对应的观测序列。
$S=(s_1,s_2,dots,s_T),;;O=(o_1,o_2,dots,o_T)$
这两个序列都是整数列表，其中每个整数索引对应的状态和观测。比如：$q_{s_2}$表示模型在时刻2时的状态，$v_{o_5}$表示模型在时刻5时的观测，而$s_2,o_5$并不代表那时模型的状态和观测本身（这里先说清楚，不然后面容易混淆）。
那么这些状态是如何初始化并相互转换，观测又是如何进行的呢？下面定义三个分布：
1、初始概率分布。时间开始时，各个状态出现的概率。定义为向量$pi$：
$egin{aligned}&pi = [pi_1,pi_2,dots,pi_N]^T\ &pi_i=P(s_1=i),;;i = 1,2,dots,N end{aligned}$
2、状态转移分布。上一时刻到下一时刻不同状态之间转换的概率。定义为矩阵$A$：
$egin{aligned} &A = [a_{ij}]_{N imes N}\& a_{ij} = P(s_{t+1}=j|s_t=i)end{aligned}$
3、观测概率分布。定义某个状态下各种观测出现的概率。定义为矩阵$B$：
$egin{aligned} &B = [b_{ij}]_{N imes M}\& b_{ij} = P(o_t=j|s_t=i) end{aligned}$
隐马尔可夫模型由以上三个分布决定，因此可以用一个三元符号表示：
$lambda = (A,B,pi)$
HMM基本假设
从定义可知，HMM做了两个基本假设：
1、齐次马尔可夫性假设。任意时刻的状态只依赖于前一时刻的状态，与其它时刻的状态及观测无关：
$P(s_t|s_{t-1},o_{t-1},dots,s_1,o_1) = P(s_t|s_{t-1})$
注意，以上条件概率中将除$s_{t-1}$以外的条件去掉，是因为$s_{t-1}$已知，并且没有之后时刻的状态或观测作为条件。如果$s_{t-1}$未知，则可以去掉$t$时刻之前的条件中，离$t$最近的$t^-$之前的状态和观测（包含$t^-$时刻的观测）。如：
$P(s_t|s_{t-2},o_{t-2},s_{t-3},s_{t-4}) = P(s_t|s_{t-2})$
另外，假如有之后时刻的状态，计算的概率就是后验概率了，那么之后时刻的状态作为条件来说也不能去掉。但是可以去掉$t$时刻之后的条件中，离$t$最近的$t^+$之后的状态和观测（包含$t^+$时刻的观测）如：
$egin{gather}P(s_t|s_{t-2},s_{t-1},o_{t-1},o_{t+1},s_{t+2},o_{t+2},o_{t+3}) = P(s_t|s_{t-1},o_{t+1},s_{t+2})label{}end{gather}$
总之，就是近的状态已知，远的状态对于计算条件概率来说就没有意义了。
2、观测独立性假设。任意时刻的观测只依赖于此刻的状态，与其它无关：
$P(o_t|s_t,o_{t-1},dots,s_1,o_1) = P(o_t|s_t)$
这个条件概率和上面也一样，也是近的状态已知，远处的状态作为条件就无意义。
以上两个假设和马尔可夫链的“链”相符，状态条件的已知就好像把链条对应结点的外侧链条砍断（靠近$s_t$的是内侧），我觉得把这种性质称作“就近原则”也不错（留意后面很多计算都用到）。下面是$(1)$式的示意图：
其中黄色结点表示已知条件，圈在虚线中的表示会影响$(1)$式概率值的状态与观测。
三个基本问题
弄懂HMM状态之间的转换规则后，现在提出HMM的三个基本问题：
1、概率计算。给定模型$lambda$和观测序列$O$，计算在该模型下观测序列$O$出现的概率$P(O|lambda)$ 。
2、学习。已知观测序列$O$，估计模型$lambda$的参数，使得在该模型下观测到这个序列的概率$P(O|lambda)$最大。
3、预测，或是解码。已知模型$lambda$，给定观测序列$O$，求与之对应的状态序列$S$，使得概率$P(S|O,lambda)$最大。
这三个基本问题分别阐述了HMM的基本推算方式和主要用途。下面依次介绍这三个基本问题的解法。
概率计算
对于给定的模型参数$lambda =(A,B,pi)$，计算观测序列$O= (o_1,o_2,dots,o_T)$出现的概率，最简单的就是把所有可能的状态序列的概率都计算出来，然后选择最大的那个。但是这个方法计算复杂度是极大的，高达$O(TN^T)$阶，所以不可行。下面介绍两种计算可行的算法：前向和后向算法。
前向算法
定义到达时刻$t$，观测序列为$(o_1,o_2,dots,o_t)$，且此时状态为$q_i$的概率
$alpha_t(i) = P(o_1,o_2,dots,o_t,s_t=i|lambda)$
为前向概率。前向算法就是用$alpha_t(i)$前向递推计算观测序列的概率。下面是算法的推算流程：
1、计算初值。即当$t=1$时所有状态的$alpha$，一共$N$个：
$egin{aligned}alpha_1(i)=P(o_1,s_1=i|lambda)=pi_ib_{io_1},;;i=1,2,dots,Nend{aligned}$
2、递推。对$t=2,3,dots,T$，计算所有状态的$alpha$，每个$t$同样都要计算$N$个：
$ egin{aligned} &alpha_t(i) \ =& P(o_1,o_2,...,o_t,s_t=i|lambda)\ =& sumlimits_{j=1}^NP(o_1,o_2,...,o_{t},s_{t-1}=j,s_t=i|lambda) \ =& sumlimits_{j=1}^N P(o_t|o_1,o_2,...,o_{t-1},s_{t-1}=j,s_t=i,lambda)P(o_1,o_2,...,o_{t-1},s_{t-1}=j,s_t=i|lambda) \ =& sumlimits_{j=1}^N P(o_t|s_t=i,lambda)P(o_1,o_2,...,o_{t-1},s_{t-1}=j|lambda)P(s_t=i|o_1,o_2,...,o_{t-1},s_{t-1}=j,lambda) \ =& P(o_t|s_t=i,lambda)sumlimits_{j=1}^N P(o_1,o_2,...,o_{t-1},s_{t-1}=j|lambda)P(s_t=i|s_{t-1}=j,lambda) \ =& b_{io_i}sumlimits_{j=1}^Nalpha_{t-1}(j)a_{ji},;;i=1,2,...,N end{aligned} $
3、计算最后的概率：
$ egin{aligned} P(O|lambda) =sumlimits_{i=1}^NP(o_1,o_2,...,o_T,s_T=i|lambda) =sumlimits_{i=1}^Nalpha_T(i) end{aligned} $
下面是书上的计算图，便于理解：
后向算法
定义在$t$时刻的状态为$q_i$的条件下，之后观测序列为$(o_{t+1},o_{t+2},dots,o_T)$的条件概率
$eta_t(i) = P(o_{t+1},o_{t+2},dots,o_T|s_t=i,lambda)$
为后向概率。与前向算法类似，后向算法是用$eta_t(i)$后向推算概率。下面是递推流程：
1、定义初值：
$egin{aligned}&eta_T(i)=1,;;i=1,2,...,N\end{aligned}$
2、对$t=T-1,T-2,...,1$，计算：
$ egin{aligned} &P(o_{t+1},...,o_T|s_t=i,lambda)\ =&sumlimits_{j=1}^NP(o_{t+1},...,o_T,s_{t+1}=j|s_t=i,lambda)\ =&sumlimits_{j=1}^NP(o_{t+1}|o_{t+2},...,o_T,s_{t+1}=j,s_t=i,lambda)P(o_{t+2},...,o_T,s_{t+1}=j|s_t=i,lambda)\ =&sumlimits_{j=1}^NP(o_{t+1}|s_{t+1}=j,lambda)P(o_{t+2},...,o_T|,s_{t+1}=j,s_t=i,lambda)P(s_{t+1}=j|s_t=i,lambda)\ =&sumlimits_{j=1}^NP(o_{t+1}|s_{t+1}=j,lambda)P(o_{t+2},...,o_T|,s_{t+1}=j,lambda)P(s_{t+1}=j|s_t=i,lambda)\ =&sumlimits_{j=1}^N a_{ij}b_{jo_{t+1}}eta_{t+1}(j),;;i=1,2,...,N\ end{aligned} $
3、计算结果：
$ egin{aligned} P(O|lambda) = sumlimits_{i=1}^Npi_ib_{io_1}eta_1(i) end{aligned} $
利用前后向算法计算某些概率
我们可以用$alpha$和$eta$计算某些概率。
1、给定$lambda$，计算以观测$O$为条件，$t$时刻状态为$q_i$的条件概率：
$ egin{aligned} &P(s_t=i|O,lambda)\ =&frac{ P(o_1,...,o_T,s_t=i|lambda) }{P(O|lambda)}\ =&frac{ P(o_1,...,o_t,s_t=i|lambda)P(o_{t+1},...,o_T|o_1,...,o_t,s_t=i,lambda) }{P(O|lambda)}\ =&frac{ P(o_1,...,o_t,s_t=i|lambda)P(o_{t+1},...,o_T| s_t=i,lambda) }{sumlimits_{j=1}^NP(o_1,...,o_T,s_t=j|lambda)}\ =&frac{ alpha_t(i)eta_t(i) }{sumlimits_{j=1}^Nalpha_t(j)eta_t(j)}\ end{aligned} $
2、给定$lambda$，计算以观测$O$为条件，$t$时刻状态为$q_i$，且$t+1$时刻状态为$q_j$的条件概率：
$egin{aligned}&P(s_t=i,s_{t+1}=j|O,lambda)\=&frac{alpha_t(i)a_{ij}b_{jo_{t+1}}eta_{t+1}(j)}{P(O|lambda)}\=&frac{alpha_t(i)a_{ij}b_{jo_{t+1}}eta_{t+1}(j)}{sumlimits_{k=1}^Nsumlimits_{l=1}^NP(s_t=k,s_{t+1}=l,O|lambda)}\=&frac{alpha_t(i)a_{ij}b_{jo_{t+1}}eta_{t+1}(j)}{sumlimits_{k=1}^Nsumlimits_{l=1}^Nalpha_t(k)a_{kl}b_{lo_{t+1}}eta_{t+1}(l)}\end{aligned}$
具体推算过程和1类似，也用到了“链”的性质把多余的条件去掉。
3、给定$lambda$，以观测$O$为条件，求某个状态或某两个连续状态出现的期望（就是出现次数的期望）。直接对所有时刻的以上条件概率进行求和即可，不再详写。
学习
如果给你观测序列和对应的状态序列，那么直接用各种转移出现的频率来估计$lambda=(A,B,pi)$即可，这很简单，就不讲了。而如果只给你观测序列，正如上面基本问题中定义的那样，就不能直接计算了，需要使用EM算法（点击链接）来迭代计算，就是所谓的Baum-Welch算法（实际上就是EM算法的应用）。
Q函数
因为只有观测序列，其未知的状态序列就可以看做隐变量。假设用于训练的观测序列集合为：
$mathbb{O} = {O_1,O_2,dots,O_{|mathbb{O}|}}$
其中任意观测序列$O_i$的长度都为$T$。
然后假定隐变量状态序列的所有可能性集合为：
$mathbb{S} = {left.S=(s_1,s_2,dots,s_T) ight|s_t = 1,2,dots,N,;;t = 1,2,dots,T}$
于是可以定义第k次迭代的$Q$函数：
$displaystyle Q(lambda,lambda^k) = sumlimits_{mathbb{O}}sumlimits_{mathbb{S}}P(S|O,lambda^k)log P(S,O|lambda)$
这里$Q$函数的定义与《统计学习方法》上定义的不同，为了更严谨，我把所有观测数据都算上了。
代入待优化参数，然后进行分解：
$ egin{aligned} Qleft(lambda,lambda^k ight) &= sumlimits_{mathbb{O}}sumlimits_{mathbb{S}}Pleft(S|O,lambda^k ight) log left(pi_{s_1}b_{s_1o_1}a_{s_1s_2}b_{s_2o_2}dots a_{s_{T-1}s_T}b_{s_To_T} ight)\ &= sumlimits_{mathbb{O}}sumlimits_{mathbb{S}}Pleft(S|O,lambda^k ight) left(log pi_{s_1}+sumlimits_{t = 1}^{T-1}log a_{s_ts_{t+1}}+sumlimits_{t = 1}^Tlog b_{s_to_t} ight) end{aligned} $
可以看到待优化式子分解为三个不同待优化变量的求和，因此可以将它们分开分别进行优化。
优化π
首先对$pi$的优化式进行整理，将每个优化变量$pi_i$的系数合并：
$ egin{gather} egin{aligned} & sumlimits_{mathbb{O}}sumlimits_{mathbb{S}}Pleft(S|O,lambda^k ight) log pi_{s_1}\ =& sumlimits_{mathbb{O}}sumlimits_{i=1}^N sumlimits_{{S|Sin mathbb{S},s_1=i}}P(S|O,lambda^k) log pi_i\ =& sumlimits_{mathbb{O}}sumlimits_{i=1}^Nsumlimits_{1le s_2,dots,s_T le N}Pleft(s_1=i,s_2,s_3,dots,s_T|O,lambda^k ight) log pi_i\ =& sumlimits_{mathbb{O}}sumlimits_{i=1}^N Pleft(s_1=i|O,lambda^k ight) log pi_i\ end{aligned}label{} end{gather} $
上式第一个等号成立是因为$sumlimits_{mathbb{S}}$和$sumlimits_{i=1}^N sumlimits_{{S|Sin mathbb{S},s_1=i}}$是等价的。
注意到等式约束$sumlimits_{i=1}^Npi_i=1$为仿射函数，并且上式为凸函数，因此这是一个带等式约束的凸优化（满足KKT条件（点击链接）即可）。写出拉格朗日函数：
$ egin{aligned} L(pi,gamma) = sumlimits_{mathbb{O}}sumlimits_{i=1}^N Pleft(s_1=i|O,lambda^k ight) log pi_i + gammaleft(sumlimits_{i=1}^Npi_i-1 ight) end{aligned} $
对每个$pi_i$求导并等于0、化简，再联立等式约束，得KKT条件：
$ left{ egin{aligned} & frac{partial L}{partialpi_i}=sumlimits_{mathbb{O}}frac{1}{pi_i} Pleft(s_1=i|O,lambda^k ight) +gamma=0,;;i=1,2,dots,N\ &sumlimits_{i=1}^Npi_i=1 end{aligned} ight. $
将分母的$pi_i$乘过去然后求和计算$gamma$：
$displaystylesumlimits_{mathbb{O}} sumlimits_{i=1}^N Pleft(s_1=i|O,lambda^k ight) +sumlimits_{i=1}^Ngammapi_i=0$
$|mathbb{O}| +gamma=0$
$gamma=-|mathbb{O}| $
其中$|mathbb{O}|$是用于训练的观测序列数。再计算$pi_i$：
$ egin{gather}egin{aligned} & pi_i = frac{1}{|mathbb{O}|}sumlimits_{mathbb{O}} Pleft(s_1=i|O,lambda^k ight) end{aligned} label{}end{gather} $
优化a
$a$的优化式：
$ egin{aligned} & sumlimits_{mathbb{O}}sumlimits_{mathbb{S}}Pleft(S|O,lambda^k ight)sumlimits_{t = 1}^{T-1}log a_{s_to_{t+1}}\ end{aligned} $
和$pi_i$一样，也要先将每个$a_{ij}$的系数合并。为了容易理解，先对单个的$a_{s_ts_{t+1}}$进行整理：
$ egin{aligned} & sumlimits_{mathbb{O}}sumlimits_{mathbb{S}}Pleft(S|O,lambda^k ight)log a_{s_to_{t+1}}\ =& sumlimits_{mathbb{O}}sumlimits_{i=1}^Nsumlimits_{j=1}^Nsumlimits_{{S|Sinmathbb{S},s_t=i,s_{t+1}=j}}Pleft(S|O,lambda^k ight)log a_{ij}\ =& sumlimits_{mathbb{O}}sumlimits_{i=1}^Nsumlimits_{j=1}^NPleft(s_t=i,s_{t+1}=j|O,lambda^k ight)log a_{ij}\ end{aligned} $
再加上$t$的求和，得到整理后的$a$的优化式：
$ egin{aligned} & sumlimits_{mathbb{O}}sumlimits_{i=1}^Nsumlimits_{j=1}^Nsumlimits_{t=1}^{T-1}Pleft(s_t=i,s_{t+1}=j|O,lambda^k ight)log a_{ij}\ end{aligned} $
类似于优化$pi$，列出具有$N$个等式约束$sumlimits_{j=1}^Na_{ij}=1$的拉格朗日函数：
$ egin{aligned} L(a,gamma) =sumlimits_{mathbb{O}}sumlimits_{i=1}^Nsumlimits_{j=1}^Nsumlimits_{t=1}^{T-1}Pleft(s_t=i,s_{t+1}=j|O,lambda^k ight)log a_{ij} +sumlimits_{i=1}^Ngamma_ileft(sumlimits_{j=1}^Na_{ij}-1 ight) end{aligned} $
列出KKT条件：
$ left{ egin{aligned} &frac{partial L}{partial a_{ij}} =frac{1}{a_{ij}}sumlimits_{mathbb{O}}sumlimits_{t=1}^{T-1}Pleft(s_t=i,s_{t+1}=j|O,lambda^k ight) +gamma_i =0,;;i,j=1,2,dots,N\ &sumlimits_{j=1}^Na_{ij}=1,;;i=1,2,dots,N end{aligned} ight. $
求和计算得$gamma$：
$displaystyle gamma_i = -sumlimits_{mathbb{O}}sumlimits_{t=1}^{T-1}Pleft(s_t=i|O,lambda^k ight) $
再算$a_{ij}$：
$ egin{gather}egin{aligned} a_{ij} =frac{sumlimits_{mathbb{O}}sumlimits_{t=1}^{T-1}Pleft(s_t=i,s_{t+1}=j|O,lambda^k ight)}{sumlimits_{mathbb{O}}sumlimits_{t=1}^{T-1}Pleft(s_t=i|O,lambda^k ight)}\ end{aligned} label{}end{gather} $
优化b
同样先将每个$b_{ij}$的系数进行合并：
$ egin{gather} egin{aligned} &sumlimits_{mathbb{O}}sumlimits_{mathbb{S}}Pleft(S|O,lambda^k ight) sumlimits_{t = 1}^Tlog b_{s_to_t}\ =&sumlimits_{mathbb{O}}sumlimits_{mathbb{S}}sumlimits_{t = 1}^TPleft(S|O,lambda^k ight) log b_{s_to_t}\ =&sumlimits_{j=1}^Msumlimits_{i=1}^Nsumlimits_{t = 1}^T sumlimits_{{O|Oin mathbb{O},o_t=j}} sumlimits_{{S|Sin mathbb{S},s_t=i}} Pleft(S|O,lambda^k ight) log b_{ij}\ =&sumlimits_{j=1}^Msumlimits_{i=1}^Nsumlimits_{t = 1}^T sumlimits_{{O|Oin mathbb{O},o_t=j}} Pleft(s_t=i|O,lambda^k ight) log b_{ij}\ end{aligned} label{} end{gather} $
列出包含约束的拉格朗日函数：
$ egin{gather} egin{aligned} L(b,gamma)=&sumlimits_{j=1}^Msumlimits_{i=1}^Nsumlimits_{t = 1}^T sumlimits_{{O|Oin mathbb{O},o_t=j}} Pleft(s_t=i|O,lambda^k ight) log b_{ij} +sumlimits_{i=1}^Ngamma_ileft(sumlimits_{j=1}^Mb_{ij}-1 ight) \ end{aligned}label{} end{gather} $
这里说明一下，$b$和$pi,a$不同，$b_{ij}$是否为0是与训练数据集$mathbb{O}$直接相关的。比如在$mathbb{O}$没有出现观测$v_3$，那我们可以直接将所有$b_{i3}$估计为0，这是符合常理的。而在$(6)$式中，因为不存在$o_t=3$的$O$，所以$log b_{i3}$也不存在，于是对$log b_{i3}$求导也没法求，就不能参与优化了，直接估计为0即可。另外，$pi,a$是与隐变量状态序列相关的，而隐变量我们是要考虑所有可能的情况的，所以所有$pi_i,a_{ij}$都可以参与优化。
下面对非0的$b_{ij}$进行求导并令其为0，联立约束列出KKT条件：
$ left{ egin{aligned} &frac{partial L}{partial b_{ij}} =frac{1}{b_{ij}}sumlimits_{t = 1}^T sumlimits_{{O|Oin mathbb{O},o_t=j}} Pleft(s_t=i|O,lambda^k ight) + gamma_i =0,;; i=1,2,dots,N,;;j=1,2,dots,M \ &sumlimits_{j=1}^Mb_{ij}= 1,;;i=1,2,dots,N end{aligned} ight. $
计算$gamma$：
$ displaystylesumlimits_{j = 1}^M sumlimits_{t = 1}^T sumlimits_{{O|Oin mathbb{O},o_t=j}} Pleft(s_t=i|O,lambda^k ight) + gamma_i = 0 $
$ displaystyle gamma_i = -sumlimits_{mathbb{O}} sumlimits_{t = 1}^T Pleft(s_t=i|O,lambda^k ight) $
计算优化$b$：
$ displaystyle b_{ij}= frac{sumlimits_{t = 1}^T sumlimits_{{O|Oin mathbb{O},o_t=j}} Pleft(s_t=i|O,lambda^k ight) } {sumlimits_{mathbb{O}} sumlimits_{t = 1}^T Pleft(s_t=i|O,lambda^k ight)} $
算上训练数据$mathbb{O}$中不包含某观测$v_j$，使得$b_{ij}=0$的情况（其实上式已经包含了，但是我不确定是看做0还是看做不存在，为了严谨还是显式写出来好了）：
egin{gather} b_{ij}=left{ egin{aligned} &frac{sumlimits_{t = 1}^T sumlimits_{{O|Oin mathbb{O},o_t=j}} Pleft(s_t=i|O,lambda^k ight) } {sumlimits_{mathbb{O}} sumlimits_{t = 1}^T Pleft(s_t=i|O,lambda^k ight)}, &exists Oin mathbb{O}, jin O \ &0,& forall Oin mathbb{O}, j otin O end{aligned} ight. label{}end{gather}
迭代
综上，使用EM算法迭代计算HMM参数的步骤如下：
0、初始化$lambda^0 = (pi^0,A^0,B^0)$，各取为符合约束的随机值。
1、第$k$次迭代，在$lambda^{k-1}$的基础上分别使用$(3),(4),(7)$式对$pi,a,b$进行优化，得到$lambda^k$。
2、判断终止条件，满足即完成迭代，否则回到第一步执行下一次迭代。
预测
预测问题最直观的方式，就是直接选择每个时刻最有可能出现的状态$s_t^*$，从而预测整个状态序列，即：
$ egin{aligned} &s_t^* = argmaxlimits_{1le i le N} frac{alpha_t(i)eta_t(i)}{sumlimits_{j=1}^Nalpha_t(j)eta_t(j)},;;t=1,2,...,T\ &S = (s_1^*,s_2^*,...,s_T^*) end{aligned} $
但这并不是整体最优的，而且它并不能剔除因为状态转移概率为0而不可能出现的状态序列。为了整体最优并且符合状态转移概率矩阵，我们通常使用维特比算法来预测状态序列。
维特比算法
维特比算法实际上就是动态规划，在递推的过程中保持步步最优，从而最终达到全局最优的目的。
定义在$t$时刻以内，且$t$时刻状态为$q_i$的所有状态序列的出现概率的最大值为：
$ egin{gather} egin{aligned} delta_t(i) = max_{s_1,...,s_{t-1}}P(s_t=i,s_{t-1},...,s_1,o_t,...,o_1|lambda),;;i=1,2,...,N end{aligned} label{}end{gather} $
由定义可以获得$delta$的递推公式：
$ egin{aligned} &delta_{t+1}(i) \ =& max_{s_1,...,s_t}P(s_{t+1}=i,s_t,...,s_1,o_{t+1},...,o_1|lambda)\ =& max_{s_t}max_{s_1,...,s_{t-1}}P(s_{t+1}=i,s_t,...,s_1,o_{t+1},...,o_1|lambda)\ =& max_{1le jle N}max_{s_1,...,s_{t-1}}P(s_{t+1}=i,s_t=j,s_{t-1},...,s_1,o_{t+1},...,o_1|lambda)\ =& max_{1le jle N}max_{s_1,...,s_{t-1}}left[P(o_{t+1}|s_{t+1}=i,lambda)P(s_{t+1}=i,s_t=j,s_{t-1},...,s_1,o_t,...,o_1|lambda) ight]\ =& P(o_{t+1}|s_{t+1}=i,lambda)max_{1le jle N}max_{s_1,...,s_{t-1}}left[P(s_{t+1}=i|s_t=j,lambda)P(s_t=j,s_{t-1},...,s_1,o_t,...,o_1|lambda) ight]\ =& P(o_{t+1}|s_{t+1}=i,lambda)max_{1le jle N}left[P(s_{t+1}=i|s_t=j,lambda)max_{s_1,...,s_{t-1}}P(s_t=j,s_{t-1},...,s_1,o_t,...,o_1|lambda) ight]\ =& b_{io_{t+1}}max_{1le jle N}a_{ji}delta_t(j) end{aligned} $
以上算的是概率，并没有获取状态编号。所以再定义：
$ egin{gather} egin{aligned} Psi_t(i) = argmax_{1le jle N} left[delta_{t-1}(j)a_{ji} ight],;;i=1,2,...,N end{aligned} label{}end{gather} $
用来回溯获取状态序列。
维特比算法用已知观测序列预测状态序列过程如下：
1、初始化。直接计算$delta_1(i)$。
2、对$t=2,3,...,T$，使用递推公式计算$delta_t(i)$，使用$(9)$式计算$Psi_t(i)$。
3、使用$Psi_t(i)$回溯，获取最大概率状态序列。
参考资料
李航《统计学习方法》