隱馬爾可夫模型HMM推導

機器學習-白板推導系列(十四)-隱馬爾可夫模型HMM（Hidden Markov Model） 課程筆記

背景介紹

介紹一下頻率派和貝葉斯派兩大流派發展出的建模方式。

頻率派

頻率派逐漸發展成了統計機器學習，該流派通常將任務建模為一個優化問題，其中有三個關鍵：

• 模型（model），建立用于解決問題的模型，如 $f(w)=w^T+b$ ；
• 策略（strategy），指導模型優化的策略，比如計算損失函數
• 算法（algorithm），根據策略優化模型的算法，比如 GD、SGD、牛頓法、擬牛頓法

著作有李航老師的統計學習方法。

貝葉斯派

貝葉斯派逐漸發展出了概率圖模型，最終是進行推斷（inference），也就是求后驗概率 $p(z|x)$ 及其期望、方差，最終是一個積分問題。這可以通過一系列數值積分方法（如 MCMC）來計算。

概率圖模型

• 有向 —> 貝葉斯網絡
• 無向 —> 馬爾克夫隨機場（馬爾可夫網絡）
• 概率圖模型 + 時間序列 —> 動態模型（Dynamic Model）
  • 隱馬兒可夫模型（HMM），隱狀態是離散的；
  • 卡爾曼濾波（Kalman Filter），隱狀態連續、線性；
  • 粒子濾波（Particle Filter），隱狀態連續、非線性。

動態模型：在 GMM 等模型中，假設每個樣本是獨立同分布的，即不同樣本之間沒有任何關系。而在動態模型中，在此基礎上有了時間序列。這里所謂的時間序列是廣義的，可以是真正的時間序列，如語音信號，也可以是其他有序序列，如自然語言中的一個句子。在動態模型中，樣本點之間是不滿足獨立同分布假設的。動態模型中，通常有隱藏狀態和觀測變量兩種隨機變量。如 HMM 示意圖中，從橫向來看，是時間序列（time），從縱向來看，是混合變量（mixture）。

HMM模型表示

HMM 模型的概率圖表示如下所示，其中空白節點表示隱藏狀態，用 $I=i_1,i_2,\dots ,i_t,\dots$ 來表示，其取值的可能為 $Q=q_1,q_2,\dots ,q_N$，陰影節點表示觀測變量，用 $O=o_1,o_2,\dots ,o_t,\dots$ 來表示，其取值的可能為 $V=v_1,v_2,\dots ,v_M$。概率圖模型的參數為 $\lambda =(\pi ,A,B)$ ，其中

• $\pi$ 是初始概率分布，即第一個隱狀態 $i_1$ 的概率分布，作為隱狀態，它的其取值可能有 $N$ 種
• $A=[a_{ij}]$ 是狀態轉移矩陣（transition matrix），矩陣中的元素 $a_{ij}$ 表示從隱狀態 $i$ 到隱狀態 $j$ 的轉移概率，即 $a_{ij}=P(i_{t+1}=q_j|i_t=q_i)$
• $B=[b_j(k)]$ 是發射矩陣（emission matrix），矩陣中的元素 $b_j(k)$ 表示同一時間點從隱狀態 $j$ 到觀測變量 $k$ 的發射概率，即 $b_j(k)=P(o_t=v_k|i_t=q_j)$

兩個假設

HMM 模型中有兩個假設，分別是齊次馬爾可夫假設和觀測獨立假設。

齊次馬爾可夫假設，即無后效性，未來狀態與過去狀態無關，只與當前狀態有關。即：
$$P(i_{t+1}|i_t,i_{t?1},\dots ,i_1,o_t,o_{t?1},\dots ,o_1)=P(i_{t+1}|i_t)$$
觀測獨立假設，某時刻的觀測變量只與其自身對應的隱狀態有關，與之前的其他變量都無關。即：
$$P(o_t|i_t,i_{t?1},\dots ,i_1,o_{t?1},\dots ,o_1)=P(o_t|i_t)$$

三個問題

HMM 模型中，通常有三種問題：evaluation、learning、decoding。

• evaluation：求值（觀測變量）

  已知 λ={π,A,B}\lambda=\{\pi, A, B\}λ={π,A,B} ，求觀測變量序列出現的概率，即求 $P(o|\lambda )$ 。通常用前向后向算法（forward-backward algorithm）求解。

• learning：求參數，參數估計

  已知觀測序列，求 $\lambda$ ，通常用 EM 算法（baum welch algorithm）求解極大似然估計 $\hat{\lambda }=\arg {\max }_{\lambda }P(O|\lambda )$ 。

• decoding：求狀態序列

  已知觀測序列和模型參數，求使得該觀測序列出現概率最大的隱狀態序列，通常用維特比算法（viterbi algorithm）估計 $\hat{I}=\arg {\max }_{i}P(O|I)$ 。

  decoding 問題又可引申出兩種問題：prediction 和 filtering。

  • prediction：已知過去和現在的觀測序列，求下一時刻的隱狀態，即求 $P(i_{t+1}|o_1,o_2,\dots ,o_t)$
  • filtering：已知過去和現在的觀測序列，求當前時刻的隱狀態，即求 $P(i_t|o_1,o_2,\dots ,o_t)$ ，HMM 一般不做這個問題。

evaluation問題

樸素解法

給定 HMM 模型的參數，求觀測變量，即求 $P(O|\lambda )$ 。
$$P(O|\lambda )=\sum _{S}P(S,O|\lambda )=\sum _{S}P(O|I,\lambda )?P(I|\lambda )$$
其中，后一項：
$$\begin{array}{rl}P(S|\lambda )& =P(i_1,i_2,\dots ,i_T|\lambda )\\ & =P(i_T|i_1,i_2,\dots ,i_{T?1},\lambda )?P(i_1,i_2,\dots ,i_{T?1},\lambda )\\ & =P(i_T|i_{T?1})?P(i_1,i_2,\dots ,i_{T?1},\lambda )齊次馬爾可夫假設\\ & =P(i_T|i_{T?1})P(i_{T?1}|i_{T?2})\dots P(i_2|i_1)遞歸展開\\ & =a_{i_{T?1},i_T}{a}_{i_{T?2},i_{T?1}}\dots a_{12}\\ & =\prod _{t=2}^T{a}_{i_{t?1}{i}_t}\end{array}\text{\hspace{1em}?}$$
同理，根據觀測獨立假設，前一項：
$$P(O|I,\lambda )=\prod _{t=1}^T{b}_{i_t}(o_t)$$
從而：
$$\begin{array}{rl}P(O|\lambda )& =\sum _I\prod _{t=2}^T{a}_{i_{t?1}{i}_t}\prod _{t=1}^T{b}_{i_t}(o_t)\\ & =\sum _{i_1}\sum _{i_2}?\sum _{i_T}\prod _{t=2}^T{a}_{i_{t?1}{i}_t}\prod _{t=1}^T{b}_{i_t}(o_t)\end{array}\text{\hspace{1em}?}$$
上式的時間復雜度為 $\mathcal{O}(T{N}^T)$ ，無疑是不能接受的。我們需要尋找時間復雜度更低的算法。

前向算法

推導前向后向算法的過程中用到的變換主要就是 HMM 的兩個假設，以及條件概率的公式。

記 ${\alpha }_t(i)=P(o_1,\dots ,o_t,i_t=q_i|\lambda )$ ，則
$$P(O|\lambda )=\sum _{i=1}^{N}P(O,i_t=q_i|\lambda )=\sum _{i=1}^N{\alpha }_T(i)$$
要求觀測變量 $P(O|\lambda )$， 現在我們只要求出 ${\alpha }_T(i)$ ，然后再 $i=1\dots N$ 上求和即可。

一個自然的想法就是找 ${\alpha }_t(i)$ 的遞推式，即其與 ${\alpha }_{t+1}(i)$ 的關系。有：
$$\begin{array}{rl}{\alpha }_{t+1}(j)& =P(o_1,\dots ,o_t,o_{t+1},i_{t+1}=q_j|\lambda )\\ & =\sum _{i=1}^{N}P(o_1,\dots ,o_t,o_{t+1},i_t=q_i,i_{t+1}=q_j|\lambda )\\ & =\sum _{i=1}^{N}P(o_{t+1}|o_1,\dots ,o_t,i_t=q_i,i_{t+1}=q_j,\lambda )?P(o_1,\dots ,o_t,i_t=q_i,i_{t+1}=q_j|\lambda )\\ & =\sum _{i=1}^{N}P(o_{i+1}|i_{t+1}=q_j)?P(o_1,\dots ,o_t,i_t=q_i,i_{t+1}=q_j|\lambda )觀測獨立假設\\ & =\sum _{i=1}^{N}P(o_{i+1}|i_{t+1}=q_j)?P(i_{t+1}=q_j|o_1,\dots ,o_t,i_t=q_i,\lambda )?P(o_1,\dots ,o_t,i_t=q_i|\lambda )\\ & =\sum _{i=1}^{N}P(o_{i+1}|i_{t+1}=q_j)?P(i_{t+1}=q_j|i_t=q_i,\lambda )?P(o_1,\dots ,o_t,i_t=q_i|\lambda )齊次馬爾可夫假設\\ & =\sum _{i=1}^{N}P(o_{i+1}|i_{t+1}=q_j)?P(i_{t+1}=q_j|i_t=q_i,\lambda )?{\alpha }_t(i)\\ & =\sum _{i=1}^N{b}_j(o_{t+1})?a_{ij}?{\alpha }_t(i)\end{array}\text{\hspace{1em}?}$$
由此，我們就得到了遞推式：
$${\alpha }_{t+1}(j)=\sum _{i=1}^N{b}_j(o_{t+1})?a_{ij}?{\alpha }_t(i)$$
從而可以前向（從前往后，從 $1$ 到 $T$）求得 ${\alpha }_T(i)$ ，從而求得 $P(O|\lambda )={\sum }_{i=1}^N{\alpha }_T(i)$ 。

后向算法

與前向算法類似，我們定義一個記號 ${\beta }_t(i)$ ：
$${\beta }_t(i)=P(o_t,\dots ,o_T|i_t=q_i,\lambda )$$
則：
$$\begin{array}{rl}P(O|\lambda )& =P(o_1,\dots ,o_T|\lambda )\\ & =\sum _{i=1}^{N}P(o_1,\dots ,o_T,i_1=q_i|\lambda )\\ & =\sum _{i=1}^{N}P(o_1,\dots ,o_T|i_1=q_i,\lambda )?P(i_1=q_i)\\ & =\sum _{i=1}^{N}P(o_1,\dots ,o_T|i_1=q_i,\lambda )?{\pi }_i\\ & =\sum _{i=1}^{N}P(o_1|o_2,\dots ,o_T,i_1=q_i,\lambda )?P(o_2,\dots ,o_T|i_1=q_i,\lambda )?{\pi }_i\\ & =\sum _{i=1}^{N}P(o_1|i_1=q_i,\lambda )?{\beta }_1(i)?{\pi }_i\\ & =\sum _{i=1}^N{b}_i(o_1)?{\pi }_i?{\beta }_1(i)\end{array}\text{\hspace{1em}?}$$
$\lambda$ 為已知，下面就省略不寫了。同樣我們試圖得到遞推式，即 ${\beta }_{t+1}(j)$ 與 ${\beta }_t(i)$ 之間的關系。有：
$$\begin{array}{rl}{\beta }_t(i)& =P(o_{t+1},\dots ,o_T|i_t=q_i)\\ & =\sum _{j=1}^{N}P(o_{t+1},\dots ,o_T,i_{t+1}=q_j|i_t=q_i)\\ & =\sum _{j=1}^{N}P(o_{t+1},\dots ,o_T|i_{t+1}=q_j,i_t=q_i)?P(i_{t+1}=q_j|i_t=q_i)\\ & =\sum _{j=1}^{N}P(o_{t+1},\dots ,o_T|i_{t+1}=q_j)?a_{ij}前面拿掉i_t=q_i的證明暫略\\ & =\sum _{j=1}^{N}P(o_{t+1}|o_{t+2},\dots ,o_T,i_{t+1}=q_j)?P(o_{t+2},\dots ,o_T|i_{t+1}=q_j)?a_{ij}\\ & =\sum _{j=1}^{N}P(o_{t+1}|i_{t+1}=q_j)?{\beta }_{t+1}(j)?a_{ij}獨立觀測假設\\ & =\sum _{j=1}^N{b}_j(o_{t+1})?a_{ij}?{\beta }_{t+1}(j)\end{array}\text{\hspace{1em}?}$$
由此，就得到了遞推式：
$${\beta }_t(i)=\sum _{j=1}^N{b}_j(o_{t+1})?a_{ij}?{\beta }_{t+1}(j)$$
從而，可以后向（從后往前，從 $T$ 到 $1$）求得 ${\beta }_1(i)$ ，從而求得 $P(O|\lambda )={\sum }_{i=1}^N{b}_i(o_1)?{\pi }_i?{\beta }_1(i)$ 。

總結

	遞推方向	設定	求解公式	復雜度
樸素解法	-	-	$P(O	\lambda)=\sum_{i_1}\sum_{i_2}\dots\sum_{i_T}\prod_{t=2}Ta_{i_{t-1}i_t}\prod_{t=1}Tb_{i_t}(o_t)$
前向算法	從前向后，從$1$到$T$	$\alpha_t(i)=P(o_1,\dots,o_t,i_t=q_i	\lambda)$	$P(O
后向算法	從后向前，從$T$到$1$	$\beta_t(i)=P(o_t,\dots,o_T	i_t=q_i,\lambda)$	$P(O

learning問題

HMM 的 learning 問題是要估計出其參數 $\lambda =(\pi ,A,B)$ 。如果直接用極大似然估計，即：
$${\lambda }_{MLE}=\arg \underset{\lambda }{\max }P(O|\lambda )$$
是無法直接解析的。因此，我們這里用迭代的 EM 算法來估計參數。EM 算法的迭代公式如下：
$${\theta }^{(t+1)}=\arg \underset{\theta }{\max }{\int }_Z\log P(X,Z|\theta )?P(Z|X,{\theta }^{(t)})dZ$$
其中 $x$ 是觀測變量，$z$ 是隱變量，$\theta$ 是參數。之前提到，HMM 除了可以橫向看做是時間（time）序列之外，還可以縱向看做是混合（mixture）變量。將觀測變量、隱變量和參數分別對應到 HMM 中符號之后，順便做一步簡化，把常數拿掉，有：
$$\begin{array}{rl}{\lambda }^{(t+1)}& =\arg \underset{\lambda }{\max }\sum _I\log P(O,I|\lambda )?P(I|O,{\lambda }^{(t)})\\ & =\arg \underset{\lambda }{\max }\sum _I\log P(O,I|\lambda )?\frac{P(O,I|{\lambda }^{(t)})}{P(O|{\lambda }^{(t)})}\\ & =\arg \underset{\lambda }{\max }\sum _I\log P(O,I|\lambda )?P(O,I|{\lambda }^{(t)})\end{array}\text{\hspace{1em}?}$$
記為 Q 函數 $Q(\lambda ,{\lambda }^{(t)})$ ，記 ${\lambda }^{(t)}=({\pi }^{(t)},A^{(t)},B^{(t)})$ ，又有 $P(O|\lambda )={\sum }_{I}P(O,I|\lambda )={\sum }_{i_1}?{\sum }_{i_T}{\pi }_{i_1}{\prod }_{t=2}^T{a}_{i_{t?1},i_t}{\prod }_{t=1}^T{b}_{i_t}(o_t)$，則：
$$\begin{array}{rl}Q(\lambda ,{\lambda }^{(t+1)})& =\sum _I\log P(O,I|\lambda )?P(O,I|{\lambda }^{(t)})\\ & =\sum _I[(\log {\pi }_{i_1}+\sum _{t=2}^T\log a_{i_{t?1},i_t}+\sum _{t=1}^T\log b_{i_t}(o_t))?P(O,I|{\lambda }^{(t)})]\end{array}\text{\hspace{1em}?}$$
這里就只介紹 $\pi$ 的參數估計，$A,B$ 類似。觀察上面 Q 函數的展開結果，只有第一項與 $\pi$ 相關，有：
$$\begin{array}{rl}{\pi }^{(t+1)}& =\arg \underset{\pi }{\max }Q(\lambda ,{\lambda }^{(t)})\\ & =\arg \underset{\pi }{\max }\sum _I\log {\pi }_{i_1}?P(O,I|{\lambda }^{(t)})\\ & =\arg \underset{\pi }{\max }\sum _{i_1}?\sum _{i_T}[\log {\pi }_{i_1}?P(O,i_1,\dots ,i_T|{\lambda }^{(t)})]\\ & =\arg \underset{\pi }{\max }\sum _{i=1}^N[\log {\pi }_{i}P(O,i_1=q_i|{\lambda }^{(t)})]\end{array}\text{\hspace{1em}?}$$
由于 $\pi$ 是初始狀態的概率分布，因此上述優化問題有約束：${\sum }_{i=1}^N{\pi }_i=1$ 。解決帶約束的優化問題，這里用拉格朗日乘子法。

令：
$$\mathcal{L}(\pi ,\eta )=\sum _{i=1}^N\log {\pi }_{i}P(O,i_1=q_i|{\lambda }^{(t)})+\eta (\sum _{i=1}^N{\pi }_i?1)$$
求其對 ${\pi }_i$ 的偏導，并令其為 0：
$$\begin{gathered} \frac{?\mathcal{L}}{?{\pi }_i}=\frac{1}{{\pi }_i}P(O,i_1=q_i|{\lambda }^{(t)})+\eta =0 \\ \sum _{i=1}^N[P(O,i_1=q_i|{\lambda }^{(t)})+{\pi }_i\eta ]=0 \\ P(O|{\lambda }^{(t)})+\eta =0 \\ \eta =?P(O|{\lambda }^{(t)}) \end{gathered}$$
將其帶回，有：
$${\pi }_i^{(t+1)}=?\frac{1}{\eta }P(O|i_1=q_i|{\lambda }^{(t)})=\frac{P(O,i_1=q_i|{\lambda }^{(t)})}{P(O|{\lambda }^{(t)})}$$
至此，我們就得到了由 ${\pi }_i^{(t)}$ 迭代計算 ${\pi }_i^{(t+1)}$ 的公式，從而，就可以估計 ${\pi }^{(t+1)}$：
$${\pi }^{(t+1)}=({\pi }_1^{(t+1)},\dots ,{\pi }_N^{(t+1)})$$
這就完成了 HMM 參數中 $\pi$ 的 learning 問題，對于另外兩個參數 $A,B$ ，方法類似。這就是 baum-welch 算法，可以看做是一種特殊形式的 EM 算法。

decoding問題

decoding 問題是已知觀測序列和模型參數的情況下，求使得該觀測序列出現概率最大的隱狀態序列，通常用維特比算法（viterbi algorithm）估計 $\hat{I}=\arg {\max }_{i}P(O|I)$ 。

在之前的題目設定中有：隱狀態 $i_t$ 取值的集合為 $Q=q_1,q_2,\dots ,q_N$ ，每個時間點的隱狀態 $i_t$ 都有可能取這些值，我們就是要每一步 $i_t$ 選取一個 $q_i$ ，最終找到一組 $I=i_1,\dots ,i_T$ 的取值（圖中藍色虛線），使得觀測變量序列 $O=o_1,\dots ,o_T$ 出現的概率最大。

實際中，我們使用動態規劃的思想來解決該問題。首先定義：${\delta }_t(i)={\max }_{i_1,i_2,\dots ,i_{t?1}}P(o_1,\dots ,o_t.i_1,\dots ,i_{t?1},i_t=q_i)$ 表示到達第 $i_t$ 個位置取值為 $i_t=q_i$ 的最大概率值。接下來我們要找到動態規劃中的狀態轉移方程，即 ${\delta }_{t+1}(j)$ 與 ${\delta }_t(i)$ 的關系。其實就是乘上從 $i$ 到 $j$ 對應的轉移概率 $a_{ij}$ 和 $t+1$ 時刻的發射概率 $b_j(o_{t+1})$ 和 ${\delta }_t(i)$ 乘積的最大者，即有：
$$\begin{array}{rl}{\delta }_{t+1}(j)& =\underset{i_1,\dots ,i_t}{\max }P(o_1,\dots ,o_{t+1},i_1,\dots ,i_t,i_{t+1}=q_j)\\ & =\underset{1\le i\le N}{\max }{\delta }_t(i)a_{ij}{b}_j(o_{t+1})\end{array}\text{\hspace{1em}?}$$
這樣，我們就求出了每一步的最大概率值，注意現在只是求得了最大概率值，沒有記錄路徑，而我們的目標是要求最大概率值對應的隱狀態序列路徑。

令
$${?}_{t+1}(j)=\arg \underset{1\le i\le N}{\max }{\delta }_t(i)?a_{ij}$$
即可將每一步所選的 $q_i$ 記錄下來。這就是求解 HMM decoding 問題的維特比算法。

總結

HMM 是一種動態模型，即含有時間序列這一維度。更廣義地來說，它是一種狀態空間模型（State Space Model），狀態空間模型可以不包含時間維度，如卡爾曼濾波、粒子濾波等，它們的概率圖都類似（如圖 1 所示）。

記隱變量為 $Z$，觀測變量為 $X$，模型參數為 $\theta$ ，這一類模型的問題可以歸納為下面幾種：

• Learning，根據觀測變量估計模型參數，如由極大似然來估計參數 $\hat{\theta }=\arg {\max }_{\theta }P(X|\theta )$，可由 Baum Welch（EM）算法求解；
• Inference，根據模型參數進行分析、求值、預測等，以下 $\theta$ 均已知
  • Decoding，根據觀測變量解碼出隱變量，即求 $\hat{Z}=\arg {\max }_{Z}P(Z|X)$，可由維特比算法（動態規劃）求解；
  • Evaluation，已知模型參數計算觀測變量，即求 prob of evidence $P(X)$ 。通常如果使用樸素方法直接計算，會有較高的復雜度，因此對該問題的求解也有對應的算法，如前向算法和后向算法；
  • Filtering，根據歷史觀測變量求當前時刻的隱變量，即求 $P(z_t|x_1,x_2,\dots ,x_t)$ 。比起只根據當前時刻觀測變量求隱變量 $P(z_t|x_t)$ ，濾波能夠根據更多歷史信息，濾除噪聲，故稱為 Filtering。常用于在線（online）分析。也可由前向算法求解。
  • Smoothing，根據全部觀測變量，求某點的隱變量，即求 $P(z_t|x_1,x_2,\dots ,x_T)$ ，常用于離線（offline）分析。可由前向后向算法求解。
  • Prediction，根據歷史觀測變量，求下一個或下兩個時刻的隱變量或觀測變量，即求 $P(z_{t+1}|x_1,x_2,\dots ,x_t)$ 或 $P(x_{t+1}|x_1,x_2,\dots ,x_t)$ 。比如 NLP 中根據 ”我愛中國共“ 這 5 個歷史觀測變量，預測出下兩個觀測變量為 ”產黨“。可用前向算法求解。