Python实现HMM（隐马尔可夫模型）

HMM的代码，这次用Python写了一遍，依据仍然是李航博士的《统计学习方法》

由于第一次用python，所以代码可能会有许多缺陷，但是所有代码都用书中的例题进行了测试，结果正确。

这里想说一下python，在编写HMM过程中参看了之前写的MATLAB程序，发现他们有太多相似的地方，用到了numpy库，在python代码过程中最让我头疼的是数组角标，和MATLAB矩阵角标从1开始不同，numpy库数组角标都是从0开始，而且数组的维数也需要谨慎，一不小心就会出现too many indices for array的错误。程序中最后是维特比算法，在运行过程中出现了__main__:190: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future的警告，还没有去掉这个警告，查了一下说不影响结果，后面会去解决这个问题，下面贴出我的代码

[python] view plain copy

1. # -*- coding: utf-8 -*-  
2. """ 
3. Created on Thu Feb 16 19:28:39 2017 
4. 2017-4-2 
5.     ForwardBackwardAlg函数功能：实现前向算法 
6.     理论依据：李航《统计学习方法》 
7. 2017-4-5 
8.     修改了ForwardBackwardAlg函数名称为ForwardAlgo以及输出的alpha数组形式 
9.     完成了BackwardAlgo函数功能：后向算法 
10.     以及函数FBAlgoAppli：计算在观测序列和模型参数确定的情况下， 
11.     某一个隐含状态对应相应的观测状态的概率 
12. 2017-4-6 
13.     完成BaumWelchAlgo函数一次迭代 
14. 2017-4-7 
15.     实现维特比算法 
16. @author: sgp 
17. """  
18.   
19. import numpy as np  
20.   
21. #输入格式如下：  
22. #A = np.array([[.5,.2,.3],[.3,.5,.2],[.2,.3,.5]])  
23. #B = np.array([[.5,.5],[.4,.6],[.7,.3]])  
24. #Pi = np.array([[.2,.4,.4]])  
25. #O = np.array([[1,2,1]])  
26.   
27. #应用ndarray在数组之间进行相互运算时，一定要确保数组维数相同！  
28. #比如：  
29. #In[93]:m = np.array([1,2,3,4])  
30. #In[94]:m  
31. #Out[94]: array([1, 2, 3, 4])  
32. #In[95]:m.shape  
33. #Out[95]: (4,)  
34. #这里表示的是一维数组  
35. #In[96]:m = np.array([[1,2,3,4]])  
36. #In[97]:m  
37. #Out[97]: array([[1, 2, 3, 4]])  
38. #In[98]:m.shape  
39. #Out[98]: (1, 4)  
40. #而这里表示的就是二维数组  
41. #注意In[93]和In[96]的区别,多一对中括号！！  
42.   
43. #N = A.shape[0]为数组A的行数， H = O.shape[1]为数组O的列数  
44. #在下列各函数中，alpha数组和beta数组均为N*H二维数组，也就是横向坐标是时间，纵向是状态  
45.   
46. def ForwardAlgo(A,B,Pi,O):  
47.     N = A.shape[0]#数组A的行数  
48.     M = A.shape[1]#数组A的列数  
49.     H = O.shape[1]#数组O的列数  
50.   
51.     sum_alpha_1 = np.zeros((M,N))  
52.     alpha = np.zeros((N,H))  
53.     r = np.zeros((1,N))  
54.     alpha_1 = np.multiply(Pi[0,:], B[:,O[0,0]-1])  
55.       
56.     alpha[:,0] = np.array(alpha_1).reshape(1,N)#alpha_1是一维数组，在使用np.multiply的时候需要升级到二维数组。#错误是IndexError: too many indices for array  

[python] view plain copy

1.     for h in range(1,H):  
2.         for i in range(N):  
3.             for j in range(M):  
4.                 sum_alpha_1[i,j] = alpha[j,h-1] * A[j,i]  
5.             r = sum_alpha_1.sum(1).reshape(1,N)#同理，将数组升级为二维数组  
6.             alpha[i,h] = r[0,i] * B[i,O[0,h]-1]  
7.     #print("alpha矩阵: \n %r" % alpha)      
8.     p = alpha.sum(0).reshape(1,H)  
9.     P = p[0,H-1]  
10.     #print("观测概率: \n %r" % P)  
11.     #return alpha  
12.     return alpha, P  
13.       
14. def BackwardAlgo(A,B,Pi,O):  
15.     N = A.shape[0]#数组A的行数  
16.     M = A.shape[1]#数组A的列数  
17.     H = O.shape[1]#数组O的列数  
18.       
19.     #beta = np.zeros((N,H))  
20.     sum_beta = np.zeros((1,N))  
21.     beta = np.zeros((N,H))  
22.     beta[:,H-1] = 1  
23.     p_beta = np.zeros((1,N))  
24.       
25.     for h in range(H-1,0,-1):  
26.         for i in range(N):  
27.             for j in range(M):  
28.                 sum_beta[0,j] = A[i,j] * B[j,O[0,h]-1] * beta[j,h]  
29.             beta[i,h-1] = sum_beta.sum(1)  
30.     #print("beta矩阵: \n %r" % beta)  
31.     for i in range(N):  
32.         p_beta[0,i] = Pi[0,i] * B[i,O[0,0]-1] * beta[i,0]  
33.     p = p_beta.sum(1).reshape(1,1)  
34.     #print("观测概率: \n %r" % p[0,0])  
35.     return beta, p[0,0]  
36.     
37. def FBAlgoAppli(A,B,Pi,O,I):  
38.     #计算在观测序列和模型参数确定的情况下，某一个隐含状态对应相应的观测状态的概率  
39.     #例题参考李航《统计学习方法》P189习题10.2  
40.     #输入格式：  
41.     #I为二维数组，存放所求概率P(it = qi,O|lambda)中it和qi的角标t和i，即P=[t,i]  
42.     alpha,p1 = ForwardAlgo(A,B,Pi,O)  
43.     beta,p2 = BackwardAlgo(A,B,Pi,O)  
44.     p = alpha[I[0,1]-1,I[0,0]-1] * beta[I[0,1]-1,I[0,0]-1] / p1  
45.     return p  
46.       
47. def GetGamma(A,B,Pi,O):  
48.     N = A.shape[0]#数组A的行数  
49.     H = O.shape[1]#数组O的列数  
50.     Gamma = np.zeros((N,H))  
51.     alpha,p1 = ForwardAlgo(A,B,Pi,O)  
52.     beta,p2 = BackwardAlgo(A,B,Pi,O)  
53.     for h in range(H):  
54.         for i in range(N):  
55.             Gamma[i,h] = alpha[i,h] * beta[i,h] / p1  
56.     return Gamma  
57.       
58. def GetXi(A,B,Pi,O):  
59.     N = A.shape[0]#数组A的行数  
60.     M = A.shape[1]#数组A的列数  
61.     H = O.shape[1]#数组O的列数  
62.     Xi = np.zeros((H-1,N,M))  
63.     alpha,p1 = ForwardAlgo(A,B,Pi,O)  
64.     beta,p2 = BackwardAlgo(A,B,Pi,O)  
65.     for h in range(H-1):  
66.         for i in range(N):  
67.             for j in range(M):  
68.                 Xi[h,i,j] = alpha[i,h] * A[i,j] * B[j,O[0,h+1]-1] * beta[j,h+1] / p1  
69.     #print("Xi矩阵: \n %r" % Xi)  
70.     return Xi  
71.       
72. def BaumWelchAlgo(A,B,Pi,O):  
73.     N = A.shape[0]#数组A的行数  
74.     M = A.shape[1]#数组A的列数  
75.     Y = B.shape[1]#数组B的列数  
76.     H = O.shape[1]#数组O的列数  
77.     c = 0  
78.     Gamma = GetGamma(A,B,Pi,O)  
79.     Xi = GetXi(A,B,Pi,O)  
80.     Xi_1 = Xi.sum(0)  
81.     a = np.zeros((N,M))  
82.     b = np.zeros((M,Y))  
83.     pi = np.zeros((1,N))  
84.     a_1 = np.subtract(Gamma.sum(1),Gamma[:,H-1]).reshape(1,N)  
85.     for i in range(N):  
86.         for j in range(M):  
87.             a[i,j] = Xi_1[i,j] / a_1[0,i]  
88.     #print(a)  
89.     for y in range(Y):  
90.         for j in range(M):  
91.             for h in range(H):  
92.                 if O[0,h]-1 == y:  
93.                     c = c + Gamma[j,h]  
94.             gamma = Gamma.sum(1).reshape(1,N)  
95.             b[j,y] = c / gamma[0,j]  
96.             c = 0  
97.     #print(b)  
98.     for i in range(N):  
99.         pi[0,i] = Gamma[i,0]  
100.     #print(pi)  
101.     return a,b,pi  
102.       
103. def BaumWelchAlgo_n(A,B,Pi,O,n):#计算迭代次数为n的BaumWelch算法  
104.     for i in range(n):  
105.         A,B,Pi = BaumWelchAlgo(A,B,Pi,O)  
106.     return A,B,Pi  
107.       
108. def viterbi(A,B,Pi,O):  
109.     N = A.shape[0]#数组A的行数  
110.     M = A.shape[1]#数组A的列数  
111.     H = O.shape[1]#数组O的列数  
112.     Delta = np.zeros((M,H))  
113.     Psi = np.zeros((M,H))  
114.     Delta_1 = np.zeros((N,1))  
115.     I = np.zeros((1,H))  
116.       
117.     for i in range(N):  
118.         Delta[i,0] = Pi[0,i] * B[i,O[0,0]-1]  
119.           
120.     for h in range(1,H):  
121.         for j in range(M):  
122.             for i in range(N):  
123.                 Delta_1[i,0] = Delta[i,h-1] * A[i,j] * B[j,O[0,h]-1]  
124.             Delta[j,h] = np.amax(Delta_1)  
125.             Psi[j,h] = np.argmax(Delta_1)+1  
126.     print("Delta矩阵: \n %r" % Delta)  
127.     print("Psi矩阵: \n %r" % Psi)  
128.     P_best = np.amax(Delta[:,H-1])  
129.     psi = np.argmax(Delta[:,H-1])  
130.     I[0,H-1] = psi + 1  
131.     for h in range(H-1,0,-1):  
132.         I[0,h-1] = Psi[I[0,h]-1,h]  
133.     print("最优路径概率: \n %r" % P_best)  
134.     print("最优路径: \n %r" % I)  

其实代码就是翻译的公式，李航博士的书中有比较详细的推理过程，或者去找一些专业的论文文献进一步了解，这里仅仅是实现了最简单的应用，其应用实例如下：

输入数据格式：

In[117]:A
Out[117]: 
array([[ 0.5,  0.2,  0.3],
       [ 0.3,  0.5,  0.2],
       [ 0.2,  0.3,  0.5]])

In[118]:B
Out[118]: 
array([[ 0.5,  0.5],
       [ 0.4,  0.6],
       [ 0.7,  0.3]])

In[119]:Pi
Out[119]: array([[ 0.2,  0.4,  0.4]])

In[120]:O
Out[120]: array([[1, 2, 1]])

输出结果为：

In[101]:alpha,p = ForwardAlgo(A,B,Pi,O)

In[102]:alpha
Out[102]: 
array([[ 0.1     ,  0.077   ,  0.04187 ],
       [ 0.16    ,  0.1104  ,  0.035512],
       [ 0.28    ,  0.0606  ,  0.052836]])

In[103]:p
Out[103]: 0.130218

In[104]:beta,p1 = BackwardAlgo(A,B,Pi,O)

In[105]:beta
Out[105]: 
array([[ 0.2451,  0.54  ,  1.    ],
       [ 0.2622,  0.49  ,  1.    ],
       [ 0.2277,  0.57  ,  1.    ]])

In[106]:p1
Out[106]: 0.130218

In[107]:gamma = GetGamma(A,B,Pi,O)

In[108]:gamma
Out[108]: 
array([[ 0.18822283,  0.31931069,  0.32153773],
       [ 0.32216744,  0.41542644,  0.27271191],
       [ 0.48960973,  0.26526287,  0.40575036]])

In[109]:xi = GetXi(A,B,Pi,O)

In[110]:xi
Out[110]: 
array([[[ 0.1036723 ,  0.04515505,  0.03939548],
        [ 0.09952541,  0.18062019,  0.04202184],
        [ 0.11611298,  0.1896512 ,  0.18384555]],


       [[ 0.14782903,  0.04730529,  0.12417638],
        [ 0.12717136,  0.16956181,  0.11869327],
        [ 0.04653735,  0.05584481,  0.16288071]]])

In[111]:a,b,pi = BaumWelchAlgo_n(A,B,Pi,O,5)

In[112]:a
Out[112]: 
array([[ 0.43972438,  0.15395857,  0.40631705],
       [ 0.309058  ,  0.45055446,  0.24038754],
       [ 0.3757005 ,  0.50361975,  0.12067975]])

In[113]:b
Out[113]: 
array([[ 0.50277235,  0.49722765],
       [ 0.49524289,  0.50475711],
       [ 0.91925551,  0.08074449]])

In[114]:pi
Out[114]: array([[ 0.08435301,  0.18040718,  0.73523981]])

In[115]:viterbi(A,B,Pi,O)
Delta矩阵: 
 array([[ 0.1    ,  0.028  ,  0.00756],
       [ 0.16   ,  0.0504 ,  0.01008],
       [ 0.28   ,  0.042  ,  0.0147 ]])
Psi矩阵: 
 array([[ 0.,  3.,  2.],
       [ 0.,  3.,  2.],
       [ 0.,  3.,  3.]])
最优路径概率: 
 0.014699999999999998
最优路径: 
 array([[ 3.,  3.,  3.]])
__main__:192: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future