快速了解并掌握HMM（隐马尔可夫模型）

最新推荐文章于 2024-07-15 17:08:02 发布

Briwisdom

最新推荐文章于 2024-07-15 17:08:02 发布

阅读量2.3k

点赞数 5

分类专栏：机器学习&深度学习文章标签：前向后向算法维特比算法鲍姆-韦尔奇算法 HMM 状态转移概率

本文链接：https://blog.csdn.net/u010420283/article/details/84075842

版权

机器学习&深度学习专栏收录该内容

38 篇文章 8 订阅

订阅专栏

导读

本文的结构安排如下目录。首先介绍了HMM的基本定义；然后引出HMM的3个基本问题中；在前两个环节中通过一个浅显易懂的案例帮助理解；其次是针对3个问题中用到算法的理论介绍；最后介绍了预测天气的HMM经典案例，并附有python代码。

隐马尔可夫模型的定义

隐马尔可夫模型是关于时序的概率模型，描述由一个隐藏的马尔科夫链随机生成不可观测的状态随机序列，再由各个状态生成一个观测而产生观测随机序列的过程。隐藏的马尔科夫链随机生成的状态序列，成为状态序列（state sequence）；每个状态生成一个观测，而由此产生的观测的随机序列，成为观测序列（observation sequence）。序列的每一个位置又可以看作是一个时刻。隐马尔可夫模型由初始概率分布、状态转移概率分布以及观测概率分布确定。

HMM的形式定义

对于HMM模型，首先我们假设Q是所有可能的隐藏状态的集合，V是所有可能的观测状态的集合，即：

其中，N是可能的隐藏状态数，M是所有的可能的观察状态数。

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

A是状态转移概率矩阵：

其中，

是在时刻 t 处于状态 $q{_{i}}^{}$ 的条件下载时刻 t+1 转移到状态 $q{_{j}}^{}$ 的概率。

B是观测概率矩阵：

其中，

是在时刻 t 处于状态 $q{_{j}}^{}$ 条件下生成观测 $v{_{k}}^{}$ 的概率。

$\pi$ 是初始状态概率向量：

其中，

是时刻 t=1处于状态 $q{_{i}}^{}$ 的概率。

隐马尔可夫模型模型，由隐藏状态初始概率分布 $\pi$ , 状态转移概率矩阵A和观测状态概率矩阵B决定。 $\pi$ 和A决定状态序列，B决定观测序列。因此，隐马尔可夫模型 $\lambda$ 可以用三元符号表示，即

直观理解HMM的各个状态

假设我手里有三个不同的骰子。第一个骰子是我们平常见的骰子（称这个骰子为D6），6个面，每个面（1，2，3，4，5，6）出现的概率是1/6。第二个骰子是个四面体（称这个骰子为D4），每个面（1，2，3，4）出现的概率是1/4。第三个骰子有八个面（称这个骰子为D8），每个面（1，2，3，4，5，6，7，8）出现的概率是1/8。

假设我们开始掷骰子，我们先从三个骰子里挑一个，挑到每一个骰子的概率都是1/3。然后我们掷骰子，得到一个数字，1，2，3，4，5，6，7，8中的一个。不停的重复上述过程，我们会得到一串数字，每个数字都是1，2，3，4，5，6，7，8中的一个。例如我们可能得到这么一串数字（掷骰子10次）：1 6 3 5 2 7 3 5 2 4

这串数字叫做可见状态链。但是在隐马尔可夫模型中，我们不仅仅有这么一串可见状态链，还有一串隐含状态链。在这个例子里，这串隐含状态链就是你用的骰子的序列。比如，隐含状态链有可能是：D6 D8 D8 D6 D4 D8 D6 D6 D4 D8

一般来说，HMM中说到的马尔可夫链其实是指隐含状态链，因为隐含状态（骰子）之间存在转换概率（transition probability）。在我们这个例子里，D6的下一个状态是D4，D6，D8的概率都是1/3。D4，D8的下一个状态是D4，D6，D8的转换概率也都一样是1/3。这样设定是为了最开始容易说清楚，但是我们其实是可以随意设定转换概率的。比如，我们可以这样定义，D6后面不能接D4，D6后面是D6的概率是0.9，是D8的概率是0.1。这样就是一个新的HMM。

同样的，尽管可见状态之间没有转换概率，但是隐含状态和可见状态之间有一个概率叫做输出概率（emission probability）。就我们的例子来说，六面骰（D6）产生1的输出概率是1/6。产生2，3，4，5，6的概率也都是1/6。我们同样可以对输出概率进行其他定义。比如，我有一个被赌场动过手脚的六面骰子，掷出来是1的概率更大，是1/2，掷出来是2，3，4，5，6的概率是1/10。

隐含状态转换关系示意图：

这里，

转换概率即为状态转移矩阵A：

这里我们假设D4，D6，D8后面跟随任意一种筛子的概率相同，都是1/3。（实际问题中可以定义为不同概率，如D4后面出现D4，D6，D8的概率分别为1/,6，2/6，3/6）

	D4	D6	D8
D4	1/3	1/3	1/3
D6	1/3	1/3	1/3
D8	1/3	1/3	1/3

输出概率即为观测概率矩阵B：

	1	2	3	4	5	6	7	8
D4	1/4	1/4	1/4	1/4	0	0	0	0
D6	1/6	1/6	1/6	1/6	1/6	1/6	0	0
D8	1/8	1/8	1/8	1/8	1/8	1/8	1/8	1/8

初始状态概率向量 $\pi$ ：

D4	D6	D8
1/3	1/3	1/3

隐马尔可夫模型的3个基本问题

概率计算问题：评估观察序列概率。即给定模型λ=(A,B, $\pi$ )和观测序列O={o1,o2,...oT}，计算在模型 λ 下观测序列O出现的概率P(O|λ)。这个问题的求解需要用到前向后向算法，是三个基本问题中最简单的算法。
预测问题，也称为解码问题。即给定模型λ=(A,B, $\pi$ )和观测序列O={o1,o2,...oT}，求给定观测序列条件概率P(I |O)下，最可能出现的对应的状态序列 I =( $i{_{1}}^{}$ , $i{_{2}}^{}$ ,... $i{_{T}}^{}$ )。即给定观测序列，求最有可能的对应的状态序列。这个问题的求解需要用到基于动态规划的维特比算法，这个问题是HMM模型三个问题中复杂度居中的算法。
模型参数学习问题。即给定观测序列O={o1,o2,...oT}，估计模型λ=(A,B, $\pi$ )的参数，使该模型下观测序列的条件概率P(O | λ)最大。这个问题的求解需要用到基于EM算法的鲍姆-韦尔奇算法，是HMM模型三个问题中最复杂的。

简单的概率计算例子

在进行上面3个问题用到的方法讲解之前，也讲一个简单的例子。以上文的骰子为例，我们已知骰子有几种，每种骰子是什么，每次掷的都是什么骰子，根据掷骰子的结果，求产生这个结果的概率。

解法无非就是概率相乘：

1. 概率计算问题的求解

概率计算问题用掷骰子为例，就是如果你在赌场中总是输钱，那么你可以产生怀疑：是不是赌场的骰子有猫腻？现在已知一个观察到的骰子序列O，依据已知的模型参数（A，B， $\pi$ ），计算出现这个序列的概率，看是否和你在赌场掷的实际概率相符，如果相差太大，那么你的骰子很大可能有问题了。比如说你怀疑自己的六面骰被赌场动过手脚了，有可能被换成另一种六面骰，这种六面骰掷出来是1的概率更大，是1/2，掷出来是2，3，4，5，6的概率是1/10。你怎么办么？

答案很简单，算一算正常的三个骰子掷出一段序列的概率，再算一算不正常的六面骰和另外两个正常骰子掷出这段序列的概率。如果前者比后者小，你就要小心了。

首先，如果我们只掷一次骰子：

看到结果为1。产生这个结果的总概率可以按照如下计算，总概率为0.18：

把这个情况拓展，我们掷两次骰子：

看到结果为1，6。产生这个结果的总概率可以按照如下计算，总概率为0.05：

继续拓展，我们掷三次骰子：

看到结果为1，6，3.产生这个结果的总概率可以按照如下计算，总概率为0.03：

同样的，我们一步一步的算，有多长算多长，再长的马尔可夫链总能算出来的。用同样的方法，也可以算出不正常的六面骰和另外两个正常骰子掷出这段序列的概率，然后我们比较一下这两个概率大小，就能知道你的骰子是不是被人换了。

解决这个问题的算法叫做前向后向算法（forward algorithm）。

2. 预测问题的求解

预测问题也就是解最大似然路径问题。举例来说，已知我有三个骰子：六面骰，四面骰，八面骰。并且知道我掷了十次的结果（1 6 3 5 2 7 3 5 2 4），但我不知道每次用了那种骰子，我想知道最有可能的骰子序列。

其实最简单而暴力的方法就是穷举所有可能的骰子序列，然后依照简单的概率计算例子的解法把每个序列对应的概率算出来。然后我们从里面把对应最大概率的序列挑出来就行了。如果马尔可夫链不长，当然可行。如果长的话，穷举的数量太大，就很难完成了。

另外一种很有名的算法叫做Viterbi algorithm. 要理解这个算法，我们先看几个简单的列子。
首先，如果我们只掷一次骰子：

看到结果为1.对应的最大概率骰子序列就是D4，因为D4产生1的概率是1/4，高于1/6和1/8.
把这个情况拓展，我们掷两次骰子：

结果为1，6.这时问题变得复杂起来，我们要计算三个值，分别是第二个骰子是D6，D4，D8的最大概率。显然，要取到最大概率，第一个骰子必须为D4。这时，第二个骰子取到D6的最大概率是

同样的，我们可以计算第二个骰子是D4或D8时的最大概率。我们发现，第二个骰子取到D6的概率最大。而使这个概率最大时，第一个骰子为D4。所以最大概率骰子序列就是D4 D6。
继续拓展，我们掷三次骰子：

同样，我们计算第三个骰子分别是D6，D4，D8的最大概率。我们再次发现，要取到最大概率，第二个骰子必须为D6。这时，第三个骰子取到D4的最大概率是：

同上，我们可以计算第三个骰子是D6或D8时的最大概率。我们发现，第三个骰子取到D4的概率最大。而使这个概率最大时，第二个骰子为D6，第一个骰子为D4。所以最大概率骰子序列就是D4 D6 D4。

写到这里，大家应该看出点规律了。既然掷骰子一二三次可以算，掷多少次都可以以此类推。我们发现，我们要求最大概率骰子序列时要做这么几件事情。首先，不管序列多长，要从序列长度为1算起，算序列长度为1时取到每个骰子的最大概率。然后，逐渐增加长度，每增加一次长度，重新算一遍在这个长度下最后一个位置取到每个骰子的最大概率。因为上一个长度下的取到每个骰子的最大概率都算过了。当我们算到最后一位时，就知道最后一位是哪个骰子的概率最大了。然后，我们要把对应这个最大概率的序列从后往前推出来。

这就是基于动态规划的维特比算法的求解最大可能隐藏序列的思路。

前向后向算法

前向后向算法是前向算法和后向算法的统称，这两个算法都可以用来求HMM观测序列的概率。我们先来看看前向算法是如何求解这个问题的。

前向算法本质上属于动态规划的算法，也就是我们要通过找到局部状态递推的公式，这样一步步的从子问题的最优解拓展到整个问题的最优解。

在前向算法中，通过定义“前向概率”来定义动态规划的这个局部状态。什么是前向概率呢, 其实定义很简单：定义时刻 t 时隐藏状态为qi, 观测状态的序列为o1,o2,...ot 的概率为前向概率。记为：

前向传播算法流程

输入：隐马尔可夫模型 $\lambda$ ，观测序列 O；

输出：观测序列概率 P（O| $\lambda$ ）。

1）初值

2）递推。对于 t =1,2,...T-1,

3) 终止

前向传播算法，步骤1）初始化前向概率，是初始时刻的状态和观测O1的联合概率。步骤2）是前向概率的递推公式，计算到时刻 t+1部分观测序列为且在时刻 t+1 处于状态q1的前向概率，如下图10.1所示。在步骤2）的方括号里，既然是到时刻 t 观测到并在时刻 t 处于状态的前向概率，那么乘积就是时刻 t 观测到并在时刻 t 处于状态而在时刻 t+1 到达状态的联合概率。对这个乘积在时刻 t 的所有可能的 N 个状态求和，其结果就是到时刻 t 观测为并在时刻 t+1 处于状态的联合概率。方括号里的值与观测概率的乘积恰好是到时刻 t+1 观测到并在时刻 t+1 处于状态的前向概率。步骤3）给出了的计算公式。

因为

所以

注：关于后向算法请参看博客最后的参考文献。

维特比算法概述

维特比算法实际是用动态规划解隐马尔可夫模型预测问题，即用动态规划（dynamic programming）求概率最大路径（最优路径）。这时一条路径对应着一个状态序列。

既然是动态规划算法，那么就需要找到合适的局部状态，以及局部状态的递推公式。在HMM中，维特比算法定义了两个局部状态用于递推。

第一个局部状态是在时刻tt隐藏状态为ii所有可能的状态转移路径i1,i2,...it中的概率最大值。记为 $\delta ^{_{t}}(i)$ :

由定义可得变量 $\delta$ 的递推公式：

第二个局部状态由第一个局部状态递推得到。定义在时刻 t 隐藏状态为 i 的所有单个状态转移路径（i1,i2,...,it−1,i）中概率最大的转移路径中第 t-1个节点的隐藏状态为 $\psi ^{_{t}}(i)$ ，其递推表达式为：

有了这两个局部状态，我们就可以从时刻0一直递推到时刻T，然后利用 $\psi ^{_{t}}(i)$ 记录的前一个最可能的状态节点回溯，直到找到最优的隐藏状态序列。

维特比算法流程：

输入：模型和观测：

输出：最优路径。

1)初始化

2）递推。对 t =2,3,...,T

3）终止

4）最优路径回溯。对 t = T-1,T-2,...,1

求得最优路径。

鲍姆-韦尔奇算法

鲍姆-韦尔奇算法原理既然使用的就是EM算法的原理，那么我们需要在E步求出联合分布P(O,I|λ)基于条件概率P(I|O, $\bar{\lambda }$ )的期望，其中 $\bar{\lambda }$ 为当前的模型参数，然后再M步最大化这个期望，得到更新的模型参数λ。接着不停的进行EM迭代，直到模型参数的值收敛为止。

首先来看看E步，当前模型参数为 $\bar{\lambda }$ , 联合分布P(O,I|λ)基于条件概率P(I|O, $\bar{\lambda }$ )的期望表达式为：

在M步，我们极大化上式，然后得到更新后的模型参数如下：

通过不断的E步和M步的迭代，直到 $\bar{\lambda }$ 收敛。

算法流程：

这里我们概括总结下鲍姆-韦尔奇算法的流程。

输入： D个观测序列样本{(O1),(O2),...(OD)}

输出：HMM模型参数

1)随机初始化所有的πi,aij,bj(k)

2) 对于每个样本d=1,2,...D用前向后向算法计算γ(d)t(i)，ξ(d)t(i,j),t=1,2...T

3) 更新模型参数：

4) 如果πi,aij,bj(k)的值已经收敛，则算法结束，否则回到第2）步继续迭代。

注：求HMM参数的算法有点难度，暂时没有完全理解。先放一个流程。有需要的可以最后参考文献部分。

案例和代码

任何一个HMM问题都可以通过下列五元组描述

:param obs:观测序列
:param states:隐状态
:param start_p:初始概率（隐状态）
:param trans_p:转移概率（隐状态）
:param emit_p: 发射概率（隐状态表现为显状态的概率）

预测隐藏状态问题：求解最可能的天气

背景：一个东京的朋友每天根据天气{下雨，天晴}决定当天的活动{公园散步,购物,清理房间}中的一种，我每天只能在twitter上看到她发的推“啊，我前天公园散步、昨天购物、今天清理房间了！”，那么我可以根据她发的推特推断东京这三天的天气。在这个例子里，显状态是活动，隐状态是天气。

这个例子用HMM来描述：

states = ('Rainy', 'Sunny')
 
observations = ('walk', 'shop', 'clean')
 
start_probability = {'Rainy': 0.6, 'Sunny': 0.4}
 
transition_probability = {
    'Rainy' : {'Rainy': 0.7, 'Sunny': 0.3},
    'Sunny' : {'Rainy': 0.4, 'Sunny': 0.6},
    }
 
emission_probability = {
    'Rainy' : {'walk': 0.1, 'shop': 0.4, 'clean': 0.5},
    'Sunny' : {'walk': 0.6, 'shop': 0.3, 'clean': 0.1},
}

求解最可能的隐状态序列是HMM的三个典型问题之一，通常用维特比算法解决。维特比算法就是求解HMM上的最短路径（log(prob)，也即是最大概率）的算法。

稍微用中文讲讲思路，很明显，第一天天晴还是下雨可以算出来：

定义V[时间][今天天气] = 概率，注意今天天气指的是，前几天的天气都确定下来了（概率最大）今天天气是X的概率，这里的概率就是一个累乘的概率了。
因为第一天我的朋友去散步了，所以第一天下雨的概率V[第一天][下雨] = 初始概率[下雨] * 发射概率[下雨][散步] = 0.6 * 0.1 = 0.06，同理可得V[第一天][天晴] = 0.24 。从直觉上来看，因为第一天朋友出门了，她一般喜欢在天晴的时候散步，所以第一天天晴的概率比较大，数字与直觉统一了。
从第二天开始，对于每种天气Y，都有前一天天气是X的概率 * X转移到Y的概率 * Y天气下朋友进行这天这种活动的概率。因为前一天天气X有两种可能，所以Y的概率有两个，选取其中较大一个作为V[第二天][天气Y]的概率，同时将今天的天气加入到结果序列中
比较V[最后一天][下雨]和[最后一天][天晴]的概率，找出较大的哪一个对应的序列，就是最终结果。

这个例子的Python代码：

# -*- coding: utf-8 -*-
states = ('Rainy', 'Sunny')

observations = ('walk', 'shop', 'clean')

start_probability = {'Rainy': 0.6, 'Sunny': 0.4}

transition_probability = {
    'Rainy': {'Rainy': 0.7, 'Sunny': 0.3},
    'Sunny': {'Rainy': 0.4, 'Sunny': 0.6},
}

emission_probability = {
    'Rainy': {'walk': 0.1, 'shop': 0.4, 'clean': 0.5},
    'Sunny': {'walk': 0.6, 'shop': 0.3, 'clean': 0.1},
}

# 打印路径概率表
def print_dptable(V):
    print("    ")
    for i in range(len(V)):
        print("%7d" % i)

    for y in V[0].keys():
        print("%.5s: " % y)
        for t in range(len(V)):
            print("%.7s" % ("%f" % V[t][y]))
        print()


def viterbi(obs, states, start_p, trans_p, emit_p):
    """
    :param obs:观测序列
    :param states:隐状态
    :param start_p:初始概率（隐状态）
    :param trans_p:转移概率（隐状态）
    :param emit_p: 发射概率 （隐状态表现为显状态的概率）
    :return:
    """
    # 路径概率表 V[时间][隐状态] = 概率
    V = [{}]
    # 一个中间变量，代表当前状态是哪个隐状态
    path = {}

    # 初始化初始状态 (t == 0)
    for y in states:
        V[0][y] = start_p[y] * emit_p[y][obs[0]]
        path[y] = [y]

    # 对 t > 0 跑一遍维特比算法
    for t in range(1, len(obs)):
        V.append({})
        newpath = {}

        for y in states:
            # 概率 隐状态 =    前状态是y0的概率 * y0转移到y的概率 * y表现为当前状态的概率
            (prob, state) = max([(V[t - 1][y0] * trans_p[y0][y] * emit_p[y][obs[t]], y0) for y0 in states])
            # 记录最大概率
            V[t][y] = prob
            # 记录路径
            newpath[y] = path[state] + [y]

        # 不需要保留旧路径
        path = newpath

    print_dptable(V)
    (prob, state) = max([(V[len(obs) - 1][y], y) for y in states])
    return (prob, path[state])

def example():
    return viterbi(observations,
                   states,
                   start_probability,
                   transition_probability,
                   emission_probability)
if __name__=='__main__':
    print(example())

输出：