本文是李航老師的《統計學習方法》[1]一書的代碼複現。
作者：黃海廣[2]
備注：代碼都可以在github[3]中下載下傳。
我将陸續将代碼釋出在公衆号“機器學習初學者”，敬請關注。
代碼目錄
代碼參考：wzyonggege[4],WenDesi[5],火燙火燙的[6]
- 第 1 章 統計學習方法概論
- 第 2 章 感覺機
- 第 3 章 k 近鄰法
- 第 4 章 樸素貝葉斯
- 第 5 章 決策樹
- 第 6 章 邏輯斯谛回歸
- 第 7 章 支援向量機
- 第 8 章 提升方法
- 第 9 章 EM 算法及其推廣
- 第 10 章 隐馬爾可夫模型
- 第 11 章 條件随機場
- 第 12 章 監督學習方法總結
第 10 章 隐馬爾可夫模型
1．隐馬爾可夫模型是關于時序的機率模型，描述由一個隐藏的馬爾可夫鍊随機生成不可觀測的狀态的序列，再由各個狀态随機生成一個觀測而産生觀測的序列的過程。
隐馬爾可夫模型由初始狀态機率向
、狀态轉移機率矩陣
和觀測機率矩陣
決定。是以，隐馬爾可夫模型可以寫成
。
隐馬爾可夫模型是一個生成模型，表示狀态序列和觀測序列的聯合分布，但是狀态序列是隐藏的，不可觀測的。
隐馬爾可夫模型可以用于标注，這時狀态對應着标記。标注問題是給定觀測序列預測其對應的标記序列。
2．機率計算問題。給定模型
和觀測序列
＝
，
，計算在模型
下觀測序列
出現的機率
。前向-後向算法是通過遞推地計算前向-後向機率可以高效地進行隐馬爾可夫模型的機率計算。
3．學習問題。已知觀測序列
＝
，
，估計模型
參數，使得在該模型下觀測序列機率
最大。即用極大似然估計的方法估計參數。Baum-Welch 算法，也就是 EM 算法可以高效地對隐馬爾可夫模型進行訓練。它是一種非監督學習算法。
4．預測問題。已知模型
和觀測序列
＝
，
，求對給定觀測序列條件機率
最大的狀态序列
＝
，
。維特比算法應用動态規劃高效地求解最優路徑，即機率最大的狀态序列。
import numpy as np class HiddenMarkov:
def forward(self, Q, V, A, B, O, PI): # 使用前向算法
N = len(Q) #可能存在的狀态數量
M = len(O) # 觀測序列的大小
alphas = np.zeros((N, M)) # alpha值
T = M # 有幾個時刻，有幾個觀測序列，就有幾個時刻
for t in range(T): # 周遊每一時刻，算出alpha值
indexOfO = V.index(O[t]) # 找出序列對應的索引
for i in range(N):
if t == 0: # 計算初值
alphas[i][t] = PI[t][i] * B[i][indexOfO] # P176（10.15）
print(
'alpha1(%d)=p%db%db(o1)=%f' % (i, i, i, alphas[i][t]))
else:
alphas[i][t] = np.dot(
[alpha[t - 1] for alpha in alphas],
[a[i] for a in A]) * B[i][indexOfO] # 對應P176（10.16）
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]))
# print(alphas)
P = np.sum([alpha[M - 1] for alpha in alphas]) # P176(10.17)
# alpha11 = pi[0][0] * B[0][0] #代表a1(1)
# alpha12 = pi[0][1] * B[1][0] #代表a1(2)
# alpha13 = pi[0][2] * B[2][0] #代表a1(3)
def backward(self, Q, V, A, B, O, PI): # 後向算法
N = len(Q) # 可能存在的狀态數量
M = len(O) # 觀測序列的大小
betas = np.ones((N, M)) # beta
for i in range(N):
print('beta%d(%d)=1' % (M, i))
for t in range(M - 2, -1, -1):
indexOfO = V.index(O[t + 1]) # 找出序列對應的索引
for i in range(N):
betas[i][t] = np.dot(
np.multiply(A[i], [b[indexOfO] 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][indexOfO],
betas[j][t + 1]),
end='')
print("0)=%.3f" % betas[i][t])
# print(betas)
indexOfO = V.index(O[0])
P = np.dot(
np.multiply(PI, [b[indexOfO] 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][indexOfO], betas[i][0]),
end="")
print("0=%f" % P)
def viterbi(self, Q, V, A, B, O, PI):
N = len(Q) #可能存在的狀态數量
M = len(O) # 觀測序列的大小
deltas = np.zeros((N, M))
psis = np.zeros((N, M))
I = np.zeros((1, M))
for t in range(M):
realT = t + 1
indexOfO = V.index(O[t]) # 找出序列對應的索引
for i in range(N):
realI = i + 1
if t == 0:
deltas[i][t] = PI[0][i] * B[i][indexOfO]
psis[i][t] = 0
print('delta1(%d)=pi%d * b%d(o1)=%.2f * %.2f=%.2f' %
(realI, realI, realI, PI[0][i], B[i][indexOfO],
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][indexOfO]
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][indexOfO],
deltas[i][t]))
psis[i][t] = np.argmax(
np.multiply(
[delta[t - 1] for delta in deltas],
[a[i]
for a in A])) + 1 #由于其傳回的是索引，是以應+1才能和正常的下标值相符合。
print('psis%d(%d)=argmax[delta%d(j)aj%d]=%d' %
(realT, realI, realT - 1, realI, psis[i][t]))
print(deltas)
print(psis)
I[0][M - 1] = np.argmax([delta[M - 1] for delta in deltas
]) + 1 #由于其傳回的是索引，是以應+1才能和正常的下标值相符合。
print('i%d=argmax[deltaT(i)]=%d' % (M, I[0][M - 1]))
for t in range(M - 2, -1, -1):
I[0][t] = psis[int(I[0][t + 1]) - 1][t + 1]
print('i%d=psis%d(i%d)=%d' % (t + 1, t + 2, t + 2, I[0][t]))
print("狀态序列I：", I) 習題 10.1
#習題10.1
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 = HiddenMarkov()
# HMM.forward(Q, V, A, B, O, PI)
# HMM.backward(Q, V, A, B, O, PI)
HMM.viterbi(Q, V, A, B, O, PI) delta1(1)=pi1 * b1(o1)=0.20 * 0.50=0.10
psis1(1)=0
delta1(2)=pi2 * b2(o1)=0.40 * 0.40=0.16
psis1(2)=0
delta1(3)=pi3 * b3(o1)=0.40 * 0.70=0.28
psis1(3)=0
delta2(1)=max[delta1(j)aj1]b1(o2)=0.06*0.50=0.02800
psis2(1)=argmax[delta1(j)aj1]=3
delta2(2)=max[delta1(j)aj2]b2(o2)=0.08*0.60=0.05040
psis2(2)=argmax[delta1(j)aj2]=3
delta2(3)=max[delta1(j)aj3]b3(o2)=0.14*0.30=0.04200
psis2(3)=argmax[delta1(j)aj3]=3
delta3(1)=max[delta2(j)aj1]b1(o3)=0.02*0.50=0.00756
psis3(1)=argmax[delta2(j)aj1]=2
delta3(2)=max[delta2(j)aj2]b2(o3)=0.03*0.40=0.01008
psis3(2)=argmax[delta2(j)aj2]=2
delta3(3)=max[delta2(j)aj3]b3(o3)=0.02*0.70=0.01470
psis3(3)=argmax[delta2(j)aj3]=3
delta4(1)=max[delta3(j)aj1]b1(o4)=0.00*0.50=0.00189
psis4(1)=argmax[delta3(j)aj1]=1
delta4(2)=max[delta3(j)aj2]b2(o4)=0.01*0.60=0.00302
psis4(2)=argmax[delta3(j)aj2]=2
delta4(3)=max[delta3(j)aj3]b3(o4)=0.01*0.30=0.00220
psis4(3)=argmax[delta3(j)aj3]=3
[[0.1 0.028 0.00756 0.00189 ]
[0.16 0.0504 0.01008 0.003024]
[0.28 0.042 0.0147 0.002205]]
[[0. 3. 2. 1.]
[0. 3. 2. 2.]
[0. 3. 3. 3.]]
i4=argmax[deltaT(i)]=2
i3=psis4(i4)=2
i2=psis3(i3)=2
i1=psis2(i2)=3
狀态序列I：[[3. 2. 2. 2.]] 習題 10.2
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 = ['紅', '白', '紅', '紅', '白', '紅', '白', '白']
PI = [[0.2, 0.3, 0.5]] HMM.forward(Q, V, A, B, O, PI)
HMM.backward(Q, V, A, B, O, PI) alpha1(0)=p0b0b(o1)=0.100000
alpha1(1)=p1b1b(o1)=0.120000
alpha1(2)=p2b2b(o1)=0.350000
alpha1(0)=[sigma alpha0(i)ai0]b0(o1)=0.078000
alpha1(1)=[sigma alpha0(i)ai1]b1(o1)=0.111000
alpha1(2)=[sigma alpha0(i)ai2]b2(o1)=0.068700
alpha2(0)=[sigma alpha1(i)ai0]b0(o2)=0.043020
alpha2(1)=[sigma alpha1(i)ai1]b1(o2)=0.036684
alpha2(2)=[sigma alpha1(i)ai2]b2(o2)=0.055965
alpha3(0)=[sigma alpha2(i)ai0]b0(o3)=0.021854
alpha3(1)=[sigma alpha2(i)ai1]b1(o3)=0.017494
alpha3(2)=[sigma alpha2(i)ai2]b2(o3)=0.033758
alpha4(0)=[sigma alpha3(i)ai0]b0(o4)=0.011463
alpha4(1)=[sigma alpha3(i)ai1]b1(o4)=0.013947
alpha4(2)=[sigma alpha3(i)ai2]b2(o4)=0.008080
alpha5(0)=[sigma alpha4(i)ai0]b0(o5)=0.005766
alpha5(1)=[sigma alpha4(i)ai1]b1(o5)=0.004676
alpha5(2)=[sigma alpha4(i)ai2]b2(o5)=0.007188
alpha6(0)=[sigma alpha5(i)ai0]b0(o6)=0.002862
alpha6(1)=[sigma alpha5(i)ai1]b1(o6)=0.003389
alpha6(2)=[sigma alpha5(i)ai2]b2(o6)=0.001878
alpha7(0)=[sigma alpha6(i)ai0]b0(o7)=0.001411
alpha7(1)=[sigma alpha6(i)ai1]b1(o7)=0.001698
alpha7(2)=[sigma alpha6(i)ai2]b2(o7)=0.000743
beta8(0)=1
beta8(1)=1
beta8(2)=1
beta7(1)=[sigma a1jbj(o8)]beta8(j)=(0.50*0.50*1.00+0.20*0.60*1.00+0.30*0.30*1.00+0)=0.460
beta7(2)=[sigma a2jbj(o8)]beta8(j)=(0.30*0.50*1.00+0.50*0.60*1.00+0.20*0.30*1.00+0)=0.510
beta7(3)=[sigma a3jbj(o8)]beta8(j)=(0.20*0.50*1.00+0.30*0.60*1.00+0.50*0.30*1.00+0)=0.430
beta6(1)=[sigma a1jbj(o7)]beta7(j)=(0.50*0.50*0.46+0.20*0.60*0.51+0.30*0.30*0.43+0)=0.215
beta6(2)=[sigma a2jbj(o7)]beta7(j)=(0.30*0.50*0.46+0.50*0.60*0.51+0.20*0.30*0.43+0)=0.248
beta6(3)=[sigma a3jbj(o7)]beta7(j)=(0.20*0.50*0.46+0.30*0.60*0.51+0.50*0.30*0.43+0)=0.202
beta5(1)=[sigma a1jbj(o6)]beta6(j)=(0.50*0.50*0.21+0.20*0.40*0.25+0.30*0.70*0.20+0)=0.116
beta5(2)=[sigma a2jbj(o6)]beta6(j)=(0.30*0.50*0.21+0.50*0.40*0.25+0.20*0.70*0.20+0)=0.110
beta5(3)=[sigma a3jbj(o6)]beta6(j)=(0.20*0.50*0.21+0.30*0.40*0.25+0.50*0.70*0.20+0)=0.122
beta4(1)=[sigma a1jbj(o5)]beta5(j)=(0.50*0.50*0.12+0.20*0.60*0.11+0.30*0.30*0.12+0)=0.053
beta4(2)=[sigma a2jbj(o5)]beta5(j)=(0.30*0.50*0.12+0.50*0.60*0.11+0.20*0.30*0.12+0)=0.058
beta4(3)=[sigma a3jbj(o5)]beta5(j)=(0.20*0.50*0.12+0.30*0.60*0.11+0.50*0.30*0.12+0)=0.050
beta3(1)=[sigma a1jbj(o4)]beta4(j)=(0.50*0.50*0.05+0.20*0.40*0.06+0.30*0.70*0.05+0)=0.028
beta3(2)=[sigma a2jbj(o4)]beta4(j)=(0.30*0.50*0.05+0.50*0.40*0.06+0.20*0.70*0.05+0)=0.026
beta3(3)=[sigma a3jbj(o4)]beta4(j)=(0.20*0.50*0.05+0.30*0.40*0.06+0.50*0.70*0.05+0)=0.030
beta2(1)=[sigma a1jbj(o3)]beta3(j)=(0.50*0.50*0.03+0.20*0.40*0.03+0.30*0.70*0.03+0)=0.015
beta2(2)=[sigma a2jbj(o3)]beta3(j)=(0.30*0.50*0.03+0.50*0.40*0.03+0.20*0.70*0.03+0)=0.014
beta2(3)=[sigma a3jbj(o3)]beta3(j)=(0.20*0.50*0.03+0.30*0.40*0.03+0.50*0.70*0.03+0)=0.016
beta1(1)=[sigma a1jbj(o2)]beta2(j)=(0.50*0.50*0.02+0.20*0.60*0.01+0.30*0.30*0.02+0)=0.007
beta1(2)=[sigma a2jbj(o2)]beta2(j)=(0.30*0.50*0.02+0.50*0.60*0.01+0.20*0.30*0.02+0)=0.007
beta1(3)=[sigma a3jbj(o2)]beta2(j)=(0.20*0.50*0.02+0.30*0.60*0.01+0.50*0.30*0.02+0)=0.006
P(O|lambda)=0.2*0.5*0.00698+0.3*0.4*0.00741+0.5*0.7*0.00647+0=0.003852 參考資料
![監督學習 | 內建學習 之Bagging、随機森林及Sklearn實作內建學習1. 投票分類器2. Bagging & Pasting3. 随機森林參考資料[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)