HMM理论及代码实现

作者: 第一个读书笔记 | 来源:发表于2020-11-29 12:55 被阅读0次

Ref：

备注：文中的例子和符号基本来自https://web.stanford.edu/~jurafsky/slp3/A.pdf。

基本概念：

Markov chain：随机变量组成的序列，下一个序列状态仅和当前状态有关，而和过去的状态无关。好比预测明天的天气，仅考虑今天的气候，而不用管昨天或以前的天气情况。

Markov chain有助于计算可观测序列状态的概率，但很多时候，我们还关心那些无法直接观测到的状态序列。好比文本序列的pos tags，我们阅读文字，然后再判断这个文字是动词，名词还是其他，这些不可直接观测的状态，又称为隐藏状态（hidden state），HMM就是来帮助处理这种情况的。

Hidden Markov Model：有两种随机状态序列，一种状态可观测，另一种状态不可直接观测，这个不可观测的序列（下文用隐藏序列）服从Markov chain，且观测序列依赖隐藏序列，HMM希望通过这些观测序列来学习隐藏序列。

数学表达：

Y：可观测状态序列，下文的O
X：隐藏状态序列，hidden state，下文的Q
假设：Y依赖X，在每个时间t，仅依赖 $X_t$ ，和 $X_{<t}$ , $Y_{<t}$ 无关，HMM的状态之间的时序变化如下图：

wiki

关键点：

Markov assumption：对一系列的状态（对应上图的X，实际应用中，表示隐藏序列hidden state） $Q = q_1,q_2,...,q_N$ ，有 $P(q_i=a|q_1,q_2,...,q_{i-1}) = P(q_i=a|q_{i-1})$
Transition probability matrix ：转移矩阵， $A = a_{11}a_{12}...a_{n1}...a_{nn}$ ， $a_{ij}$ 表示从状态i转移到状态j的概率， $\sum_j a_{ij} =1$ ；
Initial probability distribution：初始概率分布， $\pi = \pi_1,\pi_2,...,\pi_N$ ，表示Markov chain从状态i开始， $\sum_i \pi_{i} =1$ ，注意，有些状态可能为0；
观测序列T： $O = o_1o_2,...o_T$ ；
Emission probability：发射概率， $B = b_i(o_t)$ ，表示从状态i产生的可观测状态 $o_t$ 的概率。

HMM:

2个假设：

Markov assumption： $P(q_i=a|q_1,q_2,...,q_{i-1}) = P(q_i=a|q_{i-1})$
观测状态输出仅依赖于当时的隐层状态，而和之前的观测状态or之前的隐层状态无关： $P(o_i|q_1,...q_i,...q_T,o_1,...,o_i,...,o_T) = P(o_i|q_i)$

HMM可以看成3个问题：

Likelihood：已知观测序列O，参数 $\lambda=(\pi,A,B)$ ，计算likelihood $P(O|\lambda)$ ；
Decoding：已知观测序列O，参数 $\lambda=(\pi,A,B)$ ，找到最佳的隐藏状态序列Q；
Learning：已知观测序列O，一系列的HMM状态，学习参数 $\lambda=(\pi,A,B)$ 。

案例分析：吃冰淇淋与气温的关系，如下图：

https://web.stanford.edu/~jurafsky/slp3/A.pdf

2种气候状态（不可观测）：HOT，COLD，对应的初始概率为0.8，0.2，即热天的概率为0.8，气温变冷的概率为0.2，从HOT转到HOT的概率为0.6，HOT转到COLD的概率为0.4，反之0.5，COLD转到COLD的概率为0.5；

观测状态：Jason可能吃的冰淇淋数量，O = {1,2,3}。

1. Likelihood computation：Forward algorithm

动态规划算法，计算观测的概率。
即，知道了一条观测序列，所有可能的隐藏序列，观测和隐藏序列之间的关系（转移概率A和发射概率B），初始隐藏序列的概率分布（π），求这个观测序列的likelihood。

输入：观测序列O，模型参数 $\lambda=(\pi,A,B)$
输出： $P(O|Q) = \Pi_i P(o_i|q_i)$

案例分析：

输入：

观测序列：吃了冰淇淋数量为[3,1,3]；

气候序列：[hot,hot,cold]；

冰淇淋数量O和气温高低的关系： $p(3|hot) = 0.4$ ， $p(1|hot) = 0.2$ ， $p(3|cold) = 0.1$
https://web.stanford.edu/~jurafsky/slp3/A.pdf
输出：计算Jason在给定输入情况下的likelihood。
$P(O|Q) = P(3,1,3|hot,hot,cold) = p(3|hot)p(1|hot)p(3|cold) = 0.4*0.2*0.1 = 0.008$

但是实际上，我们并不知道真正的气候序列是怎样的，因此还要考虑实际的气候序列（即隐藏序列）：
$P(Q) = p(hot|start)*p(hot|hot)*p(cold|hot) = 0.8*0.6*0.4 = 0.192$

因此，给定一条观测序列和一条隐藏序列的联合概率为：
$P(O,Q) = P(O|Q)P(Q) = \Pi_iP(o_i|q_i) × \Pi_iP(q_i|q_{i-1}) = 0.008*0.192 = 0.001536$

而所有可能的隐藏序列对应的所有观测序列的总概率为： $P(O) = \sum_QP(O,Q) = \sum_QP(O|Q)P(Q)$

对于有N个隐藏状态和T个观测序列的HMM来说，一共有 $N^T$ 个可能的隐藏序列。实际应用中， $N^T$ 太大，因此无法分别计算每一个隐藏状态(N)下的观测(T)的likelihood。因此，使用Forward algorithm。

forward algorithm 的3个关键点

previous forward path probability： $\alpha_{t-1}(i)$
$\alpha_t(j) = p(o_1,o_2,...,o_t,q_t = j|\lambda )$ ：表示观测序列为t时，隐藏状态为j的概率。
$\alpha_t(j) =\sum_i^N\alpha_{t-1}(i)a_{ij}b_j(o_t)$ ：在第t个状态下，隐藏状态为j的所有可能路径概率和。
transition probability: $a_{ij}$ ，从上一个隐藏状态 $q_i$ 转到当前状态 $q_j$ 的转移概率;
state observation likelihood: $b_j(o_t)$ ，基于当前状态j，观测到的现象 $o_t$ 的likelihood。

算法

输入：观测状态序列长度为T，隐藏状态N
输出：forward-prob

初始化： $\alpha_1(j) = \pi_jb_j(o_1)$ , given $1<=j<=N$
迭代： $\alpha_t(j) = \sum_i^N \alpha_{t-1}a_{ij}b_j(o_t)$ , given $1<=j<=N,1<=t<=T$
终止： $P(O|\lambda) = \sum_i^N\alpha_T(i)$

案例分析

https://web.stanford.edu/~jurafsky/slp3/A.pdf
如上图， $\alpha_2(2)$ ，在时间t为2时，处于状态2，产生的观测序列为[3,1]，有两条路径通向这个点，分别为:
$\alpha_1(1) ×P(H|C)×P(1|H) = 0.02*0.5*0.2 = 0.002$
$\alpha_1(2) ×P(H|H)×P(1|H) = 0.32*0.6*0.2 = 0.0384$
$\alpha_2(2) = \alpha_1(1)+\alpha_1(2)=0.0404$

代码

def forward(obs,states, start_p, emission_p,trans_p):
    """
    obs = ('ice_3','ice_1','ice_3')
    states = ('HOT','COLD')
    start_p = {'HOT':0.8,'COLD':0.2}
    emission_p = {'HOT':{'ice_1':0.2,'ice_2':0.4,'ice_3':0.4},
                    'COLD':{'ice_1':0.5,'ice_2':0.4,'ice_3':0.1}}
    trans_p = {'HOT': {'HOT':0.6,'COLD':0.4},
                'COLD': {'HOT':0.5,'COLD':0.5}}

    """
    fwd = [{}]
    for state in states:
        fwd[0][state] = start_p[state] * emission_p[state][obs[0]]
    
    for t in range(1,len(obs)):
        fwd.append({})
        for state in states:
            fwd[t][state] = sum((fwd[t-1][s]*trans_p[s][state]*emission_p[state][obs[t]]) for s in states)
    prob = sum(fwd[len(obs)-1][s] for s in states)
    return prob

2. Decoding：Viterbi algorithm

给定观测序列O，参数 $\lambda=(\pi,A,B)$ ，找到最佳的隐藏状态序列Q，也就是找到最大的likelihood。一般来说，可以通过forward algorithm来计算某个隐藏状态序列下的观测序列的likelihood，然后再选择likelihood最大的那条隐藏序列。

简而言之，就是每个时间步下，都选择最大的likelihood。

但是，正如之前提到的，序列状态多种多样，因此使用Viterbi algorithm。同前向算法，它也是个动态规划问题。

3个关键点

previous Viterbi path probability: $v_{t-1}(i)$ ，从上一个时间步的Viterbi路径概率；
transition probability: $a_{ij}$ ，从状态 $q_i$ 转移到状态 $q_j$ 的转移概率；
state obervation likelihood: $b_j(o_t)$ ，基于当前状态j观测到的状态 $o_t$ 的似然。

那么，时间步t时，隐藏状态为j的Viterbi值为：
$v_t(j) =max_{q_1,...,q_{t-1}} p(o_1,o_2,...,o_t,q_t = j|\lambda ) = max_i^N v_{t-1}(i)a_{ij}b_j(o_t)$

算法

输入：观测状态序列长度为T，隐藏状态N，构造一条viterbi路径概率矩阵[N,T]
输出：最佳路径，路径概率

初始化：
given $1<=j<=N$ :
$v_1(j) = \pi_jb_j(o_1)$ , $bt_1(j) = 0$
迭代：
given $1<=j<=N,1<=t<=T$ :
$v_t(j) = \max_i^N v_{t-1}a_{ij}b_j(o_t)$
$b_t(j) = argmax_i^N v_{t-1}a_{ij}b_j(o_t)$
终止：
best score: $P* = \max_i^Nv_T(i)$
start of backtrace: $qT* = argmax_i^Nv_T(i)$

案例分析

https://web.stanford.edu/~jurafsky/slp3/A.pdf

如上图，对观测序列[3,1,3]，可能存在的隐藏序列路径为：
START->H->H->H, START->H->C>H, START->H->C>C, START->H->H>C
SATRT ->C->C>C, START->C->H>C, START->C->H>J, START->C->C>H
Step t=1 时，观测状态为3，最大概率为0.32，最可能的路径为START->H；
Step t=2 时，观测状态为1，最大概率为0.064，最可能的路径为START->H->C；
Step t=3 时，观测状态为3，最大概率为0.0128，最可能的路径为START->H->C>H。
因此，viterbi的输出为：0.0128，START->H->C>H。
同时，viterbi还可以追溯，即H->C->H->START，如图上蓝色箭头。

代码

def viterbi(obs, states, start_p, emission_p, trans_p):
    V = [{}]
    path = {}
    # 建立t0时刻个状态概率
    for state in states:
        V[0][state] = start_p[state] * emission_p[state][obs[0]]
        path[state] = [state]

    # 沿着时间1, 2..t进行计算
    for t in range(1, len(obs)):
        V.append({})
        newpath = {}

        # 根据t-1时刻状态概率，观测概率矩阵和转移概率矩阵计算t时刻最大概率转态 记录路径
        for state in states:
            # 看前一个状态中, 那个状态转移到当前状态并且当前状态喷射出当前观测值的概率大， 谁大谁就做前一个状态
            (prob, next_state) = max([(V[t-1][s] * trans_p[s][state] * emission_p[state][obs[t]], s) for s in states])
            V[t][state] = prob
            newpath[state] = path[next_state] + [state]  # state状态概率最大: 前一个状态 和 当前能喷射出最大概率的状态
            # print(V)
            # print(newpath)
        path = newpath
    # 最后结果：
    # (0.0128, ['HOT', 'COLD', 'HOT'])
    (prob, next_state) = max([(V[len(obs) - 1][s], s) for s in states])
    return (prob, path[next_state])

3. Learning：EM algorithm

已知观测序列集O，一系列的HMM状态，学习参数A,B。
注意，这里的观测序列是无标签的，也就是说，我们知道有哪些观测值，但是观测值的排列未知。同样的，隐藏状态已知，分布未知。
因此，对于前文提到的案例，我们的观测序列，即吃的冰淇淋数量可能为： $O ={1,3,2,...}$ ，以及知道存在的隐藏状态为{HOT,COLD}。但不知道状态转移概率A和发射概率B。

EM（expectation-Maximization）algorithm：学习HMM的参数A和B。
EM是一种迭代算法，计算初始概率估计，然后使用这个估计去计算更好的估计，然后继续根据这个估计再计算更好的估计。。。。。。

举个简单例子：

假如我们已经知道某几天的气温以及冰淇淋的数量了，也就是说，已知：
ice cream count: temperature：

3->hot, 3->hot, 2->cold

1->cold, 1->cold, 2->cold

1->cold, 2->hot, 3->hot

则，初步判断气温冷热的初始化概率分布为： $\pi_{hot}=1/3,\pi_{cold}=2/3$
然后，计算状态转移矩阵A：
$p(hot|hot) =2/3 ,p(cold|hot)=1/3$
$p(cold|cold) = 2/3,p(hot|cold) = 1/3$
发射概率B：
$p(1|hot) = 0, p(2|hot) = 1/4, p(3|hot) = 3/4$
$p(1|cold) = 3/5, p(2|cold) = 1/5, p(3|cold) = 0$

如果我们并不清楚冷热气温背后的冰淇淋数量呢？
Forward-backward or Baum-Welch，EM中的一种特例算法，可以不断迭代来计算冰淇淋的数量。

Baum-Welch algorithm

$\beta$ :backward probability，在时间t对应状态为i时，时间t+1到最后观测的概率：
$\beta_t(i) = P(o_{t+1},o_{t+2},...,o_{T}|q_t=i,\lambda)$

初始化：given $1<=i<=N$ , $\beta_t(i) =1$
迭代：given $1<=i<=N,1<=t<=T$ ,
$\beta_t(i) = \sum_j^N a_{ij}b_j(o_{t+1})\beta_{t+1}(j)$
终止： $P(O|\lambda) = \sum_j^N\pi_jb_j(o_1)\beta_1(j)$

数学推理

估计转移概率：
$\hat a_{ij}=$ (从状态i转移到状态j的期望次数)/( 从状态i转移出去的期望次数)
分子可用联合概率表示： $\xi_t(i,j) = P(q_t=i,q_{t+1}=j|O,\lambda)$
根据 $P(X|Y,Z) =\frac {p(X,Y|Z)}{P(Y|Z)}$
又： $P(O|\lambda) = \sum_i^N \alpha_t(i)\beta_t(j)$
$\xi_t(i,j) = P(q_t=i,q_{t+1}=j,O|\lambda)$
因此： $\xi_t(i,j) = P(q_t=i,q_{t+1}=j|O, \lambda) = \frac {P(q_t=i,q_{t+1}=j,O|\lambda)} { \sum_i^N \alpha_t(i)\beta_t(j)}$
$\hat a_{ij} = \frac{\sum_{t=1}^{T} \xi_t(i,j)}{ \sum_{t=1}^T \sum_{k=1}^{N} \xi_t(i,k)}$

同理，估计发射概率：
$\hat b_j(v_k) =$ (隐藏状态为j时观测为k的期望次数)/(隐藏状态为j的期望次数)
首先，要知道在时间步t时，状态为j的概率：
$\gamma_t(j) = P(q_t=j|O,\lambda) = \frac {P(q_t=j,O|\lambda)}{P(O|\lambda)} = \frac {\alpha_t(j)\beta_t(j)}{P(O|\lambda)}$

因此： $\hat b_j(v_k) = \frac {\sum_{t=1,s.t.o_t = v_k}^T\gamma_t(j)}{ \sum_{t=1}^T\gamma_t(j)}$
其中： $\sum_{t=1,s.t.o_t = v_k}^T$ 表示：对所有的时间步t，每个时间步t，观测到的状态为k的总和。

EM算法

已经知道转移概率和发射概率的估计方式，通过EM算法来求解HMM的A和B吧。
输入：观测状态，序列长度为T，隐藏状态集Q
输出：HMM的 $\lambda=(\pi,A,B)$

初始化： $\lambda=(\pi,A,B)$
迭代直至收敛
E-step：根据 $\lambda=(\pi,A,B)$ ，计算期望
$\gamma_t(j) = \frac {\alpha_t(j)\beta_t(j)}{\alpha_T(q_F)}$
$\xi_t(i,j) =\frac {P(q_t=i,q_{t+1}=j,O|\lambda)} { \alpha_T(q_F)}$

M-step：根据期望，进行最大值估计，重新估计参数
$\hat a_{ij} = \frac {\sum_{t=1}^{T} \xi_t(i,j)}{ \sum_{t=1}^T \sum_{k=1}^{N} \xi_t(i,k)}$
$\hat b_j(v_k) = \frac {\sum_{t=1,s.t. o_t = v_k}^T\gamma_t(j)}{ \sum_{t=1}^T\gamma_t(j)}$
返回估计的参数。

EM算法的初始化很重要，设置的不好，收敛不好。一般在实际应用时候，会根据经验，手动初始化。

补充

知识点

HMM的参数：初始概率 $\pi$ ，转移概率矩阵A，发射概率B
如果参数已知，还知道观测序列，隐藏序列：
求某个时间步下观测序列的likelihood，用forward算法；
求某个时间步下观测序列对应的最可能的隐藏序列，用Viterbi算法；
如果参数未知，就EM走起。

Viterbi 算法和Foward 算法异同：

Viterbi计算的最可能的路径，求的是max，输出最可能的路径和得分；而Foward计算的likelihood，求的是sum；
Viterbi还多了个backpointers（Viterbi模块图示的蓝色箭头）。

HMM理论及代码实现

基本概念：

数学表达：

关键点：

HMM:

案例分析：吃冰淇淋与气温的关系，如下图：

1. Likelihood computation：Forward algorithm

案例分析：

forward algorithm 的3个关键点

算法

案例分析

代码

2. Decoding：Viterbi algorithm

3个关键点

算法

案例分析

代码

3. Learning：EM algorithm

Baum-Welch algorithm

数学推理

EM算法

补充

知识点

Viterbi 算法和Foward 算法异同：

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读