隱馬爾可夫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)