技術文章 » 李航老師《統計學習方法》第二版第十章課後題答案

李航老師《統計學習方法》第二版第十章課後題答案

2020-10-14 11:00:17

10.1使用後向演演算法計算 P ( O ∣ λ ) P(O|\lambda) P(O∣λ)

課後題10.1
解：
話不多說，上程式碼，自己寫的，通用版，看到網上有多個不同版本的答案，還都是手算的，很佩服這些大佬的耐心，那麼多的小數，一點點的算。第一個函數Bw_Recurrent(A,B,start_p,list_observation)為後向演演算法，名字裡面的Bw表示BackWards

# -*- coding: utf-8 -*-
import numpy as np

def Bw_Recurrent(A,B,start_p,list_observation):
    '''
    Parameters
    ----------
    A : np.float32 matrix
        是隱馬模型的狀態轉移矩陣.
    B : np.float32 matrix
        觀測概率矩陣.
    start_p : np.float32 matrix
        初始狀態概率分佈.
    list_observation:list
        用於計算特定序列的概率
        
    Returns
    -------
    float32 .
    是P(O|lamda)

    '''
    #後向的遞迴計算需要初始化一個全一的大小為N*1的矩陣
    N=np.shape(A)[0]
    p=np.ones((N,1),dtype=np.float32)
    
    T=len(list_observation)#觀測序列的長度
    for i in range(T-1):
        #將觀測矩陣裡面的第list_observation[T-i-1]列取出來
        v=np.transpose(np.array([B[:,list_observation[T-i-1]]],dtype=np.float32))
        p=np.matmul(A,v*p)#這行程式碼既有矩陣乘法，也有矩陣點乘
        
    o1=np.transpose(np.array([B[:,list_observation[0]]],dtype=np.float32))
    return np.sum(start_p*o1*p)#注意裡面的矩陣乘法是點乘操作，也就是將對應位置的元素相乘

def Fw_Recurrent(A,B,start_p,list_observation):
    """
    Parameters
    ----------
    A : np.float32 matrix
        是隱馬模型的狀態轉移矩陣.
    B : np.float32 matrix
        觀測概率矩陣.
    start_p : np.float32 matrix
        初始狀態概率分佈.
    list_observation:list
        用於計算特定序列的概率
    Returns
    -------
    float32 .
    是P(O|lamda)
    """
    o1=np.transpose(np.array([B[:,list_observation[0]]],dtype=np.float32))
    p=start_p*o1
    T=len(list_observation)#觀測序列的長度
    for i in range(T-1):
        p=np.matmul(np.transpose(A),p)*np.transpose(np.array([B[:,list_observation[i+1]]],dtype=np.float32))
        #要注意上面這行程式碼裡面的矩陣乘法和矩陣點乘
        print(p,'\n')
    return np.sum(p)

if __name__=='__main__':
    
    A=np.array([[0.5,0.2,0.3],
                [0.3,0.5,0.2],
                [0.2,0.3,0.5]],dtype=np.float32)
    B=np.array([[0.5,0.5],
                [0.4,0.6],
                [0.7,0.3]],dtype=np.float32)
    
    start_p=np.array([[0.2],[0.4],[0.4]],dtype=np.float32)
    list_observation=[0,1,0,1]

    s_fw=Fw_Recurrent(A,B,start_p,list_observation)  
    s_bw=Bw_Recurrent(A,B,start_p,list_observation)
    print('前向演演算法的結果是： ',s_fw)
    print('後向演演算法的結果是： ',s_bw)

執行結果是：

前向演演算法的結果是：  0.060090814
後向演演算法的結果是：  0.060090803

由上圖可以看到，這裡輸出的結果有十分微小的差距，這是因為計算機的截斷誤差導致的

10.2 使用前後向演演算法計算 P ( i 4 = q 3 ∣ O , λ ) P(i_{4}=q_{3}|O,\lambda ) P(i4​=q3​∣O,λ)

課後題10.2
解：
如果想要計算這個題的答案的話，還需要修改下上題的程式。李航老師書上有計算公式
P ( i t = q i ∣ O , λ ) = P ( i t = q i , O ∣ λ ) P ( O ∣ λ ) = α t ( i ) ∗ β t ( i ) P ( O ∣ λ ) P(i_{t}=q_{i}|O,\lambda )=\frac{P(i_{t}=q_{i},O|\lambda )}{P(O|\lambda )}=\frac{\alpha _{t}(i)*\beta _{t}(i)}{P(O|\lambda )} P(it​=qi​∣O,λ)=P(O∣λ)P(it​=qi​,O∣λ)​=P(O∣λ)αt​(i)∗βt​(i)​
公式分析：上面公式的第一個等號成立使用了條件概率和聯合概率之間的計算關係----貝葉斯公式
關鍵在於如何獲取計算過程中的迭代向量，然後個別得到上面公式要求的，再按照公式計算就行了。

import numpy as np

def Bw_Recurrent(A,B,start_p,list_observation):
    '''
    Parameters
    ----------
    A : np.float32 matrix
        是隱馬模型的狀態轉移矩陣.
    B : np.float32 matrix
        觀測概率矩陣.
    start_p : np.float32 matrix
        初始狀態概率分佈.
    list_observation:list
        用於計算特定序列的概率

    Returns
    -------
    float32 .
    是P(O|lamda)

    '''
    P=[]
    #後向的遞迴計算需要初始化一個全一的大小為N*1的矩陣
    N=np.shape(A)[0]
    p=np.ones((N,1),dtype=np.float32)
    P.append(p)
    T=len(list_observation)#觀測序列的長度
    for i in range(T-1):
        #將觀測矩陣裡面的第list_observation[T-i-1]列取出來
        v=np.transpose(np.array([B[:,list_observation[T-i-1]]],dtype=np.float32))
        p=np.matmul(A,v*p)#這行程式碼既有矩陣乘法，也有矩陣點乘
        P.append(p)
        
    o1=np.transpose(np.array([B[:,list_observation[0]]],dtype=np.float32))
    return np.sum(start_p*o1*p),P#注意裡面的矩陣乘法是點乘操作，也就是將對應位置的元素相乘

def Fw_Recurrent(A,B,start_p,list_observation):
    """
    Parameters
    ----------
    A : np.float32 matrix
        是隱馬模型的狀態轉移矩陣.
    B : np.float32 matrix
        觀測概率矩陣.
    start_p : np.float32 matrix
        初始狀態概率分佈.
    list_observation:list
        用於計算特定序列的概率

    Returns
    -------
    float32 .
    是P(O|lamda)

    """
    P=[]
    o1=np.transpose(np.array([B[:,list_observation[0]]],dtype=np.float32))
    p=start_p*o1
    P.append(p)
    T=len(list_observation)#觀測序列的長度
    for i in range(T-1):
        p=np.matmul(np.transpose(A),p)*np.transpose(np.array([B[:,list_observation[i+1]]],dtype=np.float32))
        #要注意上面這行程式碼裡面的矩陣乘法和矩陣點乘
        P.append(p)
    return np.sum(p),P

if __name__=='__main__':
    
    A=np.array([[0.5,0.1,0.4],
                [0.3,0.5,0.2],
                [0.2,0.2,0.6]],dtype=np.float32)
    B=np.array([[0.5,0.5],
                [0.4,0.6],
                [0.7,0.3]],dtype=np.float32)
    
    start_p=np.array([[0.2],[0.3],[0.5]],dtype=np.float32)
    list_observation=[0,1,0,0,1,0,1,1]

    s_fw,P_fw=Fw_Recurrent(A,B,start_p,list_observation)  
    s_bw,P_bw=Bw_Recurrent(A,B,start_p,list_observation)
    print('前向演演算法的結果是： ',s_fw)
    print('後向演演算法的結果是： ',s_bw)
    print('result is :',P_fw[3][2]*P_bw[4][2]/s_bw)

計算結果是：

前向演演算法的結果是：  0.00347671
後向演演算法的結果是：  0.0034767103
result is : [0.5369518]

10.3 在習題10.1中，試用維特比(Viterbi)演演算法求最優的路徑

解：
看到維特比演演算法在很多地方都有應用，但是在不同的地方應用的方式是不同的，如果給了柵欄網圖，並且再給每一條邊都附上權值，我相信很多都會計算，但是，不同的應用場景中，如何計算這個權值是個關鍵，根據書上列好的演演算法無論是寫程式還是計算起來都很簡單，但是，我覺得關鍵是如何做到活學活用。在不同的場景裡面應用自如。好藍呀！！！
一個演演算法學了差不多一個下午，但是感覺還是差點感性認識。就這樣吧，上程式碼

# -*- coding: utf-8 -*-
"""
Created on Mon Oct 12 15:52:10 2020

@author: DELL
"""


import numpy as np

def viterbi_decode(A,B,start_p,list_observation):
    '''
    Parameters
    ----------
    A : 2D matrix float
        狀態轉移矩陣.
    B : 2D matrix float
        觀測概率矩陣.
    start_p : float matrix
        初始狀態概率分佈.
    observation_list : list
        觀測序列.

    Returns
    -------
    最優路徑.
    '''
    best_path=[]
    T=len(list_observation)
    #初始化
    N=A.shape[0]#狀態個數
    O1=np.transpose(np.array([B[:,list_observation[0]]],dtype=np.float32))
    delta=start_p*O1
    print(delta,'\n')
    all_psi=[]
    psi=np.zeros([N,1],dtype=np.float32)
    
    for i in range(T-1):
        psi=np.argmax(delta*A,axis=0)
        delta=np.transpose(np.array([np.max(delta*A,axis=0)]))*np.transpose(np.array([B[:,list_observation[i+1]]],dtype=np.float32))
        
        all_psi.append(psi)
        
    final_score=np.max(delta)
    best_path.append(np.argmax(delta))
    for i in range(T-1):
        best_path.insert(0,all_psi[T-2-i][best_path[0]])
        
    return final_score,best_path
    
        
        
if __name__=='__main__':
    
    A=np.array([[0.5,0.2,0.3],
                [0.3,0.5,0.2],
                [0.2,0.3,0.5]],dtype=np.float32)
    B=np.array([[0.5,0.5],
                [0.4,0.6],
                [0.7,0.3]],dtype=np.float32)
    
    start_p=np.array([[0.2],[0.4],[0.4]],dtype=np.float32)
    list_observation=[0,1,0,1]
    final_score,best_path=viterbi_decode(A,B,start_p,list_observation)
    print('最優路徑的概率是：',final_score)
    print('最優路徑為：',best_path)

執行結果是：

最優路徑的概率是： 0.0030240004
最優路徑為： [2, 1, 1, 1]

分析：因為我的狀態程式碼是從0開始的，因此對應起來的正確路徑就是[3,2,2,2]
書上的例子執行結果是：

最優路徑的概率是： 0.014700001
最優路徑為： [2, 2, 2]

對應起來就是[3,3,3]

10.4、試用前向概率和後向概率推導 P ( O ∣ λ ) = ∑ i = 1 N ∑ j = 1 N α t ( i ) a i j b j ( o t + 1 ) β t + 1 ( j ) P(O|\lambda )=\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha _{t}(i)a_{ij}b_{j}(o_{t+1})\beta _{t+1}(j) P(O∣λ)=i=1∑N​j=1∑N​αt​(i)aij​bj​(ot+1​)βt+1​(j)

t = 1 , 2 , . . . , T − 1 t=1,2,...,T-1 t=1,2,...,T−1
Proof:
根據後向演演算法的遞推公式計算可得： β t ( i ) = ∑ j = 1 N a i j b j ( o t + 1 ) β t + 1 ( j ) i = 1 , 2 , . . . , N (1) \beta _{t}(i)=\sum_{j=1}^{N}a_{ij}b_{j}(o_{t+1})\beta _{t+1}(j)\quad i=1,2,...,N\tag{1} βt​(i)=j=1∑N​aij​bj​(ot+1​)βt+1​(j)i=1,2,...,N(1)
又根據李航老師書上的公式可得：
P ( O ∣ λ ) = ∑ j = 1 N α t ( j ) β t ( j ) (2) P(O|\lambda )=\sum_{j=1}^{N}\alpha _{t}(j)\beta _{t}(j) \tag{2} P(O∣λ)=j=1∑N​αt​(j)βt​(j)(2)
就算不使用李航老師書上的公式，也可以很簡單的推出公式（2）
下面推導一下：
P ( O ∣ λ ) = ∑ i = 1 N ∑ j = 1 N α t ( i ) a i j b j ( o t + 1 ) β t + 1 ( j ) = ∑ i = 1 N α t ( i ) ∑ j = 1 N a i j b j ( o t + 1 ) β t + 1 ( j ) P(O|\lambda )=\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha _{t}(i)a_{ij}b_{j}(o_{t+1})\beta _{t+1}(j)=\sum_{i=1}^{N}\alpha _{t}(i)\sum_{j=1}^{N}a_{ij}b_{j}(o_{t+1})\beta _{t+1}(j) P(O∣λ)=i=1∑N​j=1∑N​αt​(i)aij​bj​(ot+1​)βt+1​(j)=i=1∑N​αt​(i)j=1∑N​aij​bj​(ot+1​)βt+1​(j)
= ∑ i = 1 N α t ( i ) β t ( i ) =\sum_{i=1}^{N}\alpha _{t}(i)\beta _{t}(i) =i=1∑N​αt​(i)βt​(i)

下面給出詳細的證明。（2）式的推導需要用到聯合分佈與邊際分佈之間的關係。還有最重要的就是一個假設：觀測獨立性假設，也就是任意時刻的觀測值只依賴於該時刻馬爾可夫鏈的狀態值，與其他狀態和觀測值無關。下面的公式（3）成立。
P ( o t ∣ i T , o T , i T − 1 , o T − 1 , . . . , i 1 , o 1 ) = p ( o t ∣ i t ) (3) P(o_{t}|i_{T},o_{T},i_{T-1},o_{T-1},...,i_{1},o_{1})=p(o_{t}|i_{t})\tag{3} P(ot​∣iT​,oT​,iT−1​,oT−1​,...,i1​,o1​)=p(ot​∣it​)(3)

根據聯合分佈於邊際分佈之間的關係。那麼就有：
P ( O ∣ λ ) = ∑ i = 1 N P ( O , i t = q i ∣ λ ) (4) P(O|\lambda )=\sum_{i=1}^{N}P(O,i_{t}=q_{i}|\lambda )\tag{4} P(O∣λ)=i=1∑N​P(O,it​=qi​∣λ)(4)
上式的意義也就是對所有的在 i t i_{t} it​時刻的狀態求和得到邊際分佈 P ( O ∣ λ ) P(O|\lambda ) P(O∣λ)
再根據前向演演算法後向演演算法的定義，所以有：
P ( O ∣ λ ) = ∑ i = 1 N P ( O , i t = q i ∣ λ ) = ∑ i = 1 N P ( o 1 , o 2 , . . . , o t , . . . , o T , i t = q t ∣ λ ) P(O|\lambda)=\sum_{i=1}^{N}P(O,i_{t}=q_{i}|\lambda )=\sum_{i=1}^{N}P(o_{1},o_{2},...,o_{t},...,o_{T},i_{t}=q_{t}|\lambda ) P(O∣λ)=i=1∑N​P(O,it​=qi​∣λ)=i=1∑N​P(o1​,o2​,...,ot​,...,oT​,it​=qt​∣λ)
= ∑ i = 1 N P ( o 1 , o 2 , . . . , o t , i t = q t ∣ λ ) ∗ P ( o t + 1 , . . . , o T ∣ o 1 , o 2 , . . . , o t , i t = q t , λ ) =\sum_{i=1}^{N}P(o_{1},o_{2},...,o_{t},i_{t}=q_{t}|\lambda )*P(o_{t+1},...,o_{T}|o_{1},o_{2},...,o_{t},i_{t}=q_{t},\lambda) =i=1∑N​P(o1​,o2​,...,ot​,it​=qt​∣λ)∗P(ot+1​,...,oT​∣o1​,o2​,...,ot​,it​=qt​,λ)
= ∑ i = 1 N P ( o 1 , o 2 , . . . , o t , i t = q t ∣ λ ) ∗ P ( o t + 1 , . . . , o T ∣ i t = q t , λ ) =\sum_{i=1}^{N}P(o_{1},o_{2},...,o_{t},i_{t}=q_{t}|\lambda )*P(o_{t+1},...,o_{T}|i_{t}=q_{t},\lambda) =i=1∑N​P(o1​,o2​,...,ot​,it​=qt​∣λ)∗P(ot+1​,...,oT​∣it​=qt​,λ)
= ∑ j = 1 N α t ( j ) β t ( j ) t = 1 , 2 , . . . , T − 1 =\sum_{j=1}^{N}\alpha _{t}(j)\beta _{t}(j) \quad t=1,2,...,T-1 =j=1∑N​αt​(j)βt​(j)t=1,2,...,T−1
得證。

獨立性的假設是一個比較強的條件，其實現實往往很複雜，這裡只是簡化了條件，不然估計很多理論難以實現到實際中去。

10.5、比較維特比演演算法中變數 δ \delta δ的計算和前向演演算法中變數 α \alpha α的計算的主要區別

解：
維特比演演算法計算的是時刻t是給定觀測序列 o 1 , o 2 , . . . , o t o_{1},o_{2},...,o_{t} o1​,o2​,...,ot​計算到目前時刻t的所有狀態取值的的概率最大值。
而前向演演算法的計算卻是給定初始狀態分佈，計算觀測序列的概率值。