from numpy import *
class HMM:
def __init__(self):
self.A=array([(0.5,0.2,0.3),(0.3,0.5,0.2),(0.2,0.3,0.5)])
self.B=array([(0.5,0.5),(0.4,0.6),(0.7,0.3)])
self.pi=array([(0.2),(0.4),(0.4)])
self.o=[0,1,0]
self.t=len(self.o)
self.m=len(self.A)
self.n=len(self.B[0])
def qianxiang(self):
self.x=array(zeros((self.t,self.m)))
for i in range(self.m):
self.x[0][i]=self.pi[i]*self.B[i][self.o[0]]
for j in range(1,self.t):
for i in range(self.m):
temp=0
for k in range(self.m):
temp=temp+self.x[j-1][k]*self.A[k][i]
self.x[j][i]=temp*self.B[i][self.o[j]]
result=0
for i in range(self.m):
result=result+self.x[self.t-1][i]
print u"前向概率矩阵及当前观测序列概率如下:"
print self.x
print result
def houxiang(self):
self.y=array(zeros((self.t,self.m)))
for i in range(self.m):
self.y[self.t-1][i]=1
j=self.t-2
while(j>=0):
for i in range(self.m):
for k in range(self.m):
self.y[j][i]+=self.A[i][k]*self.B[k][self.o[j+1]]*self.y[j+1][k]
j=j-1
result=0
for i in range(self.m):
result=result+self.pi[i]*self.B[i][self.o[0]]*self.y[0][i]
print u'后向概率矩阵及当前观测序列概率如下:'
print self.y
print result
def get_stateprobability(self,t,p):
if(t>self.t or p>self.m):
print u'输入数据超过范围'
return
print u'在时刻'+str(t)+u'处于状态'+str(p)+u'的概率是:'
temp=self.x[t-1][p-1]*self.y[t-1][p-1]
total=0
for i in range(self.m):
total=total+self.x[t-1][i]*self.y[t-1][i]
print temp/total
def viterbi(self):
self.z=array(zeros((self.t,self.m)))
self.d=array(zeros((self.t,self.m)))
for i in range(self.m):
self.z[0][i]=self.pi[i]*self.B[i][self.o[0]]
self.d[0][i]=0
for j in range(1,self.t):
for i in range(self.m):
maxnum=self.z[j-1][0]*self.A[0][i]
node=1
for k in range(1,self.m):
temp=self.z[j-1][k]*self.A[k][i]
if(maxnum<temp):
maxnum=temp
node=k+1
self.z[j][i]=maxnum*self.B[i][self.o[j]]
self.d[j][i]=node
max_probability=self.z[self.t-1][0]
last_node=[1]
temp=0
for i in range(1,self.m):
if(max_probability<self.z[self.t-1][i]):
max_probability=self.z[self.t-1][i]
last_node[0]=i+1
temp=i
i=self.t-1
while(i>=1):
last_node.append(self.d[i][temp])
i=i-1
temp=['o']
temp[0]=int(last_node[len(last_node)-1])
j=len(last_node)-2
while j>=0:
temp.append(int(last_node[j]))
j=j-1
print u'路径节点分别为'
print temp
print u'该路径概率为'+str(max_probability)
if __name__ == '__main__':
test=HMM()
test.qianxiang()
test.houxiang()
test.get_stateprobability(3,3)
test.viterbi()