HMM的前向算法和维特比算法已将在上一个博客讲过了,这节给出python的实现,例子是以李航老师《统计学习方法》上的算法和例题进行实现。
import numpy as np
class HMM:
def forward(self,Q,V,A,B,O,PI):
"""
前向算法
Q:状态序列
V:观测集合
A:状态转移概率矩阵
B:生成概率矩阵
O:观测序列
PI:状态概率向量
"""
N = len(Q)
M = len(O)
# alphas是前向概率矩阵,每列是某个时刻每个状态的前向概率
alphas = np.zeros((N,M))
T = M # 时刻
for t in range(T):
# 当前时刻观测值对应的索引
index_O = V.index(O[t])
for i in range(N):
if t == 0:
alphas[i][t] = PI[0][i] * B[i][index_O]
#print('alpha1(%d)=p%db%db(o1)=%f' % (i, i, i, alphas[i][t]))
else:
# 这里使用点积和使用矩阵相乘后再sum没区别
alphas[i][t] = np.dot([alpha[t-1] for alpha in alphas],[a[i] for a in A]) * B[i][index_O]
#print('alpha%d(%d)=[sigma alpha%d(i)ai%d]b%d(o%d)=%f' % (t, i, t - 1, i, i, t, alphas[i][t]))
# 将前向概率的最后一列的概率相加就是观测序列生成的概率
P = np.sum([alpha[M-1] for alpha in alphas])
print('alphas :')
print(alphas)
print('P=%5f' % P)
def backward(self, Q, V, A, B, O, PI):
"""后向算法"""
N = len(Q)
M = len(O)
# betas 是后向概率矩阵,第M列的概率都为1
betas = np.zeros((N,M))
for i in range(N):
#print('beta%d(%d)=1' % (M, i))
betas[i][M-1] = 1
# 倒着,从后往前递推
for t in range(M-2,-1,-1):
index_O = V.index(O[t])
for i in range(N):
betas[i][t] = np.dot(np.multiply(A[i],[b[index_O] for b in B]),[beta[t+1] for beta in betas])
realT = t + 1
realI = i + 1
#print('beta%d(%d)=[sigma a%djbj(o%d)]beta%d(j)=(' % (realT, realI, realI, realT + 1, realT + 1),end='')
#for j in range(N):
#print("%.2f*%.2f*%.2f+" % (A[i][j], B[j][index_O], betas[j][t + 1]), end='')
#print("0)=%.3f" % betas[i][t])
index_O = V.index(O[0])
P = np.dot(np.multiply(PI,[b[index_O] for b in B]),[beta[0] for beta in betas])
print("P(O|lambda)=", end="")
for i in range(N):
print("%.1f*%.1f*%.5f+" % (PI[0][i], B[i][index_O], betas[i][0]), end="")
print("0=%f" % P)
print(betas)
def viterbi(self,Q,V,A,B,O,PI):
N = len(Q)
M = len(O)-1
# deltas是记录每个时刻,每个状态的最优路径的概率
deltas = np.zeros((N,M))
# psis是记录每个时刻每个状态的前一个时刻的概率最大的路径的节点(psia的记录是比deltas晚一个时刻的)
psis = np.zeros((N,M))
# I是记录最优路径的
I = np.zeros((1,M))
for t in range(M):
realT = t+1
index_O = V.index(O[t])
for i in range(N):
realI = i+1
if t == 0 :
deltas[i][t] = PI[0][i] * B[i][index_O]
psis[i][t] = 0
#print('delta1(%d)=pi%d * b%d(o1)=%.2f * %.2f=%.2f'%(realI, realI, realI, PI[0][i], B[i][index_O], deltas[i][t]))
#print('psis1(%d)=0' % (realI))
else:
deltas[i][t] = np.max(np.multiply([delta[t-1] for delta in deltas],[a[i] for a in A])) * B[i][index_O]
psis[i][t] = np.argmax(np.multiply([delta[t-1] for delta in deltas],[a[i] for a in A]))
#print('delta%d(%d)=max[delta%d(j)aj%d]b%d(o%d)=%.2f*%.2f=%.5f'%(realT, realI, realT-1, realI, realI,realT,
# np.max(np.multiply([delta[t-1] for delta in deltas], [a[i] for a in A])),
# B[i][index_O], deltas[i][t]))
#print('psis%d(%d)=argmax[delta%d(j)aj%d]=%d' % (realT, realI, realT-1, realI, psis[i][t]))
# 直接取最后一列最大概率对应的状态为最优路径最后的节点
I[0][M-1] = np.argmax([delta[M-1] for delta in deltas])
for t in range(M-2,-1,-1):
# psis要晚一个时间步,起始将最后那个状态对应在psis那行直接取出就是最后的结果
# 但是那样体现不出回溯,下面这种每次取上一个最优路径点对应的上一个最优路径点
I[0][t] = psis[int(I[0][t+1])][t+1]
print(I)
if __name__ == '__main__':
Q = [1,2,3]
V = ['红', '白']
A = [[0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]]
B = [[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]]
# O = ['红', '白', '红', '红', '白', '红', '白', '白']
O = ['红', '白', '红','白'] #习题10.1的例子
PI = [[0.2, 0.4, 0.4]]
Hmm = HMM()
Hmm.forward(Q, V, A, B, O, PI)
#Hmm.backward(Q, V, A, B, O, PI)
#Hmm.viterbi(Q, V, A, B, O, PI)
自己的代码也是参考了别人的代码,只是将一些输出注释掉和添加了一些自己理解的注释,参考那位博主写的还是很好的。自己贴到这里只是为了以后自己查看方便,再次感谢这位博主。
这里还有另一个博客写的很好,可以在实现以上的代码后,自己进行练手。