统计学习方法第十章作业：HMM模型—概率计算问题、Baum-Welch学习算法、维特比预测算法 代码实现

HMM模型

import numpy as np

class HMM:
    def __init__(self,A=None,B=None,Pi=None,O = None):
        if A:
            self.A = np.array(A)
        else:
            self.A = None
        if Pi:
            self.Pi = np.array(Pi)
            self.i_num = len(Pi)
        else:
            self.Pi = None
            self.i_num = None
        if B:
            self.B = np.array(B)
            self.o_num = np.array(B).shape[-1]
        else:
            self.B = None
            self.o_num = None
        self.O = O
        if O:
            self.T = len(O)
        else:
            self.T = None
        self.ati = None
        self.bti = None

    def update_ati(self):
        self.ati = np.zeros((self.T,self.i_num))
        for i in range(self.i_num):
            self.ati[0][i] = self.Pi[i]*self.B[i][self.O[0]]
        for t in range(self.T-1):
            for i in range(self.i_num):
                self.ati[t+1][i] = np.sum(self.ati[t]*self.A[:,i])*self.B[i][self.O[t+1]]

    def update_bti(self):
        self.bti = np.zeros((self.T,self.i_num))
        for i in range(self.i_num):
            self.bti[self.T-1][i] = 1
        for t in range(self.T-2,-1,-1):
            for i in range(self.i_num):
                self.bti[t][i] = np.sum(self.A[i,:]*self.B[:,self.O[t+1]]*self.bti[t+1])

    def comput_ri(self, t, i):
        return self.ati[t,i]*self.bti[t,i]/np.sum(self.ati[t]*self.bti[t])

    def comput_ei(self, t, i, j):
        return self.ati[t,i]*self.A[i,j]*self.B[j][self.O[t+1]]*self.bti[t+1][j] /\
    np.sum(self.ati[t]*self.A*self.B[:,self.O[t+1]]*self.bti[t+1])
    
#--------------概率计算问题--------------------
    def predict_whole_p_o(self):
        self.update_ati()
        return np.sum(self.ati[self.T-1])

    def predict_q_give_i_t(self,t,i):
        self.update_ati()
        self.update_bti()
        return self.comput_ri(t-1,i-1)

    def predict_p_t_t_1_q(self,t,i,j):
        self.update_ati()
        self.update_bti()
        return self.comput_ei(t-1,i-1,j-1)

#-------Baum-Welch学习算法--------------------
    def Baum_welch(self,O,i_num,max_iter=20):
        self.T = len(O)
        self.O = np.array(O)
        self.i_num = i_num
        self.A = np.zeros((i_num,i_num))+ 1/i_num
        self.B = np.zeros((i_num,len(set(O)))) + 1/len(set(O))
        self.Pi = np.zeros(i_num) + 1/i_num
        for n in range(max_iter):
            self.update_ati()
            self.update_bti()
            for i in range(i_num):
                self.Pi[i] = self.comput_ri(0,i)
                for j in range(i_num):
                    sum_et = 0
                    sum_rt = 0
                    for t in range(self.T-1):
                        sum_et += self.comput_ei(t,i,j)
                        sum_rt += self.comput_ri(t,i)
                    self.A[i][j] = sum_et/sum_rt

            for j in range(i_num):
                for k in range(len(set(O))):
                    sum_rt1 = 0
                    sum_rt2 = 0
                    for t in range(self.T):
                        rtj = self.comput_ri(t,j)
                        if self.O[t] == k:
                            sum_rt1 += rtj
                        sum_rt2 += rtj
                    self.B[j][k] = sum_rt1/sum_rt2
                    
#-------维特比预测算法--------------------
    def Viterbi(self):
        p_marix = np.zeros((self.T,self.i_num))
        rout = [[0] for i in range(self.i_num)]
        for i in range(self.i_num):
            p_marix[0][i] = self.Pi[i]*self.B[i][self.O[0]]
        for t in range(1,self.T):
            for i in range(self.i_num):
                max_p = 0
                max_rout = 0
                index_rout = 0
                for j in range(self.i_num):
                    max_p = max(max_p,p_marix[t-1][j]*self.A[j][i]*self.B[i][self.O[t]])
                    if p_marix[t-1][j]*self.A[j][i] > max_rout:
                        max_rout = p_marix[t-1][j]*self.A[j][i]
                        index_rout = j
                p_marix[t][i] = max_p
                rout[i].append(index_rout)
        P_max = max(p_marix[self.T-1])
        max_index = np.argmax(p_marix[self.T-1])
        rout[max_index].append(max_index)
        return P_max , [x+1 for x in rout[max_index][1:]]

def main():
    A = [[0.5,0.2,0.3],[0.3,0.5,0.2],[0.2,0.3,0.5]]
    B = [[0.5,0.5],[0.4,0.6],[0.7,0.3]]
    Pi = [0.2,0.4,0.4]
    O = [0,1,0,1,1,0,1,0,0]
    HMM_test = HMM(A=A,B=B,Pi=Pi,O=O)
    print(HMM_test.predict_whole_p_o())
    print(HMM_test.predict_q_give_i_t(4,3))
    print(HMM_test.predict_p_t_t_1_q(3, 1, 2))
    print(HMM_test.ati)
    print(HMM_test.bti)
    P_max, rout = HMM_test.Viterbi()
    print(P_max, rout)

    HMM_test_predict = HMM()
    HMM_test_predict.Baum_welch(O=O,i_num=3,max_iter=50)
    P_max, rout = HMM_test_predict.Viterbi()
    print(P_max, rout)

if __name__ == '__main__':
    main()

###########result###################
/usr/bin/python3 /Users/zhengyanzhao/PycharmProjects/tongjixuexi/shixian2/HMM.py

0.0019524748085880934
0.21041489951972903
0.09239969242123368

[[0.1        0.16       0.28      ]
 [0.077      0.1104     0.0606    ]
 [0.04187    0.035512   0.052836  ]
 [0.0210779  0.02518848 0.01382442]
 [0.01043019 0.01257429 0.00548198]
 [0.00504189 0.00400711 0.00586943]
 [0.00244848 0.00286366 0.00157461]
 [0.00119913 0.00095756 0.00146621]
 [0.00059004 0.00046339 0.00089905]]
 
[[0.00364823 0.00395888 0.00340797]
 [0.0079628  0.00754294 0.00835972]
 [0.01479809 0.01691621 0.01385703]
 [0.03194863 0.03446953 0.02971769]
 [0.07018853 0.06518296 0.07310601]
 [0.132638   0.140463   0.122819  ]
 [0.2939     0.2588     0.3123    ]
 [0.54       0.49       0.57      ]
 [1.         1.         1.        ]]
 
2.667167999999999e-06 [3, 3, 3, 3, 2, 3, 2, 3, 3]

1.0490981222700444e-07 [3, 1, 3, 1, 1, 3, 1, 3, 3]

Process finished with exit code 0