維特比演算法詳解 Python 版¶

這是一個範例。 讓我們生成一個14天的序列，其中每1代表高溫，0代表寒冷。 我們將使用該算法查找這未來兩周中最有可能的天氣預報。

In [4]:

# 匯入需要的模組
import numpy as np
import pandas as pd

# 觀察對應表
obs_map = {'Cold':0, 'Hot':1}
# 觀察序列
obs = np.array([1,1,0,1,0,0,1,0,1,1,0,0,0,1])

# 建立觀察對應表的逆字典
inv_obs_map = dict((v,k) for k, v in obs_map.items())
# 建立觀察序列的文字列表
obs_seq = [inv_obs_map[v] for v in list(obs)]

# 建立模擬觀察序列的 DataFrame
print("Simulated Observations:\n",
      pd.DataFrame(np.column_stack([obs, obs_seq]),columns=['Obs_code', 'Obs_seq']) )

# 初始化矩陣
pi = [0.6,0.4] # initial probabilities vector
# 狀態文字清單
states = ['Cold', 'Hot']
# 隱藏天氣狀態清單
hidden_states = ['Snow', 'Rain', 'Sunshine']
# 初始化矩陣
pi = [0, 0.2, 0.8]
# 建立狀態空間，建立一個序列，資料內容為初始矩陣，索引為隱藏天氣狀態，名稱為"states"
state_space = pd.Series(pi, index=hidden_states, name='states')
# 建立 轉移矩陣 a_df ，欄位為 隱藏天氣狀態清單， 索引也為 隱藏天氣狀態清單
a_df = pd.DataFrame(columns=hidden_states, index=hidden_states)
# 建立 索引 0 下雪轉移天氣機率清單
a_df.loc[hidden_states[0]] = [0.3, 0.3, 0.4]
# 建立 索引 1 下雨轉移天氣機率清單
a_df.loc[hidden_states[1]] = [0.1, 0.45, 0.45]
# 建立 索引 2 晴天轉移天氣機率清單
a_df.loc[hidden_states[2]] = [0.2, 0.3, 0.5]
# 列印 轉移矩陣 a_df
print("\n HMM matrix:\n", a_df)
# 轉換成 ndarray 資料型態，a 即為 轉移矩陣 。
a = a_df.values

# 觀察狀趟 等於 狀態文字清單
observable_states = states
# 建立 對應變換矩陣 b_df，欄位名稱為 觀察狀態清單，索引為 隱藏天氣狀態清單
b_df = pd.DataFrame(columns=observable_states, index=hidden_states)
# 建立 索引 0 下雪對應變換矩陣（室內溫度）機率清單
b_df.loc[hidden_states[0]] = [1,0]
# 建立 索引 1 下雨對應變換矩陣（室內溫度）機率清單
b_df.loc[hidden_states[1]] = [0.8,0.2]
# 建立 索引 2 晴天對應變換矩陣（室內溫度）機率清單
b_df.loc[hidden_states[2]] = [0.3,0.7]
# 列印觀察層的矩陣 b_df
print("\n Observable layer  matrix:\n",b_df)
# 將 對應矩陣轉換為 ndarray 資料型態，即為對應變換矩陣 b
b = b_df.values

Simulated Observations:
    Obs_code Obs_seq
0         1     Hot
1         1     Hot
2         0    Cold
3         1     Hot
4         0    Cold
5         0    Cold
6         1     Hot
7         0    Cold
8         1     Hot
9         1     Hot
10        0    Cold
11        0    Cold
12        0    Cold
13        1     Hot

 HMM matrix:
          Snow  Rain Sunshine
Snow      0.3   0.3      0.4
Rain      0.1  0.45     0.45
Sunshine  0.2   0.3      0.5

 Observable layer  matrix:
          Cold  Hot
Snow        1    0
Rain      0.8  0.2
Sunshine  0.3  0.7

Now let’s use the algorithm:

現在讓我們使用以下算法：

In [5]:

'''
We used the following implementation, based on [2]:

[1] https://cse.buffalo.edu/~jcorso/t/CSE555/files/lecture_hmm.pdf
[2] http://www.blackarbs.com/blog/introduction-hidden-markov-models-python-networkx-sklearn/2/9/2017
'''
# 維特比演算法

# 輸入參數

# pi: 初始矩陣 (initial probabilities vector)
#     pi = [0, 0.2, 0.8]

# a: 轉移矩陣 (HMM matrix)
#              Snow  Rain Sunshine
#     Snow      0.3   0.3      0.4
#     Rain      0.1  0.45     0.45
#     Sunshine  0.2   0.3      0.5

# b: 對應變換矩陣 (Observable layer matrix)
#               Cold  Hot
#     Snow         1    0
#     Rain       0.8  0.2
#     Sunshine   0.3  0.7

# obs: 觀察序列 (Simulated Observations)
#     obs = np.array([1,1,0,1,0,0,1,0,1,1,0,0,0,1])

def viterbi(pi, a, b, obs):
    
    # 計算有多少狀態 np.shape(b) = (3, 2)
    # nStates = 3
    nStates = np.shape(b)[0] 
    # 計算觀察序列長度 np.shape(obs) = (14,)
    # T = 14
    T = np.shape(obs)[0]
    
    # 初始化路徑，矩陣歸零
    # init blank path
    path = path = np.zeros(T,dtype=int)
    # 到達狀態i的任何路徑之最高概率
    # delta --> highest probability of any path that reaches state i
    # 初始化 delta 矩陣為 3 x 14 ， 並且歸零
    delta = np.zeros((nStates, T))
    # 每個狀態按時間步長排列的最大值索引
    # phi --> argmax by time step for each state
    # 初始化 phi 矩陣為 3 x 14 ， 並且歸零
    phi = np.zeros((nStates, T))
    
    # 初始化 delta 和 phi
    # init delta and phi 
    # 初始化 delta 在 t = 0 時的機率， 初始機率 乘上 對應機率的第一個觀察序列。
    # delta[:, 0] = [0, 0.2, 0.8] * [0, 0.2, 0.7] = [0., 0.04, 0.56]
    delta[:, 0] = pi * b[:, obs[0]]
    # 初始化 phi 所有列的第一欄為零
    # phi[:, 0] = [0., 0., 0.]
    phi[:, 0] = 0
    
    # 開始向前計算最佳路徑表，尋找機率值最大的存入 delta，機率值最大的索引存入 phi，
    # 建立 delta 和 phi 表。
    print('\nStart Walk Forward\n')    
    # 向前延伸演算法
    # the forward algorithm extension
    # t = {1 ~ 13}，因為是使用 range ，T = 14 必須減 1
    for t in range(1, T):
        # s = {0, 1, 2}
        for s in range(nStates):
            # 計算 delta，t 初始值為 1，delta[:,t-1] = delta[:,0] 
            # 取得 delta 索引零值為 [0., 0.04, 0.56]
            # delta 乘上 轉移矩陣的 s 欄，並且使用 np.max 找出最大機率。
            # 在這裡是要找出 delta 上一個時間步的所有天氣轉移的最高機率。
            # 對應矩陣 b ，去乘上不同天氣狀態 s ，當下的觀察序列 obs[t]。
            # 這裡可以解釋成因為馬可夫性質只與前一天有關，因此 delta[:, t-1] 代表前天的機率，
            # 去乘上 轉移機率，
            # 就變成今天的天氣狀態機率，由於隱藏式馬可夫的關係，無法直接得知天氣狀態，
            # 因此再乘上 對應矩陣 於 當下狀天氣狀態及當下的觀察序列。
            # 目的是將所有天氣狀態與所有觀察序列時步都填上天氣狀態的機率值。
            delta[s, t] = np.max(delta[:, t-1] * a[:, s]) * b[s, obs[t]] 
            # 計算 phi，t 初始值為 1，delta[:,t-1] = delta[:,0] 
            # 取得 delta 索引零值為 [0., 0.04, 0.56]
            # 因為 delta 先算，所以此時的 delta 不再是初始值，而是觀測值（因與對應矩陣 b 相乘）的機率。
            # 這時候 delta 與 轉移矩陣相乘，在 t-1 的情況下，今天天氣狀況只與前天有關。
            # 就可以找到最大機率，其索引也可稱為最短路徑。而最短路徑將可求得最佳路徑。
            phi[s, t] = np.argmax(delta[:, t-1] * a[:, s])
            print('s={s} and t={t}: phi[{s}, {t}] = {phi}'.format(s=s, t=t, phi=phi[s, t]))

# delta 表格填完，就得到了每一個狀態在每一個時步的機率            
#    [[0.         0.         0.0392     0.         0.00296352 0.00088906
#      0.         0.00008066 0.         0.         0.00000213 0.00000064
#      0.00000019 0.        ]
#     [0.04       0.0336     0.04704    0.0042336  0.00355622 0.00128024
#      0.00011522 0.00009679 0.00000871 0.00000183 0.00000256 0.00000092
#      0.00000033 0.00000003]
#     [0.56       0.196      0.0294     0.0148176  0.00222264 0.00048009
#      0.00040328 0.00006049 0.00003049 0.00001067 0.0000016  0.00000035
#      0.00000012 0.0000001 ]]

# phi 表格填完，就得到每個狀態在每一時步上的索引值
#    [[0. 2. 2. 0. 2. 0. 0. 2. 0. 2. 2. 0. 0. 0.]
#     [0. 2. 2. 1. 2. 1. 1. 2. 1. 2. 2. 1. 1. 1.]
#     [0. 2. 2. 1. 2. 1. 1. 2. 1. 2. 2. 1. 1. 1.]]
    
    # 尋找最佳路徑
    # find optimal path
    print('-'*50)
    # 開始向後追蹤
    print('Start Backtrace\n')
    # 一開始 path = [0 0 0 0 0 0 0 0 0 0 0 0 0 0] 有 13 個零
    # T = 14 ， T-1 = 13
    # delta[:, T-1] = [0.         0.00000003 0.0000001 ]
    # 因此 path[13] = 2， 因為 np.argmax 會回傳機率最大值的索引
    # 所以 path[t-1] = [0 0 0 0 0 0 0 0 0 0 0 0 0 2]
    path[T-1] = np.argmax(delta[:, T-1])
    # range(T-2), -1, -1) = [12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0]
    # phi 表格填完，就得到每個狀態在每一時步上的索引值
    #      0  1  2  3  4  5  6  7  8  9 10 11 12 13
    #    [[0. 2. 2. 0. 2. 0. 0. 2. 0. 2. 2. 0. 0. 0.]
    #     [0. 2. 2. 1. 2. 1. 1. 2. 1. 2. 2. 1. 1. 1.]
    #     [0. 2. 2. 1. 2. 1. 1. 2. 1. 2. 2. 1. 1. 1.]]
    # 假設 T-2=12，path[t+1] = path[13] = 2， [t+1] = 13
    # 經由查表後，phi[2, 13] = 1，並寫入 path[12] = 1
    # 假設 T-2=11，path[t+1] = path[12] = 1， [t+1] = 12
    # 經由查表後，phi[1, 12] = 1，並寫如 path[11] = 1
    # 以此類推，查表後的結果為 [2 2 1 2 1 1 2 1 2 2 1 1 1 2] 此為最佳路徑
    for t in range(T-2, -1, -1):
        # 當下 t 只與上一個狀態有關 t+1
        path[t] = phi[path[t+1], [t+1]]
        # 列印出時步，及最佳路徑
        print('path[{}] = {}'.format(t, path[t]))
    # 回傳參數
    # path: 最佳路徑
    # delta: 到達狀態i的任何路徑之最高概率矩陣
    # phi: 每個狀態按時間步長排列的最大值索引矩陣
    return path, delta, phi

In [6]:

# 維特比演算法
# 輸入參數：
# pi: 初始矩陣 (initial probabilities vector)
# a: 轉移矩陣 (HMM matrix)
# b: 對應變換矩陣 (Observable layer matrix)
# obs: 觀察序列 (Simulated Observations)
# 輸出參數：
# path: 最佳路徑
# delta: 到達狀態i的任何路徑之最高概率矩陣
# phi: 每個狀態按時間步長排列的最大值索引矩陣
path, delta, phi = viterbi(pi, a, b, obs)
# 設定狀態對應字典
state_map = {0:'Snow', 1:'Rain', 2:'Sunshine'}
# 將最佳路徑轉換為文字
state_path = [state_map[v] for v in path]
# 建立觀察序列與最佳路徑（天氣狀態）的 DataFrame
pd.DataFrame().assign(Observation=obs_seq).assign(Best_Path=state_path)

Start Walk Forward

s=0 and t=1: phi[0, 1] = 2.0
s=1 and t=1: phi[1, 1] = 2.0
s=2 and t=1: phi[2, 1] = 2.0
s=0 and t=2: phi[0, 2] = 2.0
s=1 and t=2: phi[1, 2] = 2.0
s=2 and t=2: phi[2, 2] = 2.0
s=0 and t=3: phi[0, 3] = 0.0
s=1 and t=3: phi[1, 3] = 1.0
s=2 and t=3: phi[2, 3] = 1.0
s=0 and t=4: phi[0, 4] = 2.0
s=1 and t=4: phi[1, 4] = 2.0
s=2 and t=4: phi[2, 4] = 2.0
s=0 and t=5: phi[0, 5] = 0.0
s=1 and t=5: phi[1, 5] = 1.0
s=2 and t=5: phi[2, 5] = 1.0
s=0 and t=6: phi[0, 6] = 0.0
s=1 and t=6: phi[1, 6] = 1.0
s=2 and t=6: phi[2, 6] = 1.0
s=0 and t=7: phi[0, 7] = 2.0
s=1 and t=7: phi[1, 7] = 2.0
s=2 and t=7: phi[2, 7] = 2.0
s=0 and t=8: phi[0, 8] = 0.0
s=1 and t=8: phi[1, 8] = 1.0
s=2 and t=8: phi[2, 8] = 1.0
s=0 and t=9: phi[0, 9] = 2.0
s=1 and t=9: phi[1, 9] = 2.0
s=2 and t=9: phi[2, 9] = 2.0
s=0 and t=10: phi[0, 10] = 2.0
s=1 and t=10: phi[1, 10] = 2.0
s=2 and t=10: phi[2, 10] = 2.0
s=0 and t=11: phi[0, 11] = 0.0
s=1 and t=11: phi[1, 11] = 1.0
s=2 and t=11: phi[2, 11] = 1.0
s=0 and t=12: phi[0, 12] = 0.0
s=1 and t=12: phi[1, 12] = 1.0
s=2 and t=12: phi[2, 12] = 1.0
s=0 and t=13: phi[0, 13] = 0.0
s=1 and t=13: phi[1, 13] = 1.0
s=2 and t=13: phi[2, 13] = 1.0
--------------------------------------------------
Start Backtrace

path[12] = 1
path[11] = 1
path[10] = 1
path[9] = 2
path[8] = 2
path[7] = 1
path[6] = 2
path[5] = 1
path[4] = 1
path[3] = 2
path[2] = 1
path[1] = 2
path[0] = 2

Out[6]:

	Observation	Best_Path
0	Hot	Sunshine
1	Hot	Sunshine
2	Cold	Rain
3	Hot	Sunshine
4	Cold	Rain
5	Cold	Rain
6	Hot	Sunshine
7	Cold	Rain
8	Hot	Sunshine
9	Hot	Sunshine
10	Cold	Rain
11	Cold	Rain
12	Cold	Rain
13	Hot	Sunshine

賴岱佑（Lai Tai-Yu）