Hidden Markov Model, HMM 隐马尔可夫模型，是一种描述隐性变量(状态)和显性变量(观测状态)之间关系的模型。该模型遵循两个假设，隐性状态i只取决于前一个隐性状态i-1，而与其他先前的隐形状态无关。观测状态也只取决于当前的隐形状态。因此我们常常将隐马尔科夫模型表现为一种如下图所示链式的模型。

其中，$x_t$代表某一时刻的隐形状态链，其N个状态取值集合为{s1,s2,s3,...,N}\{s_1,s_2,s_3,...,_N\}{s1​,s2​,s3​,...,N​}。$y_t$表示为对应的该时刻的显性状态(观测状态)，其M个状态取值集合为{o1,o2,o3,...,ok,...,oM}\{o_1,o_2,o_3,...,o_k,...,o_M\}{o1​,o2​,o3​,...,ok​,...,oM​}。隐马尔科夫模型$\theta$可以由三个矩阵来进行描述$\theta =(A,B,\Pi )$。

1. 大小为 N*N (N代表N种隐性的状态)的过渡矩阵 A:

A={aij}={P(xt+1=sj∣xt=si)}={P(sj∣si)}A = \{a_{ij}\} = \{P(x_{t+1}=s_j|x_t=s_i)\} = \{P(s_j|s_i)\}A={aij​}={P(xt+1​=sj​∣xt​=si​)}={P(sj​∣si​)}
过渡矩阵A中的每一个元素表示由上一个隐性状态$s_i$变为下一个隐性状态的条件概率。

2. 大小为 1*N 的初始概率向量 $\Pi$ :

$$\Pi ={\pi }_i=P(x_1=s_i)=P(s_i)$$
初始概率向量$\Pi$中的每一个元素，表示初始隐性状态为$s_i$的概率，该向量的长度N与隐性状态的可能取值个数相同。

3. 大小为 M*N 的观测矩阵 B :

B={bki}={P(yt=ok∣xt=si)}={P(ok∣si)}B = \{b_{ki}\} = \{P(y_t=o_k|x_t=s_i)\} = \{P(o_k|s_i)\} B={bki​}={P(yt​=ok​∣xt​=si​)}={P(ok​∣si​)}

观测矩阵B中的每个元素，是用来描述N个隐形状态对应M个观测状态的概率。即在隐形状态为$s_i$ 的条件下，观测状态为$o_k$的概率。

上述三个矩阵构成了一个完整的隐马尔可夫模型。

掷骰子问题可以帮助我们更好地理解显性状态链和隐性状态链。例如我们有三个面数不一样的骰子可供选择投掷，三个骰子一个面数为4，一个面数为6，一个面数为8。每次选择的骰子是随机的且满足继续选到同一个骰子的概率是选到其他骰子概率的两倍。此时，隐性状态链$x_t$就是我们每次选择的骰子，取值集合就是骰子1，骰子2，骰子3。显性状态链就是我们掷出的一系列数值，取值集合为{1,2,3,4,5,6,7,8}\{1,2,3,4,5,6,7,8\}{1,2,3,4,5,6,7,8}。
根据上述的信息，我们不难整理出这个骰子问题HMM模型的三个核心矩阵 :

过渡矩阵A

previous state \ current state	D4	D6	D8
D4	$\frac{2}{3}$	$\frac{1}{6}$	$\frac{1}{6}$
D6	$\frac{1}{6}$	$\frac{2}{3}$	$\frac{1}{6}$
D8	$\frac{1}{6}$	$\frac{1}{6}$	$\frac{2}{3}$

初始概率向量$\Pi$
由于一开始是随机选取骰子，因此初始抽到三个骰子的概率是相同的$\frac{1}{3}$。

D4	D6	D8
$\frac{1}{3}$	$\frac{1}{3}$	$\frac{1}{3}$

观测矩阵B :
观测矩阵存放了每种隐性状态下各观测状态的条件概率

observed state \ hidden state	D4	D6	D8
1	$\frac{1}{4}$	$\frac{1}{6}$	$\frac{1}{8}$
2	$\frac{1}{4}$	$\frac{1}{6}$	$\frac{1}{8}$
3	$\frac{1}{4}$	$\frac{1}{6}$	$\frac{1}{8}$
4	$\frac{1}{4}$	$\frac{1}{6}$	$\frac{1}{8}$
5	0	$\frac{1}{6}$	$\frac{1}{8}$
6	0	$\frac{1}{6}$	$\frac{1}{8}$
7	0	0	$\frac{1}{8}$
8	0	0	$\frac{1}{8}$

使用HMM模型会衍生出三类基础问题 :

第一类问题. 在HMM模型$\theta$确定的情况下，针对一组显性状态链$y_{1:T}$，计算出该模型生成这组显性状态链的概率。

还是沿用上述的骰子为例，这类问题实际就是在确定了骰子的种类，D4D6D8(观测矩阵B)，随机选择骰子的方式(过渡矩阵A和初始概率向量$\Pi$)的情况下，用于计算使用这三个骰子掷出一组数的可能性的问题。计算显性状态链的生成概率也是对模型是否和实际情况吻合的检验。如果多次实验得出的显性状态链由该模型生成的概率很低，可以认为构成模型的三个矩阵要素可能与实际情况不符，即在该例子中使用的骰子有被掉包的可能或实际随机选取骰子的方式与模型不一致。下面给出计算该问题的其中一个算法 forward algorithm 前向算法，与之对应的还有backward algorithm 后向算法，由于两者的思路很相似，本文中不再对backward algorithm作过多阐述。

forward algorithm 前向算法

计算一组显性状态链 $y_1,y_2,y_3,...,y_T$ (下文中统一称$y_{1:T}$)生成的概率，即计算 𝔏 = $P(y_{1:T})$。我们先定义 ${\alpha }_i(t)=P(y_{1:t},x_t=s_i)$, 其中$s_i$是该显性状态链最后一个值对应的隐性状态，$t$是这个显性状态链的总长度。 ${\alpha }_i(t)$即是掷出该显性状态链，且最后一个隐形状态为$s_i$的概率，则最终要求的生成概率 𝔏 = $$\begin{gathered} P(y_{1:T})=P(y_{1:T},x_t=s_1)+P(y_{1:T},x_t=s_2)+...+P(y_{1:T},x_t=s_N) \\ =\sum _{i=1}^N{\alpha }_i(T) \end{gathered}$$
这样我们就把问题转化为了计算${\alpha }_i(T)=P(y_{1:T},x_t=s_i)$。使用前向算法forward algorithm如下图所示，通过去掉最后一个状态的状态链的生成概率${\alpha }_i(t-1)$推导至${\alpha }_i(t)$，可以归纳为一个递推模型。
首先我们计算生成$y_{1:t}$这个显性状态链且最后一个隐形状态为$s_1$的概率${\alpha }_1(t)$，通过这个概率我们就不难推导出同样生成该显性状态链，且最后一个隐性状态为$s_2,s_3,...,s_N$的概率${\alpha }_i(t)$。

$$\begin{gathered} {\alpha }_1(t)=P(y_{1:t},x_t=s_i)=P(y_{1:t-1},x_{t-1}=s_1)\ast P(s_1|s_1)\ast P(y_t|s_1) \\ +P(y_{1:t-1},x_{t-1}=s_2)\ast P(s_1|s_2)\ast P(y_t|s_1) \\ + \\ ... \\ +P(y_{1:t-1},x_{t-1}=x_N)\ast P(s_1|s_N)\ast P(y_t|s_1) \end{gathered}$$

对于每个前置状态$x_{t-1}$,计算${\alpha }_1(t)$为三个概率的乘积 : 前一时刻$t-1$隐形状态为$s_i$且状态链为$y_{1:t-1}$的概率，由前一时刻的隐形状态$s_i$过渡到当前时刻$s_1$的过渡概率，以及在隐性状态为$s_1$的情况下，显性状态呈现$y_t$的概率。
整理可得 :
$$\begin{gathered} {\alpha }_1(t)={\alpha }_1(t-1)\ast a_{11}\ast b_{yt,1} \\ + \\ {\alpha }_2(t-1)\ast a_{21}\ast b_{yt,1} \\ + \\ ... \\ + \\ {\alpha }_N(t-1)\ast a_{N1}\ast b_{yt,1} \\ =b_{yt,1}\sum _{k=1}^N{\alpha }_k(t-1)\ast a_{k1} \end{gathered}$$
tips : $a_{ij}$是过渡矩阵A中第i行第j列元素，$a_{ij}=P(x_{t-1}=s_j|x_t=s_i)$即由前一个隐形状态$s_i$过渡到现状态$s_j$的概率。

同理可得$${\alpha }_2(t)=b_{yt,2}\sum _{k=1}^N{\alpha }_k(t-1)\ast a_{k2}$$
不难推出$${\alpha }_i(t)=b_{yt,i}\sum _{k=1}^N{\alpha }_k(t-1)\ast a_{ki}for(t>=2)$$
我们就得到了一个计算${\alpha }_i(t)$的递推模型，初始值$${\alpha }_i(1)={\pi }_i\ast b_{y_1,i}$$ tips : ${\pi }_i$是初始隐性状态为i的概率(见初始概率向量$\Pi$)。

且𝔏 = $$\begin{gathered} P(y_{1:T})=P(y_{1:T},x_t=s_1)+P(y_{1:T},x_t=s_2)+...+P(y_{1:T},x_t=s_N) \\ =\sum _{i=1}^N{\alpha }_i(T) \end{gathered}$$易得$$P(y_{1:T})=\sum _{i=1}^N(b_{y_{T,1}}\sum _{k=1}^N{\alpha }_k(T-1)\ast a_{k1})for(T>=2)$$

第二类问题. 在HMM模型$\theta$确定的情况下，给出一组显性状态链$y_{1:T}$，找出与其对应的可能性最大 的隐性状态序列$x_{1:T}$。

这个问题广泛应用于自然语言处理NLP中。例如在我们获取了一组手写符号时，找出其最大概率对应的字符。解决这个问题我们就要用到viterbi algorithm 维特比算法。维特比算法仍然利用递推关系,由$t$时刻的状态推出$t+1$时刻的状态。定义${\delta }_i(t)$为$t$时刻且显性状态为$s_i$的最优路径的概率。
初始情况非常简单，在$t=1$时刻，${\delta }_i(t)={\pi }_i\ast b_{y_1,i}$，每个隐性状态对应的概率均是初始状态概率与对应的显性状态观测概率的乘积。紧接着有递推关系 :
$${\delta }_i(t+1)=[\underset{j}{\max }{\delta }_j(t)a_{ji}]\ast b_{y_{t+1},i}$$在这个递推关系中，下一时刻$t+1$每个隐性状态$i$的最大概率${\delta }_i(t+1)$总是对应上一时刻$t$中最优路径概率${\delta }_j(t)$与过渡概率$a_{ji}$乘积最大的隐性状态$j$。因此，最优的隐性状态路径序列$\Psi_i(t) = [\argmax_j\delta_i(t)a_{ji}]。下面引用索邦大学课件上的例题作为使用维特比算法的例子 :

给定一个HMM模型$\theta$，其中
$$A=\left|\begin{array}{ccc}0.3& 0.5& 0.2\\ 0& 0.3& 0.7\\ 0& 0& 1\end{array}\right|$$

$$\Pi =\left|\begin{array}{c}0.6\\ 0.4\\ 0\end{array}\right|$$

$$B=\left|\begin{array}{ccc}1& 0.5& 0\\ 0& 0.5& 1\end{array}\right|$$
计算显性序列为 y = [1,1,2,2] 时，对应可能性最大的隐性序列。
结果如下 :
最终的最大概率隐性状态链为[1 2 3 3]，对应的概率为0.105。

第三类问题. Baum-Welch算法，给定一个序列$y_{1:T}$，估算出能最大化该序列生成概率的HMM模型的参数。

这类问题是非常常见的，HMM建模的问题。通过一组给定的显性状态链，推测出可能性最大的HMM模型。 首先我们需要用到前向算法中提到的${\alpha }_i(t)$以及后向算法中的${\beta }_i(t)$。前向算法的推导过程可以参考之前的部分，类似的后向算法如果有兴趣可以自行再查找资料，这里不过多赘述只列出使用的公式。
根据之前前向算法部分，我们定义了${\alpha }_i(t)=P(y_{1:t},x_t=s_i)$这个可以递推出的概率，${\alpha }_i(t)$代表了观测状态为$y_{1:t}$且第$t$个隐性状态$s_i$的概率。同样在后向算法中我们定义了${\beta }_i(t)=P(y_{t+1},...,y_T,x_t=s_i)$代表了观测状态为$y_{t+1:T}$且第$t$个隐性状态为$s_i$的概率。

显然我们有${\alpha }_i(t){\beta }_i(t)=P(y_{1:T},x_t=s_i)$ 且 𝔏 = $P(y_{1:T})$。

根据贝叶斯公式$P(A|B)=\frac{P(A,B)}{P(B)}$, [2] 我们不难得出在给定观测序列$y_{1:T}$以及HMM模型$\theta$的情况下，在$t$时刻状态是$s_i$的概率 :

$${\gamma }_i(t)=P(x_t=s_i|y_{1:T})=\frac{P(x_t=s_i,y_{1:T})}{P(y_{1:T})}=\frac{{\alpha }_i(t){\beta }_i(t)}{{\sum }_{j=1}^N{\alpha }_j(t){\beta }_j(t)}$$

在$t$时刻状态是$s_i$，在$t+1$时刻状态时$s_j$的概率 :
$${\xi }_{ij}=P(x_t=s_i,x_{t+1}=s_j|y_{1:T})=\frac{P(s_j|s_i){\alpha }_i(t){\beta }_j(t+1)P(y_{t+1}|s_j)}{P(y_{1:T})}=\frac{a_{ij}{\beta }_j(t+1)b_{y_{t+1},j}}{P(y_{1:T})}$$

通过这些基础的概率，我们更新$\theta$中的参数。

初始状态是$s_i$的概率向量:
$${\Pi }_i^{\ast }={\gamma }_i(1)$$

更新过渡矩阵:
$$a_{ij}^{\ast }=P(s_j|s_i)=\frac{{\sum }_{t=1}^{T-1}{\xi }_{ij}(t)}{{\sum }_{t=1}^{T-1}{\gamma }_i(t)}$$

更新观测矩阵:
$$b_{o_k,i}^{\ast }=\frac{{\sum }_{t=1}^T{\alpha }_i(t)P(y_t|s_i)1_{y_t=o_k}}{{\sum }_{t=1}^T{\alpha }_i(t){\beta }_i(t)}$$
其中
$1_{y_t=o_k}=1when{y}_t=o_k$
$1_{y_t=o_k}=0whenelse$

不断重复上述更新操作直至收敛，就得到了新的HMM模型。

\newline
\newline
\newline
\newline
\newline

References :

[1] Viterbi AJ (April 1967). “Error bounds for convolutional codes and an asymptotically optimum decoding algorithm”. IEEE Transactions on Information Theory. 13

[2] Baum-Welch algorithm