代码经书上正文中例子测试,通过
10.1本题要求用后向计算,这里同时提供正向算法代码
-
import numpy as np def back_prob(): beta = np.ones(N) for t in range(T - 2, -1, -1): beta_old = beta.copy() for i in range(N): beta[i] = 0 for j in range(N): beta[i] += A[i, j] * B[j, O[t + 1] ] * beta_old[j] p = 0 for i in range(N): p += PI[i] * B[i, O[0]] * beta[i] return p def forward_prob(): alpha = PI*B[:,O[0]] for t in range(T-1): alpha_old = alpha.copy() for i in range(N): alpha[i] = 0 for j in range(N): alpha[i] += alpha_old[j]*A[j,i] alpha[i] *= B[ i,O[t+1] ] p = sum(alpha) return p 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] ]) PI = np.array([0.2, 0.4, 0.4]) O = [0, 1, 0, 1] # 0表示红,1表示白 T = 4 N = 3 M = 2 p = back_prob() print("后向 P(O|lambda) = ", p) p = forward_prob() print("前向 P(O|lambda) = ", p)
10.2 这里的概率,其实就是γ_t4(q3)
按照书上公式(10.24),前向后向过程记录alpha,beta
import numpy as np
def back():
BETA = []
beta = np.ones(N)
BETA.append(beta)
for t in range(T - 2, -1, -1):
beta_old = beta.copy()
for i in range(N):
beta[i] = 0
for j in range(N):
beta[i] += A[i, j] * B[j, O[t + 1] ] * beta_old[j]
BETA.append(beta)
return np.array(BETA).reshape(T, N)
def forward():
ALPHA = []
alpha = PI*B[:,O[0]]
ALPHA.append(alpha)
for t in range(T-1):
alpha_old = alpha.copy()
for i in range(N):
alpha[i] = 0
for j in range(N):
alpha[i] += alpha_old[j]*A[j,i]
alpha[i] *= B[ i,O[t+1] ]
ALPHA.append(alpha)
return np.array(ALPHA).reshape(T,N)
A = np.array([
[0.5,0.1,0.4],
[0.3,0.5,0.2],
[0.2,0.2,0.6]
])
B = np.array([
[0.5,0.5],
[0.4,0.6],
[0.7,0.3]
])
PI = np.array([0.2,0.3,0.5])
T = 8
N = 3
M = 2
O = np.array([0,1,0,0,1,0,1,1])
BETA = back()
ALPHA = forward()
p = ALPHA[3,2]*BETA[3,2] / sum(ALPHA[3,:]*BETA[3,:])
print("P(i4=q3|O,lambda) = ",p)
10.3
最优路径:
import numpy as np
def viterbi():
delta = PI * B[:, O[0]]
FAI = []
fai = np.zeros(N)
FAI.append(fai)
for t in range(1, T):
delta_old = delta.copy()
for i in range(N):
delta[i] = np.max(delta_old * A[:, i] * B[i, O[t]])
fai[i] = np.argmax(delta_old * A[:, i])
FAI.append(fai.copy()) #注意,应为copy
FAI = np.array(FAI).reshape(T, N)
path = []
I = np.argmax(delta)
path.append(I)
for t in range(T - 2, -1, -1):
I = int( FAI[t + 1, I] )
path.append(I)
return (np.array(path)+1)[::-1] # 使状态从1开始计数 状态正向顺序
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]
])
PI = np.array([0.2, 0.4, 0.4])
O = np.array([0, 1, 0, 1]) # 0表示红,1表示白
T = 4
# O = np.array([0, 1, 0]) # 0表示红,1表示白
# T = 3
N = 3
M = 2
path = viterbi()
print(path)