隐马尔可夫模型(HMM)

隐马尔科夫模型(Hidden Markov Model,HMM)是离散状态的动态模型，用于对时间序列的分析研究，HMM属于生成模型在语音识别等领域有着广泛的应用。目录1.隐马尔可夫模型的基本概念1.1 隐马尔科夫模型定义及三要素1.1.1 HMM定义1.1.2 HMM三要素1.2隐马尔可夫模型的3个基本问题2.求值问题2.1 直接求值法2.2 前向算法2.3...

Tonywu2018

1767人浏览 · 2019-01-15 15:03:37

Tonywu2018 · 2019-01-15 15:03:37 发布

隐马尔科夫模型(Hidden Markov Model,HMM)是离散状态的动态模型，用于对时间序列的分析研究，HMM属于生成模型在语音识别等领域有着广泛的应用。

1.隐马尔可夫模型的基本概念

1.1 隐马尔科夫模型定义及三要素

说到一个技术之前总是先说这个技术的概念，关于HMM的定义，引用李航博士的《统计学习方法》中的定义。

1.1.1 HMM定义

定义：隐马尔可夫模型是关于时序的概率模型，描述由一个隐藏的马尔可夫链随机生成不可观测的状态随机序列，再由各个状态生成一个观测而产生观测随机序列的过程。状态序列(state sequence)用S表示，观测序列(observation sequence)用O表示。

隐马尔科夫模型，我们从名字上就能看出它的特点，即：符合马尔可夫过程、观测序列由隐状态产生。

有关马尔科夫过程有两个假设：

（1）、齐次马尔科夫性假设。即马尔可夫链在任意时刻t的状态只与前一时刻t-1的状态有关，与其他时刻的状态和观测都无关，即有: $P(s_{t}|s_{t-1},o_{t-1},...,s_{1},o_{1}) = P(s_{t}|s_{t-1}), t = 1,2,...,T$ (1.1)

（2）、观测独立性假设。即任意时刻的观测只依赖于该时刻的马尔可夫链的状态，与其他观测及状态无关。

$P(o_{t}|s_{T},o_{T},...,s_{t+1},o_{t+1},s_{t},s_{t-1},o_{t-1},...,s_{1},o_{1}) = P(o_{t}|s_{t})$ (1.2)

1.1.2 HMM三要素

如上图所示，是一个股市的简易图。其中Bull、Bear和Even分别表示牛市、熊市和平仓，它们表示股市的三种状态。从图中可以看到当前状态为Bull时，下一状态还是为Bull的概率为0.6，下一状态为Even的概率为0.2，下一概率为Bear的概率为0.2；同样的，当前状态为Bear时，下一状态还是Bear的概率为0.3，下一状态为Bull的概率为0.5，下一转态为Even的概率为0.2；但前状态为Even时，下一状态还是Even的概率为0.5，下一状态为Bull的概率为0.4，下一状态为Bear的概率为0.1。我们称这些概率为转移概率，用状态转移概率矩阵A表示为：

矩阵中的元素表示概率，以 $a_{11} = P(s_{t+1} = Bull|s_{t} = Bull) = 0.6$ 为例，0.6就表示在时刻t处于状态Bull的条件下时刻t+1转移到Bull的概率为0.6；同样的 $a_{23} = P(s_{t+1} = Even|s_{t} = Bear) = 0.2$ 表示在时刻t处于状态Bear的条件下时刻t+1转移到Even的概率为0.2。可以看到在状态转移概率矩阵A中，矩阵的每一行上的元素相加的和为1，这也很容易理解。以A矩阵的第一行为例，行和为1就是表示在t时刻的状态为Bull，在t+1时刻状态要么是Bull要么是Bear要么是Even，所以这三个转移概率的和为1。还需要注意的一点是状态转移概率矩阵A是k*k的方阵，其中k表示状态数。

对于每一个状态都有三个虚线箭头指向的三个观测值，分别是：up、down和unchanged。以Bull状态为例，0.7表示在此状态下，股市上升的概率为0.7；0.1表示在此状态下，股市下降的概率为0.1；0.2表示在此状态下，股市不变的概率为0.2。我们把这些概率称为观测概率，把所有的观测概率用观测概率矩阵B表示为：

以矩阵B中的元素 $b_{32} = P(o_{t} = down|s_{t} = Even) = 0.3$ 为例，它表示在t时刻处于状态Even的条件下生成观测down的概率。同样的矩阵B的每一行上的元素和也为1；观测概率矩阵B是k*l的矩阵，k表示状态数，l表示观测数。

初始状态概率向量用 $\pi$ 表示， $\large \pi = (\pi_{i})$ ,其中， $\large \pi_{i} = P(s_{1} = q_{i}),i = 1,2,3,...,k$ 是时刻t=1处于状态 $q_{i}$ 的概率,初始状态概率向量的各元素之和为1。

初始状态概率向量 $\pi$ 、状态转移概率矩阵A以及观测概率矩阵B称为隐马尔科夫模型的三要素，用 $\lambda = (A,B,\pi)$ 表示。

状态转移概率矩阵A与初始状态概率向量 $\pi$ 确定了隐藏的马尔可夫链，生成不可观测的状态序列。观测概率矩阵B确定了如何从状态生成观测，与状态序列综合确定了如何产生观测序列。

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

(1) 概率计算问题，又称求值问题。给定模型 $\lambda = (A,B,\pi)$ 和观测序列 $O = (o_{1},o_{2},...,o_{T})$ ,计算在模型 $\lambda$ 下观测序列O出现的概率 $P(O|\lambda )$ 。

(2) 学习问题，又称模型的训练。已知观测序列 $O = (o_{1},o_{2},...,o_{T})$ ，估计模型 $\lambda = (A,B,\pi)$ 参数，使得在该模型下观测序列概率 $P(O|\lambda )$ 最大，可用极大似然估计的方法估计参数。 $\large \lambda_{MLE} = argmaxP(O|\lambda)$ 。

(3) 预测问题，又称解码问题(decoding)。已知模型 $\lambda = (A,B,\pi)$ 和观测序列 $O = (o_{1},o_{2},...,o_{T})$ ，求对给定观测序列条件概率 P(S|O) 最大的状态序列 $S = (s_{i},s_{2},...,s_{T})$ 。即给定观测序列，求最有可能对应的状态序列。

2.求值问题

求值问题就是在给定模型参数 $\lambda = (A,B,\pi)$ 和观测序列 $O = (o_{1},o_{2},...,o_{T})$ 的情况下，计算概率 $P(O|\lambda )$ 。

2.1 直接求值法

我们可以直接对其进行计算。

$P(O|\lambda) & = \sum_{S}P(O,S|\lambda)\\ = \sum_{s_{1} = 1}^{k},...,\sum_{s_{T} = 1}^{k}P(o_{1},o_{2},...,o_{T};s_{1},s_{2},...,s_{T}|\lambda)\\ = \sum_{s_{1} = 1}^{k},...,\sum_{s_{T} = 1}^{k}P(s_{1}|\lambda)P(o_{1}|s_{1},\lambda)P(s_{2}|s_{1},\lambda),...,P(s_{T}|s_{T-1},\lambda)P(o_{T}|s_{T},\lambda)\\ =\sum_{s_{1} = 1}^{k},...,\sum_{s_{T} = 1}^{k}\pi_{s_{i}}\prod_{t = 2}^{T}a_{s_{t-1},s_{t}}b_{s_{t}}(o_{t})$ (2.1)

其中 $p(s_{t} = j|s_{t-1} = i) = a_{ij}$ , $p(o_{t}|s_{t} = j) = b_{j}(o_{t})$ ,由上式可知。我们可以通过 $\lambda = (A,B,\pi)$ 和 $O = (o_{1},o_{2},...,o_{T})$ 求出 $P(O|\lambda )$ 接下来对上面的式子进行简要的说明。

这种直接求解的方法在理论上是可行的，但是它的计算量很大，一共需要做 $k^{T}$ 次加法运算，T表示时间序列的时间点数，k表示状态数。因此引出了前向-后向算法。

2.2 前向算法

前向概率：给定HMM参数 $\lambda$ ，定义到时刻t部分观测序列为 $o_{1},o_{2},...,o_{t}$ 且状态为 $q_{i}$ 的概率为前向概率，记作： $\alpha_{i}(t) = P(o_{1},o_{2},...,o_{t},s_{t} = q_{i}|\lambda)$ (2.2)

取t = 1,得到： $\alpha_{i}(1) = P(o_{1},s_{1} = q_{i}|\lambda) = P(s_{1} =q_{i})P(o_{1}|s_{1} =q_{i},\lambda) = \pi_{i}b_{i}(o_{1})$ (2.3)

取t = 2,得到：

递归得：

$\alpha_{j}(T) = b_{j}(o_{T})[\sum_{i = 1}^{k}\alpha_{i}(T-1)a_{ij}]\\ =P(o_{1},o_{2},...,o_{T},s_{T} =q_{j} |\lambda)\\$ (2.4)

$P(o_{1},o_{2},...,o_{T},s_{T} = q_{j}|\lambda)$ 是 $o_{1},o_{2},...,o_{T}$ 和 $s_{T}$ 的联合概率，所以对 $s_{T}$ 累加得到：

$P(O|\lambda) = P(o_{1},o_{2},...,o_{T}|\lambda)\\ = \sum_{j = 1}^{k}P(o_{1},o_{2},...,o_{T},s_{T} = q_{j}|\lambda)\\ = \sum_{j = 1}^{k}\alpha_j(T)$ (2.5)

前向计算要比直接计算高效很多，其关键是局部计算前向概率，然后递推到全局。可以看到每一个时刻的前向概率只需要做k次加法运算，而下一时刻的计算直接引用前一时刻的计算结果，一共有T个时间点的话，那么用前向算法计算概率 $P(O|\lambda )$ 需要kT次加法运算，这显然比 $k^{T}$ 小很多很多。

2.3 后向算法

后向概率：给定HMM参数 $\lambda$ ,定义在时刻t状态为 $q_{i}$ 的条件下，从t+1到T的部分观测序列为 $o_{t+1},o_{t+2},...,o_{T}$ 的概率为后向概率，记作： $\beta_{i}(t) = P(o_{t+1},o_{t+2},...,o_{T}|s_{t} = q_{i},\lambda)$ (2.6)

规定 $\beta_{i}(T) = 1$

取t = T-1得：

取t = T-2得

递归得： $\beta_{i}(t) = \sum_{j = 1}^{k}a_{ij}b_{j}(o_{t+1})\beta_{j}(t+1)$ (2.7)

最后得到： $P(O|\lambda) = \sum_{i = 1}^{k}\beta_{i}(1)\pi_{i}b_{i}(o_{1})$ (2.8)

2.4 前向-后向算法

$\alpha_{i}(t) = P(o_{1},o_{2},...,o_{t},s_{t} = q_{i}|\lambda)$ $\beta_{i}(t) = P(o_{t+1},o_{t+2},...,o_{T}|s_{t} = q_{i},\lambda)$

在求值时，可以用前向算法也可以用后向算法，还可以同时使用两种算法，下面举例说明一下。

2.4.1 单个状态概率计算公式

给定模型参数 $\lambda = (A,B,\pi)$ 和观测序列O，求在t时刻处于状态 $q_{i}$ 的概率： $P(s_{t} = q_{i}|O,\lambda)$ ,对此概率用贝叶斯公式转换一下得到：

$\large P(s_{t} = q_{i}|O,\lambda) = \frac{P(s_{t} = q_{i},O|\lambda)}{P(O|\lambda)} = \frac{P(s_{t} = q_{i},O|\lambda)}{\sum_{j= 1}^{k}P(O,s_{t} = q_{j}|\lambda)}$ (2.9) 接下来求解概率 $P(O,s_{t} = q_{i}|\lambda)$ ：

同样的得到 $\sum_{j = 1}^{k}P(O,s_{t} = q_{j}|\lambda) = \sum_{j = 1}^{k}\alpha_{j}(t)\beta_{j}(t)$ ,所以t时刻处于状态 $q_{i}$ 的概率为：

$\large P(s_{t} = q_{i}|O,\lambda) = \frac{\alpha_{i}(t)\beta_{i}(t)}{\sum_{j = 1}^N\alpha_{j}(t)\beta_{j}(t)}$ (2.10)

3.学习问题

HMM的学习问题就是求得参数 $\lambda = (A,B,\pi)$ 的问题，HMM的学习，根据训练数据是包括观测序列和对应的状态序列还是只有观测序列，可以分为监督学习和非监督学习。由于监督学习的方法需要人工标注数据，代价较高，所以使用的比较少，且监督学习的算法比较简单，所以下面只分析非监督学习的情况，即对只给出观测序列的数据进行训练。

3.1 EM算法

训练数据只包含观测序列 $O_{1},O_{2},...,O_{T}$ ，而没有状态序列，这时我们可以吧状态序列看成隐变量S,那么HMM是一个含有隐变量的概率模型： $P(O|\lambda) = \sum_{S}P(O|S,\lambda)P(S|\lambda)$ (3.1)

因为模型中含有隐变量，所以训练HMM时可使用EM算法迭代求解参数。

EM算法迭代求解参数步骤：

1、选择参数的初始值，得到模型 $\lambda^{(0)} = (A^{(0)},B^{(0)},\pi^{(0)})$

2、E步：E步是求期望，那么求得是谁的期望呢。求得是联合分布 $P(O,S|\lambda)$ 的对数似然函数以 $P(S|O,\lambda^{(i)})$ 为分布的期望，即 $E_{P(S|O,\lambda^{(i)})}(logP(O,S|\lambda)) = \sum_S logP(O,S|\lambda))P(S|O,\lambda^{(i)})$ ，这里需要对期望进行简单的修改：

修改后的期望定义为 $Q(\lambda,\lambda^{(i)})$ ，其中 $\large \lambda$ 表示要极大化的HMM参数， $\lambda^{(i)}$ 表示当前HMM参数的估计值。由(2.1)式得

$P(O,S|\lambda) =\pi_{s_{i}}\prod_{t = 2}^{T}a_{s_{t-1},s_{t}}b_{s_{t}}(o_{t})$ （3.2）

将其带入函数 $Q(\lambda,\lambda^{(i)})$ 得：

$Q(\lambda,\lambda^{(i)}) = \sum_{S}log(\pi_{s_{i}}\prod_{t = 2}^{T}a_{s_{t-1},s_{t}}b_{s_{t}}(o_{t}))P(O,S|\lambda^{(i)})\\ = \sum_{S}log\pi_{s_{i}}P(O,S|\lambda^{(i)}) + \sum_{S}(\sum_{t = 2}^{T}loga_{s_{t-1},s_{t}})P(O,S|\lambda^{(i)}) + \sum_{S}(\sum_{t = 2}^{T}logb_{s_{t}}(t))P(O,S|\lambda^{(i)})$ （3.3）

3、M步：最大化 $Q(\lambda,\lambda^{(i)})$ 求模型参数 $\lambda = (A,B,\pi)$

在式(3.3)中可以看到， $Q(\lambda,\lambda^{(i)})$ 函数分成了三个部分，可以分别求三个部分的最大值从而得到模型的参数 $\lambda = (A,B,\pi)$ 。最大化的过程中会用到拉格朗日乘子法。

4、终止程序。得到模型参数 $\lambda^{(n+1)} = (A^{(n+1)},B^{(n+1)},\pi^{(n+1)})$

4.预测问题

预测问题就是给定模型 $\lambda = (A,B,\pi)$ 和观测序列 $O = (o_{1},o_{2},...,o_{T})$ ，求对给定观测序列条件概率 P(S|O) 最大的状态序列 $S = (s_{i},s_{2},...,s_{T})$ 。即给定观测序列，求最有可能对应的状态序列。

4.1 维特比算法

维特比算法是一种动态规划的算法，有关动态规划的解释可以查看参考文献3，文章用漫画的形式解释了动态规划，简单易懂。以维基百科的感冒发烧的例子进行说明。

假设村民的身体特质是理想的，要么健康要么发烧。假设某诊所的医生通过询问病人的感觉来诊断病人是否发烧，而村民只需要告诉医生感觉正常、头晕或冷。病人的状态有两种：‘健康’和‘发烧’，这个是医生无法直接观察到的，因此把它看做是隐状态；而村民告诉医生的感觉：正常、头晕和发烧是由病人状态决定的，这些看做是观测结果。那么这就是一个HMM的例子。

现假设有如下的HMM状态链。用参数表示为：

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},
   }

Graphical representation of the given HMM

假设病人连续三天看医生，医生发现第一天他感觉正常，第二天感觉冷，第三天感觉头晕。于是医生产生了一个问题：怎样的健康状态序列最能够解释这些观察结果。这个问题就可以用维特比算法解决。

def viterbi(obs, states, start_p, trans_p, emit_p):
    V = [{}]
    path = {}

    # Initialize base cases (t == 0)
    for y in states:
        V[0][y] = start_p[y] * emit_p[y][obs[0]]
        path[y] = [y]

    # Run Viterbi for t > 0
    for t in range(1,len(obs)):
        V.append({})
        newpath = {}

        for y in states:
            (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]

        # Don't need to remember the old paths
        path = newpath

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

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

最后结果为：(0.01512, ['Healthy', 'Healthy', 'Fever'])。也就是说病人第一天健康，第二天也健康，第三天发烧产生上面的观测结果的概率最大，为0.01512。

Viterbi算法
输入：	模型 $\large \lambda = (A,B,\pi)$ ；观测数据 $\large O = o_{1},o_{2},...,o_{T}$
输出：	最优路径 $\large S^{\star} = (s_{1}^{\star},...,s_{T}^{\star})$
(1)、初始化 $\large \delta_{i}(1) = \pi_{i}b_{i}(o_{1})$ $\large i = 1,2,...,N$ $\large \psi_{i}(1) = 0$ $\large i = 1,2,...,N$ (2)、递推。对 $\large t = 2,3,...,T$ $\large \delta_{i}(t) =\mathop{\max}_{1\leq j\leq N}[\delta_{j}(t-1)a_{ji}]b_{i}(o_{t})$ $\large i = 1,2,...,N$ $\large \psi_{i}(t) = \mathop{\arg\max}_{1 \leq j \leq N}[\delta_{j}(t-1)a_{ji}]$ $\large i = 1,2,...,N$ (3)、终止 $\large P^{\star} = \mathop{max}_{1 \leq i \leq N}\delta_{i}(T)$ $\large s_{T}^{\star} = \mathop{argmax}_{1 \leq i \leq N}[\delta_{i}(T)]$ (4)、最优路径回溯。对 $\large t = T-1,T-2,...,1$ $\large s^{\star}_{t} = \psi_{s^{\star}_{t+1}}(t+1)$ 求得最优路径 $\large S^{\star} = (s_{1}^{\star},...,s_{T}^{\star})$