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的警告，还没有去掉这个警告，查了一下说不影响结果，后面会去解决这个问题，下面贴出我的代码
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. 0]#数组A的行数
48. 1]#数组A的列数
49. 1]#数组O的列数
50.
51. sum_alpha_1 = np.zeros((M,N))
52. alpha = np.zeros((N,H))
53. 1,N))
54. 0,:], B[:,O[0,0]-1])
55.
56. 0] = np.array(alpha_1).reshape(1,N)#alpha_1是一维数组，在使用np.multiply的时候需要升级到二维数组。#错误是IndexError: too many indices for array
1. for h in range(1,H):
2. for i in range(N):
3. for j in range(M):
4. 1] * A[j,i]
5. 1).reshape(1,N)#同理，将数组升级为二维数组
6. 0,i] * B[i,O[0,h]-1]
7. #print("alpha矩阵: \n %r" % alpha)
8. 0).reshape(1,H)
9. 0,H-1]
10. #print("观测概率: \n %r" % P)
11. #return alpha
12. return alpha, P
13.
14. def BackwardAlgo(A,B,Pi,O):
15. 0]#数组A的行数
16. 1]#数组A的列数
17. 1]#数组O的列数
18.
19. #beta = np.zeros((N,H))
20. 1,N))
21. beta = np.zeros((N,H))
22. 1] = 1
23. 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. 0,j] = A[i,j] * B[j,O[0,h]-1] * beta[j,h]
29. 1] = sum_beta.sum(1)
30. #print("beta矩阵: \n %r" % beta)
31. for i in range(N):
32. 0,i] = Pi[0,i] * B[i,O[0,0]-1] * beta[i,0]
33. 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. 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. 0]#数组A的行数
49. 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. 0]#数组A的行数
60. 1]#数组A的列数
61. 1]#数组O的列数
62. 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. 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. 0]#数组A的行数
74. 1]#数组A的列数
75. 1]#数组B的列数
76. 1]#数组O的列数
77. 0
78. Gamma = GetGamma(A,B,Pi,O)
79. Xi = GetXi(A,B,Pi,O)
80. 0)
81. a = np.zeros((N,M))
82. b = np.zeros((M,Y))
83. 1,N))
84. 1),Gamma[:,H-1]).reshape(1,N)
85. for i in range(N):
86. for j in range(M):
87. 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. 1).reshape(1,N)
95. 0,j]
96. 0
97. #print(b)
98. for i in range(N):
99. 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. 0]#数组A的行数
110. 1]#数组A的列数
111. 1]#数组O的列数
112. Delta = np.zeros((M,H))
113. Psi = np.zeros((M,H))
114. 1))
115. 1,H))
116.
117. for i in range(N):
118. 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. 0] = Delta[i,h-1] * A[i,j] * B[j,O[0,h]-1]
124. Delta[j,h] = np.amax(Delta_1)
125. 1
126. print("Delta矩阵: \n %r" % Delta)
127. print("Psi矩阵: \n %r" % Psi)
128. 1])
129. 1])
130. 0,H-1] = psi + 1
131. for h in range(H-1,0,-1):
132. 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
由于第一次用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的警告，还没有去掉这个警告，查了一下说不影响结果，后面会去解决这个问题，下面贴出我的代码
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. 0]#数组A的行数
48. 1]#数组A的列数
49. 1]#数组O的列数
50.
51. sum_alpha_1 = np.zeros((M,N))
52. alpha = np.zeros((N,H))
53. 1,N))
54. 0,:], B[:,O[0,0]-1])
55.
56. 0] = np.array(alpha_1).reshape(1,N)#alpha_1是一维数组，在使用np.multiply的时候需要升级到二维数组。#错误是IndexError: too many indices for array
1. for h in range(1,H):
2. for i in range(N):
3. for j in range(M):
4. 1] * A[j,i]
5. 1).reshape(1,N)#同理，将数组升级为二维数组
6. 0,i] * B[i,O[0,h]-1]
7. #print("alpha矩阵: \n %r" % alpha)
8. 0).reshape(1,H)
9. 0,H-1]
10. #print("观测概率: \n %r" % P)
11. #return alpha
12. return alpha, P
13.
14. def BackwardAlgo(A,B,Pi,O):
15. 0]#数组A的行数
16. 1]#数组A的列数
17. 1]#数组O的列数
18.
19. #beta = np.zeros((N,H))
20. 1,N))
21. beta = np.zeros((N,H))
22. 1] = 1
23. 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. 0,j] = A[i,j] * B[j,O[0,h]-1] * beta[j,h]
29. 1] = sum_beta.sum(1)
30. #print("beta矩阵: \n %r" % beta)
31. for i in range(N):
32. 0,i] = Pi[0,i] * B[i,O[0,0]-1] * beta[i,0]
33. 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. 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. 0]#数组A的行数
49. 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. 0]#数组A的行数
60. 1]#数组A的列数
61. 1]#数组O的列数
62. 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. 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. 0]#数组A的行数
74. 1]#数组A的列数
75. 1]#数组B的列数
76. 1]#数组O的列数
77. 0
78. Gamma = GetGamma(A,B,Pi,O)
79. Xi = GetXi(A,B,Pi,O)
80. 0)
81. a = np.zeros((N,M))
82. b = np.zeros((M,Y))
83. 1,N))
84. 1),Gamma[:,H-1]).reshape(1,N)
85. for i in range(N):
86. for j in range(M):
87. 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. 1).reshape(1,N)
95. 0,j]
96. 0
97. #print(b)
98. for i in range(N):
99. 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. 0]#数组A的行数
110. 1]#数组A的列数
111. 1]#数组O的列数
112. Delta = np.zeros((M,H))
113. Psi = np.zeros((M,H))
114. 1))
115. 1,H))
116.
117. for i in range(N):
118. 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. 0] = Delta[i,h-1] * A[i,j] * B[j,O[0,h]-1]
124. Delta[j,h] = np.amax(Delta_1)
125. 1
126. print("Delta矩阵: \n %r" % Delta)
127. print("Psi矩阵: \n %r" % Psi)
128. 1])
129. 1])
130. 0,H-1] = psi + 1
131. for h in range(H-1,0,-1):
132. 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
