本文是李航老師的《統計學習方法》^[1]一書的代碼復現。

作者：黃海廣^[2]

備注：代碼都可以在github^[3]中下載。

我將陸續將代碼發布在公眾號“機器學習初學者”，敬請關注。

代碼目錄

• 第 1 章 統計學習方法概論

• 第 2 章 感知機

• 第 3 章 k 近鄰法

• 第 4 章 樸素貝葉斯

• 第 5 章 決策樹

• 第 6 章 邏輯斯諦回歸

• 第 7 章 支持向量機

• 第 8 章 提升方法

• 第 9 章 EM 算法及其推廣

• 第 10 章 隱馬爾可夫模型

• 第 11 章 條件隨機場

• 第 12 章 監督學習方法總結

代碼參考：wzyonggege^[4],WenDesi^[5],火燙火燙的^[6]

第 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)) # betafor 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 + 1realI = i + 1print('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 + 1indexOfO = V.index(O[t]) # 找出序列對應的索引for i in range(N):realI = i + 1if t == 0:deltas[i][t] = PI[0][i] * B[i][indexOfO]psis[i][t] = 0print('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

參考資料

[1] 《統計學習方法》:?https://baike.baidu.com/item/統計學習方法/10430179
[2] 黃海廣:?https://github.com/fengdu78
[3] github:?https://github.com/fengdu78/lihang-code
[4] wzyonggege:?https://github.com/wzyonggege/statistical-learning-method
[5] WenDesi:?https://github.com/WenDesi/lihang_book_algorithm
[6] 火燙火燙的:?https://blog.csdn.net/tudaodiaozhale

往期精彩回顧

• 那些年做的學術公益-你不是一個人在戰斗

• 適合初學者入門人工智能的路線及資料下載

• 吳恩達機器學習課程筆記及資源（github標星12000+，提供百度云鏡像）
  

• 吳恩達深度學習筆記及視頻等資源（github標星8500+，提供百度云鏡像）

• 《統計學習方法》的python代碼實現（github標星7200+）

• 機器學習的數學精華（在線閱讀版）
  

備注：加入本站微信群或者qq群，請回復“加群”

加入知識星球（4300+用戶，ID：92416895），請回復“知識星球”

總結

以上是生活随笔為你收集整理的复现经典：《统计学习方法》第 10 章 隐马尔可夫模型的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。