機器學習之隱馬爾科夫模型(HMM)原理及Python實現 (大章節)

HMM

隱馬爾可夫模型（hidden Markov model, HMM）是可用於標註問題的統計學模型，是生成模型。

本章節內容參考李航博士的《統計學習方法》
本章節添加了一些結論性結果的推導過程。

1. 從一個天然語言處理的例子開始

例若有三個個句子：
句子一：我/名詞 看見/動詞 貓/名詞
句子二：貓/名詞 是/動詞 可愛的/形容詞
句子三：我/名詞 是/動詞 可愛的/形容詞
通常只能觀察到具體的詞，因此像"我 看見 貓 …"是觀測集合，而詞性如"名詞 動詞 形容詞 …"是狀態序列

設 $Q$是全部可能的狀態集合， $V$是全部可能的觀測集合：

$Q = \{q_1, q_2, ..., q_N\}, V=\{v_1, v_2, ..., v_M\}$

其中， N是可能的狀態數，M是可能的觀測數。

例如： $Q=\{名詞，動詞，形容詞 \}，V=\{我， 看見， 貓， 是，可愛的\}，N=3, M=5$

$I$是長度爲 $T$的狀態序列， $O$是對應的觀測序列：

$I = \{i_1, i_2,..., i_T \}, O=\{o_1, o_2,..., o_T\}$

例如： $I=(名詞，動詞，名詞)， O=(我，看見，貓)$

$A$是狀態轉移矩陣：

$A=[a_{ij}]_{N*N} \tag1$

其中，

$a_{ij} = p(i_{t+1}=q_j|i_t=q_i), i=1,2,...,N; j=1,2,...,N \tag2$

例如：

轉態轉移機率	名詞	動詞	形容詞
名詞	0	1	0
動詞	1/3	0	2/3
形容詞	1/3	1/3	1/3

$B$是觀測機率矩陣，也就是發射矩陣：

$B=[b_j(k)]_{N*M} \tag3$

其中，

$b_j(k) = p(o_t=v_k|i_t=q_j), k=1,2,...,M; j=1,2,...,N \tag4$

例如：

觀測矩陣機率	我	看見	貓	是	可愛的
名詞	1	0	1	0	0
動詞	0	1	0	1	0
形容詞	0	0	0	0	1

$\pi$是初始狀態機率向量：

$\pi = (\pi_i) \tag5$

其中，

$\pi_i = p(i_1 = q_i), i = 1,2,...,N \tag6$

$A,B$和 $\pi$是HMM的參數，用 $\lambda$表示：

$\lambda = (A,B,\pi) \tag7$

例如：

名詞	動詞	形容詞
1	0	0

隱馬爾可夫的三個基本問題
1.機率計算問題。給定模型 $\lambda=(A,B,\pi)$和觀測序列 $O=(o_1,o_2,...,o_T)$,計算在已知模型參數的狀況下，觀測序列的機率，即 $p(O|\lambda)$。
2.學習問題。已知觀測序列 $O=(o_1,o_2,...,o_T)$,估計模型參數 $\lambda=(A,B,\pi)$，使 $p(O|\lambda)$最大。
3.預測問題，也稱解碼問題。已知模型 $\lambda=(A,B,\pi)$和 $O=(o_1,o_2,...,o_T)$，求條件機率最大 $p(I|O)$最大的狀態序列 $I=(i_1,i_2,...,i_T)$。

2. 機率預測問題

機率問題預測用直接計算法，計算複雜度高，能夠採用動態規劃形式的前向和後向算法下降計算複雜度。
爲了表示方便，記：

$(o_{1:t} )= (o_1,o_2,...,o_n); (o_{t_:T})=(o_t,o_{t+1},...,o_T)$

2.1 前向算法

接下來就是解前向機率 $p(i_t,o_{1:t}|\lambda)$：

$\begin{aligned} p(i_t,o_{1:t}|\lambda) &=\sum_{i_{t-1}} p(i_{t-1},i_t,o_{1:t-1},o_t|\lambda) \\ &=\sum_{i_{t-1}} p(o_t|i_{t-1},i_t,o_{1:t-1},\lambda)p(i_t|i_{t-1},o_{1:t-1},\lambda)p(i_{t-1},o_{1:t-1}|\lambda) \end{aligned}$

由隱馬爾科夫的條件獨立性假設可得：

$p(o_t|i_{t-1},i_t,o_{1:t-1},\lambda) = p(o_t|i_t,\lambda)$

$p(i_t|i_{t-1},o_{1:t-1},\lambda)=p(i_t|i_{t-1},\lambda)$

故

$p(i_t,o_{1:t}|\lambda)=\sum_{i_{t-1}} p(o_t|i_t,\lambda) p(i_t|i_{t-1},\lambda)p(i_{t-1},o_{1:t-1}|\lambda)=[\sum_{i_{t-1} } p(i_{t-1},o_{1:t-1}|\lambda) p(i_t|i_{t-1},\lambda)] p(o_t|i_t,\lambda)$

設：

$\alpha_{t+1}(i) = p(o_{1:t+1},i_{t+1}=q_i|\lambda) \tag8$

且：

$p(i_{t+1}=q_i|i_t=q_j,\lambda)] = a_{ji}$

$p(o_{t+1}|i_{t+1},\lambda)=b_i(o_{t+1})$

則：

$\alpha_{t+1}(i)=[\sum_{j=1}^N \alpha_t(j)a_{ji}]b_i(o_{t+1}) \tag9$

因此前向算法就可迭代進行。

前向算法：
1.初值

$\alpha_1(i) = \pi_ib_i(o_1)$

2.遞推 $t=1,2,...,T-1$

$\alpha_{t+1}(i)=[\sum_{j=1}^N \alpha_t(j)a_{ji}]b_i(o_{t+1})$

3.終止
$p(O|\lambda) = \sum_{i=1}^N \alpha_T(i)$

2.2 後向算法

後向算法解決後向機率 $p(o_{t+1:T}|i_t, \lambda)$:

$\begin{aligned} p(o_{t+1:T}|i_t, \lambda) &= \sum_{i_{t+1}} p(i_{t+1},o_{t+1},o_{t+2:T} | i_t, \lambda) \\ &= \sum_{i_{t+1}} p(o_{t+2:T}|i_{t+1}, i_t, o_{t+1}, \lambda) p(o_{t+1}|i_{t+1}, i_t, \lambda) p(i_{t+1}|i_t,\lambda)\\ \end{aligned}$

由隱馬爾科夫的條件獨立假設得：

$p(o_{t+2:T}|i_{t+1}, i_t, o_{t+1}, \lambda)=p(o_{t+2:T}|i_{t+1}, \lambda)$

$p(o_{t+1}|i_{t+1}, i_t, \lambda) = p(o_{t+1}|i_{t+1}, \lambda)$

設：

$\beta_t(i) = p(o_{t+1:T}|i_t=q_i, \lambda) \tag{10}$

又：

$p(i_{t+1}=q_j|i_t=q_i,\lambda) = a_{ij}$

$p(o_{t+1}|i_{t+1}=q_j, \lambda) = b_j(o_{t+1})$

則：

$\beta_t(i) = \sum_{j=1}^N a_{ij} b_j(o_{t+1}) \beta_{t+1}(i) \tag{11}$

後向算法：
(1)

$\beta_T (i) = 1$

(2) 對t=T-1,T-2,…,1

$\beta_t(i) = \sum_{j=1}^N a_{ij} b_j(o_{t+1}) \beta_{t+1}(i)$

(3)

$p(O|\lambda) = \sum_{i=1}^N \pi_i b_i(o_1) \beta_1(i)$

2.3 一些機率與指望值

這兩個指望值都是後面EM算法用到的中間參量
1.計算 $t$時刻處於狀態 $q_i$的機率。
機率計算問題是計算 $p(O|\lambda)$，則有：

$p(O|\lambda)=\sum_{i_t}p(O,i_t|\lambda)$

依據隱馬爾科夫的獨立性假設：

$p(o_{t+1:T}|i_t,o_{1:t}, \lambda) = p(o_{t+1:T}|i_t, \lambda)$

因此：

$\begin{aligned} p(O|\lambda) &=\sum_{i_t}p(O,i_t|\lambda) \\ &=\sum_{i_t} p(o_{t+1:T}|i_t,o_{1:t}, \lambda) p(i_t,o_{1:t}|\lambda) \\ &=\sum_{i_t} p(o_{t+1:T}|i_t, \lambda) p(i_t,o_{1:t}|\lambda) \\ \end{aligned}$

又有：

$\alpha_t(i) = p(o_{1:t},i_t=q_i|\lambda) \tag{12}$

$\beta_t(i) = p(o_{t+1:T}|i_t=q_i, \lambda) \tag{13}$

故：

$p(O,i_t=q_i|\lambda) = p(o_{t+1:T}|i_t=q_i, \lambda) p(i_t=q_i,o_{1:t}|\lambda) = \alpha_t(i) \beta_t(i)$

$p(O|\lambda) = \sum_{i_t} \alpha_t(i) \beta_t(i)$

設：

$\gamma_t(i) = p(i_t=q_i|O,\lambda)$

因而能夠獲得：

$\gamma_t(i) = p(i_t=q_i|O,\lambda) = \frac {p(i_t=q_i,O|\lambda)}{p(O|\lambda)} = \frac {\alpha_t(i) \beta_t(i)}{\sum_{j=1}^N \alpha_t(j) \beta_t(j)} \tag{14}$

2.計算計算 $t$時刻處於狀態 $q_i$且計算 $t+1$時刻處於狀態 $q_j$的機率

$\begin{aligned} p(O|\lambda) &=\sum_{i_t} \sum_{i_{t+1}} p(O,i_t, i_{t+1}|\lambda) \\ &=\sum_{i_t} \sum_{i_{t+1}} p(o_{1:t},o_{t+1},o_{t+2:T},i_t, i_{t+1}|\lambda) \\ &=\sum_{i_t} \sum_{i_{t+1}} p(o_{t+2:T}|o_{1:t},o_{t+1},i_t, i_{t+1},\lambda)p(o_{t+1}|o_{1:t},i_t,i_{t+1},\lambda) p(i_{t+1}|i_t,o_{1:t},\lambda) p(i_t,o_{1:t}|\lambda) \\ \end{aligned}$

由隱馬爾科夫的獨立性假設可得：

$p(O|\lambda) = \sum_{i_t} \sum_{i_{t+1}} p(o_{t+2:T}| i_{t+1},\lambda)p(o_{t+1}|i_{t+1},\lambda) p(i_{t+1}|i_t,\lambda) p(i_t,o_{1:t}|\lambda)$

設：

$\xi_t(i,j)=p(i_t=q_i,i_{t+1}=q_j|O,\lambda)$

又有公式(2)(4)(12)(13)

得：

$\xi_t(i,j) = \frac {p(i_t=q_i,i_{t+1}=q_j|O,\lambda)}{p(O|\lambda)} =\frac {\alpha_t(i) a_{ij} b_j(o_{t+1}) \beta_{t+1}(j)} {\sum_{i=1}^N \sum_{j=1}^N \alpha_t(i) a_{ij} b_j(o_{t+1}) \beta_{t+1}(j)} \tag{15}$

3. 學習問題

3.1 監督學習

若是有標記好狀態序列的樣本，那就太好辦了，直接將接個矩陣統計的各個維度定義後進行統計就能夠了。統計過程當中注意機率之和爲一的約束。

3.2 無監督學習

若是沒有標記狀態序列的樣本，能夠用Baum-Welch算法(EM算法)實現。

已知：包含 $S$個長度爲 $T$的觀測序列的觀測序列 $\{O_1,O_2,...,O_S \}$
目標：學習隱馬爾可夫模型的參數 $\lambda=(A,B,\pi)$

記觀測數據 $O$,隱數據 $I$，那麼隱馬爾可夫模型能夠表示爲：

$p(O|\lambda) = \sum_I p(O|I,\lambda) p(I|\lambda)$

E步：

由於對 $\lambda$而言， $1/p(O| \overline \lambda)$是常數項，因此

$\begin{aligned} Q(\lambda,\overline \lambda) &= E_I[\log p(O,I|\lambda)|O, \overline \lambda] \\ &= \sum_I \log p(O,I|\lambda) p(I|O,\overline \lambda) \\ &= \sum_I \log p(O,I|\lambda) \frac {p(I,O|\overline \lambda)}{p(O| \overline \lambda)} \\ &= \sum_I \log p(O,I|\lambda) p(I,O|\overline \lambda) \\ \end{aligned}$

將機率計算問題2.1小姐中前向算法的遞歸公式展開就能夠獲得：

$p(O,I|\lambda) = \pi_{i_1} b_{i_1}(o_1) a_{i_1i_2} b_{i_2}(o_2) ... a_{i_{T-1}i_T} b_{iT}(o_T) = \pi_{i_1} [\prod_{t=1}^{T-1} a_{i_ti_{t+1}}][\prod_{t=1}^T b_{it}(o_t)]$

因而：

$Q(\lambda, \overline \lambda) = \sum_I \log \pi_{i_1} p(O, I| \overline \lambda) + \sum_I (\sum_{t=1}^{T-1} a_{i_ti_{t+1}}) p(O, I| \overline \lambda) + \sum_I (\sum_{t=1}^T b_{it}(o_t)) p(O, I| \overline \lambda) \tag{16}$

特此說明隱變量
隱馬爾可夫模型的隱變量就是觀測序列對應的狀態序列，因此隱變量能夠用（14）式的變量表示
後面在M步中更新模型參數的時候也用到了（15）式，是否是就說明隱變量是兩個，其實不是的，這兒只是爲了表示的方便和算法的方便。
也就是在E步中，用 $\gamma$和 $\xi$表示隱變量，只是爲了編程和表示的便利，這兩個變量在E步中信息是重複的。

M步：

1.求解 $\pi_i$
由(15)式可得：

$L(\pi_{i_1}) = \sum_I \log \pi_{i_1} p(O, I| \overline \lambda) = \sum_{i}^N \log \pi_{i_1} p(O, i_1=i| \overline \lambda)$

又由於 $\pi_i$知足約束條件 $\sum_{i=1}^N \pi_{i_1}=1$，利用拉格朗日乘子法，寫出拉格朗日函數：

$\sum_{i=1}^N \log \pi_{i} p(O, i_1=i| \overline \lambda) + \gamma(\sum_{i=1}^N \pi_{i} - 1)$

對其求偏導而且令其結果爲0得：

$\frac {\partial} {\partial \pi_i} [\sum_{i=1}^N \log \pi_{i} p(O, i=i| \overline \lambda) + \gamma(\sum_{i_1=1}^N \pi_{i} - 1)]=0 \tag{17}$

得：

$p(O, i_1=i| \overline \lambda) + \gamma \pi_i=0$

獲得：

$\pi_i = \frac {p(O, i_1=i| \overline \lambda)} {-\lambda}$

帶入 $\sum_{i=1}^N \pi_{i_1}=1$的：

$-\lambda = \sum_{i=1}^N p(O, i_1=i| \overline \lambda) = p(o|\overline \lambda)$

求得並有公式(14)：

$\pi_i = \frac {p(O, i_1=i| \overline \lambda)} {p(o|\overline \lambda)} = \gamma_1(i) \tag{18}$

2.求解 $a_{ij}$:

$L(a_{ij})=\sum_I (\sum_{t=1}^{T-1} a_{i_ti_{t+1}}) p(O, I| \overline \lambda) = \sum_{i=1}^N (\sum_{t=1}^{T-1} a_{i_ti_{t+1}}) ( \sum_{j=1}^N p(O, i_t=i, i_{t+1}=j| \overline \lambda) ) \\ = \sum_{i=1}^N \sum_{j=1}^N \sum_{t=1}^{T-1} a_{ij} p(O, i_t=i, i_{t+1}=j| \overline \lambda)$

應用約束條件 $\sum_{j=1}^N a_{ij} = 1$,用拉格朗日乘子法能夠求出：

$\sum_{i=1}^N \sum_{j=1}^N \sum_{t=1}^{T-1} a_{ij} p(O, i_t=i, i_{t+1}=j| \overline \lambda) + \lambda(\sum_{j=1}^N a_{ij} - 1)$

對上式求騙到並等於0獲得：

$\frac {\partial}{\partial a_{ij}} [\sum_{i=1}^N \sum_{j=1}^N \sum_{t=1}^{T-1} a_{ij} p(O, i_t=i, i_{t+1}=j| \overline \lambda) + \lambda(\sum_{j=1}^N a_{ij} - 1)] = 0$

獲得：

$\sum_{t=1}^{T-1} p(O, i_t=i, i_{t+1}=j| \overline \lambda) + \lambda a_{ij} = 0$

因此：

$a_{ij} = \frac {\sum_{t=1}^{T-1} p(O, i_t=i, i_{t+1}=j| \overline \lambda)}{- \lambda}$

將上式帶入 $\sum_{j=1}^N a_{ij} = 1$：

$- \lambda = \sum_{j=1}^N \sum_{t=1}^{T-1} p(O, i_t=i, i_{t+1}=j| \overline \lambda) = \sum_{t=1}^{T-1} p(O, i_t=i| \overline \lambda)$

故得：

$a_{ij} = \frac {\sum_{t=1}^{T-1} p(O, i_t=i, i_{t+1}=j| \overline \lambda)}{\sum_{t=1}^{T-1} p(O, i_t=i| \overline \lambda) } = \frac {\sum_{t=1}^{T-1} p(O, i_t=i, i_{t+1}=j| \overline \lambda) / p(o|\overline \lambda)} {\sum_{t=1}^{T-1} p(O, i_t=i| \overline \lambda) / p(o|\overline \lambda) }$

將（14）和（15）帶入的：

$a_{ij} = \frac {\sum_{t=1}^{T-1} \xi_t(i,j)} {\sum_{t=1}^{T-1} \gamma_t(i) } \tag{19}$

3.求解 $b_j{k}$:

$L(b_j{k}) = \sum_I (\sum_{t=1}^T b_{it}(o_t)) p(O, I| \overline \lambda) = \sum_{j=1}^N \sum_{t=1}^T b_{j}(o_t) p(O, i_t=j| \overline \lambda)$

在約束條件 $\sum_{k=1}^M b_j(k) = 1$的拉格朗日乘子法：

$\sum_{j=1}^N \sum_{t=1}^T b_{j}(o_t) p(O, i_t=j| \overline \lambda) + \lambda(\sum_{k=1}^M b_j(k) - 1)$

對其求偏導得：

$\frac {\partial}{\partial b_j(k)} [\sum_{j=1}^N \sum_{t=1}^T b_{j}(o_t) p(O, i_t=j| \overline \lambda) + \lambda(\sum_{k=1}^M b_j(k) - 1)] = 0$

由於只有在 $o_t=v_k$時偏導纔不會等於0，以 $I(o_t=v_k)$表示，則：

$\sum_{t=1}^T p(O, i_t=j| \overline \lambda) I(o_t=v_k) + \lambda b_{j}(o_t)I(o_t=v_k) = 0$

$b_{j}(o_t)I(o_t=v_k)$能夠寫做 $b_{j}(k)$，故：

$b_{j}(k) = \frac {\sum_{t=1}^T p(O, i_t=j| \overline \lambda) I(o_t=v_k)} {- \lambda}$

將上式帶入 $\sum_{k=1}^M b_j(k) = 1$得：

$- \lambda = \sum_{k=1}^M \sum_{t=1}^T p(O, i_t=j| \overline \lambda) I(o_t=v_k) = \sum_{t=1}^T p(O, i_t=j| \overline \lambda)$

獲得：

$b_{j}(k) = \frac {\sum_{t=1}^T p(O, i_t=j| \overline \lambda) I(o_t=v_k)} {\sum_{t=1}^T p(O, i_t=j| \overline \lambda)}$

又有（14）式可得：

$b_{j}(k) = \frac {\sum_{t=1,o_t=v_k}^T \gamma_t(j)} {\sum_{t=1}^T \gamma_t(j)} \tag{20}$

EM算法總結：
E步：

$\gamma_t(i) = p(i_t=q_i|O,\lambda) = \frac {p(i_t=q_i,O|\lambda)}{p(O|\lambda)} = \frac {\alpha_t(i) \beta_t(i)}{\sum_{j=1}^N \alpha_t(j) \beta_t(j)}$

$\xi_t(i,j) = \frac {p(i_t=q_i,i_{t+1}=q_j|O,\lambda)}{p(O|\lambda)} =\frac {\alpha_t(i) a_{ij} b_j(o_{t+1}) \beta_{t+1}(j)} {\sum_{i=1}^N \sum_{j=1}^N \alpha_t(i) a_{ij} b_j(o_{t+1}) \beta_{t+1}(j)}$

M步：
$\pi_i = \frac {p(O, i_1=i| \overline \lambda)} {p(o|\overline \lambda)} = \gamma_1(i)$

$a_{ij} = \frac {\sum_{t=1}^{T-1} \xi_t(i,j)} {\sum_{t=1}^{T-1} \gamma_t(i) }$

$b_{j}(k) = \frac {\sum_{t=1,o_t=v_k}^T \gamma_t(j)} {\sum_{t=1}^T \gamma_t(j)}$

4. 預測問題（解碼問題）

用維特比算法進行求解：
已知：模型 $\lambda=(A,B,\pi)$和 $O=(o_1,o_2,...,o_T)$
求：條件機率最大 $p(I|O,\lambda)$最大的狀態序列 $I=(i_1,i_2,...,i_T)$
由於 $p(O)$是一個定值，因此：

$\max_I p(I|O,\lambda) = \max_I p(I, O|\lambda) / p(O|\lambda) = \max_I p(I, O|\lambda)$

定義在時刻 $t$狀態爲 $i$的全部單個路徑 $(i_1,i_2,...,i_t)$中機率最大值爲：

$\delta_t(i) = \max_{i_1,i_2,...,i_{t-1}} p(i_t=i, i_{t-1:i_1},o_{t:1}|\lambda)$

遞推推導：

$\begin{aligned} &p(i_{t+1}=i,i_{t:1},o_{t+1:1}| \lambda) \\ &=p(i_{t+1}=i,i_t,i_{t-1:1},o_{t+1},o_{t:1}| \lambda) \\ &= p(o_{t+1}|i_{t+1}=i,i_t,o_{t:1},\lambda) p(i_{t+1}=i|i_t,i_{t-1:1},o_{t:1}, \lambda) p(i_t,i_{t-1:1},o_{t:1}|\lambda) \\ &= p(o_{t+1}|i_{t+1}=i,\lambda) p(i_{t+1}=i|i_t,\lambda) p(i_t,i_{t-1:1},o_{t:1}|\lambda) \\ \end{aligned}$

故：

$\delta_{t+1}(i) = \max_{i_1,i_2,...,i_{t-1}} p(i_{t+1}=i,i_{t:1},o_{t+1:1}| \lambda) = \max_{1 \le j \le N} [\delta _t(j) a_{ji}] b_i(o_{t+1}) \tag{21}$

定義在時刻 $t$狀態爲 $i$的全部單個路徑 $(i_1,i_2,...,i_{t-1})$中機率最大的第 $t-1$