HMM隐马尔科夫模型(附维特比代码）

背景知识：马尔科夫模型

1 马尔科夫的局限性

在一些情况下，我们并不能直接得到观测的结果，比如在天气系统中，我们不能直接得到天气的状态，但是我们有一堆蚂蚁，可以从蚂蚁的行为状态找到天气变化的关系规律。此时我们就有两组状态：

• 观察的状态（即蚂蚁的行为状态）
• 隐藏的状态（天气状态)

在这种情况下，我们希望设计出一个模型，在不能直接观测到天气状态的情况下，通过观测蚂蚁的状态，预测出下一天气状态。这就是隐马尔可夫模型（hidden Markov Model）的思想。

2 隐马尔可夫模型

1 介绍

隐马尔可夫模型（Hidden Markov Models, HMM）描述由一个隐藏的马尔可夫链随机生成的不可观测的状态序列，再由各个状态生成一个观测而产生观测随机序列的过程，属于生成模型（同属于生成模型的还有朴素贝叶斯）关于生成/判别模型，点这里查看更多介绍。

HMM(隐马尔可夫模型)是结构最简单的动态贝叶斯网，是一种著名的有向(无环)图模型，主要用于时序数据建模（语音识别、自然语言处理等）。

在马尔可夫模型中，状态对于观察者来说是直接可见的。唯一的参数就是状态转移概率。而在隐马尔可夫模型中，状态并不是直接可见的，但受状态影响的某些变量是可见的。每一个状态在可能输出的符号上都有一概率分布，因此输出符号的序列能够透露出状态序列的一些信息，由此形成了双重随机过程：

• 不仅仅是状态之间转移的随机
• 还是状态到观测变量之间的随机。

其难点是从可观察的参数中确定该过程的隐含参数，然后利用这些参数来作进一步的分析。

2 概率图表示

该图分为上下两行，上面那行就是一个隐含的马尔科夫转移过程，称作隐藏状态（如词性、天气）；下面这一行则是可以观察到的输出值，称为观察状态（词、蚂蚁行为情况）。

3 两个假设

1、齐次马尔可夫假设
隐藏的马尔可夫链在任意时刻t的状态只和上一时刻的状态有关，与其他状态和观测无关，也与时刻t无关。
2、观测独立性假设
任意时刻的观测只和当前时刻的状态有关。

4 举例

有三个不同的骰子，对应的面数和数值如下图：

接下来随机掷骰子，从三个骰子随机取出一个并掷骰，记录正面数值，10次之后，可能得到一个这样的序列（1 6 3 5 2 7 3 5 2 4），这个序列就是可观测序列，而相对的隐含的状态链就是选取骰子的顺序，是由D6，D4，D8的组合。

由于是随机选取，所以$D_i到D_j$的概率都是1/3，当然如果加入一定的规则，这个概率也是可以改变的。

同样的，我们知道各个状态到观测变量之间的概率，例如D6到1，2，3，4，5，6观测的概率都是1/6。同样也可以对这个概率进行定义。

这样一来，我们就得到了：
初始概率分布$\pi =(1/8,1/8,...,1/8)$
状态转移概率分布：$A=\left[\begin{array}{ccc}1/3& 1/3& 1/3\\ 1/3& 1/3& 1/3\\ 1/3& 1/3& 1/3\end{array}\right]$
概念转移概率分布：$B=\left[\begin{array}{cccccccc}1/6& 1/6& 1/6& 1/6& 1/6& 1/6& 0& 0\\ 1/4& 1/4& 1/4& 1/4& 0& 0& 0& 0\\ 1/8& 1/8& 1/8& 1/8& 1/8& 1/8& 1/8& 1/8\end{array}\right]$

以上就是一个简单的隐马尔可夫模型$\lambda =(A,B,\pi )$。

5 观测序列的生成

输入:一个HMM :$\lambda =(A,B,\pi )$
输出：观测序列 $O=(o_1,o_2,...,o_T)$

(1)按照初始分布$\pi ，产生i_1，t=1$
(2)按照状态$i_t$的观测概率分布$b_{i_t}(k)生成o_t$
(3)按照状态$i_t$的状态转移概率生成状态$i_{t+1}$
(4)令t = t+1，当t<T时，返回step2，否则，终止。

3 HMM的三个基本问题

1 概率计算问题

已知：

• HMM模型$\lambda =(A,B,\pi )$
• 观测序列 $O=(o_1,o_2,...,o_T)$

计算：

• 观测序列出现的概率，即：$P(O|\lambda )$

这个问题的意义更多在于验证我们的模型是否正确，通常使用前向算法和后向算法。

前向算法

定义前向概率：给定HMM模型$\lambda$,定义到时刻t的部分观测序列为$o_1,o_2,..,o_t$,且状态为$q_i$的概率为前向概率，记作：
$${\alpha }_t(i)=P(o_1,o_2,..,o_t,i_t=q_i|\lambda )$$
通过递推即可以得到观测序列$P(O|\lambda )$

---

算法：
Step1：初值：
$${\alpha }_1(i)={\pi }_i{b}_i(o_1)$$
Step2：递推，对于t=1，2，…,T
$${\alpha }_{t+1}(i)=[\sum {\alpha }_t(j)a_{ji}]b_i(o_{t+1}),i=1,2,3...N$$

Step3:终止
$$P(O|\lambda )=\sum {\alpha }_T(i)$$

后向算法

定义后概率：给定HMM模型$\lambda$,定义t到时刻T的部分观测序列为$o_{t+1},o_{t+2},..,o_T$,且状态为$q_i$的概率为前向概率，记作：
$${\beta }_t(i)=P(o_{t+1},o_{t+2},..,o_T|i_t=q_i,\lambda )$$
通过递推即可以得到观测序列$P(O|\lambda )$

---

后向计算算法

2 学习问题

使用极大似然方法进行参数估计，这是最常见的问题也是很重要的问题。
已知：

• 观测序列 $O=(o_1,o_2,...,o_T)$

估计：

• $\lambda =(A,B,\pi )$参数

监督学习算法

假设给定了S个长度相同的观测序列和对应的状态序列：$(O_1,I_1),(O_2,I_2),...,(O_S,I_S)$,可以使用伯努利大数定律来进行参数估计下,即频数等价于概率。

伯努利大数定理：设μ是n次试验中A发生的次数，一次试验中A发生的概率为 p ，则对任意正数 $\epsilon$有:

1）状态转移概率估计：根据频数统计概率
$${\alpha }_{ij}=A_{ij}/\sum _{j=1}^N{A}_{ij}$$

2）观测概率$b_j(k)$的估计
同样根据已知数据的频数来统计.
$${\beta }_j(k)=B_{jk}/\sum _{k=1}^M{B}_{jk}$$
3)初始状态概率${\pi }_i$的估计，就是初始状态的频率。

这就直接估计出参数。

无监督学习算法（Baum-Welch算法）

当给定的个数据只有观测数据$O=(o_1,o_2,...,o_S)$,而没有对应的状态序列时，这也是最常见的情况，那么隐马尔可夫模型就是一个含隐变量的概率模型。
$$P(O|\lambda )=\sum _{I}P(O|I,\lambda )P(I|\lambda )$$

Baum-Welch也就是EM算法，可以见这篇EM算法,我写的非常好理解EM。

3 预测（解码）问题

给定观测序列，求最可能对应的状态序列
已知：

• HMM模型$\lambda =(A,B,\pi )$
• 观测序列 $O=(o_1,o_2,...,o_T)$

计算：

• 给定观测序列O的条件概率$P(I|O)$最大的状态序列$I=(i_1,i_2,...,i_T)$

我们可以计算出所有中结果的可能性，那就是N！种预测结果，所以不切实际。

近似算法

在每个时刻都计算该时刻最可能出现的状态，从而得到一个预测状态序列作为输出结果。
优点在于：计算简单
缺点：不能保证每一时刻的最优就是全局的最优。
尽管如此，近似算法仍然是有效的。

维特比算法（Viterbi）

Viterbi就是动态规划算法用在HMM问题，本质上还是Dynamic Programming。
动态规划思想：对于一条最佳路径o1-ot，其部分路径oi-ot也是从i-t之间的最优路径。否则该部分就存在最优路径的替代，产生矛盾。
算法解读见维特比算法：

根据这一原理，可以快速高效的计算出最大的概率的路径，也就是待预测的状态转移序列。

Viterbi代码实现

Hidden_states = ('Healthy', 'Fever') # 隐状态集合
 
Observations_states = ('normal', 'cold', 'dizzy') # 观测状态集合
 
Start_probability = {'Healthy': 0.6, 'Fever': 0.4} # 表示病人第一次到访时医生认为其所处的HMM状态，他唯一知道的是病人倾向于是健康的（可以理解为这是基于一个对大众身体信息的了解得出的初始状态）
 
Hidden_transition_probability = { # 隐马尔可夫链中身体状态的状态转移概率，我们能够看到，当天健康的人，第二天有30%的概率会发烧
   'Healthy' : {'Healthy': 0.7, 'Fever': 0.3},
   'Fever' : {'Healthy': 0.4, 'Fever': 0.6},
   }
 
Hidden_observations_probability = {    # 原来叫emission_probability。这里表示病人每天的感觉的可能性。即，如果他是一个健康人，有50%的可能会感觉正常，40%觉得冷，10%觉得头晕
   'Healthy' : {'normal': 0.5, 'cold': 0.4, 'dizzy': 0.1},
   'Fever' : {'normal': 0.1, 'cold': 0.3, 'dizzy': 0.6},
   }

# Helps visualize the steps of Viterbi.
def print_dptable(V): # 打印dp矩阵
    print ("    "),
    for i in range(len(V)): 
        print("%7d" % i)
    print()

    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, h2o_p): # Viterbi算法
    V = [{}]
    path = {}

    # Initialize base cases (t == 0)
    for y in states:
        V[0][y] = start_p[y] * h2o_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] * h2o_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

    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,
                   Hidden_states,
                   Start_probability,
                   Hidden_transition_probability,
                   Hidden_observations_probability)
print(example())

Reference：

1.https://blog.csdn.net/pipisorry/article/details/50722178