一直想用隐马可夫模型做图像识别，但是python的scikit-learn组件包的hmm module已经不再支持了，需要安装hmmlearn的组件，不过hmmlearn的多项式hmm每次出来的结果都不一样，= =||，难道是我用错了？？后来又只能去参考网上C语言的组件，模仿着把向前向后算法“复制”到python里了，废了好大功夫，总算结果一样了o(╯□╰)o。。
把代码贴出来把，省的自己不小心啥时候删掉。。。
1#-*-coding:UTF-8-*-2'''3Created on 2014年9月25日4@author: Ayumi Phoenix5'''6importnumpy as np78classHMM:9def__init__(self, Ann, Bnm, pi1n):10self.A =np.array(Ann)11self.B =np.array(Bnm)12self.pi =np.array(pi1n)13self.N =self.A.shape[0]14self.M = self.B.shape[1]1516defPRinthmm(self):17print"=================================================="18print"HMM content: N =",self.N,",M =",self.M19foriinrange(self.N):20ifi==0:21print"hmm.A",self.A[i,:],"hmm.B",self.B[i,:]22else:23print"",self.A[i,:],"",self.B[i,:]24print"hmm.pi",self.pi25print"=================================================="2627#函数名称：Forward *功能：前向算法估计参数 *参数:phmm:指向HMM的指针28#T:观察值序列的长度 O:观察值序列29#alpha:运算中用到的临时数组 pprob:返回值,所要求的概率30defForward(self,T,O,alpha,pprob):31#1. Initialization 初始化32foriinrange(self.N):33alpha[0,i] = self.pi[i]*self.B[i,O[0]]3435#2. Induction 递归36fortinrange(T-1):37forjinrange(self.N):38sum = 0.039foriinrange(self.N):40sum += alpha[t,i]*self.A[i,j]41alpha[t+1,j] =sum*self.B[j,O[t+1]]42#3. Termination 终止43sum = 0.044foriinrange(self.N):45sum += alpha[T-1,i]46pprob[0] *=sum4748#带修正的前向算法49defForwardWithScale(self,T,O,alpha,scale,pprob):50scale[0] = 0.051#1. Initialization52foriinrange(self.N):53alpha[0,i] = self.pi[i]*self.B[i,O[0]]54scale[0] +=alpha[0,i]5556foriinrange(self.N):57alpha[0,i] /=scale[0]5859#2. Induction60fortinrange(T-1):61scale[t+1] = 0.062forjinrange(self.N):63sum = 0.064foriinrange(self.N):65sum += alpha[t,i]*self.A[i,j]6667alpha[t+1,j] = sum * self.B[j,O[t+1]]68scale[t+1] += alpha[t+1,j]69forjinrange(self.N):70alpha[t+1,j] /= scale[t+1]7172#3. Termination73fortinrange(T):74pprob[0] +=np.log(scale[t])7576#函数名称：Backward * 功能:后向算法估计参数 * 参数:phmm:指向HMM的指针77#T:观察值序列的长度 O:观察值序列78#beta:运算中用到的临时数组 pprob:返回值，所要求的概率79defBackWord(self,T,O,beta,pprob):80#1. Intialization81foriinrange(self.N):82beta[T-1,i] = 1.083#2. Induction84fortinrange(T-2,-1,-1):85foriinrange(self.N):86sum = 0.087forjinrange(self.N):88sum += self.A[i,j]*self.B[j,O[t+1]]*beta[t+1,j]89beta[t,i] =sum9091#3. Termination92pprob[0] = 0.093foriinrange(self.N):94pprob[0] += self.pi[i]*self.B[i,O[0]]*beta[0,i]9596#带修正的后向算法97defBackwardWithScale(self,T,O,beta,scale):98#1. Intialization99foriinrange(self.N):100beta[T-1,i] = 1.0101102#2. Induction103fortinrange(T-2,-1,-1):104foriinrange(self.N):105sum = 0.0106forjinrange(self.N):107sum += self.A[i,j]*self.B[j,O[t+1]]*beta[t+1,j]108beta[t,i] = sum / scale[t+1]109110#Viterbi算法111#输入：A，B，pi,O 输出P(O|lambda)最大时Poptimal的路径I112defviterbi(self,O):113T =len(O)114#初始化115delta =np.zeros((T,self.N),np.float)116phi =np.zeros((T,self.N),np.float)117I =np.zeros(T)118foriinrange(self.N):119delta[0,i] = self.pi[i]*self.B[i,O[0]]120phi[0,i] =0121#递推122fortinrange(1,T):123foriinrange(self.N):124delta[t,i] = self.B[i,O[t]]*np.array([delta[t-1,j]*self.A[j,i]forjinrange(self.N)]).max()125phi[t,i] = np.array([delta[t-1,j]*self.A[j,i]forjinrange(self.N)]).argmax()126#终结127prob = delta[T-1,:].max()128I[T-1] = delta[T-1,:].argmax()129#状态序列求取130fortinrange(T-2,-1,-1):131I[t] = phi[t+1,I[t+1]]132returnI,prob133134#计算gamma : 时刻t时马尔可夫链处于状态Si的概率135defComputeGamma(self, T, alpha, beta, gamma):136fortinrange(T):137denominator = 0.0138forjinrange(self.N):139gamma[t,j] = alpha[t,j]*beta[t,j]140denominator +=gamma[t,j]141foriinrange(self.N):142gamma[t,i] = gamma[t,i]/denominator143144#计算sai(i,j) 为给定训练序列O和模型lambda时：145#时刻t是马尔可夫链处于Si状态，二时刻t+1处于Sj状态的概率146defComputeXi(self,T,O,alpha,beta,gamma,xi):147fortinrange(T-1):148sum = 0.0149foriinrange(self.N):150forjinrange(self.N):151xi[t,i,j] = alpha[t,i]*beta[t+1,j]*self.A[i,j]*self.B[j,O[t+1]]152sum +=xi[t,i,j]153foriinrange(self.N):154forjinrange(self.N):155xi[t,i,j] /=sum156157#Baum-Welch算法158#输入 L个观察序列O，初始模型：HMM={A,B,pi,N,M}159defBaumWelch(self,L,T,O,alpha,beta,gamma):160print"BaumWelch"161DELTA = 0.01 ; round = 0 ; flag = 1 ; probf = [0.0]162delta = 0.0 ; deltaprev = 0.0 ; probprev = 0.0 ; ratio = 0.0 ; deltaprev = 10e-70163164xi =np.zeros((T,self.N,self.N))165pi =np.zeros((T),np.float)166denominatorA =np.zeros((self.N),np.float)167denominatorB =np.zeros((self.N),np.float)168numeratorA =np.zeros((self.N,self.N),np.float)169numeratorB =np.zeros((self.N,self.M),np.float)170scale =np.zeros((T),np.float)171172whileTrue :173probf[0] =0174#E - step175forlinrange(L):176self.ForwardWithScale(T,O[l],alpha,scale,probf)177self.BackwardWithScale(T,O[l],beta,scale)178self.ComputeGamma(T,alpha,beta,gamma)179self.ComputeXi(T,O[l],alpha,beta,gamma,xi)180foriinrange(self.N):181pi[i] +=gamma[0,i]182fortinrange(T-1):183denominatorA[i] +=gamma[t,i]184denominatorB[i] +=gamma[t,i]185denominatorB[i] += gamma[T-1,i]186187forjinrange(self.N):188fortinrange(T-1):189numeratorA[i,j] +=xi[t,i,j]190forkinrange(self.M):191fortinrange(T):192ifO[l][t] ==k:193numeratorB[i,k] +=gamma[t,i]194195#M - step196#重估状态转移矩阵 和 观察概率矩阵197foriinrange(self.N):198self.pi[i] = 0.001/self.N + 0.999*pi[i]/L199forjinrange(self.N):200self.A[i,j] = 0.001/self.N + 0.999*numeratorA[i,j]/denominatorA[i]201numeratorA[i,j] = 0.0202203forkinrange(self.M):204self.B[i,k] = 0.001/self.M + 0.999*numeratorB[i,k]/denominatorB[i]205numeratorB[i,k] = 0.0206207pi[i]=denominatorA[i]=denominatorB[i]=0.0;208209ifflag == 1:210flag =0211probprev =probf[0]212ratio = 1213continue214215delta = probf[0] -probprev216ratio = delta /deltaprev217probprev =probf[0]218deltaprev =delta219round += 1220221ifratio <=DELTA :222print"num iteration",round223
break224225226if__name__=="__main__":227print"python my HMM"228229A = [[0.8125,0.1875],[0.2,0.8]]230B = [[0.875,0.125],[0.25,0.75]]231pi = [0.5,0.5]232hmm =HMM(A,B,pi)233234O = [[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]]237L =len(O)238T = len(O[0])#T等于最长序列的长度就好了239alpha =np.zeros((T,hmm.N),np.float)240beta =np.zeros((T,hmm.N),np.float)241gamma =np.zeros((T,hmm.N),np.float)242hmm.BaumWelch(L,T,O,alpha,beta,gamma)243244hmm.printhmm()
View Code
由于为了自己理解方便，直接翻译公式。。。其实可以用numpy的函数写的简单点的O(∩_∩)O