使用python实现HMM

一直想用隐马可夫模型做图像识别，但是python的scikit-learn组件包的hmm module已经不再支持了，需要安装hmmlearn的组件，不过hmmlearn的多项式hmm每次出来的结果都不一样，= =||，难道是我用错了？？后来又只能去参考网上C语言的组件，模仿着把向前向后算法“复制”到python里了，废了好大功夫，总算结果一样了o(╯□╰)o。。
把代码贴出来把，省的自己不小心啥时候删掉。。。
  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.0
101     
102     #     2. Induction
103         for t in range(T-2,-1,-1):
104             for i in range(self.N):
105                 sum = 0.0
106                 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的路径I
112     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] = 0
121         # 递推
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,prob
133     
134     # 计算gamma : 时刻t时马尔可夫链处于状态Si的概率    
135     def ComputeGamma(self, T, alpha, beta, gamma):
136         for t in range(T):
137             denominator = 0.0
138             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]/denominator
143     
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.0
149             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] /= sum
156                     
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-70
163         
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] = 0
174             # E - step
175             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 - step
196             # 重估状态转移矩阵 和 观察概率矩阵
197             for i in range(self.N):
198                 self.pi[i] = 0.001/self.N + 0.999*pi[i]/L
199                 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.0
202                 
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.0
206                 
207                 pi[i]=denominatorA[i]=denominatorB[i]=0.0;
208             
209             if flag == 1:
210                 flag = 0
211                 probprev = probf[0]
212                 ratio = 1
213                 continue
214             
215             delta = probf[0] - probprev
216             ratio = delta / deltaprev
217             probprev = probf[0]
218             deltaprev = delta
219             round += 1
220             
221             if ratio <= DELTA :
222                 print "num iteration ",round
223                 break
224          
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()
View Code
由于为了自己理解方便，直接翻译公式。。。其实可以用numpy的函数写的简单点的O(∩_∩)O