HMM 回顧

《統計學習方法》 p.174

隱馬爾科夫模型（HMM）有三個基本的問題

• （1）概率計算問題。給定模型 λ=(A,B,Pi) 和觀測序列 O(o1,o2,...,oT),計算在模型 λ 下觀測序列 O 的概率 P(O|λ)
• （2）學習問題。已知觀測序列 O(o1,o2,...,oT),估計模型 λ=(A,B,Pi) 的引數，使得在該模型下觀測序列概率 P(O|λ) 最大。即用極大似然估計的方法估計引數。
• （3）預測問題，也稱為解碼問題。已知模型 λ=(A,B,Pi) 和觀測序列 O(o1,o2,...,oT),求對給定觀測序列條件概率 P(I|O) 最大的狀態序列 I=(i
  1,i2,...,iT)，即給定觀測序列，求最有可能的對應的狀態序列。

維特比譯碼

《統計學習方法》 p.186

問題型別：

（3）預測問題。已知模型 lambda = (A, B, Pi) 和觀測序列，求最優的狀態序列.

問題描述：

一共有三個盒子 <=> 狀態集合 Q = {1,2,3}
盒子裡邊裝著 紅色 和 白色 兩種顏色的球 <=> 觀測集合 V = {紅，白}
從盒子中有放回的取 3 次球，顏色分別為 O=(紅，白，紅)，問取球的盒子順序最可能是什麼？

注意：程式碼中所有下標從 0 開始，而課本中問題描述是從 1 開始。

import numpy as np

# 模型引數
A = np.asarray([[0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]]) # 轉移矩陣
B = np.asarray([[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]])
Pi = np.asarray([0.2, 0.4, 0.4]).transpose()

O = np.asarray([0,1,0]) 
T = O.shape[0]
N = A.shape[0]   # 狀態數

p_nodes = Pi * B[:, O[0]]     # 記錄每個節點的路徑概率
path_nodes = list()           # 記錄每個節點的路徑
 

# 計初始化路徑
for node in xrange(N):
    path_nodes.append([node])
# T 個時刻
for step in xrange(1, T):
    for this_node in xrange(N):   # 計算每個節點的新概率
        p_news = list()
        for last_node in xrange(N):
            p_trans = A[last_node, this_node]  # 轉移概率
            p_out = B[this_node, O[step]]       # 輸出概率
            p_new = p_nodes[last_node] * p_trans * p_out
            p_news.append(p_new)
        p_nodes[this_node] = np.max(p_news)    # 更新節點路徑概率
        last_index = np.argmax(p_news)         # 更新節點路徑
        temp = path_nodes[last_index][:]
        temp.append(this_node)
        path_nodes[this_node] = temp

print p_nodes     # 最有一步每個節點的概率
print path_nodes  
max_index = np.argmax(p_nodes)
max_path = path_nodes[max_index]
print max_path   # 最優路徑

[ 0.00756  0.01008  0.0147 ]
[[2, 2, 0], [2, 2, 1], [2, 2, 2]]
[2, 2, 2]

實際專案中，一般都是通過對概率取對數，將乘法轉為加法進行求解的。但是，要注意轉移概率為 0 的特殊情況，下面不考慮這種情況。

##############################
# 維特比譯碼
# 《統計學習方法》 p.186
#############################

# 模型引數
A = np.asarray([[0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]]) # 轉移矩陣
B = np.asarray([[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]])
Pi = np.asarray([0.2, 0.4, 0.4]).transpose()
A = np.log(A)
B = np.log(B)
Pi = np.log(Pi)

O = np.asarray([0,1,0])  # 觀測序列


def Viterbi_decode(A, B, Pi, O):
    """維特比解碼,所有概率為對數概率。
    輸入：
        A: N×N 的轉移矩陣
        B: N×M 的輸出矩陣
        Pi: list, 初始狀態概率分佈
        O: list, 觀測序列
    返回：
        max_path: list, 最優路徑。
    """
    T = O.shape[0]
    N = A.shape[0]   # 狀態數
    p_nodes = Pi + B[:, O[0]]     # 記錄每個節點的路徑概率
    path_nodes = list()           # 記錄每個節點的路徑
    # 計初始化路徑
    for node in xrange(N):
        path_nodes.append([node])
    # T 個時刻
    for step in xrange(1, T):
        for this_node in xrange(N):   # 計算每個節點的新概率
            p_news = list()
            for last_node in xrange(N):
                p_trans = A[last_node, this_node]  # 轉移概率
                p_out = B[this_node, O[step]]       # 輸出概率
                p_new = p_nodes[last_node] + p_trans + p_out
                p_news.append(p_new)
            p_nodes[this_node] = np.max(p_news)    # 更新節點路徑概率
            last_index = np.argmax(p_news)         # 更新節點路徑
            temp = path_nodes[last_index][:]
            temp.append(this_node)
            path_nodes[this_node] = temp
    max_index = np.argmax(p_nodes)
    max_path = path_nodes[max_index]        
    return max_path

max_path = Viterbi_decode(A, B, Pi, O)
print max_path

[2, 2, 2]

• 初始路徑概率 pnodes=Pi+B[:,O[0]], 改成使用 Bi-LSTM 模型輸出的分類概率
• 輸出概率 pout=B[thisnode,O[step]], 改成使用 Bi-LSTM 模型輸出的分類概率