一直想用隐马可夫模型做图像识别,但是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