《统计学习方法》第十章隐马尔可夫模型与python实现

一、写在前面

本文是《统计学习方法》第10章隐马尔可夫模型的读书笔记。用一段百余行的Python代码实现了隐马模型观测序列的生成、前向算法、维特比算法。公式与代码相互对照，循序渐进。
HMM算是个特别常见的模型，在自然语言处理中有很多的应用，比如基于字符序列标注的分词和词性标注了。但我的理解仅仅停留在“前向算法”“Viterbi”等层次。这次静下心来，从头到尾将这章认真看完，与自己原有的理解做一个对照，加深理解。

二、隐马尔可夫模型基本概念

隐马尔可夫模型是关于时序的概率模型，描述由一个隐藏的马尔可夫链随机生成不可观测的状态随机序列，再由各个状态生成一个观测而产生观测随机序列的过程。隐藏的马尔可夫链随机生成的状态的序列，称为状态序列（state sequence);每个状态生成一个观测，而由此产生的观测的随机序列，称为观测序列（observation sequence)。序列的每一个位置又可以看作是一个时刻。

什么叫隐马尔科夫链呢？简单说来，就是把时间线看做一条链，每个节点只取决于前N个节点。就像你打开朋友圈，发现你可以根据你的基友最近的几条状态猜测出接下来TA狗嘴里要吐什么东西一样。

接下来引入一些符号来表述这个定义——

设Q是所有可能的状态的集合，V是所有可能的观测的集合。

其中，N是可能的状态数，M是可能的观测数。

状态q是不可见的，观测v是可见的。应用到词性标注系统，词就是v，词性就是q。

I是长度为T的状态序列，O是对应的观测序列。

这可以想象为相当于给定了一个词（O）+词性（I）的训练集，于是我们手上有了一个可以用隐马尔可夫模型解决的实际问题。
定义A为状态转移概率矩阵：

其中，

是在时刻t处于状态qi的条件下在时刻t+1转移到状态qj的概率。

这实际在表述一个一阶的HMM，所作的假设是每个状态只跟前一个状态有关。

B是观测概率矩阵:

其中，

是在时刻t处于状态qj的条件下生成观测vk的概率（也就是所谓的“发射概率”）。

这实际上在作另一个假设，观测是由当前时刻的状态决定的，跟其他因素无关，这有点像Moore自动机。

π是初始状态概率向量：

其中，

是时刻t=1处于状态qj的概率。

隐马尔可夫模型由初始状态概率向量π、状态转移概率矩阵A和观测概率矩阵B决定。π和A决定状态序列，B决定观测序列。因此，隐马尔可夫模型屏幕快照 2016-06-28 下午9.59.49.png可以用三元符号表示，即:

$A,B,\pi$称为隐马尔可夫模型的三要素。如果加上一个具体的状态集合Q和观测序列V，则构成了HMM的五元组，又是一个很有中国特色的名字吧。
状态转移概率矩阵A与初始状态概率向量π确定了隐藏的马尔可夫链，生成不可观测的状态序列。观测概率矩阵B确定了如何从状态生成观测，与状态序列综合确定了如何产生观测序列。
从定义可知，隐马尔可夫模型作了两个基本假设：
(1)齐次马尔可夫性假设，即假设隐藏的马尔可夫链在任意时刻t的状态只依赖于其前一时刻的状态，与其他时刻的状态及观测无关。

从上式左右两边的复杂程度来看，齐次马尔可夫性假设简化了许多计算。
(2)观测独立性假设，即假设任意时刻的观测只依赖于该时刻的马尔可夫链的状态，与其他观测及状态无关。

三、隐马尔可夫python实例

让我们抛开教材上拗口的红球白球与盒子模型吧，来看看这样一个来自wiki的经典的HMM例子：

你本是一个城乡结合部修电脑做网站的小码农，突然春风吹来全民创业。于是跟村头诊所的老王响应总理号召合伙创业去了，有什么好的创业点子呢？对了，现在不是很流行什么“大数据”么，那就搞个“医疗大数据”吧，虽然只是个乡镇诊所……但管它呢，投资人就好这口。

数据从哪儿来呢？你把老王，哦不，是王老板的出诊记录都翻了一遍，发现从这些记录来看，村子里的人只有两种病情：要么健康，要么发烧。但村民不确定自己到底是哪种状态，只能回答你感觉正常、头晕或冷。有位村民是诊所的常客，他的病历卡上完整地记录了这三天的身体特征(正常、头晕或冷)，他想利用王老板的“医疗大数据”得出这三天的诊断结果(健康或发烧)。

这时候王老板掐指一算，说其实感冒这种病，只跟病人前一天的病情有关，并且当天的病情决定当天的身体感觉。

于是你一拍大腿，天助我也，隐马尔可夫模型的两个基本假设都满足了，于是统计了一下病历卡上的数据，撸了这么一串Python代码：

states = ('Healthy', 'Fever')
 
observations = ('normal', 'cold', 'dizzy')
 
start_probability = {'Healthy': 0.6, 'Fever': 0.4}
 
transition_probability = {
    'Healthy': {'Healthy': 0.7, 'Fever': 0.3},
    'Fever': {'Healthy': 0.4, 'Fever': 0.6},
}
 
emission_probability = {
    'Healthy': {'normal': 0.5, 'cold': 0.4, 'dizzy': 0.1},
    'Fever': {'normal': 0.1, 'cold': 0.3, 'dizzy': 0.6},
}

states代表病情，observations表示最近三天观察到的身体感受，start_probability代表病情的分布，transition_probability是病情到病情的转移概率，emission_probability则是病情表现出身体状况的发射概率。隐马的五元组都齐了，就差哪位老总投个几百万了。
有了前面的知识，我们就可以动手写一个一阶HMM模型的定义了：

class HMM:
    """
    Order 1 Hidden Markov Model
 
    Attributes
    ----------
    A : numpy.ndarray
        State transition probability matrix
    B: numpy.ndarray
        Output emission probability matrix with shape(N, number of output types)
    pi: numpy.ndarray
        Initial state probablity vector
    """
    def __init__(self,A,B,pi):
        self.A=A
        self.B=B
        self.pi=pi

此外，我们还需要一些utility来辅助我们实现该马尔可夫模型：

def generate_index_map(lables):
    index_label = {}
    label_index = {}
    i = 0
    for l in lables:
        index_label[i] = l
        label_index[l] = i
        i += 1
    return label_index, index_label
 
 
states_label_index, states_index_label = generate_index_map(states)
observations_label_index, observations_index_label = generate_index_map(observations)
 
 
def convert_observations_to_index(observations, label_index):
    list = []
    for o in observations:
        list.append(label_index[o])
    return list
 
 
def convert_map_to_vector(map, label_index):
    v = np.empty(len(map), dtype=float)
    for e in map:
        v[label_index[e]] = map[e]
    return v
 
 
def convert_map_to_matrix(map, label_index1, label_index2):
    m = np.empty((len(label_index1), len(label_index2)), dtype=float)
    for line in map:
        for col in map[line]:
            m[label_index1[line]][label_index2[col]] = map[line][col]
    return m
 
 
A = convert_map_to_matrix(transition_probability, states_label_index, states_label_index)
print(A)
B = convert_map_to_matrix(emission_probability, states_label_index, observations_label_index)
print(B)
observations_index = convert_observations_to_index(observations, observations_label_index)
pi = convert_map_to_vector(start_probability, states_label_index)
print(pi)
h = hmm.HMM(A, B, pi)

上述代码的输出为矩阵$A,B,\pi$:

[[0.7 0.3]
 [0.4 0.6]]
[[0.5 0.4 0.1]
 [0.1 0.3 0.6]]
[0.6 0.4]

观测序列生成Python实现:
为什么要生成观测序列呢？虽然我很想剧透并没有什么卵用，但王老板说他基于三年乡镇诊所病历卡的“医疗大数据库”里面并没有足够的数据给你跑什么算法。你只好自己动手，写算法生成一些了：

def simulate(self, T):
 
    def draw_from(probs):
        return np.where(np.random.multinomial(1,probs) == 1)[0][0]
 
    observations = np.zeros(T, dtype=int)
    states = np.zeros(T, dtype=int)
    states[0] = draw_from(self.pi)
    observations[0] = draw_from(self.B[states[0],:])
    for t in range(1, T):
        states[t] = draw_from(self.A[states[t-1],:])
        observations[t] = draw_from(self.B[states[t],:])
    return observations,states

这段代码其实就是用隐马尔可夫模型生成一个长度为T的观测和状态序列，比较短小精悍。具体的调用方法是：

observations_data,states_data=hmm.simulate(10)
print(observations_data)
print(states_data)

输出是随机的状态和序列。这里只打印了index，如果要糊弄投资人的话，还应该转换成相应的状态和观测。

[0 2 2 2 2 0 0 2 2 2]
[1 1 1 1 0 0 1 1 1 0]

四、隐马尔可夫模型的三个基本问题

隐马尔可夫模型有3个基本问题：
(1) 概率计算问题。给定模型$\lambda =(\pi ,A,B)$和观测序列$O$计,算$\lambda$模型下观测序列$O$出现的概率$p$.
(2) 学习问题。己知观测序列$O$,估计模型$\lambda =(\pi ,A,B)$的参数，使得在该模型下观测序列概率$P(O|\lambda )$最大。即用极大似然估计的方法估计参数。
(3)预测问题，也称为解码（decoding)问题。已知模型$\lambda =(\pi ,A,B)$和观测序列$O$，求对给定观测序列条件概率$P(O|\lambda )$最大的状态序列$I$。即给定观测序列，求最有可能的对应的状态序列。
下面各节将逐一介绍这些基本问题的解法。

1.概率计算算法

这节介绍计算观测序列概率的前向（forward)与后向（backward)算法，以及概念上可行但计算上不可行的直接计算法（枚举）。前向后向算法无非就是求第一个状态的前向概率或最后一个状态的后向概率，然后向后或向前递推即可。
直接计算法:
给定模型，求给定长度为T的观测序列的概率，直接计算法的思路是枚举所有的长度T的状态序列，计算该状态序列与观测序列的联合概率（隐状态发射到观测），对所有的枚举项求和即可。在状态种类为N的情况下，一共有N^T种排列组合，每种组合计算联合概率的计算量为T。计算量太大，不可取。
前向算法：
首先定义前向概率。定义(前向概率）给定隐马尔可夫模型$\lambda$，定义到时刻t为止的观测序列为$o_1,o_2......o_t$，且状态为$q_i$的前向概率，记作：

可以递推地求得前向概率${\alpha }_i$和观测序列概率$p(O|\lambda )$。
前向算法流程如下：
(1)初值：

前向概率的定义中一共限定了两个条件，一是到当前为止的观测序列，另一个是当前的状态。所以初值的计算也有两项（观测和状态），一项是初始状态概率，另一项是发射到当前观测的概率。
(2)递推对$t=1,2,3......T-1$:

每次递推同样由两部分构成，大括号中是当前状态为i且观测序列的前t个符合要求的概率，括号外的是状态i发射观测t+1的概率。
(3)终止:

由于到了时间T，一共有N种状态发射了最后那个观测，所以最终的结果要将这些概率加起来。

由于每次递推都是在前一次的基础上进行的，所以降低了复杂度。准确来说，其计算过程如下图所示：

下方标号表示时间节点，每个时间点都有N种状态，所以相邻两个时间之间的递推消耗$N^2$次计算。而每次递推都是在前一次的基础上做的，所以只需累加O(T)次，所以总体复杂度是O(T)个N^2。
前向代码python实现如下：

    def _forward(self,obs_seq):
        N=self.A.shape[0]
        T=len(obs_seq)
        F=np.zeros((N,T))
        F[:,0]=self.pi*self.B[:,obs_seq[0]]
        for t in range(1,T):
            for n in range(N):
                F[n,t]=np.dot(F[:,t-1],(self.A[:,n]))*self.B[n,obs_seq[t]]
        return F
 hmm._forward([1,2,1,2,0])

输出如下：

array([[0.24      , 0.0216    , 0.019872  , 0.00209088, 0.00270691],
       [0.12      , 0.0864    , 0.017496  , 0.00987552, 0.00065526]])

其中最后一列的数求和即为给定观测序列和模型参数后，模型输出此观测的概率。

2.学习算法

隐马尔可夫模型的学习，根据训练数据是包括观测序列和对应的状态序列还是只有观测序列，可以分别由监督学习与非监督学习实现。本文只介绍监督学习算法。
监督学习方法;
假设已给训练数据包含S个长度相同的观测序列和对应的状态序列:

那么可以利用极大似然估计法来估计隐马尔可夫模型的参数。具体方法如下。
1.转移概率${\alpha }_{ij}$的估计
设样本中时刻t处于状态i时刻t+1转移到状态j的频数为$A_{ij}$，那么状态转移概率${\alpha }_{ij}$估计为：

很简单的最大似然估计。
2. 观测概率矩阵的估计
设样本中状态为j并观测为k的频数是$B_{jk}$，那么状态为j观测为k的概率的估计是:

3. 初始状态概率$\pi$，这个就对所有样本的初始状态进行极大似然估计，也很简单。
学习算法较为简单，就是极大似然估计的反复运用。此处不再赘述。

3.预测算法

下面介绍隐马尔可夫模型预测的两种算法：近似算法与维特比算法（Viterbi algorithm)。其中维特比算法高效准确，是最为常用的预测、解码算法。
近似算法:
近似算法的想法是，在每个时刻t选择在该时刻最有可能出现的状态$i^{\ast }$，从而得到一个状态序列：

将它作为预测的结果。
给定隐马尔可夫模型$\lambda$和观测序列$O$，在时刻t处于状态$q_i$的概率为：

在每一时刻t最优可能的状态$i^{\ast }$是:

从而得到目标状态序列。
似算法的优点是计算简单，其缺点是不能保证预测的状态序列整体是最有可能的状态序列，因为预测的状态序列可能有实际不发生的部分。事实上，上述方法得到的状态序列中有可能存在转移概率为0的相邻状态，即对某些${\alpha }_{ij}=0$时。尽管如此，近似算法仍然是有用的。（没看出来有什么用，嗷呜！）
维特比算法：
维特比算法实际是用动态规划解隐马尔可夫模型预测问题，即用动态规划(dynamic programming)求概率最大路径（最优路径）。这时一条路径对应着一个状态序列。
维特比算法shiyizhong 动态规划算法用于最可能产生观测时间序列的-维特比路径-隐含状态序列，特别是在马尔可夫信息源上下文和隐马尔科夫模型中。术语“维特比路径”和“维特比算法”也被用于寻找观察结果最有可能解释的相关dongtai 规划算法。例如在统计句法分析中动态规划可以被用于发现最有可能的上下文无关的派生的字符串，有时被称为“维特比分析”。

利用动态规划寻找最短路径

动态规划是运筹学的一个分支，是求解决策过程最优化的数学方法，通常情况下应用于最优化的问题，这类问题一般有很多可行的解，每个解有一个值，而我们希望从中找到最优的答案。

在计算机科学领域，应用动态规划的思想解决的最基本的一个问题就是：寻找有向无环图（篱笆网络）当中两个点之间的最短路径（实际应用于地图导航、语音识别、分词、机器翻译等等）

下面举一个比较简单的例子做说明：求S到E的最短路径，如下图（各点之间距离不相同）：

我们知道，要找到S到E之间最短路径，最容易想到的方法就是穷举法。也就是把所有可能的路径都例举出来。从S走向A层共有4种走法，从A层走向B层又有4种走法，从B层走向C层又有4种走法，然后C层走向E点只有一种选择。所以最终我们穷举出了444=64种可能。显然，这种方法必定可行，但在实际的应用当中，对于数量及其庞大的节点数和边数的图，其计算复杂度也将会变得非常大，而计算效率也会随之降低。

因此，这里选择适用一种基于动态规划的方式来寻找最佳路径。

所谓动态规划。其核心就是“动态”的概念，把大的问题细分为多个小的问题，基于每一步的结果再去寻找下一步的策略，通过每一步走过之后的局部最优去寻找全局最优，这样解释比较抽象，下面直接用回刚刚的例子说明。如下图：

首先，我们假设S到E之间存在一条最短路径（红色），且这条路径经过C2点，那么我们便一定能够确定从S到C2的64条（444=64）子路经当中，该子路经一定最短。（证明：反证法。如果S到C2之间存在一条更短的子路经，那么便可以用它来替代原先的路径，而原来的路径显然就不是最短了，这与原假设自相矛盾）。

同理，我们也可以得出从S到B2点为两点间最短子路经的结论。这时候，真相便慢慢浮出水面：既然如此，我们计算从S点出发到点C2的最短路径，是不是只要考虑从S出发到B层所有节点的最短路径就可以了？答案是肯定的！因为，从S到E的“全局最短”路径必定经过这些“局部最短”子路经。没错！这就是上面提及到的通过局部最优的最优去寻找全局最优，问题的规模被不断缩小！

接下来，要揭晓答案了！继续看下图：

回顾之前的分析：我们计算从S到C2点的最短路径时候只需要考虑从S出发到B层所有节点的最短路径，B层也是。对B2来说，一共有4条路线可达，分别是A1->B2,A2->B2,A3->B2, A4->B2。我们需要做的就是A2->B2这条路线保留，而其他3条删掉（因为根据以上的分析，他们不可能构成全程的最短路线）。Ok，来到这里，我们会发现一个和小“漏洞”，这段S->A2->B2->C2->E的路线只是我一厢情愿的假设，最短路径下不一定是经过以上这些点。所以，我们要把每层的每个节点都考虑进来。
以下是具体做法：

step1：从点S出发。对于第一层的4个节点，算出他们的距离d(S,A1),d(S,A2),d(S,A3),d(S,A4),因为只有一步，所以这些距离都是S到它们各自的最短距离

step2:对于B层的所有节点(B1,B2,B3,B4),要计算出S到他们的最短距离。我们知道，对于特定的节点B2，从S到它的路径可以经过A层的任何一个节点(A1,A2,A3,A4)。对应的路径长就是d(S,B2)=d(S,Ai)+d(Ai,B2)（其中i=1，2，3，4）。由于A层有4个节点（即i有4个取值），我们要一一计算，然后找到最小值。这样，对于B层的每个节点，都需要进行4次运算，而B层有4个节点，所以共有4*4=16次运算。

step3:这一步是该算法的核心。我们从step2计算得出的阶段结果只保留4个最短路径值（每个节点保留一个）。那么，若从B层走向C层来说，该步骤的级数已经不再是44，而是变成4！也就是说，从B层到C层的最短路径只需要基于B层得出的4个结果来计算。这种方法一直持续到最后一个状态，每一步计算的复杂度为相邻两层的计算复杂度为44乘积的正比！再通俗点说，连接着两两相邻层的计算符合变成了“+”号，取代了原先的“”号。用这种方法，只需要进行44*2=32次计算！

其实上述的算法就是著名的维特比算法，事实上非常简单！

若假设整个网格的宽度为D，网格长度为N，那么弱适用穷举法整个最短路径的算法复杂度为O(D^{N),而适用这种算法的计算复杂度为O（ND}2）.试想一下，弱D与N都非常大，适用维特比算法的效率将会提高几个数量级！

五、维特比算法Python实现

还是利用我们刚刚建立的马尔可夫模型，在诊所来了一位病人，他最近三天的病状是：

observations = ('normal', 'cold', 'dizzy')

如何用Viterbi算法计算他的病情以及相应的概率呢？

    def viterbi(self,obs_seq):
        N=self.A.shape[0]
        T=len(obs_seq)
        prev=np.zeros((T-1,N),dtype=int)
        V=np.zeros((N,T))
        V[:,0]=self.pi*self.B[:,obs_seq[0]]
        for t in range(1,T):
            for n in range(N):
                seq_probs=V[:,t-1]*self.A[:,n]*self.B[n,obs_seq[t]]#前一层的所有可能状态的概率分别乘上转移到新的各个状态的概率，再乘上该状态发射出当前观测的概率
                prev[t-1,n]=np.argmax(seq_probs)
                V[n,t]=np.max(seq_probs)#只保留到这一点概率最大的路径对应概率
        return V,prev

这应该是目前最简洁最优雅的实现了，调用方法如下：

V, p = hmm.viterbi(observations_index)
print(" " * 7, " ".join(("%10s" % observations_index_label[i]) for i in observations_index))
for s in range(0, 2):
    print("%7s: " % states_index_label[s] + " ".join("%10s" % ("%f" % v) for v in V[s]))

输出结果如下：

            normal       cold      dizzy
Healthy:   0.300000   0.084000   0.005880
  Fever:   0.040000   0.027000   0.015120

解释：每一列概率最大的值对应状态为当前层预测状态，如，上式显示预测的状态序列为：Healthy,Healthy,Fever。该观测序列对应该状态序列的概率为：0.015120，是所有可能状态序列中概率最大的那种。
对比维基百科的结果：

         0          1          2
  Healthy: 0.30000 0.08400 0.005884 
    Fever: 0.04000 0.02700 0.015125 
The steps of states are Healthy Healthy Fever with highest probability of 0.01512

两者是完全一致的，小数点后最后一位有差异是保留位数不同引起的差异，但计算方式和结果以及预测结果是完全正确且相同的。