《马尔可夫》隐马尔可夫链模型的预测问题

Ddcc 2018年8月22日 14:12 1021916684@qq.com
HMM 预测问题 维特比算法

如需转载请注明出处:http://zczzxz.top,整理不易请谅解。

一、 预测问题

      给定模型λ=(A,B,π)和观测序列O={o1,o2,...,oT},计算在模型λ下观测序列O出现时,最有可能的隐状态序列S={i1,i2....iT}。在实际的应用上,我们更加的关注隐状态的变化情况,这使得预测问题变得非常重要了。

二、篱笆网络(Lattice)的最短路径问题

      寻找有向无环图(篱笆网络,各列的节点必经过一个)当中首尾节点之间的最短路径。

      首先,我们假设 s 到 e 之间存在一条最短路径(红色),且这条路径经过 c2 点,那么我们便一定能够确定从 s到 c2 的子路径当中找到最短路径。(证明:反证法。如果 s 到 c2 之间存在一条更短的子路径,那么便可以用它来代替原先的路径,而原先的路径显然就不是最短了,这与原假设自相矛盾)。同理,我们也可以得出从 s 到 b2 点为两点间最短子路径的结论。既然如此,我们计算从 s 点出发到点 c2 的最短路径,只要考虑从 s 出发到 b 列所有节点的最短路径(不确定时 b 列中的哪一个,所以记录 b 列中可能的节点); 我们计算从 s 点出发到点 b 的最短路径,只要考虑从 s 出发到 a 列所有节点的最短路径。

      然后, 标记所有可能的最短路径,以次假设 b1、 b2 和 b3 为最短路径中的节点。 计算路径的距离,当从 s->a->b->c-e 得到最短路径之后,通过回溯发找到每一列的确定节点,即可找到最短路径。

三、 维特比算法
      维特比算法是用动态规划解HMM的预测问题,根据动态规划原理,最优路径具有这样的特性:如果最优路径在时刻t通过节点it,那么这一路径从节点it到终点iT的部分路径,对于从it到iT的所有可能的部分路径来说,必须是最优的。假如不是这样,那么就有另一条更好的部分路径存在,形成一条比原来的路径更优的路径,这是矛盾的。依此,我们只需从时刻t=1开始,递推地计算在时刻t状态为i的各部分路径的最大概率,直至得到时刻t=T状态为i的各条路径的最大概率。时刻t=T的最大概率即为最优路径的概率P,最优路径的终点iT也同时得到。之后,为了找出最优路径的各个节点,从终点iT开始,由后向前逐步求得节点,得到最优路径I=(i1,i2,…,iT)。
      定义一个:
 
      代表t时刻,所有的路径中概率最大的路径(代表t时刻这一列的概率)。那么在t+1时刻有:
 
      定义在时刻t状态为i的所有单个路径中概率最大的路径的第t-1个节点为:
 
       最优路径回溯:
 
       为了便于理解:
 
       在概率计算中求的是所有指向q1的红线和作为q1出现的概率,现在不再求红线的和,而是求红线中概率最大值作为q1出现的概率,此时需要记录是那一条红线的概率最大,再计算所有指向q2的蓝线,再计算所有指向q3的棕线,最后这一列对应位置乘以O2

最优路径回溯的意思为找t+1这一列中最大的概率,假设对应q1,那么还需要再找哪根红线最大,这个作为最优路径。


class HMMForward(object):

    def __int__(self):
        pass

    def ForWardViterbi(self, A, B, O, pi):
        N, T = np.shape(A)[0], np.shape(O)[0]
        alpha = np.zeros((N, T))
        result_postion = np.zeros((1, T))
        result_postion_temp = np.zeros((N, T))
        for t in range(T):
            if 0 == t:
                alpha[:, t] = np.multiply(pi.T,  B[:, O[t]])
                result_postion[0, 0] = np.where(alpha[:, t] == np.max(alpha[:, t]))[0][0]  # 第一次的直接存进去,后面不再存了
            else:
                for k in range(N):
                    alpha_temp = np.multiply(alpha[:, t-1], A[:, k])
                    result_postion_temp[k, t] = np.where(alpha_temp == np.max(alpha_temp))[0][0]  #存放当前节点是从哪里来的(概率最大的路径),存放在第k列
                    alpha[k, t] = np.max(alpha_temp) * B[k, O[t]]    #存放当前节点发生的概率,存放在第k列
                    max = np.where(alpha[:, t] == np.max(alpha[:, t]))[0][0]
        print('alpha is:\n', alpha)
        print('result_postion_temp is:\n', result_postion_temp)
        # 回溯
        for t in range(T):
            n = T-t-1
            if t==0:
                result_postion[0, n] = max
            else:
                if 0 == n:
                    break
                else:
                    result_postion[0, n] = result_postion_temp[int(result_postion[0, n+1]), n]
        return result_postion

if __name__ == '__main__':
    A = np.array([[0.5, 0.2, 0.3],
                  [0.3, 0.5, 0.2],
                  [0.2, 0.3, 0.5]])
    # 隐马尔可夫模型B参数
    B = np.array([[0.5, 0.5],
                  [0.4, 0.6],
                  [0.7, 0.3]])
    # 初始状态概率分布pi
    pi = np.array([[0.2],
                   [0.4],
                   [0.4]])
    # 所有可能的观测序列
    O = [0, 1, 0]
    hmm = HMMForward()
    prob = hmm.ForWardViterbi(A, B, O, pi)
    print('\n观测序列:', O, '\n出现的最大概率隐状态为:', prob) 

四、代码验证 

       将李航的《统计学习方法》一个练习题带入程序验证程序的维特比算法。

调用函数

程序输出

统计学习方法的例题