隱馬爾科夫模型HMM（三）鮑姆-韋爾奇演算法求解HMM引數
在本篇我們會討論HMM模型引數求解的問題，這個問題在HMM三個問題裡算是最複雜的。在研究這個問題之前，建議先閱讀這個系列的前兩篇以熟悉HMM模型和HMM的前向後向演算法，以及EM演算法原理總結，這些在本篇裡會用到。在李航的《統計學習方法》中，這個演算法的講解只考慮了單個觀測序列的求解，因此無法用於實際多樣本觀測序列的模型求解，本文關注於如何使用多個觀測序列來求解HMM模型引數。
1. HMM模型引數求解概述
HMM模型引數求解根據已知的條件可以分為兩種情況。
第一種情況較為簡單，就是我們已知$D$個長度為$T$的觀測序列和對應的隱藏狀態序列，即$\{(O_1, I_1), (O_2, I_2), ...(O_D, I_D)\}$是已知的，此時我們可以很容易的用最大似然來求解模型引數。
假設樣本從隱藏狀態$q_i$轉移到$q_j$的頻率計數是$A_{ij}$,那麼狀態轉移矩陣求得為：$$A = \Big[a_{ij}\Big], \;其中a_{ij} = \frac{A_{ij}}{\sum\limits_{s=1}^{N}A_{is}}$$
假設樣本隱藏狀態為$q_j$且觀測狀態為$v_k$的頻率計數是$B_{jk}$,那麼觀測狀態概率矩陣為：$$B= \Big[b_{j}(k)\Big], \;其中b_{j}(k) = \frac{B_{jk}}{\sum\limits_{s=1}^{M}B_{js}}$$
假設所有樣本中初始隱藏狀態為$q_i$的頻率計數為$C(i)$,那麼初始概率分佈為：$$\Pi = \pi(i) = \frac{C(i)}{\sum\limits_{s=1}^{N}C(s)}$$
可見第一種情況下求解模型還是很簡單的。但是在很多時候，我們無法得到HMM樣本觀察序列對應的隱藏序列，只有$D$個長度為$T$的觀測序列，即$\{(O_1), (O_2), ...(O_D)\}$是已知的，此時我們能不能求出合適的HMM模型引數呢？這就是我們的第二種情況，也是我們本文要討論的重點。它的解法最常用的是鮑姆-韋爾奇演算法，其實就是基於EM演算法的求解，只不過鮑姆-韋爾奇演算法出現的時代，EM演算法還沒有被抽象出來，所以我們本文還是說鮑姆-韋爾奇演算法。
2. 鮑姆-韋爾奇演算法原理
鮑姆-韋爾奇演算法原理既然使用的就是EM演算法的原理，那麼我們需要在E步求出聯合分佈$P(O,I|\lambda)$基於條件概率$P(I|O,\overline{\lambda})$的期望，其中$\overline{\lambda}$為當前的模型引數，然後再M步最大化這個期望，得到更新的模型引數$\lambda$。接著不停的進行EM迭代，直到模型引數的值收斂為止。
首先來看看E步，當前模型引數為$\overline{\lambda}$, 聯合分佈$P(O,I|\lambda)$基於條件概率$P(I|O,\overline{\lambda})$的期望表示式為：$$L(\lambda, \overline{\lambda}) = \sum\limits_{I}P(I|O,\overline{\lambda})logP(O,I|\lambda)$$
在M步，我們極大化上式，然後得到更新後的模型引數如下： $$\overline{\lambda} = arg\;\max_{\lambda}\sum\limits_{I}P(I|O,\overline{\lambda})logP(O,I|\lambda)$$
通過不斷的E步和M步的迭代，直到$\overline{\lambda}$收斂。下面我們來看看鮑姆-韋爾奇演算法的推導過程。
3. 鮑姆-韋爾奇演算法的推導
我們的訓練資料為$\{(O_1, I_1), (O_2, I_2), ...(O_D, I_D)\}$，其中任意一個觀測序列$O_d = \{o_1^{(d)}, o_2^{(d)}, ... o_T^{(d)}\}$,其對應的未知的隱藏狀態序列表示為：$I_d = \{i_1^{(d)}, i_2^{(d)}, ... i_T^{(d)}\}$
首先看鮑姆-韋爾奇演算法的E步，我們需要先計算聯合分佈$P(O,I|\lambda)$的表示式如下：$$P(O,I|\lambda) = \prod_{d=1}^D\pi_{i_1^{(d)}}b_{i_1^{(d)}}(o_1^{(d)})a_{i_1^{(d)}i_2^{(d)}}b_{i_2^{(d)}}(o_2^{(d)})...a_{i_{T-1}^{(d)}i_T^{(d)}}b_{i_T^{(d)}}(o_T^{(d)})$$
我們的E步得到的期望表示式為：$$L(\lambda, \overline{\lambda}) = \sum\limits_{I}P(I|O,\overline{\lambda})logP(O,I|\lambda)$$
在M步我們要極大化上式。由於$P(I|O,\overline{\lambda}) = P(I,O|\overline{\lambda})/P(O|\overline{\lambda})$,而$P(O|\overline{\lambda})$是常數，因此我們要極大化的式子等價於：$$\overline{\lambda} = arg\;\max_{\lambda}\sum\limits_{I}P(O,I|\overline{\lambda})logP(O,I|\lambda)$$
我們將上面$P(O,I|\lambda)$的表示式帶入我們的極大化式子，得到的表示式如下：$$\overline{\lambda} = arg\;\max_{\lambda}\sum\limits_{d=1}^D\sum\limits_{I}P(O,I|\overline{\lambda})(log\pi_{i_1} + \sum\limits_{t=1}^{T-1}log\;a_{i_t,i_{t+1}} + \sum\limits_{t=1}^Tb_{i_t}(o_t))$$
我們的隱藏模型引數$\lambda =(A,B,\Pi)$,因此下面我們只需要對上式分別對$A,B,\Pi$求導即可得到我們更新的模型引數$\overline{\lambda}$
首先我們看看對模型引數$\Pi$的求導。由於$\Pi$只在上式中括號裡的第一部分出現，因此我們對於$\Pi$的極大化式子為：$$\overline{\pi_i} = arg\;\max_{\pi_{i_1}} \sum\limits_{d=1}^D\sum\limits_{I}P(O,I|\overline{\lambda})log\pi_{i_1} = arg\;\max_{\pi_{i}} \sum\limits_{d=1}^D\sum\limits_{i=1}^NP(O,i_1^{(d)} =i|\overline{\lambda})log\pi_{i}$$
由於$\pi_i$還滿足$\sum\limits_{i=1}^N\pi_i =1$，因此根據拉格朗日子乘法，我們得到$\pi_i$要極大化的拉格朗日函式為：$$arg\;\max_{\pi_{i}}\sum\limits_{d=1}^D\sum\limits_{i=1}^NP(O,i_1^{(d)} =i|\overline{\lambda})log\pi_{i} + \gamma(\sum\limits_{i=1}^N\pi_i -1)$$
其中，$\gamma$為拉格朗日系數。上式對$\pi_i$求偏導數並令結果為0， 我們得到：$$\sum\limits_{d=1}^DP(O,i_1^{(d)} =i|\overline{\lambda}) + \gamma\pi_i = 0$$
令$i$分別等於從1到$N$，從上式可以得到$N$個式子，對這$N$個式子求和可得：$$\sum\limits_{d=1}^DP(O|\overline{\lambda}) + \gamma = 0 $$
從上兩式消去$\gamma$,得到$\pi_i$的表示式為：$$\pi_i =\frac{\sum\limits_{d=1}^DP(O,i_1^{(d)} =i|\overline{\lambda})}{\sum\limits_{d=1}^DP(O|\overline{\lambda})} = \frac{\sum\limits_{d=1}^DP(O,i_1^{(d)} =i|\overline{\lambda})}{DP(O|\overline{\lambda})} = \frac{\sum\limits_{d=1}^DP(i_1^{(d)} =i|O, \overline{\lambda})}{D} = \frac{\sum\limits_{d=1}^DP(i_1^{(d)} =i|O^{(d)}, \overline{\lambda})}{D}$$
利用我們在隱馬爾科夫模型HMM（二）前向後向演算法評估觀察序列概率裡第二節中前向概率的定義可得：$$P(i_1^{(d)} =i|O^{(d)}, \overline{\lambda}) = \gamma_1^{(d)}(i)$$
因此最終我們在M步$\pi_i$的迭代公式為：$$\pi_i = \frac{\sum\limits_{d=1}^D\gamma_1^{(d)}(i)}{D}$$
現在我們來看看$A$的迭代公式求法。方法和$\Pi$的類似。由於$A$只在最大化函式式中括號裡的第二部分出現，而這部分式子可以整理為：$$\sum\limits_{d=1}^D\sum\limits_{I}\sum\limits_{t=1}^{T-1}P(O,I|\overline{\lambda})log\;a_{i_t,i_{t+1}} = \sum\limits_{d=1}^D\sum\limits_{i=1}^N\sum\limits_{j=1}^N\sum\limits_{t=1}^{T-1}P(O,i_t^{(d)} = i, i_{t+1}^{(d)} = j|\overline{\lambda})log\;a_{ij}$$
由於$a_{ij}$還滿足$\sum\limits_{j=1}^Na_{ij} =1$。和求解$\pi_i$類似，我們可以用拉格朗日子乘法並對$a_{ij}$求導，並令結果為0，可以得到$a_{ij}$的迭代表示式為：$$a_{ij} = \frac{\sum\limits_{d=1}^D\sum\limits_{t=1}^{T-1}P(O^{(d)}, i_t^{(d)} = i, i_{t+1}^{(d)} = j|\overline{\lambda})}{\sum\limits_{d=1}^D\sum\limits_{t=1}^{T-1}P(O^{(d)}, i_t^{(d)} = i|\overline{\lambda})}$$
利用隱馬爾科夫模型HMM（二）前向後向演算法評估觀察序列概率裡第二節中前向概率的定義和第五節$\xi_t(i,j)$的定義可得們在M步$a_{ij}$的迭代公式為：$$a_{ij} = \frac{\sum\limits_{d=1}^D\sum\limits_{t=1}^{T-1}\xi_t^{(d)}(i,j)}{\sum\limits_{d=1}^D\sum\limits_{t=1}^{T-1}\gamma_t^{(d)}(i)}$$
現在我們來看看$B$的迭代公式求法。方法和$\Pi$的類似。由於$B$只在最大化函式式中括號裡的第三部分出現，而這部分式子可以整理為：$$\sum\limits_{d=1}^D\sum\limits_{I}\sum\limits_{t=1}^{T}P(O,I|\overline{\lambda})log\;b_{i_t}(o_t) = \sum\limits_{d=1}^D\sum\limits_{j=1}^N\sum\limits_{t=1}^{T}P(O,i_t^{(d)} = j|\overline{\lambda})log\;b_{j}(o_t)$$
由於$b_{j}(o_t)$還滿足$\sum\limits_{k=1}^Mb_{j}(o_t =v_k) =1$。和求解$\pi_i$類似，我們可以用拉格朗日子乘法並對$b_{j}(k)$求導，並令結果為0，得到$b_{j}(k)$的迭代表示式為：$$b_{j}(k) = \frac{\sum\limits_{d=1}^D\sum\limits_{t=1}^{T}P(O,i_t^{(d)} = j|\overline{\lambda})I(o_t^{(d)}=v_k)}{\sum\limits_{d=1}^D\sum\limits_{t=1}^{T}P(O,i_t^{(d)} = j|\overline{\lambda})}$$
其中$I(o_t^{(d)}=v_k)$當且僅當$o_t^{(d)}=v_k$時為1，否則為0. 利用隱馬爾科夫模型HMM（二）前向後向演算法評估觀察序列概率裡第二節中前向概率的定義可得$b_{j}(o_t)$的最終表示式為：$$b_{j}(k) = \frac{\sum\limits_{d=1}^D\sum\limits_{t=1, o_t^{(d)}=v_k}^{T}\gamma_t^{(d)}(j)}{\sum\limits_{d=1}^D\sum\limits_{t=1}^{T}\gamma_t^{(d)}(j)}$$
有了$\pi_i, a_{ij},b_{j}(k)$的迭代公式，我們就可以迭代求解HMM模型引數了。
4. 鮑姆-韋爾奇演算法流程總結
這裡我們概括總結下鮑姆-韋爾奇演算法的流程。
輸入： $D$個觀測序列樣本$\{(O_1), (O_2), ...(O_D)\}$
輸出：HMM模型引數
1)隨機初始化所有的$\pi_i, a_{ij},b_{j}(k)$
2) 對於每個樣本$d = 1,2,...D$，用前向後向演算法計算$\gamma_t^{(d)}(i)，\xi_t^{(d)}(i,j), t =1,2...T$
3) 更新模型引數：
$$\pi_i = \frac{\sum\limits_{d=1}^D\gamma_1^{(d)}(i)}{D}$$
$$a_{ij} = \frac{\sum\limits_{d=1}^D\sum\limits_{t=1}^{T-1}\xi_t^{(d)}(i,j)}{\sum\limits_{d=1}^D\sum\limits_{t=1}^{T-1}\gamma_t^{(d)}(i)}$$
$$b_{j}(k) = \frac{\sum\limits_{d=1}^D\sum\limits_{t=1, o_t^{(d)}=v_k}^{T}\gamma_t^{(d)}(j)}{\sum\limits_{d=1}^D\sum\limits_{t=1}^{T}\gamma_t^{(d)}(j)}$$
4) 如果$\pi_i, a_{ij},b_{j}(k)$的值已經收斂，則演算法結束，否則回到第2）步繼續迭代。
以上就是鮑姆-韋爾奇演算法的整個過程。
（歡迎轉載，轉載請註明出處。歡迎溝通交流： liujianping-ok@163.com）