import  numpy as np

#you must install the numpy  

A=np.array([[0.5,0.2,0.3],[0.3,0.5,0.2],[0.2,0.3,0.5]])

B=np.array([[0.5,0.5],[0.4,0.6],[0.7,0.3]])

O=np.array([0 ,1, 0])#T=3

pi=np.array([0.2,0.4,0.4])  

N=3#N kind state

M=2#M kind of observation

T=3  

#initialize:  

delta=np.zeros((T,N),np.float)

w=np.zeros((T,3),np.float)

I=np.zeros(T)  

for i in range(N):

    delta[0,i]=pi[i]*B[i,O[0]]

    w[0,i]=0  

#Recursion:

for t in range(1,T):

    for i in range(N):

        delta[t,i]=B[i,O[t]]*np.array([delta[t-1,j]*A[j,i]  for j in range(N)]).max()

        w[t,i]=np.array([delta[t-1,j]*A[j,i]  for j in range(N)]).argmax()  

#termination

Poptimal=delta[T-1,:].max()

I[T-1]=delta[T-1,:].argmax()

for t in range(T-2,-1,-1):

        I[t]=w[t+1,I[t+1]]  

print 'the optimal I is ' I