隐马尔可夫模型(Hidden Markov Model)

隐马尔可夫模型(Hidden Markov Model)

隐马尔可夫模型(Hidden Markov Model, HMM)是一个重要的机器学习模型。直观地说，它可以解决一类这样的问题：有某样事物存在一定的状态，但我们无法得知某个时刻（或位置）它所处在的状态，但是我们有一个参照事物，我们知道这个参照事物在某个时刻（或位置）的状态并认为参照事物的状态和原事物的状态存在联系，那么我们可以使用机器学习来推测原事物最有可能在一个时刻（或位置）处在什么样的状态。也就是说，这是一个基于概率统计的模型。

举一个直观的例子吧。比如你要预测某市连续(n)天内的天气是冷还是热，但是你不知道这个城市的气温，你只知道某市某天冰淇淋的总销量，那么我们可以训练隐马尔可夫模型来用上述数据推测城市的某天的天气最有可能是冷还是热。这里，我们要预测的事物是天气，有 {冷，热} 两个状态，参照事物是冰淇淋销量，有({1, 2, …})等若干种有限的状态，要预测的结果是(t)时刻最有可能的天气状态。

下面我们给出隐马尔可夫模型的严格数学定义。

一、马尔可夫链 (Markov Chain)

要介绍隐马尔可夫模型，必须先介绍马尔可夫链的基本概念。不严格地说，马尔可夫链是一个描述随机过程的数学模型。随机过程是指一个事物有一组状态，会随着时间的流逝按照一定概率随机变化的一个过程。马尔可夫链有两个重要假设，一是在(t)时刻事物状态变化的概率只与(t)时刻本身事物的状态有关，不与(t-1)时刻及之前的事物状态有关，二是从(i)状态改变到(j)状态的概率只与(i)与(j)状态本身有关，与时间(t)的值无关。这些假设可以大大简化我们对这个随机过程的研究，而且依然可以解释现实世界中的很多现象。

比如说，假如一个呼叫中心要接听电话，有两个状态：{没电话，有电话}。没电话时，下一秒有电话来的概率为(p)，没电话来的概率为(1-p)；正在打电话时，电话接着打的概率为(q)，电话挂断的概率为(1-q)。这就构成了一个最简单的具有两个状态的以秒为离散时间单位的马尔可夫链。

下面我们给出马尔可夫链的一个严格数学定义。

考察一个离散时间的随机过程(X_n)，这个(X_n)在一组有限状态集合(S={1, 2,…, N })上取值，那么对每个时刻(t)，如果(X_1=i_1, …, X_{t-1}=i_{t-1},X_t=i, mbox{ and } i_1, i_2,…,i_{t-1},iin S)，我们定义(t+1)时刻(X_{t+1}=j)的概率

[P_{ij}=P{X_{t+1}|X_t=i, X_{t-1}=i_{t-1},…,X_1=i_1} ]

我们假设(P_{ij}=P{X_{t+1}|X_t=i})。满足这个条件，就称此随机过程满足马尔可夫性。

我们再假设在任意(t)时刻，任意状态(i)到状态(j)的变化概率(P_{ij}=p(i, j))（与(t)无关），(p)是(S imes S ightarrow [0,1])的函数。满足这个条件，就称此随机过程满足时齐马尔可夫性。

于是我们可以将(p(i, j))写成一个(N imes N)矩阵(A={a_{ij}})，(a_{ij}=p(i,j))。这个矩阵就是马尔可夫链的转移概率矩阵。

给定初始概率分布(phi(i)=P{X_0=i})，那么就能得到任意一个具体的状态转移链的发生概率

[P{X_0=i_0,X_2=i_2,...,X_n=i_n}=phi(i_0)p(i_0,i_1)p(i_1,i_2)...p(i_{n-1},i_n) ]

我们还可以证明，(n)步转移概率(p_n(i,j))就是矩阵(A^n)在((i,j))位置上的元素，即(P{X_{n+1}=j|X_n=i}=A^n_{ij}.)

我们用一个图模型就能描述马尔可夫链。令所有的状态(S={1,2,…, {N }})为顶点集合，转移概率矩阵(A)描述图上的边，再定义特殊的开始状态和结束状态(q_0,q_F)，矩阵第(0)行表示初始概率分布和初始的边即可。在之后的讨论中，我们将频繁地用图来描述马尔可夫链。

二、隐马尔可夫模型（HMM)

隐马尔科夫模型(HMM)是马尔可夫链的一个拓展。一个HMM含有如下元素：

1. 一组内部状态序列(Q=q_1q_2,…,q_T)，(q_iin S={1, 2, .., N})，(S)是状态集合，有(N)种状态。
2. 一个(N imes N)状态转移概率矩阵(A={a_{ij}}). (a_{ij})表示从状态(i)转移到状态(j)的概率。
3. 一组对这个状态的外部观察序列(O=o_1o_2…o_T), 每个观察从集合(V={v_1,v_2,…,v_V})中取值。
4. 一个描述产生某样观察的可能性的序列(B=b_i(o_t))，表示(t)时刻对状态(i)观察到(o_t)的概率。
5. 特殊的开始状态(q_0)和结束状态(q_F)，用(a_{01}a_{02}…a_{0n})表示初始状态到某个状态的概率分布，(a_{1F}a_{2F}…a_{nF})表示某个状态转移到结束状态的概率分布。

对HMM而言之前谈的马尔可夫链的两个假设依然是必要的。

对于一个HMM，有三个根本性的问题需要解决：

1. Likelihood：给定一个HMM的(A)和(B)，(lambda=(A,B))和观察序列(O)，判断观察到这个(O)的概率(P(O|lambda ))。

2. Decoding：给定一个HMM的(A)和(B)，(lambda=(A,B))和观察到的序列(O)，判断最有可能的内部状态序列(Q)。

3. Learning：给定观察序列(O)和状态集合(S)，学习HMM的参数(A)和(B)。

之后将会分别讨论解决这三个问题的算法，第三个会用到前两个的重要结论。

三、The Forward Algorithm

第一个问题的解决方案被称为Forward Algorithm，是一种动态规划算法。其思想是从时刻1开始计算一个时刻的所有状态下这个观察序列产出的概率，计算到时刻(T)并将所有状态的概率求和就得到了我们的结果。

定义$alpha_t(j) (为)t(时刻产生给定观察序列)o_1o_2…o_t(，且内部状态为)j$的概率。也就是

[alpha_t(j)=P{o_1,o_2,...,o_t,q_t=j|lambda} ]

我们发现这个问题(alpha_t(j))只与(t-1)时刻的状态有关，因为如果我们知道了(t-1)时刻产生观察序列(o_1o_2…o_{t-1})并且处于状态(i)的概率，那么下一时刻处在状态(j)的概率就是之前处在状态(i)的概率乘以从状态(i)转移到状态(j)的概率再乘以产生观察(o_t)的概率，再将所有的状态$求和即可。

所以有如下递推式

[alpha_t(j)=sum_{i=1}^{N}alpha_{t-1}(i)a_{ij}b_j(o_t) ]

初始状态

[alpha_1(j)=a_{0j}b_j(o_1), 1le jle N ]

最后结果

[P(O|lambda)=sum_{i=1}^Nalpha_T(i)a_{iF} ]

这就是动态规划算法的递推式。

四、The Viterbi Algorithm

第二个问题的解决方法和第一个问题非常相似，都是动态规划算法。

类似地，定义(v_t(j))表示前(t)时刻观察到了(o_1o_2…o_t)序列并且经过了最有可能的状态序列(q_0q_1…q_{t-1})，并且在(t)时刻处在状态(j)的概率，也就是说

[v_t(j)=max_{q_0,q_1,...,q_t-1}P{q_0,q_1,...,q_{t-1},o_1,o_2,...,o_t,q_t=j|lambda} ]

那么有如下递推方程

[v_t(j)=max_{i=1}^{N}v_{t-1}(i)a_{ij}b_j(o_t) ]

初始状态

[v_{1j}=a_{0j}b_j(o_1)，1le jle N ]

最终的最大概率

[P=max_{i=1}^{N}v_T(i)a_{iF} ]

递推过程中设置标记数组记录每一步选择了哪个(i)就能在递推完后恢复出整个路径。

五、HMM训练：The Forward-Backward Algorithm

回顾一下要解决的问题：给定观察序列(O)和状态集合(S)，学习HMM的参数(A)和(B)。如果我们学习到了这些参数，给定新的观察序列我们就能评估出最有可能的内部状态序列。

这个算法的得出过程比较复杂，需要一些数学推导。

我们先考虑马尔可夫链的情况。这个时候观察序列(O)直接与内部状态序列(Q)一一对应，那么我们可以直接估计其转移概率矩阵。其中从状态(i)转移到状态(j)的概率(a_{ij})，就是从状态(i)转移到状态(j)的次数，除以从状态(i)转移到某个状态的总次数，也就是

[a_{ij}=frac{C(i ightarrow j)}{sum_{qin Q}C(i ightarrow q)} ]

但是对于HMM，我们不能从观察序列直接估计内部状态序列的概率。因此，算法的思路是先估计一个(A)和(B)，然后不断重新估计…直到收敛，以完成参数的训练。而对概率的估计需要用到前向概率(alpha_t(i))和后向概率(eta_t(i))（Forward and backward probabilities）。其中

[alpha_t(i)=P{o_1,o_2,...,o_t,q_t=i|lambda}\ eta_t(i)=P{o_{t+1},o_{t+2},...,o_T|q_t=i,lambda} ]

(alpha_t(i))就是第一个算法里的(alpha_t(i))，(eta_t(i))与之相似，但注意有所不同，如此定义是为了其从后向前推导的计算过程与第一个算法类似，也就是

[eta_T(i)=a_{iF}, 1 le i le N\ eta_t(i)=sum_{j=1}^{N}a_{ij}b_j(o_{t+1})eta_{t+1}(j), 1le ile N,1le t < T ]

那么，要怎么估计(a_{ij})呢？令记号(hat{a}_{ij})表示对(a_{ij})的估计值，那么

[hat a_{ij}=frac{从状态i到状态j的期望转移次数}{从状态i转走的期望转移次数} ]

令(xi_t)为时刻(t)为状态(i)，时刻(t+1)时为状态(j)的概率，即(xi_t(i,j)=P{q_t=i,q_{t+1}=j|O,lambda}).

这个东西不太好求，但是由于概率论公式

[P{X|Y,Z}=frac{P{X,Y|Z}}{P{Y|Z}} ]

所以

[egin{align} xi_t(i,j)&=frac{P{q_t=i,q_{t+1}=j,O|lambda}}{P{O|lambda}}\ &=frac{alpha_t(i)a_{ij}b_j(o_{t+1})eta_{t+1}(j)}{P{O|lambda}} end{align} ]

分母就是之前动态规划算法的结果，分子则是用两个动态规划算法的递推矩阵里的值计算。

那么

[hat a_{ij}=frac{sum_{t=1}^{T-1}xi_t(i,j)}{sum_{t=1}^{T-1}sum_{k=1}^{N}xi_t(i,k)} ]

所以说只要给定了一个(B)的估计值，(B={b_j(o_t)})，我们就能估计一个(A={a_{ij}}).

之后我们有了一个(A)，我们还需要一个方法来估计(B)。对(B)中的每个元素，在状态(j)观测到观察(v_k)的概率的估计是(hat{b}_j(v_k))，我们将其表示为

[hat{b}_j(v_k)=frac{状态j并且观测到v_k的总数期望}{状态j出现的总数期望} ]

为了计算(hat b_j(v_k))，我们需要知道在时间(t)处在状态(j)的概率(gamma_t(j))

[gamma_t(j)=P{q_t=j|O,lambda}=frac{P{q_t=j,O|lambda}}{P{O|lambda}}=frac{alpha_t(j)eta_t(j)}{P{O|lambda}} ]

那么

[hat b_j(v_k)=frac{sum_{t=1, s.t. o_t=v_k}^Tgamma_t(j)}{sum_{t=1}^Tgamma_t(j)} ]

因此，对于算法而言，这是一个EM(Expectation, Maximization)算法，在E阶段我们计算(xi)和(gamma)，在M阶段我们估计(A)和(B)直到收敛。算法伪代码如下

Forward-Backward(长度为T的观察，观察取值集合V，隐藏状态取值集合S):
    初始化A和B
    迭代直到收敛:
        # E阶段
        for each t and j:
            计算gamma[t][j]
        for each t, i and j
            计算xi[t][i][j]
        # M阶段
        计算A[i][j]
        计算B[j][v[k]]
    return A, B

参考文献

1. Daniel Jurafsky, James H. Martin. Speech and Language Processing. Draft of August 7, 2017.
2. Gregory F. Lawler. Introduction to Stochastic Processes(Second Edition). China Machine Press.