（系列笔记）17.HMM系列（3）

HMM——三个基本问题的计算

1、概率计算问题

我们要计算的是$O=(红宝石，珍珠，珊瑚)$的出现概率，我们要求的是$P(O|\lambda )$

直接计算

用直接计算法来求$\lambda$情况下长度为T的观测序列O的概率：
$P(O|\lambda )={\sum }_{S\in S_T}P(O,S|\lambda )$
其中$S_T$表示所有长度为 $T$ 的状态序列的集合，S为其中一个状态序列。对所有长度为T的状态序列$S_t$和观测序列O求以$\lambda$为条件的联合概率，然后对所有可能的状态序列求和，就得到了$P(O|\lambda )$的值。
因为： $P(O,S|\lambda )=P(O|S,\lambda )P(S|\lambda )$
又因为：$P(O,S|\lambda )=b_{11}{b}_{22}...b_{TT}$
所以：$P(O|\lambda )={\sum }_{s_1,s_2,...,s_T}{\pi }_1{b}_{11}{a}_{12}{b}_{22}{a}_{23}...a_{T-1}{b}_{TT}$
理论上讲，我们可以把所有状态序列都按照上述公式计算一遍，但它的时间复杂度是$O(T{N}^T)$计算量太大，不可行。那该如何计算一个观测序列出现的概率呢？用前向-后向算法。

前向-后向算法

前向-后向算法是一种动态规划算法，它分两条路景来计算观测序列概率，一条从前向后（前向），另一条从后向前（后向），这两条路径，都可以分别计算观测序列出现概率。在实际应用中，选择其中之一来计算就可以。所以前向-后向算法其实也可以被看作两个算法：前向算法和后向算法，它们可以分别用来求解$P(O|\lambda )$。

前向算法

设${\alpha }_t(i)=P(o_1,o_2,...,o_t,s_t=S_i|\lambda )$为前向概率。
它表示的是：给定$\lambda$的情况下，到时刻 t时，已经出现的观测序列为$o_1,o_2,...,o_t$且此时状态值为$S_i$的概率。

换句话说，现在是t时刻，到此刻为止，我们获得的观测序列是$(o_1,o_2,...,o_t)$，这些$(o_1,o_2,...,o_t)$的下标标志是时刻（而非观测空间的元素下标），具体某个$o_i$的取值才是观测空间中的一项，因此观测序列中可能出现$o_i=o_j$的状况。

$S_i$为一个具体的状态值，当前HMM的状态空间一共有N个状态值：$(S_1,S_2,...,S_N)$，$S_i)$是其中的一项。$s_t$在此处表示t时刻的状态，具体取值是$S_i$
因此：
（1）${\alpha }_1(i)={\pi }_i{b}_i{o}_1，i=1,2,...,N$；
（2）递推，对于$t=1,2,...,T-1$，有${\alpha }_{t+1}(i)=[{\sum }_{j=1}^N{\alpha }_t(k)a_{ji}]b_i(o_{t+1})，i=1,2,...N$;
（3）最终$P(O|\lambda )={\sum }_{i=1}^N{\alpha }_T(i)，i=1,2,...,N$。
如此一来，概率计算问题的计算复杂度就变成了$O(N^{2}T)$。

后向算法

设${\beta }_t(i)=P(o_1,o_2,...,o_t,s_t=S_i|\lambda )$为后向概率。
它表示的是：给定$\lambda$的情况下，到时刻 t时，状态为$S_i$的条件下，从$t+1$到$T$时刻出现的观察序列为$o_1,o_2,...,o_t$的概率。
（1）${\beta }_T(i)=1，i=1,2,...,N$，到了最终时刻，不论状态是什么，我们都规定当时的后向概率为1；
（2）对$t=T-1,T-2,..,1$，有${\beta }_t(i)={\sum }_{i=1}^N{a}_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j)，i=1,2,...,N$
（3）$P(O|\lambda )={\sum }_{i=1}^N{\pi }_i{b}_i(o_1){\beta }_1(i)，i=1,2,...,N$

结合前向后向算法，定义$P(O|\lambda )$：

2、预测算法

直接求解

预测算法是在$\lambda$既定，观测序列已知的情况下，找出最有可能产生如此观测序列的状态序列的算法。比如上面的例子中，真的出现了$O=(红宝石，珍珠，珊瑚)$，那么，最有可能导致O出现的S是什么？

比较直接的算法是：在每个时刻t，选择在该时刻最有可能出现的状态$s_t^{\ast }$，从而得到一个状态序列$S^{\ast }=(s_1^{\ast },s_2^{\ast },...,s_T^{\ast })$，直接就把它作为预测结果。

给定$\lambda$和观测序列O的情况下，t时刻处于状态$S_i$的概率为：

因为${\gamma }_t(i)=P(s_t=S_i|O,\lambda )=\frac{P(s_t=S_i,O|\lambda )}{P(O|\lambda )}$，通过前向概率${\alpha }_t(i)$和后向概率${\beta }_t(i)$定义可知：${\alpha }_t(i){\beta }_t(i)=P=(s_t=S_i,O|\lambda )$
所以有：

有了这个公式，我们完全可以求出t时刻，状态空间中每一个状态值$S_1,S_2,...,S_N$所对应的${\gamma }_t(i)$，然后选取其中概率最大的作为t时刻的状态即可。

这种预测方法在每一个点上都求最优，然后再“拼”成整个序列。它的有点是简单明了，缺点同样明显：不能保证状态序列对于观测序列的支持整体最优。

维特比算法

如何避免掉哪些实际上不会发生的状态序列，从而求得整体上的“最有可能”呢？
有一个很有名的算法：维特比算法——用动态规划求解概率最大路径。

维特比算法是求解HMM预测问题的经典算法，笔记中有不少地方会用到动态规划算法（前面的SMO），因为动态规划是一类相对独立且不算非常“基础”的算法，对它我们不再展开介绍

3、学习算法

HMM的学习算法根据训练数据不同，可以分为有监督学习和无监督学习两种。如果训练数据既包括观测序列，又包括对应的状态序列，且两者之间的对应关系已经明确标注了出来，那么就可以用有监督学习算法。否则，如果只有观测序列而没有明确对应的状态序列，就要用无监督学习算法训练。

有监督学习

训练数据是若干观测序列和对应的状态序列的样本对（Pair）,也就上说训练数据不仅包含观测序列，同时每个观测值对应的状态值是什么，也是已知的。

那么，最简单的，可以用频数估计频率。我们经过统计训练数据中所有的状态值和观测值，可以得到状态空间$(S_1,S_2,...,S_N)$和观测空间（$O_1,O_2,...,O_M$）。

结社样本中时刻t处于的状态为$S_i$，到了t+1时刻，状态属于$S_j$的频数为$A_{ij}$，那么状态转移概率$a_{ij}$的估计是：

初始状态概率${\pi }_i$的估计${\pi }_i^{\prime }$为所有样本中初始状态为$S_i$的频率。这样通过简单的统计估计，就得到了$\lambda =(A^{\prime },B^{\prime },{\pi }^{\prime })$。

无监督学习

当训练数据仅有观测序列（设观测序列为O），而没有与其对应的状态序列时，状态序列S实则处于隐藏状态。这种情况下，不能直接用频率来估计概率。

有专门的Baum-Welch算法，针对这种情况求$\lambda =(A,B,\pi )$

Baum-Welch算法利用了前向-后向算法，同时还是EM算法的一个特例。简单点形容，大致是一个嵌套了EM算法的前向后向算法。

4、HMM实例

用代码实现小马轮盘赌的问题：

    from __future__ import division
    import numpy as np

    from hmmlearn import hmm

    def calculateLikelyHood(model, X):
        score = model.score(np.atleast_2d(X).T)

        print "\n\n[CalculateLikelyHood]: "
        print "\nobservations:"
        for observation in list(map(lambda x: observations[x], X)):
            print " ", observation

        print "\nlikelyhood:", np.exp(score)

    def optimizeStates(model, X):
        Y = model.decode(np.atleast_2d(X).T)
        print"\n\n[OptimizeStates]:"
        print "\nobservations:"
        for observation in list(map(lambda x: observations[x], X)):
            print " ", observation

        print "\nstates:"
        for state in list(map(lambda x: states[x], Y[1])):
            print " ", state


    states = ["Gold", "Silver", "Bronze"]
    n_states = len(states)

    observations = ["Ruby", "Pearl", "Coral", "Sapphire"]
    n_observations = len(observations)

    start_probability = np.array([0.3, 0.3, 0.4])

    transition_probability = np.array([
        [0.1, 0.5, 0.4],
        [0.4, 0.2, 0.4],
        [0.5, 0.3, 0.2]
    ])

    emission_probability = np.array([
        [0.4, 0.2, 0.2, 0.2],
        [0.25, 0.25, 0.25, 0.25],
        [0.33, 0.33, 0.33, 0]
    ])


    model = hmm.MultinomialHMM(n_components=3)

    # 直接指定pi: startProbability, A: transmationProbability 和B: emissionProbability

    model.startprob_ = start_probability
    model.transmat_ = transition_probability
    model.emissionprob_ = emission_probability


    X1 = [0,1,2]

    calculateLikelyHood(model, X1)
    optimizeStates(model, X1)

    X2 = [0,0,0]

    calculateLikelyHood(model, X2)
    optimizeStates(model, X2)

result

    [CalculateLikelyHood]: 

    observations:
      Ruby
      Pearl
      Coral

    likelyhood: 0.021792431999999997


    [OptimizeStates]:

    observations:
      Ruby
      Pearl
      Coral

    states:
      Gold
      Silver
      Bronze


    [CalculateLikelyHood]: 

    observations:
      Ruby
      Ruby
      Ruby

    likelyhood: 0.03437683199999999


    [OptimizeStates]:

    observations:
      Ruby
      Ruby
      Ruby

    states:
      Bronze
      Gold
      Bronze