HMM-前向後向算法與實作

目錄

• 基本要素
• HMM三大問題
• 機率計算問題
  • 前向算法
  • 後向算法
  • 前向-後向算法

基本要素

• 狀态 \(N\)個
• 狀态序列 \(S = s_1,s_2,...\)
• 觀測序列 \(O=O_1,O_2,...\)
• \(\lambda(A,B,\pi)\)
  • 狀态轉移機率 \(A = \{a_{ij}\}\)
  • 發射機率 \(B = \{b_{ik}\}\)
  • 初始機率分布 \(\pi = \{\pi_i\}\)
• 觀測序列生成過程
  • 初始狀态
  • 選擇觀測
  • 狀态轉移
  • 傳回step2

HMM三大問題

• 機率計算問題（評估問題）

給定觀測序列 \(O=O_1O_2...O_T\)，模型 \(\lambda (A,B,\pi)\)，計算 \(P(O|\lambda)\)，即計算觀測序列的機率

• 解碼問題

給定觀測序列 \(O=O_1O_2...O_T\)，模型 \(\lambda (A,B,\pi)\)，找到對應的狀态序列 \(S\)

• 學習問題

給定觀測序列 \(O=O_1O_2...O_T\)，找到模型參數 \(\lambda (A,B,\pi)\)，以最大化 \(P(O|\lambda)\)，

機率計算問題

給定模型 \(\lambda\) 和觀測序列 \(O\)，如何計算\(P(O| \lambda)\)？

暴力枚舉每一個可能的狀态序列 \(S\)

• 對每一個給定的狀态序列

  \[P(O|S,\lambda) = \prod^T_{t=1} P(O_t|s_t,\lambda) =\prod^T_{t=1} b_{s_tO_t}

  \]

• 一個狀态序列的産生機率

  \[P(S|\lambda) = P(s_1)\prod^T_{t=2}P(s_t|s_{t-1})=\pi_1\prod^T_{t=2}a_{s_{t-1}s_t}

  \]

• 聯合機率

  \[P(O,S|\lambda) = P(S|\lambda)P(O|S,\lambda) =\pi_1\prod^T_{t=2}a_{s_{t-1}s_t}\prod^T_{t=1} b_{s_tO_t}

  \]

• 考慮所有的狀态序列

  \[P(O|\lambda)=\sum_S\pi_1b_{s_1O_1}\prod^T_{t=2}a_{s_{t-1}s_t}b_{s_tO_t}

  \]

\(O\) 可能由任意一個狀态得到，是以需要将每個狀态的可能性相加。

這樣做什麼問題？時間複雜度高達 \(O(2TN^T)\)。每個序列需要計算 \(2T\) 次，一共 \(N^T\) 個序列。

前向算法

在時刻 \(t\)，狀态為 \(i\) 時，前面的時刻觀測到 \(O_1,O_2, ..., O_t\) 的機率，記為 \(\alpha _i(t)\) ：

\[\alpha_{i}(t)=P\left(O_{1}, O_{2}, \ldots O_{t}, s_{t}=i | \lambda\right)

\]

當 \(t=1\) 時，輸出為 \(O_1\)，假設有三個狀态，\(O_1\) 可能是任意一個狀态發出，即

\[P(O_1|\lambda) = \pi_1b_1(O_1)+\pi_2b_2(O_1)+\pi_2b_3(O_1) = \alpha_1(1)+\alpha_2(1)+\alpha_3(1)

\]

當 \(t=2\) 時，輸出為 \(O_1O_2\) ，\(O_2\) 可能由任一個狀态發出，同時産生 \(O_2\) 對應的狀态可以由 \(t=1\) 時刻任意一個狀态轉移得到。假設 \(O_2\) 由狀态

1

發出，如下圖

\[P(O_1O_2,s_2=q_1|\lambda) = \pi_1b_1(O_1)a_{11}b_1(O_2)+\pi_2b_2(O_1)a_{21}b_1(O_2)+\pi_2b_3(O_1)a_{31}b_1(O_2) \\

=\bold{\alpha_1(1)}a_{11}b_1(O_2)+\bold{\alpha_2(1)}a_{21}b_1(O_2)+\bold{\alpha_3(1)}a_{31}b_1(O_2) = \bold{\alpha_1(2)}

\]

同理可得 \(\alpha_2(2),\alpha_3(2)\)

\[\bold{\alpha_2(2)} = P(O_1O_2,s_2=q_2|\lambda)

=\bold{\alpha_1(1)}a_{12}b_2(O_2)+\bold{\alpha_2(1)}a_{22}b_2(O_2)+\bold{\alpha_3(1)}a_{32}b_2(O_2)

\\

\bold{\alpha_3(2)} = P(O_1O_2,s_2=q_3|\lambda)

=\bold{\alpha_1(1)}a_{13}b_3(O_2)+\bold{\alpha_2(1)}a_{23}b_3(O_2)+\bold{\alpha_3(1)}a_{33}b_3(O_2)

\]

是以

\[P(O_1O_2|\lambda) =P(O_1O_2,s_2=q_1|\lambda)+ P(O_1O_2,s_2=q_2|\lambda) +P(O_1O_2,s_2=q_3|\lambda)\\

= \alpha_1(2)+\alpha_2(2)+\alpha_3(2)

\]

是以前向算法過程如下：

​ step1：初始化 \(\alpha_i(1)= \pi_i*b_i(O_1)\)

​ step2：計算 \(\alpha_i(t) = (\sum^{N}_{j=1} \alpha_j(t-1)a_{ji})b_i(O_{t})\)

​ step3：\(P(O|\lambda) = \sum^N_{i=1}\alpha_i(T)\)

相比暴力法，時間複雜度降低了嗎？

目前時刻有 \(N\) 個狀态，每個狀态可能由前一時刻 \(N\) 個狀态中的任意一個轉移得到，是以單個時刻的時間複雜度為 \(O(N^2)\),總時間複雜度為 \(O(TN^2)\)

代碼實作

例子：

假設從三個 袋子

{1,2,3}

中 取出 4 個球

O={red,white,red,white}

，模型參數\(\lambda = (A,B,\pi)\) 如下，計算序列

O

出現的機率

#狀态 1 2 3
A = [[0.5,0.2,0.3],
	 [0.3,0.5,0.2],
	 [0.2,0.3,0.5]]

pi = [0.2,0.4,0.4]

# red white
B = [[0.5,0.5],
	 [0.4,0.6],
	 [0.7,0.3]]
           

​ step1：初始化 \(\alpha_i(1)= \pi_i*b_i(O_1)\)

​ step2：計算 \(\alpha_i(t) = (\sum^{N}_{j=1} \alpha_j(t-1)a_{ji})b_i(O_{t})\)

​ step3：\(P(O|\lambda) = \sum^N_{i=1}\alpha_i( T)\)

#前向算法
def hmm_forward(A,B,pi,O):
    T = len(O)
    N = len(A[0])
    #step1 初始化
    alpha = [[0]*T for _ in range(N)]
    for i in range(N):
        alpha[i][0] = pi[i]*B[i][O[0]]

    #step2 計算alpha(t)
    for t in range(1,T):
        for i in range(N):
            temp = 0
            for j in range(N):
                temp += alpha[j][t-1]*A[j][i]
            alpha[i][t] = temp*B[i][O[t]]
            
    #step3
    proba = 0
    for i in range(N):
        proba += alpha[i][-1]
    return proba,alpha

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]]
pi = [0.2,0.4,0.4]
O = [0,1,0,1]
hmm_forward(A,B,pi,O)  #結果為 0.06009
           

結果

後向算法

在時刻 \(t\)，狀态為 \(i\) 時，觀測到 \(O_{t+1},O_{t+2}, ..., O_T\) 的機率，記為 \(\beta _i(t)\) ：

\[\beta_{i}(t)=P\left(O_{t+1},O_{t+2}, ..., O_T | s_{t}=i, \lambda\right)

\]

當 \(t=T\) 時，由于 \(T\) 時刻之後為空，沒有觀測，是以 \(\beta_i(t)=1\)

當 \(t = T-1\) 時，觀測 \(O_T\) ，\(O_T\) 可能由任意一個狀态産生

\[\beta_i(T-1) = P(O_T|s_{t}=i,\lambda) = a_{i1}b_1(O_T)\beta_1(T)+a_{i2}b_2(O_T)\beta_2(T)+a_{i3}b_3(O_T)\beta_3(T)

\]

當 \(t=1\) 時，觀測為 \(O_{2},O_{3}, ..., O_T\)

\[\begin{aligned}

\beta_1(1)

&= P(O_{2},O_{3}, ..., O_T|s_1=1,\lambda)\\

&=a_{11}b_1(O_2)\beta_1(2)+a_{12}b_2(O_2)\beta_2(2)+a_{13}b_3(O_2)\beta_3(2)

\\

\quad

\\

\beta_2(1)

&= P(O_{2},O_{3}, ..., O_T|s_1=2,\lambda)\\

&=a_{21}b_1(O_2)\beta_1(2)+a_{22}b_2(O_2)\beta_2(2)+a_{23}b_3(O_2)\beta_3(2)

\\

\quad

\\

\beta_3(1)

&=P(O_{2},O_{3}, ..., O_T|s_1=3,\lambda)\\

&=a_{31}b_1(O_2)\beta_1(2)+a_{32}b_2(O_2)\beta_2(2)+a_{33}b_3(O_2)\beta_3(2)

\end{aligned}

\]

是以

\[P(O_{2},O_{3}, ..., O_T|\lambda) = \beta_1(1)+\beta_2(1)+\beta_3(1)

\]

後向算法過程如下：

​ step1：初始化 \(\beta_i(T)=1\)

​ step2：計算 \(\beta_i(t) = \sum^N_{j=1}a_{ij}b_j(O_{t+1})\beta_j(t+1)\)

​ step3：\(P(O|\lambda) = \sum^N_{i=1}\pi_ib_i(O_1)\beta_i(1)\)

• 時間複雜度 \(O(N^2T)\)

代碼實作

還是上面的例子

#後向算法
def hmm_backward(A,B,pi,O):
    T = len(O)
    N = len(A[0])
    #step1 初始化
    beta = [[0]*T for _ in range(N)]
    for i in range(N):
        beta[i][-1] = 1
        
    #step2 計算beta(t)
    for t in reversed(range(T-1)):
        for i in range(N):
            for j in range(N):
                beta[i][t]  += A[i][j]*B[j][O[t+1]]*beta[j][t+1]
            
    #step3
    proba = 0
    for i in range(N):
        proba += pi[i]*B[i][O[0]]*beta[i][0]
    return proba,beta

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]]
pi = [0.2,0.4,0.4]
O = [0,1,0,1]
hmm_backward(A,B,pi,O)  #結果為 0.06009
           

結果

前向-後向算法

回顧前向、後向變量：

• \(a_i(t)\) 時刻 \(t\)，狀态為 \(i\) ，觀測序列為 \(O_1,O_2, ..., O_t\) 的機率
• \(\beta_i(t)\) 時刻 \(t\)，狀态為 \(i\) ，觀測序列為 \(O_{t+1},O_{t+2}, ..., O_T\) 的機率

\[\begin{aligned}

P(O,s_t=i|\lambda)

&= P(O_1,O_2, ..., O_T,s_t=i|\lambda)\\

&= P(O_1,O_2, ..., O_t,s_t=i,O_{t+1},O_{t+2}, ..., O_T|\lambda)\\

&= P(O_1,O_2, ..., O_t,s_t=i|\lambda)*P(O_{t+1},O_{t+2}, ..., O_T|O_1,O_2, ..., O_t,s_t=i,\lambda) \\

&= P(O_1,O_2, ..., O_t,s_t=i|\lambda)*P(O_{t+1},O_{t+2}, ..., O_T,s_t=i|\lambda)\\

&= a_i(t)*\beta_i(t)

\end{aligned}

\]

即在給定的狀态序列中，\(t\) 時刻狀态為 \(i\) 的機率。

使用前後向算法可以計算隐狀态，記 \(\gamma_i(t) = P(s_t=i|O,\lambda)\) 表示時刻 \(t\) 位于隐狀态 \(i\) 的機率

\[P\left(s_{t}=i, O | \lambda\right)=\alpha_{i}(t) \beta_{i}(t)

\]

\[\begin{aligned}

\gamma_{i}(t)

&=P\left(s_{t}={i} | O, \lambda\right)=\frac{P\left(s_{t}={i}, O | \lambda\right)}{P(O | \lambda)} \\

&=\frac{\alpha_{i}(t) \beta_{i}(t)}{P(O | \lambda)}=\frac{\alpha_{i}(t) \beta_{i}(t)}{\sum_{i=1}^{N} \alpha_{i}(t) \beta_{i}(t)}

\end{aligned}

\]

references：

[1] https://www.cs.sjsu.edu/~stamp/RUA/HMM.pdf

[2]https://www.cnblogs.com/fulcra/p/11065474.html

[3] https://www.cnblogs.com/sjjsxl/p/6285629.html

[4] https://blog.csdn.net/xueyingxue001/article/details/52396494