HMM_機率計算——forwar_algorithm實現

1、前向算法（遞推）

　　1）觀測機率初值：

　　　　　

　　2）遞推：

　　  　　

　　3）終止：

　　　　　　

2、理解前向算法

哈哈，先欠着，沒時間碼字！

3、代碼實現：

 1 def getPro():
 2     iniStatePro=[]   #初始狀態機率分佈列表
 3     pi=np.array(iniStatePro)
 4     
 5     stateTransPro=[]           #狀態轉移機率矩陣  N x N 
 6     A=np.array(stateTransPro)  # A[a][b]： 狀態a -> 狀態b 的機率
 7     
 8     observationPro=[]          #觀測機率矩陣     N x M 
 9     B=np.array(observationPro) # B[a][o]：狀態a->觀測值o的機率  
10     return A,B,pi;
11     
12 def getForwardPro(A,B,pi,T,O):  #O爲觀測序列
13     alpha=np.zeros(shape=(T,len(A)))     #定義前向機率
14     for i in range(len(B)):                  #初始化前向機率分佈
15         alpha[0][i]=pi[i]*B[ i ][ O[0] ] 
16     
17     for t in range(1,T):            # 遞推
18         for i in range( len(A) ):   #  t時刻每一個狀態
19             pro=0
20             for j in range( len(B) ):        #  前一狀態到當前時刻 聯合機率
21                 pro=pro+ A[j][i]*alpha[t-1][j];
22             alpha[t][i]=pro*B[i][ O[t] ]          #預測機率
23     OSeqBaseLambdaPro=sum( alpha[T-1] );          # p( O| lambda )
24     return OSeqBaseLambdaPro;

4、驗證：統計學習方法 177 編程實現

 1 def getPro():
 2     iniStatePro=[ 0.2 ,0.4, 0.4 ]   #初始狀態機率分佈列表 3個狀態
 3     pi=np.array(iniStatePro)
 4     
 5     stateTransPro=[[ 0.5,0.2,0.3 ],[ 0.3,0.5,0.2 ],[ 0.2, 0.3, 0.5 ]]           #狀態轉移機率矩陣  N x N 
 6     A=np.array(stateTransPro)  # A[a][b]： 狀態a -> 狀態b 的機率
 7     
 8     observationPro=[[0.5,0.5],[0.4,0.6],[0.7,0.3]]          #觀測機率矩陣     N x M    [紅，白]
 9     B=np.array(observationPro) # B[a][o]：狀態a->觀測值o的機率  
10     
11     return A,B,pi;
12 
13 def getForwardPro(A,B,pi,T,O):  #O爲觀測序列
14     alpha=np.zeros(shape=(T,len(A)))     #定義前向機率
15     for i in range(len(B)):                  #初始化前向機率分佈
16         alpha[0][i]=pi[i]*B[ i ][ O[0] ] 
17     
18     for t in range(1,T):            # 遞推
19         for i in range( len(A) ):   #  t時刻每一個狀態
20             pro=0
21             for j in range( len(B) ):        #  前一狀態到當前時刻 聯合機率
22                 pro=pro+ A[j][i]*alpha[t-1][j];
23             alpha[t][i]=pro*B[i][ O[t] ]          #預測機率
24     OSeqBaseLambdaPro=sum( alpha[T-1] );          # p( O| lambda )
25     return OSeqBaseLambdaPro;
26 
27 A,B,pi=getPro();
28 T=3
29 O=[0,1,0]
30 pro=getForwardPro(A,B,pi,T,O)
31 print(pro)

結果：

書上答案：0.13022