Ref：

1. https://web.stanford.edu/~jurafsky/slp3/A.pdf
2. https://en.wikipedia.org/wiki/Hidden_Markov_model

文中的例子和符号基本来自Ref[1]

基本概念：

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

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

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

数学表达：

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

wiki

概念

1. Markov assumption
   对隐藏序列（对应上图的X），有
2. Transition probability matrix
   状态转移矩阵，
   表示从隐藏状态i转移到隐藏状态j的概率，
3. Initial probability distribution
   初始概率分布，，表示隐藏序列从哪个状态开始；
   ，注意，有些状态可能为0；
4. 观测序列T
   
5. Emission probability
   发射概率， ，表示隐藏状态为i时，对应的观测状态为t的概率。

HMM:

2个假设：

1. Markov assumption： 隐藏状态i只和上一个隐藏状态i-1有关，
2. 观测状态仅依赖于当时的隐层状态，而和之前的观测状态or之前的隐层状态无关：

HMM可以看成3个问题：

1. Likelihood：已知观测序列O，参数 ，计算likelihood ；
2. Decoding：已知观测序列O，参数 ，找到最佳的隐藏序列Q，用Viterbi；
3. Learning：已知观测序列O，一系列的HMM状态，学习参数 ，用EM。

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

1. Likelihood computation：Forward algorithm

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

输入：观测序列O，模型参数
输出：

案例分析

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

1. 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；
2. 观测状态：Jason可能吃的冰淇淋数量，O = {1,2,3}。

1. 观测序列：吃了冰淇淋数量为[3,1,3]；
2. 气候序列：[hot,hot,cold]；
3. 冰淇淋数量O和气温高低的关系： ， ，
   
   https://web.stanford.edu/~jurafsky/slp3/A.pdf

• 计算条件概率 or likelihood：
  

• 但是实际上，我们并不知道真正的气候序列是怎样的，因此还要考虑实际的气候序列（即隐藏序列）：
  

• 因此，给定一条观测序列和一条隐藏序列的联合概率为：
  

• 而所有可能的隐藏序列对应的所有观测序列的总概率为：

如果有N个隐藏状态和T个观测序列，则可能有个可能的隐藏序列。
太大，因此无法分别计算每一个隐藏状态(N)下的观测(T)的likelihood。因此，使用Forward algorithm。

算法

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

设，previous forward path probability：
则，在第t个状态下，隐藏状态为j的所有可能路径概率和：

步骤：

1. 初始化：, given
2. 迭代：, given
3. 终止：

案例分析

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

如上图，，在时间t为2时，处于状态2，产生的观测序列为[3,1]，有两条路径通向这个点，分别为:

代码

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，参数 ，找到最佳的隐藏状态序列Q，也就是找到最大的likelihood。一般来说，可以通过forward algorithm来计算某个隐藏状态序列下的观测序列的likelihood，然后再选择likelihood最大的那条隐藏序列。

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

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

算法

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

设，上一个时间步的Viterbi路径概率previous Viterbi path probability:

则，时间步t时，隐藏状态为j的Viterbi值为：

步骤：

1. 初始化：
   对: ,
2. 迭代：
   对 :
   
   
3. 终止：
   best score:
   start of backtrace:

案例分析

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。
注意，这里的观测序列是无标签的，也就是说，我们知道有哪些观测值，但是观测值的排列未知。同样的，隐藏状态已知，分布未知。
因此，对于前文提到的案例，我们的观测序列，即吃的冰淇淋数量可能为：，以及知道存在的隐藏状态为{HOT,COLD}。但不知道状态转移概率A和发射概率B。

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

案例分析

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

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

则，初步判断气温冷热的初始化概率分布为：
然后，计算状态转移矩阵A：


发射概率B：

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

Baum-Welch 算法

设，时间t状态为i时，backward probability：
则，

步骤：

1. 初始化：,
2. 迭代：对,
   
3. 终止：

数学推理

1. 估计转移概率
   (从状态i转移到状态j的期望次数)/( 从状态i转移出去的期望次数)
   其中，分子可用联合概率表示：
   
   因为： 且
   所以：
   

2. 估计发射概率
   (隐藏状态为j时观测为k的期望次数)/(隐藏状态为j的期望次数)
   首先，要知道在时间步t时，状态为j的概率：
   
   因此：
   其中：
   ：对所有的时间步t，观测到的状态为k的总和。

EM算法

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

步骤：

1. 初始化：

2. 迭代直至收敛
   E-step：根据 ，计算期望
   
   

   M-step：根据期望，进行最大值估计，重新估计参数
   
   

3. 返回估计的参数。

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

补充

知识点

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

Viterbi 算法和Foward 算法异同：

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