1 #-*-coding:UTF-8-*-  2 ‘‘‘  3 Created on 2014年9月25日  4 @author: Ayumi Phoenix  5 ‘‘‘  6 import numpy as np  7   8 class HMM:  9     def __init__(self, Ann, Bnm, pi1n): 10         self.A = np.array(Ann) 11         self.B = np.array(Bnm) 12         self.pi = np.array(pi1n) 13         self.N = self.A.shape[0] 14         self.M = self.B.shape[1] 15          16     def printhmm(self): 17         print "==================================================" 18         print "HMM content: N =",self.N,",M =",self.M 19         for i in range(self.N): 20             if i==0: 21                 print "hmm.A ",self.A[i,:]," hmm.B ",self.B[i,:] 22             else: 23                 print "      ",self.A[i,:],"       ",self.B[i,:] 24         print "hmm.pi",self.pi 25         print "==================================================" 26  27     # 函数名称：Forward *功能：前向算法估计参数 *参数:phmm:指向HMM的指针 28     # T:观察值序列的长度 O:观察值序列     29     # alpha:运算中用到的临时数组 pprob:返回值,所要求的概率 30     def Forward(self,T,O,alpha,pprob): 31     #     1. Initialization 初始化 32         for i in range(self.N): 33             alpha[0,i] = self.pi[i]*self.B[i,O[0]] 34      35     #     2. Induction 递归 36         for t in range(T-1): 37             for j in range(self.N): 38                 sum = 0.0 39                 for i in range(self.N): 40                     sum += alpha[t,i]*self.A[i,j] 41                 alpha[t+1,j] =sum*self.B[j,O[t+1]] 42     #     3. Termination 终止 43         sum = 0.0 44         for i in range(self.N): 45             sum += alpha[T-1,i] 46         pprob[0] *= sum 47      48     # 带修正的前向算法 49     def ForwardWithScale(self,T,O,alpha,scale,pprob): 50         scale[0] = 0.0 51     #     1. Initialization 52         for i in range(self.N): 53             alpha[0,i] = self.pi[i]*self.B[i,O[0]] 54             scale[0] += alpha[0,i] 55          56         for i in range(self.N): 57             alpha[0,i] /= scale[0] 58          59     #     2. Induction 60         for t in range(T-1): 61             scale[t+1] = 0.0 62             for j in range(self.N): 63                 sum = 0.0 64                 for i in range(self.N): 65                     sum += alpha[t,i]*self.A[i,j] 66                  67                 alpha[t+1,j] = sum * self.B[j,O[t+1]] 68                 scale[t+1] += alpha[t+1,j] 69             for j in range(self.N): 70                 alpha[t+1,j] /= scale[t+1] 71          72     #     3. Termination 73         for t in range(T): 74             pprob[0] += np.log(scale[t]) 75  76     # 函数名称：Backward * 功能:后向算法估计参数 * 参数:phmm:指向HMM的指针  77     # T:观察值序列的长度 O:观察值序列  78     # beta:运算中用到的临时数组 pprob:返回值，所要求的概率 79     def Backword(self,T,O,beta,pprob): 80     #     1. Intialization 81         for i in range(self.N): 82             beta[T-1,i] = 1.0 83     #     2. Induction 84         for t in range(T-2,-1,-1): 85             for i in range(self.N): 86                 sum = 0.0 87                 for j in range(self.N): 88                     sum += self.A[i,j]*self.B[j,O[t+1]]*beta[t+1,j] 89                 beta[t,i] = sum 90                  91     #     3. Termination 92         pprob[0] = 0.0 93         for i in range(self.N): 94             pprob[0] += self.pi[i]*self.B[i,O[0]]*beta[0,i] 95      96     # 带修正的后向算法 97     def BackwardWithScale(self,T,O,beta,scale): 98     #     1. Intialization 99         for i in range(self.N):100             beta[T-1,i] = 1.0101     102     #     2. Induction103         for t in range(T-2,-1,-1):104             for i in range(self.N):105                 sum = 0.0106                 for j in range(self.N):107                     sum += self.A[i,j]*self.B[j,O[t+1]]*beta[t+1,j]108                 beta[t,i] = sum / scale[t+1]109     110     # Viterbi算法111     # 输入：A，B，pi,O 输出P(O|lambda)最大时Poptimal的路径I112     def viterbi(self,O):113         T = len(O)114         # 初始化115         delta = np.zeros((T,self.N),np.float)  116         phi = np.zeros((T,self.N),np.float)  117         I = np.zeros(T)118         for i in range(self.N):  119             delta[0,i] = self.pi[i]*self.B[i,O[0]]  120             phi[0,i] = 0121         # 递推122         for t in range(1,T):  123             for i in range(self.N):                                  124                 delta[t,i] = self.B[i,O[t]]*np.array([delta[t-1,j]*self.A[j,i]  for j in range(self.N)]).max()125                 phi[t,i] = np.array([delta[t-1,j]*self.A[j,i]  for j in range(self.N)]).argmax()126         # 终结127         prob = delta[T-1,:].max()  128         I[T-1] = delta[T-1,:].argmax()129         # 状态序列求取   130         for t in range(T-2,-1,-1): 131             I[t] = phi[t+1,I[t+1]]132         return I,prob133     134     # 计算gamma : 时刻t时马尔可夫链处于状态Si的概率    135     def ComputeGamma(self, T, alpha, beta, gamma):136         for t in range(T):137             denominator = 0.0138             for j in range(self.N):139                 gamma[t,j] = alpha[t,j]*beta[t,j]140                 denominator += gamma[t,j]141             for i in range(self.N):142                 gamma[t,i] = gamma[t,i]/denominator143     144     # 计算sai(i,j) 为给定训练序列O和模型lambda时：145     # 时刻t是马尔可夫链处于Si状态，二时刻t+1处于Sj状态的概率146     def ComputeXi(self,T,O,alpha,beta,gamma,xi):147         for t in range(T-1):148             sum = 0.0149             for i in range(self.N):150                 for j in range(self.N):151                     xi[t,i,j] = alpha[t,i]*beta[t+1,j]*self.A[i,j]*self.B[j,O[t+1]]152                     sum += xi[t,i,j]153             for i in range(self.N):154                 for j in range(self.N):155                     xi[t,i,j] /= sum156                     157     # Baum-Welch算法158     # 输入 L个观察序列O，初始模型：HMM={A,B,pi,N,M}159     def BaumWelch(self,L,T,O,alpha,beta,gamma):160         print "BaumWelch"161         DELTA = 0.01 ; round = 0 ; flag = 1 ; probf = [0.0]162         delta = 0.0 ; deltaprev = 0.0 ; probprev = 0.0 ; ratio = 0.0 ; deltaprev = 10e-70163         164         xi = np.zeros((T,self.N,self.N))165         pi = np.zeros((T),np.float)166         denominatorA = np.zeros((self.N),np.float)167         denominatorB = np.zeros((self.N),np.float)168         numeratorA = np.zeros((self.N,self.N),np.float)169         numeratorB = np.zeros((self.N,self.M),np.float)170         scale = np.zeros((T),np.float)171         172         while True :173             probf[0] = 0174             # E - step175             for l in range(L):176                 self.ForwardWithScale(T,O[l],alpha,scale,probf)177                 self.BackwardWithScale(T,O[l],beta,scale)178                 self.ComputeGamma(T,alpha,beta,gamma)179                 self.ComputeXi(T,O[l],alpha,beta,gamma,xi)180                 for i in range(self.N):181                     pi[i] += gamma[0,i]182                     for t in range(T-1): 183                         denominatorA[i] += gamma[t,i]184                         denominatorB[i] += gamma[t,i]185                     denominatorB[i] += gamma[T-1,i]186                     187                     for j in range(self.N):188                         for t in range(T-1):189                             numeratorA[i,j] += xi[t,i,j]190                     for k in range(self.M):191                         for t in range(T):192                             if O[l][t] == k:193                                 numeratorB[i,k] += gamma[t,i]194                             195             # M - step196             # 重估状态转移矩阵 和 观察概率矩阵197             for i in range(self.N):198                 self.pi[i] = 0.001/self.N + 0.999*pi[i]/L199                 for j in range(self.N):200                     self.A[i,j] = 0.001/self.N + 0.999*numeratorA[i,j]/denominatorA[i]201                     numeratorA[i,j] = 0.0202                 203                 for k in range(self.M):204                     self.B[i,k] = 0.001/self.M + 0.999*numeratorB[i,k]/denominatorB[i]205                     numeratorB[i,k] = 0.0206                 207                 pi[i]=denominatorA[i]=denominatorB[i]=0.0;208             209             if flag == 1:210                 flag = 0211                 probprev = probf[0]212                 ratio = 1213                 continue214             215             delta = probf[0] - probprev216             ratio = delta / deltaprev217             probprev = probf[0]218             deltaprev = delta219             round += 1220             221             if ratio <= DELTA :222                 print "num iteration ",round223                 break224          225 226 if __name__ == "__main__":227     print "python my HMM"228   229     A = [[0.8125,0.1875],[0.2,0.8]]230     B = [[0.875,0.125],[0.25,0.75]]231     pi = [0.5,0.5]232     hmm = HMM(A,B,pi)233     234     O = [[1,0,0,1,1,0,0,0,0],235      [1,1,0,1,0,0,1,1,0],236      [0,0,1,1,0,0,1,1,1]]237     L = len(O)238     T = len(O[0])  # T等于最长序列的长度就好了239     alpha = np.zeros((T,hmm.N),np.float)240     beta = np.zeros((T,hmm.N),np.float)241     gamma = np.zeros((T,hmm.N),np.float)242     hmm.BaumWelch(L,T,O,alpha,beta,gamma)243     244     hmm.printhmm()