一直想用隐马可夫模型做图像识别，但是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()
由于为了自己理解方便，直接翻译公式。。。其实可以用numpy的函数写的简单点的O(∩_∩)O