使用python实现HMM

qaqa12345667

　　一直想用隐马可夫模型做图像识别，但是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*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*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*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*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 += gamma[0,i]
182                     for t in range(T-1):
183                         denominatorA += gamma[t,i]
184                         denominatorB += gamma[t,i]
185                     denominatorB += 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 = 0.001/self.N + 0.999*pi/L
199                 for j in range(self.N):
200                     self.A[i,j] = 0.001/self.N + 0.999*numeratorA[i,j]/denominatorA
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
205                     numeratorB[i,k] = 0.0
206                 
207                 pi=denominatorA=denominatorB=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