隐马尔可夫模型-standford

chapter6: 隐马尔科夫模型 standford tutorial, 可以说是看过的关于隐马尔科夫最好的教程了。要是看英文原版能和看中文一样easy该多好，可惜文化差异的情况下，即使是单词都认识，理解起来也会有点困难。对此，只能自己重新总结一下吧～
可以说，斯坦福从未让人失望过，太赞了！

也是无意中在google上看到这篇文章，才发现了这么好的一本书, Speech and language processing.

前言： 隐马尔可夫模型马尔可夫模型的后裔，它是一个序列模型(sequence model). 对于一个序列模型，它的工作是给序列中的每一个小单元分配标签label或是类别class.其中包括：part-of-speech tagging, named entity tagging, and speech recognition.

马尔可夫链 Markov chains

马尔可夫链和隐马尔科夫模型都是有限自动机(finite automation)的拓展。可以将他们看作是有权重的有限自动机(weighted finite automation),包括一些列状态和状态之间的转移关系，其中从一个状态到另一个状态的转移弧线是有权重的。在马尔可夫链中，这样的权重代表着概率，且从一个节点出来的概率之和为1.

上图中不仅包括了状态之间的转移概率，还包括了start和end两种特定的状态。

states	hot1	cold2	warm3
hot1	$a_{11}$	$a_{12}$	$a_{13}$
cold2	$a_{21}$	$a_{22}$	$a_{23}$
warm3	$a_{31}$	$a_{32}$	$a_{33}$

start0: {$a_{01},a_{02},a_{03}$}

end4: {$a_{14},a_{24},a_{34}$}

根据图9.1，马尔可夫链可以看做是一个概率图模型，其中包括以下几部分：

Q代表状态集合，A是状态转移矩阵， $a_{ij}$ 表示从状态 $q_i$ 到状态 $q_j$ 的概率,那么有 $\sum_{j=1}^na_{ij}=1,i=1,2...N$
$q_0$ 和 $q_F$是初始状态和终止状态。

做两个重要的假设： 以简化模型

(1) The Limited horiqon assumption 齐次假设
对于t时刻的状态，只取决于之前的一个状态。

$$Markov Assumption: P(q_i|q_1...q_{i-1})=P(q_i|q_{i-1})$$
也就是一阶马尔科夫链(a first-order Markov)。

图中的 $x_i$ 表示可观测时的 $q_i$

(2) Stationary process assumption 静态过程假设
the conditional distribution
over next state given current state does not change over time. 对于状态之间的条件概率不会随时间的变化而改变，也就是状态转移矩阵只有一个～

$$P(q_t|q_{t-1})=P(q_2|q_1);t \in 2...T$$

在很多其他的论文中，用 $\pi$ 来表示初始状态

都是一样的，不过这篇教程中用第一种表示方法。

隐马尔可夫模型

为了更形象的描述隐马尔可夫模型，文章举了这样的一个栗子：

Imagine that you are a climatologist in the year 2799 studying the history of global warming. You cannot find any records of the weather in Baltimore, Maryland, for the summer of 2007, but you do find Jason Eisner’s diary, which lists how many ice creams Jason ate every day that summer. Our goal is to use these observations to estimate the temperature every day. We’ll simplify this weather task by assuming there are only two kinds of days: cold (C) and hot (H). So the Eisner task is as follows:

Given a sequence of observations O, each observation an integer corresponding to the number of ice creams eaten on a given day, figure
out the correct ‘hidden’ sequence Q of weather states (H or C) which caused Jason to eat the ice cream.

总结下就是，观察到的序列是每天吃的冰淇淋数目2,3,1,2,3,….，从而判断每天的天气是hot or cold这两种状态中的哪一种。

定义一个隐马尔可夫模型，有以下组成部分：

两个假设

(1)一阶马尔科夫模型：

$$Markov Assumption: P(q_i|q_1...q_{i-1})=P(q_i|q_{i-1})$$
(2)条件独立，在状态 $q_i$ 的条件下，观测 $o_i$ 只取决于 $q_i$， 而与其他的状态和观测值够无关
$$Output Indepence: P(o_i|q_1...q_i,...,q_T,o_1,...,o_i,...,o_T)=P(o_i|q_i)$$

对ice cream task.问题的描述如下图：

注意到，图中所有的概率都不为零，这种HMM模型叫做 fully connected 或是 ergodic HMM. 但并不是所有状态都可以互相转移的，比如 left-to-right HMM,又称Bakis HMMs,通常用于语音处理。

左图表示Bakis HMM，右图是ergodic HMM.

关于隐马尔可夫模型的三个问题：

1. 问题1：计算似然概率。
   
   • 前向/后向算法
2. 问题2：解码问题，又称预测问题，已知模型参数和观测序列，求最有可能的状态序列
   
   • 维特比算法
3. 问题3：学习问题，已知观测序列，估计模型参数使得该模型下观测序列的概率最大。
   
   • 极大似然估计，Baum-Welch, EM算法

概率计算：前向算法

状态已知的话，是监督学习

以图9.3为例，观测序列为（3,1,3）假如我们知道隐藏状态是（hot, hot, cold）的话，在此基础上计算似然概率就很简单了。

也就是说，已知观测序列 $O=o_1,o_2,...,o_T.且隐藏状态序列已知$Q=q_0,q_1,…,q_T$,那么似然概率为：$$P(O|Q)=\prod_{i=1}^TP(o_i|q_i)$$

状态无法观测，无监督学习

但大多数情况下，状态是不知道的，因为我们需要取计算出现状态 (3,1,3) 的所有可能的隐藏状态序列。

观测序列$O=o_1,o_2,...,o_T，假定存在一个特定的隐藏状态序列$Q=q_0,q_1,…,q_T$，那么联合概率分布：$$P(O,Q)=P(O|Q)\times P(Q)=\prod_{i=1}^Tp(o_i|q_i)\times \prod_{i=1}^TP(q_i|q_{i-1})\tag{9.10}$$
所以隐马尔科夫是一个双重随机过程。

那么观察序列为（3,1,1）和状态序列为（hot, hot, cold）的联合概率为：

这样我们知道了怎么求一个特定的隐藏序列和观测序列的联合概率，那么所有可能隐藏序列的类和就是观测序列的总似然概率了。

$$P(O)=\sum_QP(O,Q)=\sum_QP(O|Q)P(Q)\tag{9.12}$$

对冰淇淋的例子，如果观测序列是（3,1,3），那么似然概率为：

如果有N中状态的话，对于长度为T的序列，其计算复杂度就是 $O(N^T)$ 这就太大了，所以得寻求更简单的解法。

forward algorithm

前向算法是一种动态规划算法，其计算复杂度是 $O(N^2T)$

The forward algorithm is a kind of dynamic programming algorithm, that is, an algorithm that uses a table to store intermediate values as it builds up the probability of the observation sequence. The forward algorithm computes the observation probability by summing over the probabilities of all possible hidden state paths that could generate the observation sequence, but it does so efficiently by implicitly folding each of these paths into a single forward trellis.

之所以高效的原因，是它将所有的路径都隐式的折叠到一个前向网格中。。

前向网格(forward trellis)中的每一个单元(cell) $\alpha_t(j)$ 表示给定模型 $\lambda$，t时刻观测序列为 $o_1,o_2,...,o_t$,状态为j的概率。

$$\alpha_t(j)=P(o_1,o_2,...,o_t,q_t=j|\lambda)\tag{9.13}$$
其中， $\alpha_t(j)$ 的计算是叠加所有可能的路径。
$$\alpha_t(j)=[\sum_{i=1}^N\alpha_{t-1}a_{ij}]b_j(o_t)$$
表示t时刻状态为j,从t-1时刻到t时刻，有N条路径，叠加～

以图9.7中的第二时间步和状态2为例，

$$\alpha_2(1)=\alpha_1(1)×P(H|H)×P(1|H) + α_1(2)×P(H|C)×P(1|H)$$

下图描述了前向网格中计算一个cell的步骤：

整个似然概率的计算过程，伪代码：

真的讲的太好了太清楚了！！！感动哭。。。说真的，要不中国的大学都改成英文教学吧。。看着那些翻译过来的书籍都头疼。。这么好的英文教学材料，为什么翻译过来之后就那么难理解了。。感觉很多老师可能自己懂了，但是讲出来的课或是写出来的书，完全就是应付任务的吧。。

好吧，回到主题。伪代码中：
概率矩阵 forward [N+2,T] 是包括了初始状态start和结束状态end,那么 forward[s,t]表示 $\alpha_t(s)$

1. initialization step:
   

   $$\alpha_1(j)=a_{0j}b_j(o_1),1\le j \le N$$
   伪代码中：forward[s,1] <– $a_{0,s}* b_s(o_1)$
   是从状态0到t=1时刻的状态s

2. Recursion (since states 0 anf F are non-emittinf):
   

   $$\alpha_t(j)=[\sum_{i=1}^N\alpha_{t-1}(i)a_{ij}]b_j(o_t)$$
   在伪代码中有两个循环，分别是对时间步2到T，以及每个时间步，其中的状态从1到N

for each time step from 2 to T do:
  for each state s from 1 to N do:
    forward[s,t] = sum_{s'=1}^N forward[s', t-1] * a_{s',s} * b_s{o_t}

$forward[s,t] = sum_{s'=1}^N forward[s', t-1] * a_{s',s}$
可以发现，概率矩阵forward中的每一个值表示的是t时刻状态为s的概率，也就是 $\alpha_t(s)$

1. Termination:
   $$P(O|\lambda)=\alpha_T(q_F)=\sum_{i=1}^N\alpha_T(i)a_{iF}$$
   伪代码中：

forward[$q_F$,T] = $\sum_{s=1}^N$ forward[s,T] * $a_{s,q_F}$

return forward[$q_F$,T] T时刻状态为 $q_F$的概率。感觉应该是T+1时刻吧。。。？？？

预测问题：维特比算法

Decoding: The Viterbi Algorithm
解码问题（预测问题）：给定HMM模型 $\lambda=(A,B)$ 和观察序列 $O=o_1,o_2,...,o_T$, 找出概率最大的隐藏状态序列 $Q=q_1q_2...q_T$

在前向算法中，我们知道了怎么计算特定隐藏状态序列和观测序列的联合概率，也就是公式（9.13），然后找出其中概率最大的序列就可以了对吧？但是我们知道状态序列有 $N^2$ 个，这样计算就太复杂了。于是，有了 Viterbi algorithm. 维特比算法是一种动态规划算法，类似于最小化编辑距离。

上图展示了，冰淇淋例子中，HMM模型参数已知，观测序列为（3,1,3）的情况下，计算最大隐藏状态序列的过程。Viterbi网格中每一个cell为 $v_T(j)$ 表示t时刻，观察序列为 $o_1,o_2,...,o_t$， 隐藏状态为j，前t-1的状态序列为 $q_1q_2...q_{t-1}$ 的概率。

$$v_t(j)=P(q_0q_1...q_{t-1},o_1,o_2,...,o_t,q_t=j|\lambda)$$

换句话说，t-1时刻的状态序列是这样 $q_1q_2...q_{t-1}$，有N个这样的序列（因为t-1时刻的状态有N个），然后计算出这N个序列中到t时刻状态为j的概率，找出其中最大值，就是从开始到t时刻状态为j的最大概率序列。

$$v_t(j)=max_{i=1}^N v_{t-1}(i)a_{ij}b_j(o_t)$$
对应到图（9.10）中，以 $v_2(2)$为例，
$$v_2(2)=max(v_1(1)* a_{12}* b_2(1),v_1(2)* a_{22}* b_2(1))$$

$$v_2(2)=max(v_1(1)* P(H|C)* P(1|H), v_1(2)* P(H|H)* P(1|H))$$

维特比算法的整个过程，伪代码如下：

我们发现Viterbi算法跟前向算法非常相似，除了前向算法是计算的sum，而Viterbi是计算的max，同时我们也发现，Viterbi相比前向算法多了一个部分：backpointers. 原因是因为前向算法只需要计算出最后的似然概率，但Viterbi不仅要计算出最大的概率，还要得到对应的状态序列。因此，在类似于前向算法计算概率的同时，记录下路径，并在最后backtracing最大概率的路径。

The Viterbi backtrace.

1. initialization:
   

   $$v_1(j)=a_{0j}b_j(o_1)\tag{9.20}$$
   $$b_{t1}(j)=0\tag{9.21}$$
   初始状态只有一个节点start，可确定为0

2. Recursion(recall that states 0 and $q_F$ are non-emitting):
   

   $$v_t(j)=max_{i=1}^Nv_{t-1}(i)a_{ij}b_j(o_t);1\le j \le N, 1\le t\le T\tag{9.22}$$
   $$b_{t_t}(j)=argmax_{i=1}^Nv_{t-1}(i)a_{ij}b_i(o_t);1\le j \le N, 1\le t\le T\tag{9.23}$$
   $b_{t_t}(j)表示t时刻状态为j的所有N个路径$(q_1,q_2,…,q_{t-1})$$ 概率最大的路径的第k-1个节点。
   以图9.12为例，对于t=2，状态为1的节点，其max()中有2项，分别是 $v_1(1)* a_{11}* b_1(1)$ 和 $v_1(2)* a_{21}* b_1(1)$,其中较大的项的节点就是 $max_{i=1}^Nv_{t-1}(i)a_{ij}b_i(o_t)$, 但t时刻也有N个节点，所以也需要记录，因此用arg

3. Termination:

The best score: $P*=v_T(q_F)=max_{i=1}^Nv_T(i)* a_{iF}\tag{9.24}$

计算出T时刻到end状态的最大概率，也就是所有路径中的最大概率。

The start of backtrace:$q_T*=b_{t_T}(q_F)=argmax_{i=1}^Nv_T(i)* a_{iF}\tag{9.25}$

表示T时刻中N个路径中概率最大的结点。

其实Viterbi算法的路径数量和前向算法的路径数量是一模一样的，只是一个是max，一个是sum，因此也可以参考图9.8.

学习问题：HMM Training: The Forward-Backward Algorithm

Learning: Given an observation sequence O and the set of possible states in the HMM, learn the HMM parameters A and B.
第三个问题，给定观测序列 $O=(o_1,o_2,...,o_T)$, 估计模型参数使得在该HMM模型下，观测序列的概率 $P(O|\lambda)$ 概率最大，即用极大似然估计的方法估计参数。

先考虑马尔可夫链

马尔可夫链其状态是可观察的，可以看作是退化的隐马尔可夫模型。即没有发射概率(emmision probablities)B.因此，我们需要学习的蚕食只有状态转移矩阵（probability matrix）A.

其中 $a_ij$ 表示从状态i转移到状态j的概率，可以用大数定律来计算。 $C(i\rightarrow)$ 表示观察到的序列中从状态i转移到状态j的数量。然后除以所有从状态i转移的总数量。

$$a_{ij}=\dfrac{C(i\rightarrow j)}{\sum_{q\in Q}C(i\rightarrow q)}\tag{9.26}$$
显然分母不包括最后T时刻出现状态i，因为end不属于Q.

隐马尔可夫模型： Baum-Welch算法

The Baum-Welch algorithm uses two neat intuitions to solve this problem. The first idea is to iteratively estimate the counts. We will start with an estimate for the transition and observation probabilities and then use these estimated probabilities to derive better and better probabilities. The second idea is that we get our estimated probabilities by computing the forward probability for an observation and then dividing that probability mass among all the different paths that contributed to this forward probability.

其实就是EM算法～

在此之前，先了解一下后向算法，可以看做反向的前向算法。

后向算法

对应的后向概率（Backward probability）$\beta_t(i)$ 表示给定hmm模型 $\lambda$, 在t时刻状态为i的条件下，t+1时刻的观测序列为 $o_{t+1},o_{t+2},...,o_T$的概率.

$$\beta_t(i)=P(o_{t+1},o_{t+2},...,0_T|q_t=i,\lambda)\tag{9.27}$$

1. initialization:
   

   $$\beta_T(i)=a_{iF},1\le i \le N$$
   在李航老师的《统计学习导论》这本书上 $\beta_T(i)=1$. $\beta_T(i)$ 的定义表示在T时刻状态为i，观测到序列为 $o_{T+1}$, 这东西不存在，可以看做是end吧，所以从T到end其概率应该是 $a_{iF}$,但是在李航老师的书上只有初始状态的概率 $\pi=a_{01}$,而没有 $a_{iF}$. 感觉跟具体在什么场景下有关系。。

2. Recursion(again since stetes 0 and $q_F$ are non-emitting)
   

   $$\beta_t(i)=\sum_{j=1}^N\beta_{t+1}(j)a_{ij}b_j(o_{t+1}),1\le j \le N, 1\le t\le T$$
   根据定义很好理解。$\beta_{t+1}(j)$ 表示给定hmm模型 $\lambda$, 在t+1时刻状态为j的条件下，t+1时刻的观测序列为 $o_{t+2},o_{t+3},...,o_T$的概率.那么就可以得到 $\beta_t(i)$ 了。

3. Termination:
   

   $$P(O|\lambda)=\alpha_T(q_F)=\beta_1(q_0)=\sum_{j=1}^\beta_1(j)a_{0j}b_j(o_1)$$

在初始模型参数下，用大数定律估计新的模型参数，也就是极大似然估计

根据公式(9.26)我们可以知道，状态ｉ到ｊ的概率：

$$\hat a_{ij}=\dfrac{expected\ number\ of\ transitions\ from\ state\ i \ to\ state\ j}{expected\ number\ of\ transitions\ from\ state\ i}$$
然而怎么计算这些numerator？试想，如果我们知道特定时刻t，从状态i转移到j的概率，那么就能计算所有时刻t的从i转移到j的数量。

定义 $\zeta_t$ 表示在给定模型参数和观察序列条件下，t时刻状态为i，t+1时刻状态为j的概率。

$$\zeta_t(i,j)=P(q_t=i,q_{t+1}=j|O,\lambda)\tag{9.32}$$
但是模型参数我们不知道呀，也是我们需要学习得到的。
这样我们先计算一个和 $\zeta_t$ 相似的概率。
$$not-quite-\zeta_t(i,j)=P(q_t=i,q_{t+1}=j,O|\lambda)\tag{9.33}$$

$\alpha_t(i)$ 和 $\beta_t(j)$ 是前向\ 后向算法中的定义。我们先看下前向算法计算的条件：

也就是problem1中的条件，给定 $\lambda$ 和 观察序列，求 $P(O|\lambda)$
我们可以用 $\alpha_t(i)$ 和 $\beta_t(j)$ ，来表示 $\zeta_t$, 是因为我们在计算 $\zeta_t$ 时是先假定有这个一个模型参数，比如初始参数～

那么：

$$not-quite-\zeta_t=\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)\tag{9.34}$$

根据bayes公式：

$$P(X|Y,Z)=\dfrac{P(X,Y,Z)}{P(Y,Z)}=\dfrac{P(X,Y,Z)}{P(Z)P(Y|Z)}=\dfrac{P(X,Y|Z)}{P(Y|Z)}\tag{9.35}$$
对应起来就是：
$$P(q_t=i,q_{t+1}=j|O,\lambda)=\dfrac{P(q_t=i,q_{t+1}=j,O|\lambda)}{P(O|\lambda)}=\dfrac{not-quite-\zeta_t}{P(O|\lambda)}\tag{9.36}$$
其中：
$$P(O|\lambda)=\alpha_T(q_F)=\beta_1(q_0)=\sum_{j=1}^N\alpha_t(j)\beta_t(j)\tag{9.37}$$
这一步最后面一个式子的理解可以看做是前向算法和后向算法在时刻t相遇。

看到这里会发现，李航老师书中179页，公式25-26的推导就有点逻辑不通了。

因此，现在就可以推导出 $\zeta_t$：

$$\zeta_t{i,j}=\dfrac{\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)}{\alpha_T(q_F)}\tag{9.37}$$

$\zeta_t$ 表示的是某一个时刻t，那么对于参数 $a_{ij}$ 的估计就是所有的时刻中i到j的总数除以i到k(k=1,2…N)的总数

$$\hat a_{ij}=\dfrac{\sum_{t=1}^{T-1}\zeta_t(i,j)}{\sum_{t=1}^{T-1}\sum_{k=1}^N\zeta_t(i,k)}\tag{9.38}$$

同样的道理，我们可以推理得到发射矩阵B的参数估计 $b_j(v_k)$
状态为j，观察得到 $v_k,v_k\in V$的概率：

在t时刻状态为j的概率，定义为 $\gamma_t(j)$

$$\gamma_t(j)=P(q_t=j|O,\lambda)\tag{9.40}$$
同样的道理：
$$\gamma_t(j)=\dfrac{P(q_t=j,O|\lambda)}{P(O|\lambda)}\tag{9.41}$$

同公式（9.37）一样，前向后向算法在t时刻相遇：

$$\gamma_t(j)=\dfrac{\alpha_t(j)\beta_t(j)}{P(O|\lambda)}$$
然后求整个时间段内的j到 $v_k$的总数，除以状态为j的总数。
$$\hat b_j(v_k)=\dfrac{\sum_{t=1,o_t=v_k}^T}{\sum_{t=1}^T}\gamma_t(j)$$

仔细回顾以下这个过程，在初始模型参数和观测序列的条件下，根据大数定律对模型参数进行更新。其实就是在初始模型参数下，根据极大似然估计求得观测序列的极大似然估计，然后在似然概率最大的条件下求得相应的模型参数。

可以看到这里直接用大数定律和李航老师书上，使用极大似然估计，然后求导得到的公式是一样的。

EM算法

E step:
- 根据初始模型参数或是M step得到的模型参数，得到后验概率。

M step:
- 根据E step中得到的概率，估计出新的模型参数。这里直接用大数定律得到的，其实其本质原理就是极大似然估计，也就是求出使得概率最大的模型参数。

那么迭代条件呢？什么情况下终止？在GMM中有log函数，这里呢。。
- 这里应该就是 $P(O|\lambda)$ 吧，在前向算法中有计算到在模型参数和观察序列条件下的极大似然估计。

根据GMM和HMM对使用EM算法进行参数估计的一点想法：

所以EM算法中的E step并不是求期望，而是在对模型参数进行估计时，在初始模型或previous模型的情况下，求得基于观测序列或是训练样本的用极大似然估计或是大数定律求得后验概率。

然后M STEP就是让这个概率最大的条件下更新模型参数。

总结：

在回顾下隐马尔可夫模型的三个问题：

1. 第一个问题，计算概率

   • 已知模型参数 $\lambda$ 和观测序列 $O$，求在该模型下，出现观测序列的概率。
   • 使用前向算法，一个动态回归的算法，把求长度为T的概率转换为t到t+1的概率sum
   • 这一问题其实主要是为后面两个问题铺垫的，因为一般的场景都是状态未知，更不可能知道模型参数了。

2. 预测问题，又称解码问题

   • 已知模型参数 $\lambda$ 和观测序列 $O$, 求概率最大的状态序列。
   • 使用Viterbi算法，类似于前向算法，不过每一步不是sum，而是max，并且需要回溯backpointers
   • 这个问题的应用场景就比较广了。

3. 学习问题：模型参数估计

   • 已知观测序列 $O$，估计模型参数 $\lambda$, 使得观测序列的概率 $P(O|\lambda)$ 最大。
   • 使用Baum-Welch(极大似然估计)或forward-backward(大数定律)算法，并使用EM算法迭代，对参数进行估计，