Viterbi 算法 Python實現 [NLP學習一]

　　最近思考了一下將來，結合老師的意見，仍是決定挑一個方向開始研究了，雖然我的更喜歡鼓搗。深思熟慮後，結合本身的興趣點，選擇了NLP方向，感受比純粹的人工智能、大數據之類的方向有趣多了，我的仍是不適合純粹理論研究 :)。發現圖書館一本語言處理方面的書也沒有後，在京東找了一本書--《NLP漢語天然語言處理原理與實踐》，到今天看了大約150頁，發現仍是很模糊，決定找點代碼來看。

---

 

　　從最簡單的分詞開始，發現分詞的庫已經不少了，選擇了比較輕巧的jieba來研究。看了一下GitHub的基本介紹，忽然感受：我次奧，這也不過如此嘛，來來來寫一個。jieba對於詞典外的詞用HMM模型進行解決，用Viterbi算法實現。網上對於HMM的解釋不少，我我的也不太可以經過數學公司解釋，其實模型比較簡單，先了解了馬爾可夫模型後就能比較容易地理解隱馬爾可夫模型。

---

 

　　這裏記錄一下對於HMM模型中，解決求隱藏狀態鏈問題的Viterbi算法的學習。

在知乎問題：https://www.zhihu.com/question/20136144 中，我看了高票答案地回答，基本理解了思想，可是這個回答稍微有一點不全面，並無強調回溯的思想，這裏讓我對算法產生了一點誤解，後面有提到。

　　大約花了半天時間閱讀網上各類資料後，我選了一個Python代碼決定仿照寫一次 （https://blog.csdn.net/youfefi/article/details/74276546）代碼用了numpy數組，因爲我對numpy的函數不是太熟悉，決定先用list寫一遍。令我驚訝的是，因爲我腦殼太笨加上當時對維特比算法不是很清晰，我沒有看懂做者的代碼（略顯尷尬），沒辦法，決定先按照本身的思路寫出來。寫的過程當中發現算法有點出奇簡單，竟然簡單3層循環就出來了，大概40分鐘寫完了。不出意外，寫完運行發現和網上代碼結果不一致。

　　打斷點開始調試，發現本身求出的機率矩陣是正確的，但最後路徑不正確，意識到可能本身的理解有問題。再次上網查詢資料。參考 https://www.2cto.com/kf/201609/544539.html 的文章，發現本身的理解是有問題的：

個人理解是依次迭代每一個時刻，找到機率最大的狀態即爲該時刻狀態，這樣理解錯誤在於求的隱藏狀態是一個鏈，根據馬爾可夫假設，下一刻時刻的狀態是依賴前一個狀態的，個人理解就將狀態之間割裂了，沒法行成鏈。

正確的算法是每次迭代過程當中，記錄每種狀態機率最大時其前驅狀態，這樣到最後一個時刻，選擇機率最大的狀態，再進行回溯。

代碼以下：直接使用List、迭代過程當中沒有對機率進行判斷，還有優化空間。

 1 # state 存放隱藏序列，sunny 0 rainy 1
 2 # obser 存放觀測序列  0 1 2 對應 walk shop clean
 3 # start_p 是初始機率，0元素對應sunny的初始機率 1元素對應rainy的機率
 4 # transition_p 轉移機率矩陣 2*2 行爲初始狀態  列爲新狀態
 5 # emission_p 發射機率矩陣 2*3 行爲隱藏狀態  列爲可觀測狀態
 6 
 7 # 迭代過程，每次只須要記錄第t個時間點 每一個節點的最大機率便可，後續計算時直接使用前序節點的最大機率便可
 8 def compute(obser, state, start_p, transition_p, emission_p):
 9     # max_p 記錄每一個時間點每一個狀態的最大機率，i行j列，（i,j）記錄第i個時間點 j隱藏狀態的最大機率
10     max_p = [[0 for col in range(len(state))] for row in range(len(obser))]
11     # path 記錄max_p 對應機率處的路徑 i 行 j列 （i,j）記錄第i個時間點 j隱藏狀態最大機率的狀況下 其前驅狀態
12     path = [[0 for col in range(len(state))] for row in range(len(obser))]
13     # 初始狀態(1狀態)
14     for i in range(len(state)):
15         # max_p[0][i]表示初始狀態第i個隱藏狀態的最大機率
16         # 機率 = start_p[i] * emission_p [state[i]][obser[0]]
17         max_p[0][i] = start_p[i] * emission_p[state[i]][obser[0]]
18         path[0][i] = i
19     # 後續循環狀態(2-t狀態)
20     # 此時max_p 中已記錄第一個狀態的兩個隱藏狀態機率
21     for i in range(1, len(obser)):  # 循環t-1次，初始已計算
22         max_item = [0 for i in range(len(state))]
23         for j in range(len(state)):  # 循環隱藏狀態數，計算當前狀態每一個隱藏狀態的機率
24             item = [0 for i in state]
25             for k in range(len(state)):  # 再次循環隱藏狀態數，計算選定隱藏狀態的前驅狀態爲各類狀態的機率
26                 p = max_p[i - 1][k] * emission_p[state[j]][obser[i]] * transition_p[state[k]][state[j]]
27                 # k即表明前驅狀態 k或state[k]均爲前驅狀態
28                 item[state[k]] = p
29             # 設置機率記錄爲最大狀況
30             max_item[state[j]] = max(item)
31             # 記錄最大狀況路徑(下面語句的做用：當前時刻下第j個狀態機率最大時，記錄其前驅節點)
32             # item.index(max(item))尋找item的最大值索引，因item記錄各類前驅狀況的機率
33             path[i][state[j]] = item.index(max(item))
34         # 將單個狀態的結果加入總列表max_p
35         max_p[i] = max_item
36     #newpath記錄最後路徑
37     newpath = []
38     #判斷最後一個時刻哪一個狀態的機率最大
39     p=max_p[len(obser)-1].index(max(max_p[len(obser)-1]))
40     newpath.append(p)
41     #從最後一個狀態開始倒着尋找前驅節點
42     for i in range(len(obser) - 1, 0, -1):
43         newpath.append(path[i][p])
44         p = path[i][p]
45     newpath.reverse()
46     return newpath
47 
48 
49 if __name__ == '__main__':
50     #   隱狀態
51     hidden_state = ['rainy', 'sunny']
52     #   觀測序列
53     obsevition = ['walk', 'shop', 'clean']
54     state_s = [0, 1]
55     obser = [0, 1, 2]
56     #   初始狀態，測試集中，0.6機率觀測序列以sunny開始
57     start_probability = [0.6, 0.4]
58     #   轉移機率，0.7：sunny下一天sunny的機率
59     transititon_probability = [[0.7, 0.3], [0.4, 0.6]]
60     #   發射機率，0.4：sunny在0.4機率下爲shop
61     emission_probability = [[0.1, 0.4, 0.5], [0.6, 0.3, 0.1]]
62     result = compute(obser, state_s, start_probability, transititon_probability, emission_probability)
63     for k in range(len(result)):
64         print(hidden_state[int(result[k])])