隐马尔可夫HMM中viterbi算法

引言

viterbi算法简化最有可能的天气序列的运算过程，forward算法简化该该观察值的概率。

问题描述

你在中国，你朋友F在美国，F的作息有walk, shop, clean，但这选择跟天气有关，我们又知道Rainy的概率比Sunny的概率大
这是初始概率

这是天气转移矩阵

这是在相应天气下发生相应事的概率分布
然后，F这三天是walk，walk，shop，问最有可能的天气序列

问题分析

同样的，我们先用穷举法来算，即

Sunny Sunny Sunny 
Sunny Sunny Rainy
Sunny Rainy Sunny
Sunny Rainy Rainy
Rainy Sunny Sunny
Rainy Sunny Rainy
Rainy Rainy Sunny
Rainy Rainy Rainy

分 别求出在各个情况下，{walk，walk，shop}的概率，然后比较概率，最大的就是我们想要的结果，我们发现这跟我们forward的应用场景相 似，只不过，forward算法是求和，而viterbi算法是求最大，forward算法是所有路径之和，而viterbi算法只是求一条路径。这也就 意味着这两算法可以写在一起。

那我们重点讲讲算法吧，这次我是看wiki上的代码，再写，再比较，发现wiki上的代码有一处甚是费解，因此，加以修改，并一目了然了，我把我的思考罗列如下

1）我的思路是把第一天的情况特殊考虑，可wiki上的代码直接整合在一起了，其意义甚是费解，所以，我把第一天特殊考虑了

2）再者，我发现wiki上的代码，明明要求只需要确定三天的天气却出现四天的天气，估计是第一天没用的

3）我修改后代码跟wiki上代码，forward算法最终结果是一样的，但过程中求得值是不一样的

4）用viterbi求得概率是不同的，但wiki上的代码，把第一天去掉之后，剩下的天气序列跟我的是一样

不管怎样，wiki上的代码整体框架还是非常好的（或许还是我改错了），思路如下：

1）我上一篇自己写forward算法，差不多用二维数组来记录整个过程，其实没必要，只需一维足以，用完一行，替代一行，所以才会有字典T和字典U

2）返回三样，分别是该观察序列的最大概率（即之前用forward求出来的），最可能出现的天气序列，该天气序列的概率，同时每次都更新这三个值

3）三重for循环，第一重是第几天，第二重是第几天晴天还是雨天，第三重是前一天是晴天还是雨天下求得概率，再具体求和，求最大值

我稍加修改的Python代码及结果

def forward_viterbi(obs, states, start_p, trans_p, emit_p):
    T = {}
    #for the first one
    for state in states:
        T[state] = (start_p[state] * emit_p[state][obs[0]], [state], start_p[state])
    for output in obs:
        U = {}
        #pass the first one
        if output == obs[0]:
            continue
        for next_state in states:
            total = 0
            argmax = None
            valmax = 0
            for source_state in states:
                #get the previous one
                (prob, v_path, v_prob) = T[source_state]
                #caculate the commone one
                p = trans_p[source_state][next_state] * emit_p[next_state][output]
                prob *= p
                v_prob *= p
                #sum up
                total += prob
                #find the max
                if v_prob > valmax:
                    argmax = v_path + [next_state]
                    valmax = v_prob
            U[next_state] = (total, argmax, valmax)
        T = U
    #sum up & find the max one
    total = 0
    argmax = None
    valmax = 0
    for state in states:
        (prob, v_path, v_prob) = T[state]
        total += prob
        if v_prob > valmax:
            valmax = v_prob
            argmax = v_path
    return (total, argmax, valmax)
states = ('Rainy', 'Sunny')
observations = ('walk', 'shop', 'clean')
start_probability = {'Rainy': 0.6, 'Sunny': 0.4}
transition_probability = {
   'Rainy' : {'Rainy': 0.7, 'Sunny': 0.3},
   'Sunny' : {'Rainy': 0.4, 'Sunny': 0.6},
   }
emission_probability = {
   'Rainy' : {'walk': 0.1, 'shop': 0.4, 'clean': 0.5},
   'Sunny' : {'walk': 0.6, 'shop': 0.3, 'clean': 0.1},
   }
#A simple example of the using algorithm
def example():
    return forward_viterbi(observations,
                           states,
                           start_probability,
                           transition_probability,
                           emission_probability)
print example()

#(0.033611999999999996, ['Rainy', 'Rainy', 'Rainy'], 0.05879999999999999)

wiki代码及结果

def forward_viterbi(obs, states, start_p, trans_p, emit_p):
    T = {}
    for state in states:
        ##          prob.           V. path  V. prob.
        T[state] = (start_p[state], [state], start_p[state])
    for output in obs:
        print T
        U = {}
        for next_state in states:
            total = 0
            argmax = None
            valmax = 0
            for source_state in states:
                (prob, v_path, v_prob) = T[source_state]
                p = emit_p[source_state][output] * trans_p[source_state][next_state]
                #caculate
                prob *= p
                v_prob *= p
                #the different way to gain
                total += prob
                if v_prob > valmax:
                    argmax = v_path + [next_state]
                    valmax = v_prob
            U[next_state] = (total, argmax, valmax)
        T = U
    print T
    ## apply sum/max to the final states:
    total = 0
    argmax = None
    valmax = 0
    for state in states:
        (prob, v_path, v_prob) = T[state]
        total += prob
        if v_prob > valmax:
            argmax = v_path
            valmax = v_prob
    return (total, argmax, valmax)
   
states = ('Rainy', 'Sunny')
observations = ('walk', 'shop', 'clean')
start_probability = {'Rainy': 0.6, 'Sunny': 0.4}
transition_probability = {
   'Rainy' : {'Rainy': 0.7, 'Sunny': 0.3},
   'Sunny' : {'Rainy': 0.4, 'Sunny': 0.6},
   }
emission_probability = {
   'Rainy' : {'walk': 0.1, 'shop': 0.4, 'clean': 0.5},
   'Sunny' : {'walk': 0.6, 'shop': 0.3, 'clean': 0.1},
   }
#A simple example of the using algorithm
def example():
    return forward_viterbi(observations,
                           states,
                           start_probability,
                           transition_probability,
                           emission_probability)
print example()

#(0.033611999999999996, ['Sunny', 'Rainy', 'Rainy', 'Rainy'], 0.009407999999999998)