最近，我们使用隐马尔可夫模型开发了一种解决方案，并被要求解释这个方案。

由Kaizong Ye，Sherry Deng撰写

HMM用于建模数据序列，无论是从连续概率分布还是从离散概率分布得出的。它们与状态空间和高斯混合模型相关，因为它们旨在估计引起观测的状态。状态是未知或“隐藏”的，并且HMM试图估计状态，类似于无监督聚类过程。

可下载资源

例子

在介绍HMM背后的基本理论之前，这里有一个示例，它将帮助您理解核心概念。有两个骰子和一罐软糖。B掷骰子，如果总数大于4，他会拿几颗软糖再掷一次。如果总数等于2，则他拿几把软糖，然后将骰子交给A。现在该轮到A掷骰子了。如果她的掷骰大于4，她会吃一些软糖，但是她不喜欢黑色的其他颜色（两极分化的看法），因此我们希望B会比A多。他们这样做直到罐子空了。

隐马尔科夫模型（Hidden Markov Model，以下简称HMM）是比较经典的机器学习模型了，它在语言识别，自然语言处理，模式识别等领域得到广泛的应用。

当然，随着目前深度学习的崛起，尤其是RNN，LSTM等神经网络序列模型的火热，HMM的地位有所下降。

但是作为一个经典的模型，学习HMM的模型和对应算法，对我们解决问题建模的能力提高以及算法思路的拓展还是很好的。

1.1 什么样的问题需要HMM模型

首先我们来看看什么样的问题解决可以用HMM模型。

使用HMM模型时我们的问题一般有这两个特征：

１）我们的问题是基于序列的，比如时间序列，或者状态序列。

２）我们的问题中有两类数据，一类序列数据是可以观测到的，即观测序列；而另一类数据是不能观察到的，即隐藏状态序列，简称状态序列。

有了这两个特征，那么这个问题一般可以用HMM模型来尝试解决。这样的问题在实际生活中是很多的。比如：我现在在打字写博客，我在键盘上敲出来的一系列字符就是观测序列，而我实际想写的一段话就是隐藏序列，输入法的任务就是从敲入的一系列字符尽可能的猜测我要写的一段话，并把最可能的词语放在最前面让我选择，这就可以看做一个HMM模型了。再举一个，我在和你说话，我发出的一串连续的声音就是观测序列，而我实际要表达的一段话就是状态序列，你大脑的任务，就是从这一串连续的声音中判断出我最可能要表达的话的内容。

从这些例子中，我们可以发现，HMM模型可以无处不在。但是上面的描述还不精确，下面我们用精确的数学符号来表述我们的HMM模型。

1.2 HMM模型的定义

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

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

对于一个长度为 $T$ 的序列， $I$ 对应的状态序列, $O$ 是对应的观察序列，即：

HMM模型做了两个很重要的假设如下：

1）齐次马尔科夫链假设。即任意时刻的隐藏状态只依赖于它前一个隐藏状态。当然这样假设有点极端，因为很多时候我们的某一个隐藏状态不仅仅只依赖于前一个隐藏状态，可能是前两个或者是前三个。但是这样假设的好处就是模型简单，便于求解。如果在时刻 $t$ 的隐藏状态是 $i_{t}=q_{i}$ ，在时刻 $t+1$ 的隐藏状态是 $i_{t+1}=q_{j}$ ，则从时刻 $t$ 到时刻 $t+1$ 的HMM状态转移概率 $a_{ij}$ 可以表示为：

这样 $a_{ij}$ 可以组成马尔科夫链的状态转移矩阵 $A$ ：

2）观测独立性假设。即任意时刻的观察状态只仅仅依赖于当前时刻的隐藏状态，这也是一个为了简化模型的假设。如果在时刻 $t$ 的隐藏状态是 $i_{t}=q_{j}$ ，而对应的观察状态为 $o_{t}=v_{k}$ ，则该时刻观察状态 $v_{k}$ 在隐藏状态 $q_{j}$ 下生成的概率 $b_{j}(k)$ 满足：

这样 $b_{j}(k)$ 可以组成观测状态生成的概率矩阵 $B$ ：

除此之外，我们需要一组在时刻 $t=1$ 的隐藏状态概率分布 $\Pi$ ：

一个HMM模型，可以由隐藏状态初始概率分布 $\Pi$ ，状态转移概率矩阵 $A$ 和观测状态概率矩阵 $B$ 决定。 $\Pi,A$ 决定状态序列， $B$ 决定观测序列。因此，HMM模型可以由一个三元组 $\lambda$ 表示如下： $\lambda=(A,B,\Pi)$

1.3 一个HMM模型实例

下面我们用一个简单的实例来描述上面抽象出的HMM模型。这是一个盒子与球的模型，例子来源于李航的《统计学习方法》。

假设我们有3个盒子，每个盒子里都有红色和白色两种球，这三个盒子里球的数量分别是：

按照下面的方法从盒子里抽球，开始的时候，从第一个盒子抽球的概率是0.2，从第二个盒子抽球的概率是0.4，从第三个盒子抽球的概率是0.4。以这个概率抽一次球后，将球放回。然后从当前盒子转移到下一个盒子进行抽球。

规则是：如果当前抽球的盒子是第一个盒子，则以0.5的概率仍然留在第一个盒子继续抽球，以0.2的概率去第二个盒子抽球，以0.3的概率去第三个盒子抽球。如果当前抽球的盒子是第二个盒子，则以0.5的概率仍然留在第二个盒子继续抽球，以0.3的概率去第一个盒子抽球，以0.2的概率去第三个盒子抽球。如果当前抽球的盒子是第三个盒子，则以0.5的概率仍然留在第三个盒子继续抽球，以0.2的概率去第一个盒子抽球，以0.3的概率去第二个盒子抽球。如此下去，直到重复三次，得到一个球的颜色的观测序列:

O={红，白，红}

注意在这个过程中，观察者只能看到球的颜色序列，却不能看到球是从哪个盒子里取出的。

那么按照我们上一节HMM模型的定义，我们的观察集合是:

V={红，白}，M=2

我们的状态集合是：

Q={盒子1，盒子2，盒子3}，N=3

而观察序列和状态序列的长度为3.

初始状态分布为：

$\Pi=(0.2,0.4,0.4)^{T}$

状态转移概率分布矩阵为：

观测状态概率矩阵为：

1.4 HMM观测序列的生成

从上一节的例子，我们也可以抽象出HMM观测序列生成的过程。

1.5 HMM模型的三个基本问题

HMM模型一共有三个经典的问题需要解决：

1）评估观察序列概率。即给定模型 $\lambda=(A,B,\Pi)$ 和观察序列 $O=\left\{ o_{1},o_{2},...,o_{T} \right\}$ ，计算在模型 $\lambda$ 下观测序列 $O$ 出现的概率 $P(O|\lambda)$ 。这个问题的求解需要用到前向后向算法，这个问题是HMM模型三个问题中最简单的。

2）模型参数学习问题。即给定观测序列 $O=\left\{ o_{1},o_{2},...,o_{T} \right\}$ ，估计模型 $\lambda=(A,B,\Pi)$ 的参数，使该模型下观测序列的条件概率 $P(O|\lambda)$ 最大。这个问题的求解需要用到基于EM算法的鲍姆-韦尔奇算法，这个问题是HMM模型三个问题中最复杂的。

3）预测问题，也称为解码问题。即给定模型 $\lambda=(A,B,\Pi)$ 和观测序列 $O=\left\{ o_{1},o_{2},...,o_{T} \right\}$ ，求给定观测序列条件下，最可能出现的对应的状态序列，这个问题的求解需要用到基于动态规划的维特比算法，这个问题是HMM模型三个问题中复杂度居中的算法。

2 前向后向算法评估观察序列概率

2.1 回顾HMM问题一：求观测序列的概率

首先我们回顾下HMM模型的问题一。这个问题是这样的。我们已知HMM模型的参数 $\lambda=(A,B,\Pi)$ 。其中A是隐藏状态转移概率的矩阵，B是观测状态生成概率的矩阵， $\Pi$ 是隐藏状态的初始概率分布。同时我们也已经得到了观测序列 $O=\left\{ o_{1},o_{2},...,o_{T} \right\}$ ，现在我们要求观测序列O在模型λ下出现的条件概率 $P(O|\lambda)$ 。

乍一看，这个问题很简单。因为我们知道所有的隐藏状态之间的转移概率和所有从隐藏状态到观测状态生成概率，那么我们是可以暴力求解的。

虽然上述方法有效，但是如果我们的隐藏状态数 $N$ 非常多的那就麻烦了，此时我们预测状态有 $N^{T}$ 种组合，算法的时间复杂度是 $O(TN^{T})$ 阶的。因此对于一些隐藏状态数极少的模型，我们可以用暴力求解法来得到观测序列出现的概率，但是如果隐藏状态多，则上述算法太耗时，我们需要寻找其他简洁的算法。

前向后向算法就是来帮助我们在较低的时间复杂度情况下求解这个问题的。

2.2 用前向算法求HMM观测序列的概率

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

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

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

既然是动态规划，我们就要递推了，假设已经找到了在时刻t时各个隐藏状态的前向概率，现在需要递推出时刻t+1时各个隐藏状态的前向概率。

我们的动态规划从时刻1开始，到时刻 $T$ 结束，由于 $\alpha_{T}(i)$ 表示在时刻 $T$ 观测序列为 $o_{1},o_{2},...,o_{T}$ ，并且时刻 $T$ 隐藏状态 $q_{i}$ 的概率，只要将所有隐藏状态对应的概率相加，即 $\sum_{i=1}^{N}{\alpha_{T}(i)}$ 就得到了在时刻 $T$ 观测序列为 $o_{1},o_{2},...,o_{T}$ 的概率。

下面总结下前向算法。

从递推公式可以看出，我们的算法时间复杂度是 $O(TN^{2})$ ，比暴力解法的时间复杂度 $O(TN^{T})$ 少了几个数量级。

2.3 HMM前向算法求解实例

这里用盒子与球的例子来显示前向概率的计算。

首先计算时刻1三个状态的前向概率：

时刻1是红色球，隐藏状态是盒子1的概率为：

隐藏状态是盒子2的概率为：

隐藏状态是盒子3的概率为：

现在可以开始递推了，首先递推时刻2三个状态的前向概率：

时刻2是白色球，隐藏状态是盒子1的概率为：

隐藏状态是盒子2的概率为：

隐藏状态是盒子3的概率为：

继续递推，现在我们递推时刻3三个状态的前向概率：

时刻3是红色球，隐藏状态是盒子1的概率为：

隐藏状态是盒子2的概率为：

隐藏状态是盒子3的概率为：

2.4 用后向算法求HMM观测序列的概率

熟悉了用前向算法求HMM观测序列的概率，现在我们再来看看怎么用后向算法求HMM观测序列的概率。

后向算法和前向算法非常类似，都是用的动态规划，唯一的区别是选择的局部状态不同，后向算法用的是“后向概率”，那么后向概率是如何定义的呢？

这样我们得到了后向概率的递推关系式如下：

现在我们总结下后向算法的流程,注意下和前向算法的相同点和不同点：

此时我们的算法时间复杂度仍然是 $o(TN^{2})$ 。

2.5 HMM常用概率的计算

利用前向概率和后向概率，可以计算出HMM中单个状态和两个状态的概率公式。

上面这些常用的概率值在求解HMM问题二，即求解HMM模型参数的时候需要用到。

3 鲍姆-韦尔奇算法求解HMM参数

HMM模型参数求解的问题，这个问题在HMM三个问题里算是最复杂的。

3.1 HMM模型参数求解概述

HMM模型参数求解根据已知的条件可以分为两种情况。

可见第一种情况下求解模型还是很简单的。但是在很多时候，我们无法得到HMM样本观察序列对应的隐藏序列，只有 $D$ 个长度为 $T$ 的观测序列，即 $\left\{ \left( O_{1} \right)，\left( O_{2} \right),...,\left( O_{D} \right)\right\}$ 是已知的，此时我们能不能求出合适的HMM模型参数呢？这就是我们的第二种情况，也是我们本文要讨论的重点。

它的解法最常用的是鲍姆-韦尔奇算法，其实就是基于EM算法的求解，只不过鲍姆-韦尔奇算法出现的时代，EM算法还没有被抽象出来，所以我们本文还是说鲍姆-韦尔奇算法。

3.2 鲍姆-韦尔奇算法原理

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

3.3 鲍姆-韦尔奇算法的推导

利用前向概率的定义可得：

现在我们来看看A的迭代公式求法。方法和 $\Pi$ 的类似。由于A只在最大化函数式中括号里的第二部分出现，而这部分式子可以整理为：

现在我们来看看B的迭代公式求法。方法和 $\Pi$ 的类似。由于B只在最大化函数式中括号里的第三部分出现，而这部分式子可以整理为：

3.4 鲍姆-韦尔奇算法流程总结

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

以上就是鲍姆-韦尔奇算法的整个过程。

4 维特比算法解码隐藏状态序列

HMM模型最后一个问题的求解，即给定模型和观测序列，求给定观测序列条件下，最可能出现的对应的隐藏状态序列。

HMM模型的解码问题最常用的算法是维特比算法，当然也有其他的算法可以求解这个问题。同时维特比算法是一个通用的求序列最短路径的动态规划算法，也可以用于很多其他问题，比如用维特比算法来做分词。

本文关注于用维特比算法来解码HMM的的最可能隐藏状态序列。

4.1 HMM最可能隐藏状态序列求解概述

一个可能的近似解法是求出观测序列 $O$ 在每个时刻 $t$ 最可能的隐藏状态 $i_{t}^{*}$ ，然后得到一个近似的隐藏状态序列 $I^{*}=\left\{ i_{1}^{*},i_{2}^{*},...,i_{T}^{*} \right\}$ 。在给定模型 $\lambda$ 和观测序列 $O$ 时，在时刻 $t$ 处于状态 $q_{i}$ 的概率 $\gamma_{t}(i)$ ，这个概率可以通过HMM的前向算法与后向算法计算。这样我们有：

近似算法很简单，但是却不能保证预测的状态序列的整体是最可能的状态序列，因为预测的状态序列中某些相邻的隐藏状态可能存在转移概率为0的情况。

而维特比算法可以将HMM的状态序列作为一个整体来考虑，避免近似算法的问题，下面我们来看看维特比算法进行HMM解码的方法。

4.2 维特比算法概述

维特比算法是一个通用的解码算法，是基于动态规划的求序列最短路径的方法。

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

4.3 维特比算法流程总结

现在我们来总结下维特比算法的流程：

4.4 HMM维特比算法求解实例

下面用盒子与球的例子来看看HMM维特比算法求解。

按照我们上一节的维特比算法，首先需要得到三个隐藏状态在时刻1时对应的各自两个局部状态，此时观测状态为1：

现在开始递推三个隐藏状态在时刻2时对应的各自两个局部状态，此时观测状态为2：

4.5 HMM模型维特比算法总结

维特比算法的核心是定义动态规划的局部状态与局部递推公式，这一点在中文分词维特比算法和HMM的维特比算法是相同的，也是维特比算法的精华所在。

维特比算法也是寻找序列最短路径的一个通用方法，和dijkstra算法有些类似，但是dijkstra算法并没有使用动态规划，而是贪心算法。同时维特比算法仅仅局限于求序列最短路径，而dijkstra算法是通用的求最短路径的方法。

现在假设A和B在不同的房间里，我们看不到谁在掷骰子。取而代之的是，我们只知道后来吃了多少软糖。我们不知道颜色，仅是从罐子中取出的软糖的最终数量。我们怎么知道谁掷骰子？HMM。

在此示例中，状态是掷骰子的人，A或B。观察结果是该回合中吃了多少软糖。如果该值小于4，骰子的掷骰和通过骰子的条件就是转移概率。由于我们组成了这个示例，我们可以准确地计算出转移概率，即1/12。没有条件说转移概率必须相同，例如A掷骰子2时可以将骰子移交给他，例如，概率为1/36。

视频

R语言中的隐马尔可夫HMM模型实例

探索见解 ➜

去bilibili观看 ➜

探索更多视频 ➜

模拟

首先，我们将模拟该示例。B平均要吃12颗软糖，而A则需要4颗。

 
 
 
 
# 设置
simulate <- function(N, dice.val = 6, jbns, switch.val = 4){
  ＃模拟变量
    ＃可以只使用一个骰子样本
    ＃不同的机制，例如只丢1个骰子，或任何其他概率分布
  
    b<- sample(1:dice.val, N, replace = T) + sample(1:dice.val, N, replace = T)
    a <- sample(1:dice.val, N, replace = T) + sample(1:dice.val, N, replace = T)
    bob.jbns <- rpois(N, jbns[1])
    alice.jbns <- rpois(N, jbns[2])
    # 状态 
    draws <- data.frame(state = rep(NA, N), obs = rep(NA, N), 
    # 返回结果
    return(cbind(roll = 1:N, draws))
# 模拟场景
draws <- simulate(N, jbns = c(12, 4), switch.val = 4)
# 观察结果
ggplot(draws, aes(x = roll, y = obs)) + geom_line()

如您所见，仅检查一系列计数来确定谁掷骰子是困难的。我们将拟合HMM。由于我们正在处理计数数据，因此观察值是从泊松分布中得出的。

fit.hmm <- function(draws){
  # HMM 
  mod <- fit(obs ~ 1, data = draws, nstates = 2, family = poisson()
  # 通过估计后验来预测状态
  est.states <- posterior(fit.mod)
  head(est.states)
  # 结果
 
  hmm.post.df <- melt(est.states, measure.vars = 
  # 输出表格
  print(table(draws[,c("state", "est.state.labels")]))

## iteration 0 logLik: -346.2084 
## iteration 5 logLik: -274.2033 
## converged at iteration 7 with logLik: -274.2033 
##        est.state.labels
## state   alice bob
##   a   49   2
##   b    3  46

模型迅速收敛。使用后验概率，我们估计过程处于哪个状态，即谁拥有骰子，A或B。要具体回答该问题，我们需要更多地了解该过程。在这种情况下，我们知道A只喜欢黑软糖。否则，我们只能说该过程处于状态1或2。下图显示了HMM很好地拟合了数据并估计了隐藏状态。

# 绘图输出
 
 
    g0 <- (ggplot(model.output$draws, aes(x = roll, y = obs)) + geom_line() +
        theme(axis.ticks = element_blank(), axis.title.y = element_blank())) %>% ggplotGrob
    g1 <- (ggplot(model.output$draws, aes(x = roll, y = state, fill = state, col = state)) + 
      
 
 
    g0$widths <- g1$widths
    return(grid.arrange(g0, g1
 
 
plot.hmm.output(hmm1)

最受欢迎的见解

1.用R语言模拟混合制排队随机服务排队系统

2.R语言中使用排队论预测等待时间

3.R语言中实现马尔可夫链蒙特卡罗MCMC模型

4.R语言中的马尔科夫机制转换(Markov regime switching)模型

5.matlab贝叶斯隐马尔可夫hmm模型

6.用R语言模拟混合制排队随机服务排队系统

7.Python基于粒子群优化的投资组合优化

8.R语言马尔可夫转换模型研究交通伤亡人数事故预测

9.用机器学习识别不断变化的股市状况——隐马尔可夫模型的应用

令人印象深刻的是，该模型拟合数据和滤除噪声以估计状态的良好程度。公平地说，可以通过忽略时间分量并使用EM算法来估计状态。但是，由于我们知道数据形成一个序列，因为观察下一次发生的概率取决于前一个即\（P（X_t | X_ {t-1}）\），其中\（X_t \ ）是软糖的数量。

用机器学习识别不断变化的股市状况—隐马尔可夫模型(HMM)股票指数预测实战

阅读文章 ➜

考虑到我们构造的问题，这可能是一个相对简单的案例。如果转移概率大得多怎么办？

 
 simulate(100, jbns = c(12, 4), switch.val = 7)

## iteration 0 logLik: -354.2707 
## iteration 5 logLik: -282.4679 
## iteration 10 logLik: -282.3879 
## iteration 15 logLik: -282.3764 
## iteration 20 logLik: -282.3748 
## iteration 25 logLik: -282.3745 
## converged at iteration 30 with logLik: -282.3745 
##        est.state.labels
## state   alice bob
##   alice    54   2
##   bob       5  39

plot(hmm2)

这有很多噪音数据，但是HMM仍然做得很好。性能的提高部分归因于我们对从罐中取出的软糖数量的选择。分布越明显，模型就越容易拾取转移。公平地讲，我们可以计算中位数，并将所有低于中位数的值都归为一个状态，而将所有高于中位数的值归为另一状态，您可以从结果中看到它们做得很好。这是因为转移概率非常高，并且预计我们会从每个状态观察到相似数量的观察结果。