Hidden Markov Model

Markov Chain

馬爾科夫鏈(Markov chain)是一個具備馬氏性的隨機過程，其時間和狀態參數都是離散的。馬爾科夫鏈可用於描述系統在狀態空間中的各類狀態之間的轉移狀況，其中下一個狀態僅依賴於當前狀態。由於系統是隨機變化的，因此不可能百分百預測出將來某個時刻的系統狀態，可是咱們能夠預測出將來時刻系統處在某個狀態的機率。 下面咱們從實際生活中的天氣預測問題入手解析馬爾科夫鏈。現將天氣的狀態粗分爲三種：1-雨雪天氣、2-多雲、3-天晴。假設明天的天氣狀況僅和今天的天氣有關，根據大量的氣象數據，咱們統計出了今明兩天之間的天氣變化規律，如圖1中的表格所示。天氣狀態變化對應的機率圖模型爲圖1(b)，其中狀態值之間的邊的權值表示轉移機率，且每一個狀態的全部指出去的邊的權值之和爲1。若是今天爲雨天，那麼明天下雨的機率爲0.4；若是今天爲多雲，明天出現多雲天氣的機率爲0.6；若是今天爲晴天，明天最有可能出現的天氣依然是晴天，其機率爲0.8。基於今天的天氣情況和明天可能出現的天氣情況，咱們能夠估算一個星期後的各類天氣情況對應的機率。根據該實例，咱們可知運用馬爾科夫模型預測將來天氣情況須要幾個要素：1)當前的天氣(初始狀態)；2)連續兩天的天氣之間轉變模式(轉移機率)。

咱們進一步用數學語言對馬爾科夫鏈相應的知識點進行抽象和歸納。馬爾科夫鏈是由隨機變量組成的序列，每一個隨機變量可能的值組成該馬爾科夫鏈的狀態空間\(\mathcal{S}\)，且該序列知足馬氏性，馬氏性的數學表述以下： \begin{equation} P(X_{t+1}=x|X_1=x_1,X_2=x_2,\cdots,X_t=x_t)=P(X_{t+1}=x|X_t=x_t) \end{equation} 其中條件機率\(a_{ij}(t)=P(X_{t+1}=j|X_t=i)\)的含義爲隨機遊動的質點在時刻\(t\)處於狀態\(i\)的前提下，下一步轉移到狀態\(j\)的機率。當\(a_{ij}(t)\)與時刻\(n\)無關時，馬爾科夫鏈具備平穩轉移機率，下面咱們只討論這種狀況。設\(A\)爲轉移機率組成的矩陣，且狀態空間\(\mathcal{S}=\{1,2,\cdots,N\}\)，則狀態轉移矩陣\(A\)的形式以下： \begin{equation} A=\left[\begin{array}{cccc} a_{11}& a_{12} &\cdots &a_{1N}\\ a_{21}& a_{22} &\cdots &a_{2N}\\ \vdots & \vdots &\ddots&\vdots\\ a_{N1}& a_{N2} &\cdots &a_{NN} \end{array} \right] \end{equation} 其中轉移機率具備以下性質： \begin{equation} \begin{cases} a_{i,j}\geq 0, &i,j\in\mathcal{S}\\ \sum_{j\in\mathcal{S}}a_{ij}=1, &i,j\in\mathcal{S} \end{cases} \end{equation} 由定義可知，從時刻0到\(n\)的馬爾科夫鏈狀態序列對應的機率爲： \begin{equation} \begin{array}{ll} &P(X_0=x_0,X_1=x_1,\cdots,X_t=x_t)\\ =&P(X_t=x_t|X_0=x_0,X_1=x_1,\cdots,X_{t-1}=x_{t-1})\\ & \cdot P(X_0=x_0,X_1=x_1,\cdots,X_{t-1}=x_{t-1})\\ =&P(X_t=x_t|X_{t-1}=x_{t-1})P(X_0=x_0,X_1=x_1,\cdots,X_{t-1}=x_{t-1})\\ =& \cdots\\ =&P(X_0=x_0)\prod_{i=0}^{t-1}P(X_{i+1}=x_{i+1}|X_i=x_i)\\ =&\pi_{x_0}\prod_{i=0}^{t-1}a_{x_ix_{i+1}} \end{array} \end{equation} 其中\(\pi_i\)表示初始狀態爲\(i\)的機率。若是今天是晴天，那麼後面一個星期的天氣情況觀察序列\(O\)是"sun-sun-rain-rain-sun-cloudy-sun"的機率是多大呢？在前面給出的天氣變化模型中，天氣狀態有三種\(s1\)(雨雪)、\(s2\)(多雲)和\(s3\)(晴)，根據這個公式，咱們能夠計算出對應的機率值： \begin{equation} \begin{array}{rl} P(O|model)=&P(s3,s3,s3,s1,s1,s3,s2,s3|model)\\ =&P(s3)P(s3|s3)P(s3|s3)P(s1|s3)\\ &\cdot P(s1|s1)P(s3|s1)P(s2|s3)P(s3|s2)\\ =&P(s3)a_{33}a_{33}a_{13}a_{11}a_{31}a_{23}a_{32}\\ =&1\cdot 0.8\cdot 0.8\cdot 0.1\cdot 0.4\cdot 0.3\cdot 0.1\cdot 0.2\\ =&1.536\times 10^{-4} \end{array} \end{equation} 如圖2所示，如何求解隨機遊動的質點在時刻\(t\)處於狀態\(s\in\mathcal{S}\)的機率\(P(X_t=s)\)呢？

1. 一個簡單可是效率很低的解法以下： \begin{equation} P(X_t=s)=\sum_{i=0}^{t-1}\sum_{x_i=1}^NP(X_0=x_0,X_1=x_1,\cdots,X_t=s) \end{equation} 結合馬爾科夫鏈狀態序列機率的求解規則，有 \begin{equation} P(X_t=s)=\sum_{i=0}^{t-1}\sum_{x_i=1}^NP(X_0=x_0)\prod_{i=0}^{t-1}P(X_{i+1}=x_{i+1}|X_i=x_i) \end{equation} 因爲前\(t\)個隨機變量\(X_i\)各有\(N\)種可能的取值，計算每一個狀態序列的時間複雜度爲\(O(t)\)，所以上式總的時間複雜度爲\(O(tN^t)\)。
2. 第一種解法複雜度過高，由於針對全部的公共子路徑都會重複計算，若是在圖2中的每一個結點都存儲機率值\(P(X_i=x_i)\)，就能夠消除重複計算，算法的複雜度也會下降。基於這個想法，咱們給出以下的計算規則： \begin{equation} \begin{array}{rl} P(X_i=x_i)&=\sum_{x_{i-1}=1}^NP(X_{i-1}=x_{i-1},X_i=x_i)\\ &=\sum_{x_{i-1}=1}^NP(X_{i-1}=x_{i-1})P(X_i=x_i|X_{i-1}=x_{i-1})\\ &=\sum_{x_{i-1}=1}^NP(X_{i-1}=x_{i-1})a_{x_{i-1}x_{i}} \end{array} \end{equation} 由上式可知，計算每一個結點的機率\(P(X_i=x_i)\)的時間複雜度爲\(O(N)\)，圖2有\(N(t-1)+1\)個這樣的結點須要計算，所以總的時間複雜度爲\(O(tN^2)\)。這種方法不管是用遞歸仍是非遞歸都很容易實現，只不過用非遞歸方式實現時，須要按時間遞增的順序計算\(P(X_i=x_i)\)且須要存儲每一個時刻的\(N\)種可能的狀態值對應的機率。

上述問題涉及到馬爾科夫鏈的\(n\)步轉移機率，表示隨機移動的質點在兩個先後相差\(n\)個時刻的狀態值分別爲\(i\)和\(j\)的機率，數學描述形式以下： \begin{equation} a_{ij}^{(n)}=P(X_{m+n}=j|X_m=i), i,j\in\mathcal{S},m\geq 0,n\geq 1 \end{equation} 則\(n\)步狀態轉移矩陣\(A^{(n)}\)的形式以下： \begin{equation} A^{(n)}=\left[\begin{array}{cccc} a_{11}^{(n)}& a_{12}^{(n)} &\cdots &a_{1N}^{(n)}\\ a_{21}^{(n)}& a_{22}^{(n)} &\cdots &a_{2N}^{(n)}\\ \vdots & \vdots &\ddots&\vdots\\ a_{N1}^{(n)}& a_{N2}^{(n)} &\cdots &a_{NN}^{(n)} \end{array} \right] \end{equation} 其中\(n\)步轉移機率具備以下性質： \begin{equation} \begin{cases} a_{i,j}^{(n)}\geq 0, &i,j\in\mathcal{S}\\ \sum_{j\in\mathcal{S}}a_{ij}^{(n)}=1, &i,j\in\mathcal{S}\\ a_{i,j}^{(n)}=\sum_{k=1}^Na_{i,k}^{(l)}a_{k,j}^{(n-l)}\\ A^{(n)}=AA^{(n-1)} \end{cases} \end{equation} 根據全機率公式及馬氏性，可證實第三條性質： \begin{equation} \begin{array}{rl} a_{i,j}^{(n)}&=P(X_{m+n}=j|X_m=i)=\frac{P(X_m=i,X_{m+n}=j)}{P(X_m=i)}\\ &=\sum_{k=1}^N\frac{P(X_m=i,X_{m+l}=k,X_{m+n}=j)}{P(X_m=i,X_{m+l}=k)}\cdot \frac{P(X_m=i,X_{m+l}=k)}{P(X_m=i)}\\ &=\sum_{k=1}^NP(X_{m+n}=j|X_{m+l}=k)P(X_{m+l}=k|X_m=i)\\ &=\sum_{k=1}^Na_{i,k}^{(l)}a_{k,j}^{(n-l)} \end{array} \end{equation} 求出了\(n\)步轉移機率矩陣後，求某個狀態在特定時刻出現的機率\(P(X_t=s)\)也就迎刃而解了 \begin{equation} P(X_t=s)=\sum_{j=1}^NP(X_0=j)a_{js}^{(n)} \end{equation} 計算\(A^{(n)}\)須要完成\(n-1\)個\(N\times N\)的矩陣相乘，計算\(P(X_t=s)\)須要兩個\(N\)維向量相乘，所以最終計算\(P(X_t=s)\)的時間爲複雜度\(O(tN^2)\)，空間複雜度爲\(O(N^2)\)。 假設現有\(L\)條馬爾可夫序列\(\{X^1,X^2,\cdots X^L\}\)，如何基於這些訓練數據學習一個最合理的Markov模型呢？咱們採用最大似然估計的方法估計Markov模型的最優參數。假設訓練數據之間相互獨立，則整個訓練集上的似然函數的對數形式以下： \begin{equation} \begin{array}{rl} &L(A,\pi)\\ =&\log P(X^1,X^2,\cdots X^L)\\ =&\sum_{s=1}^L\log P(X^s)\\ =&\sum_{s=1}^L\log\pi_{x_0^s}+\sum_{t=1}^{T^s}\log a_{x^s_{t-1}x_t^s}\\ =&\sum_{s=1}^L\sum_{i=1}^N1\{x^s_0=i\}\log \pi_i\\ \quad &+\sum_{s=1}^L\sum_{i=1}^N\sum_{j=1}^N\sum_{t=1}^{T^s}1\{x^s_{t-1}=i,x_t^s=j\}\log a_{ij} \end{array} \end{equation} 其中各符號的上標\(s\)表示該符號與第\(s\)條訓練數據\(X^s\)有關。又轉移機率\(a_{ij}\)和初始機率必須知足約束條件 \begin{equation} \begin{cases} \sum_{i=1}^N\pi_i=1\\ \sum_{j=1}^Na_{ij}=1 \end{cases} \end{equation} 咱們引入Lagrange乘子法求得使上述似然函數最大且知足約束條件的參數，構造的Lagrange函數以下： \begin{equation} \begin{array}{rl} \mathcal{L}(A,\pi,\alpha,\beta)=L(A,\pi)+\beta(1-\sum_{i=1}^N\pi_i)+\sum_{i=1}^N\alpha_i(1-\sum_{j=1}^Na_{ij}) \end{array} \end{equation} 用Lagrange函數分別對參數\(a_{ij}\)和\(\alpha_i\)求偏導 \begin{equation} \frac{\partial\mathcal{L}(A,\pi,\alpha,\beta)}{\partial a_{ij}}=\frac{1}{a_{ij}}\sum_{s=1}^L\sum_{t=1}^{T^s}1\{x^s_{t-1}=i,x_t^s=j\}-\alpha_i=0 \end{equation} \begin{equation} \frac{\partial\mathcal{L}(A,\pi,\alpha,\beta)}{\partial \alpha_{i}}=1-\sum_{j=1}^Na_{ij}=0 \end{equation} 結合上面兩個等式，可得 \begin{equation} a_{ij}=\frac{\sum_{s=1}^L\sum_{t=1}^{T^s}1\{x^s_{t-1}=i,x_t^s=j\}}{\sum_{s=1}^L\sum_{t=1}^{T^s}1\{x^s_{t-1}=i\}} \end{equation} 用Lagrange函數分別對參數\(\pi_{i}\)和\(\beta\)求偏導 \begin{equation} \frac{\partial\mathcal{L}(A,\pi,\alpha,\beta)}{\partial \pi_{i}}=\frac{\sum_{s=1}^L1\{x_0^s=i\}}{\pi_{i}}-\beta=0 \end{equation} \begin{equation} \frac{\partial\mathcal{L}(A,\pi,\alpha,\beta)}{\partial \beta}=1-\sum_{i=1}^N\pi_i=0 \end{equation} 結合上面兩個等式，可得 \begin{equation} \pi_{i}=\frac{\sum_{s=1}^L1\{x^s_{0}=i\}}{L} \end{equation} 根據上面的參數學習規則，很容易看出來這些參數的學習都是基於統計的，\(a_{ij}\)等於全部從狀態\(i\)跳轉到狀態\(j\)的次數與全部狀態爲\(i\)的比值，而\(\pi_i\)則爲\(L\)條訓練數據中以狀態\(i\)開始的比率。

Hidden Markov Model

在馬爾可夫鏈中，狀態是直接可見的。但是，在現實世界中，咱們觀察到的狀況不免被噪聲或錯誤干擾，使得"What you see is the truth"再也不是有效的假設。某種現象的背後必然存在與實際相符的內在因素，也許這些內在因素還未被揭示出來，但由內在因素驅動的外部現象每每能夠做爲研究內在因素的線索，這也就是隱馬爾可夫模型(Hidden Markov Model,HMM)的內在哲理。在HMM中，狀態是不可見的，但依賴於狀態的輸出值是可見的；每一個狀態在全部可能的輸出上都有一個機率分佈，由HMM模型產生的輸出序列爲相應的狀態序列提供了部分信息[1]。 咱們結合圖3所示的瓷缸問題具體描述HMM的思想。假設有\(N\)個瓷缸，\(M\)種不一樣顏色的球隨機分佈在瓷缸中。實驗人員根據某種隨機過程選擇一個初始的瓷缸，從該瓷缸隨機取出一個球並記錄下其顏色，而後將球放回原瓷缸；接着，實驗人員根據與當前選擇的瓷缸相關的隨機選擇過程選擇一個新的瓷缸，記錄球的顏色後放回原瓷缸；不斷重複該過程，最後獲得了一個有限的顏色觀測序列。咱們想要基於這些已知的顏色序列運用HMM建模數學模型，那麼HMM的狀態對應特定的瓷缸，取出的球的顏色則對於該狀態對應的輸出值。HMM在具備時序性質的模式識別中有着普遍的應用，包括語音識別、手寫體識別、手勢識別和基因序列分析等。

HMM能夠被視爲混合模型(mix model)的推廣，由於HMM中的隱含狀態間存在必定的上下文關係；而在Gaussian Mixture Model(GMM)等混合模型中，也是由多個隱含狀態支配全部觀測值，只是其中的隱含狀態彼此間是獨立的。所以，HMM中存在兩個隨機過程，一是隱含狀態間隨機跳轉且具備馬氏性，二是隱含狀態產生的輸出值也具備隨機性。與這兩個隨機過程相關聯的全部參數組成了HMM的參數，如圖4所示，\(x\)表示隱含狀態，各狀態間的轉移機率(transition probabilities)爲\(a\)，每一個狀態值對應的輸出變量\(o\)取值爲\(y\)的輸出機率(output probabilities)爲\(b\)。

HMM的通常性結構如圖5所示，其中每一個橢圓表示隨機變量，\(x(t)\)爲\(t\)時刻的隱含狀態，\(o(t)\)爲表示\(t\)時刻的觀察值的隨機變量，箭頭表示條件依賴性。從圖示可知，\(x(t)\)僅依賴於\(x(t-1)\)，\(o(t)\)僅依賴於\(x(t)\)。在標準的HMM中，狀態空間是離散的，但觀測變量既能夠是離散的變量也能夠是服從高斯分佈等連續型變量。

接下來，咱們對HMM進行形式化描述。HMM的模型參數爲 \begin{equation} \lambda=(A\in\mathbb{R}^{N\times N},B\in\mathbb{R}^{N\times M},\pi\in\mathbb{R}^{1\times N}) \end{equation} 參數的具體說明以下：

1. 隱含狀態的數目\(N\)。儘管HMM中的狀態是不可見的，但在不少實際應用中，狀態是具備必定物理意義的。例如，在瓷缸問題中，每一個狀態表明一個特定的瓷缸。通常而言，狀態間都是相互連通的，但也有例外。狀態集合表示爲\(\mathcal{S}=\{s_1,s_2,\cdots,s_N\}\)，\(t\)時刻的狀態標號用\(x_t\)表示，其中\(x_t\in\{1,2,\cdots,N\}\)。
2. 每一個狀態可能對應的觀測值\(M\)。觀測值對應模型的輸出，在瓷缸問題中就是球的顏色。觀察集合表示爲\(\mathcal{Y}=\{y_1,y_2,\cdots,y_M\}\)，\(t\)時刻的觀測值標號用\(o_t\)表示，其中\(o_t\in\{1,2,\cdots,M\}\)。 \item 狀態轉移機率矩陣\(A=\{a_{ij}\}\)，其中 \begin{equation} a_{ij}=P(x_{t+1}=j|x_t=i)\geq 0, 1\leq i,j\leq N \end{equation}
3. 觀測值關於狀態的輸出機率矩陣\(B=\{b_{j}(k)\}\)，其中 \begin{equation} b_{j}(k)=P(o_t=k|x_t=j),1\leq j\leq N,1\leq k\leq M \end{equation}
4. 初始狀態的機率分佈\(\pi=\{\pi_i\}\)，其中 \begin{equation} \pi_i=P(x_0=i),1\leq i\leq N \end{equation}

Inference in HMM

圍繞着HMM有三個基本問題，這也是HMM在現實應用中發揮做用的關鍵性問題[3]。

1. Evaluation:給定HMM的模型參數\(\lambda\)，如何計算產生某個特定觀測序列的機率。這個問題可理解爲評估觀察序列與模型間的匹配程度，若是咱們須要從能產生這個觀測序列的諸多模型中進行模型選擇時，該問題的解決方案使得咱們能夠選出與觀察序列最匹配的模型。
2. Decoding:給定HMM的模型參數\(\lambda\)和觀測序列，如何計算最有可能產生這個觀測序列的狀態序列。
3. Training:給定足夠的觀測序列，如何估計HMM的模型參數\(\lambda\)(狀態轉移機率和輸出機率)。

下面，咱們就對這三個關鍵問題逐個擊破。

Evaluation

給定模型參數\(\lambda\)和長度爲\(T\)的觀測序列\(O=o_0o_1\cdots o_{T-1}\)，求觀測序列對應的機率值： \begin{equation} P(O|\lambda)=\sum_XP(X,O|\lambda)=\sum_XP(O|X,\lambda)P(X|\lambda) \end{equation} 其中\(X=x_0x_1\cdots x_{T-1}\)爲觀測序列對應的未知的狀態序列，且有 \begin{equation} P(O|X,\lambda)=\prod_{t=0}^{T-1}P(o_t|x_t)=\prod_{t=0}^{T-1}b_{x_t}(o_t) \end{equation} \begin{equation} P(X|\lambda)=P(x_0)\prod_{t=1}^{T-1}P(x_t|x_{t-1})=\pi_{x_0}\prod_{t=1}^{T-1}a_{x_{t-1}x_t} \end{equation} 解決上述問題的最簡單的方法就是窮舉法，時間複雜度爲\(O(TN^T)\)。窮舉法中存在不少重複計算，更爲有效的是採用前向或後向的動態規劃算法，這裏討論前向和後向兩種動態規劃的策略。 假設\(t\)時刻的狀態\(x_t=j\)且已知的觀察序列爲\(o_0o_2\cdots o_t\)的機率爲： \begin{equation} \begin{array}{rl} \alpha_t(j)&=P(o_0o_2\cdots o_t,x_t=j)\\ &=\sum_{i=1}^NP(o_0o_2\cdots o_t,x_{t-1}=i,x_t=j)\\ &=\sum_{i=1}^NP(o_t,x_t=j|o_0\cdots o_{t-1},x_{t-1}=i)P(o_0\cdots o_{t-1},x_{t-1}=i)\\ &=\sum_{i=1}^NP(o_0\cdots o_{t-1},x_{t-1}=i)P(x_t=j|x_{t-1}=i)P(o_t|x_t=j)\\ &=\sum_{i=1}^N\alpha_{t-1}(i)a_{ij}b_j(o_t) \end{array} \end{equation} 由上式可知，只要計算出全部的\(\{\alpha_{t-1}(i)|i=1\cdots N\}\)，就能算出\(\alpha_t(j)\)，因此稱之爲前向算法(forward algorithm)。那麼在初始時刻對應的\(\alpha_0(i)\)爲 \begin{equation} \alpha_0(i)=\pi_ib_i(o_0),1\leq i\leq N \end{equation} 根據前向算法的推理規則不斷計算，最後合併\(T-1\)時刻的全部結果便可獲得觀測序列的機率： \begin{equation} P(O|\lambda)=\sum_{i=1}^N\alpha_{T-1}(i) \end{equation} 對應的算法描述參見1，算法時間複雜讀爲\(O(TN^2)\)，空間複雜度爲\(O(N)\)。

假設在\(t\)時刻的狀態\(x_t=j\)的前提下，\(t+1\)至\(T-1\)的時間段內對應的觀察序列爲\(o_0o_2\cdots o_t\)的機率爲： \begin{equation} \begin{array}{rl} \beta_t(i)&=P(o_{t+1}\cdots o_{T-1}|x_t=i)\\ &=\sum_{j=1}^NP(o_{t+1}\cdots o_{T-1},x_{t+1}=j|x_t=i)\\ &=\sum_{j=1}^NP(x_{t+1}=j|x_t=i)P(o_{t+1}|x_{t+1}=j)P(o_{t+2}\cdots o_{T-1}|x_{t+1}=j)\\ &=\sum_{j=1}^N\beta_{t+1}(j)a_{ij}b_j(o_{t+1}) \end{array} \end{equation} 由上式可知，要要計算出\(\beta_t(i)\)，必須先算出全部的\(\{\beta_{t+1}(j)|j=1\cdots N\}\)，因此稱之爲後向算法(backward algorithm)。在初始時刻對應的\(\beta_{T-1}(j)\)都爲\(1\)。根據後向算法的推理規則不斷計算，最後合併\(0\)時刻的全部結果便可獲得觀測序列的機率： \begin{equation} \begin{array}{rl} P(O|\lambda)&=\sum_{i=1}^NP(o_0\cdots o_{T-1},x_0=i)\\ &=\sum_{i=1}^NP(o_1,\cdots o_{T-1}|x_0=i,o_0)P(o_0|x_0=i)P(x_0=i)\\ &=\sum_{i=1}^NP(o_1,\cdots o_{T-1}|x_0=i)P(o_0|x_0=i)P(x_0=i)\\ &=\sum_{i=1}^N\beta_0(i)b_i(o_0)\pi_i \end{array} \end{equation} 對應的算法描述參加2，算法時間複雜讀爲\(O(TN^2)\)，空間複雜度爲\(O(N)\)。

Decoding

給定模型參數\(\lambda\)和長度爲\(T\)的觀測序列\(O=o_0o_1\cdots o_{T-1}\)，求最優的狀態序列\(X^\star=x_0x_1\cdots x_{T-1}\)： \begin{equation} X^\star=\underset{X}{arg\max}P(O,X|\lambda)=\underset{X}{arg\max}P(O|X,\lambda)P(X|\lambda) \end{equation} 基於動態規劃的Viterbi算法[2]可用於尋找一條最優的隱含狀態序列。假設已知的前\(t\)個觀測值有多條以狀態值\(i\)收尾的狀態序列，這些狀態序列與觀測序列最大的聯合機率爲 \begin{equation} \delta_t(i)=\underset{x_0x_1\cdots x_{t-1}}{\max}P(x_0x_1\cdots x_{t-1},x_t=i,o_0o_1\cdots o_t) \end{equation} 進一步推理，可知前\(t+1\)個觀察值對應的狀態序列以\(j\)收尾的最大機率爲 \begin{equation} \delta_{t+1}(j)=b_j(o_{t+1})\underset{1\leq i\leq N}{\max}\delta_t(i)a_{ij} \end{equation} 根據上式不斷向前計算\(\delta_{t}(j)\)直至最後時刻\(T-1\)，此時可知最優狀態序列\(X^\star\)對應的機率爲 \begin{equation} P(O,X^\star|\lambda)=\underset{1\leq i\leq N}{\max}\delta_{T-1}(i) \end{equation} 那麼在\(T-1\)時刻對應的最優狀態\(x^\star_{T-1}\)爲 \begin{equation} x^\star_{T-1}=\underset{1\leq i\leq N}{arg\max}\delta_{T-1}(i) \end{equation} 假設咱們已經知道最優狀態序列在\(t+1\)至\(T-1\)時時間段內對應局部狀態序列，且\(t+1\)時刻的狀態值爲\(x_{t+1}=j\),那麼\(t\)時刻的狀態\(x_t\)爲什麼值時爲最優呢？因爲\(x_{t+1}\cdots x_{T-1}\)均已肯定，意味着最優狀態序列從\(t+1\)時刻開始的狀態轉移機率和輸出機率對剩下的各狀態變量而言是常數了，又\(P(X,O)\)可經過初始狀態機率及一系列的狀態轉移機率和輸出機率相乘獲得，則 \begin{equation} \begin{array}{rl} &\underset{X}{\max} P(O,X|\lambda)\\ =&\underset{X}{\max} P(O|X,\lambda)P(X|\lambda)\\ =&\underset{X}{\max} \pi_{x_0}b_{x_0}(o_0)\prod_{t'=1}^{T-1}a_{x_{t'-1}x_{t'}}b_{x_{t'}}(o_{t'})\\ =&\underset{X}{\max} \pi_{x_0}b_{x_0}(o_0)(\prod_{t'=1}^{t}a_{x_{t'-1}x_{t'}}b_{x_{t'}}(o_{t'}))(a_{x_{t}x_{t+1}}b_{x_{t+1}}(o_{t+1}))\\ &\quad \cdot(\prod_{t'=t+2}^{T-1}a_{x_{t'-1}x_{t'}}b_{x_{t'}}(o_{t'}))\\ =&\underset{{x_0\cdots x_t}}{\max}P(o_0\cdots o_t,x_0\cdots x_t)a_{x_{t}x_{t+1}}\\ &\quad\cdot \underbrace{b_{x_{t+1}}(o_{t+1})\prod_{t'=t+2}^{T-1}a_{x_{t'-1}x_{t'}}b_{x_{t'}(o_{t'})}}_{constant}\\ =&\underset{1\leq i\leq N}{\max}\delta_t(i)a_{ij}\cdot {constant} \end{array} \end{equation} 由上式可知，當肯定\(t+1\)及之後全部時刻的局部最優狀態序列後，\(t\)時刻的最優狀態值爲 \begin{equation} x_t^\star=\Psi_{t+1}(j)=\underset{1\leq i\leq N}{arg\max}\delta_t(i)a_{ij}=\underset{1\leq i\leq N}{arg\max}\delta_t(i)a_{ij}b_{x_{t+1}}(o_{t+1}) \end{equation} 根據上述分析，在遍歷完整個觀測序列後，\(\Psi\)記錄下的狀態間的鏈接關係如圖6所示。在Viterbi算法計算過程當中，關鍵是要記錄全部時間段內的狀態鏈鏈接關係\(\{\Psi_t(j)|t=1\cdots T-1,j=1\cdots N\}\)便可，而局部路徑對應的最大機率\(\delta_t(i)\)只須要保存前一個時刻和當前時刻的便可。Viterbi算法對應的算法描述見3，時間複雜度爲\(O(TN^2)\)，空間複雜度爲\(O(TN)\)。 

 

Training

給定觀測序列的集合\(\mathcal{O}=\{O^1,O^2,\cdots,O^L|O^s=o_0^s\cdots o^s_{T^s}\}\)，HMM的模型參數\(\lambda=\{A,B,\pi\}\)爲什麼值時才能是觀測序列集合\(\mathcal{O}\)出現的機率最大？接下來，咱們基於前面的前向算法、後向算法，用Lagrange乘子處理參數的約束條件，最後用EM(Expection Maximization)算法的思想完成模型的訓練問題[4]。 根據前面知識，觀察序列和隱含狀態序列的聯合機率計算形式以下： \begin{equation} P(O,X|\lambda)=\pi_{x_0}\left(\prod_{t=0}^{T-1}a_{x_{t}x_{t+1}}\right)\left(\prod_{t=0}^{T}b_{x_t}(o_t)\right) \end{equation} EM算法的基本思想是利用Jasen不等式爲對數似然函數構造一個下界函數，E步驟中在固定模型參數的狀況下經過使隱含狀態序列知足特定的機率分佈使得下界函數與真實目標函數相等，M步驟則在固定隱含狀態序列機率分佈狀況下優化模型參數使得下界函數取得局部最優解；如此不斷重複前面兩個步驟，直至收斂爲止。雖然EM算法只能找到局部最優解，幸運的是大多時候局部最優解已經能夠知足咱們的要求。假設訓練集中的觀測序列相互獨立，則EM算法中的下界函數\(L(A,B,\pi)\)推導以下： \begin{equation} \begin{array}{rl} &\log P(O^1,O^2,\cdots,O^L|\lambda)\\ =&\log\prod_{l=1}^LP(O^l|\lambda)=\sum_{l=1}^L\log P(O^l|\lambda)\\ \geq&\sum_{l=1}^L\sum_{X^l}Q(X^l)\log\frac{P(O^l,X^l)}{Q(X^l)}\\ =&\sum_{l=1}^L\sum_{X^l}Q(X^l)\left[\log P(O^l,X^l)-\log Q(X^l)\right]\\ =&\sum_{l=1}^L\sum_{X^l}Q(X^l)\left[\log\pi_{x^l_0}+\sum_{t=0}^{T^l-1}\log a_{x^l_{t}x^l_{t+1}}\right.\\ \quad&+\left.\sum_{t=0}^{T^l} \log b_{x^l_t}(o^l_t)-\log Q(X^l)\right]\\ =&\sum_{l=1}^L\sum_{X^l}Q(X^l)\left[\sum_{i=1}^N1\{x_0^l=i\}\log\pi_{i}\right.\\ \quad&+\sum_{t=0}^{T^l-1}\sum_{i=1}^N\sum_{j=1}^N1\{x_t^l=i,x_{t+1}^l=j\}\log a_{ij}\\ \quad&+\left.\sum_{t=0}^{T^l} \sum_{i=1}^N\sum_{k=1}^M 1\{x_t^l=i,o_t^l=k\}\log b_{i}(j)-\log Q(X^l)\right]\\ =&L(A,B,\pi) \end{array} \end{equation} 在M步驟中須要求得使\(L(A,B,\pi)\)最大的參數，並且參數必須知足下列約束 \begin{equation} \begin{cases} \sum_{j=1}^Na_{ij}=1\\ \sum_{k=1}^Mb_i(k)=1\\ \sum_{i=1}^N\pi_i=1 \end{cases} \end{equation} 引入Lagrange乘子處理上述約束條件，構造以下的Lagrange函數 \begin{equation} \begin{array}{rl} \mathcal{L}(A,B,\pi)=&L(A,B,\pi)+\sum_{i=1}^N\theta_i(1-\sum_{j=1}^Na_{ij})\\ &+\sum_{i=1}^N\mu_i\left(1-\sum_{k=1}^Mb_i(k)\right)+\eta(1-\sum_{i=1}^N\pi_i) \end{array} \end{equation} 下面用最大似然的參數估計方法計算最優參數。Lagrange函數分別對狀態轉移機率\(a_{ij}\)和\(\theta_i\)求偏導並使其爲0。 \begin{equation} \frac{\partial \mathcal{L}(A,B,\pi)}{\partial a_{ij}}=\frac{1}{a_{ij}}\sum_{l=1}^L\sum_{X^l}Q(X^l)\sum_{t=0}^{T^l-1}1\{x_t^l=i,x_{t+1}^l\}=0 \end{equation} \begin{equation} \frac{\partial \mathcal{L}(A,B,\pi)}{\partial \theta_i}=1-\sum_{j=1}^Na_{ij}=0 \end{equation} 結合上述兩式，可得\(a_{ij}\)的參數求解以下 \begin{equation} a_{ij}=\frac{\sum_{l=1}^L\sum_{X^l}Q(X^l)\sum_{t=0}^{T^l-1}1\{x_t^l=i,x_{t+1}^l=j\}}{\sum_{l=1}^L\sum_{X^l}Q(X^l)\sum_{t=0}^{T^l-1}1\{x_t^l=i\}} \end{equation} 同理，Lagrange函數分別對輸出機率\(b_i(k)\)和\(\mu_i\)求偏導並使其爲0。 \begin{equation} \frac{\partial \mathcal{L}(A,B,\pi)}{\partial b_i(k)}=\frac{1}{b_i(k)}\sum_{l=1}^L\sum_{X^l}Q(X^l)\sum_{t=0}^{T^l}1\{o_t^l=k,x^l_{t}=j\}-\mu_i=0 \end{equation} \begin{equation} \frac{\partial \mathcal{L}(A,B,\pi)}{\partial \mu_i}=1-\sum_{k=1}^Mb_i(k)=0 \end{equation} 結合上述兩式，可得輸出機率\(b_i(k)\)的計算規則爲 \begin{equation} b_i(k)=\frac{\sum_{l=1}^L\sum_{X^l}Q(X^l)\sum_{t=0}^{T^l}1\{o_t^l=k,x_{t}^l=i\}}{\sum_{l=1}^L\sum_{X^l}Q(X^l)\sum_{t=0}^{T^l}1\{x_t^l=i\}} \end{equation} Lagrange函數分別對初始狀態機率\(\pi_i\)和\(\eta\)求偏導並使其爲0。 \begin{equation} \frac{\partial \mathcal{L}(A,B,\pi)}{\pi_i}=\frac{1}{\pi_i}\sum_{l=1}^L\sum_{X^l}Q(X^l)1\{x_0^l=i\}-\eta=0 \end{equation} \begin{equation} \frac{\partial \mathcal{L}(A,B,\pi)}{\eta}=1-\sum_{i=1}^N\pi_i=0 \end{equation} 結合上述兩式，可得初始狀態機率\(\pi_i\)的計算方式爲 \begin{equation} \pi_i=\frac{\sum_{l=1}^L\sum_{X^l}Q(X^l)1\{x_0^l=i\}}{\sum_{l=1}^L\sum_{i=1}^N\sum_{X^l}Q(X^l)1\{x_0^l=i\}} \end{equation} 是否經過上述的參數計算規則，咱們就能獲得了最優參數呢？由於隱含狀態序列\(X^l\)未知，全部以沒法計算\(Q(X^l)\)。咱們先來看兩個基本問題，在已知模型參數\(\lambda\)和觀測序列\(O\)的前提下，求任意單個時刻或連續兩個時刻的隱含狀態爲特定值的機率

• 任意連續兩個時刻的隱含狀態爲特定值的機率 \begin{equation} \begin{array}{rl} &P(O,x_t=i,x_{t+1}=j)\\ =&P(o_0o_1\cdots o_t,x_t=i,o_{t+1}\cdots o_T,x_{t+1}=j)\\ =&P(o_{t+1}\cdots o_T|x_{t+1}=j)P(x_{t+1}=j|x_t=i)\\ \quad&\cdot P(o_{t+1}|x_{t+1}=j)P(o_0o_1\cdots o_t,x_t=i)\\ =&\alpha_t(i)a_{ij}\beta_{t+1}(j)b_j(o_{t+1}) \end{array} \end{equation}
• 任意單個時刻的隱含狀態爲特定值的機率 \begin{equation} P(O,x_t=i)=\sum_{j=1}^NP(O,x_t=i,x_{t+1}=j)=\sum_{j=1}^N\alpha_t(i)a_{ij}\beta_{t+1}(j)b_j(o_{t+1}) \end{equation}

在E步驟中，咱們有\(Q(X)=P(X|O;\lambda)\)，則 \begin{equation} \begin{array}{rl} &\sum_XQ(X)\sum_{t=0}^{T-1}1\{x_t=i,x_{t+1}=j\}\\ =&\frac{1}{P(O)}\sum_X\sum_{t=0}^{T-1}1\{x_t=i,x_{t+1}=j\}P(O,X)\\ =&\frac{1}{P(O)}\sum_{t=0}^{T-1}P(O,x_t=i,x_{t+1}=j)\\ =&\frac{1}{P(O)}\sum_{t=0}^{T-1}\alpha_t(i)a_{ij}\beta_{t+1}(j)b_j(o_{t+1}) \end{array} \end{equation} \begin{equation} \begin{array}{rl} &\sum_XQ(X)\sum_{t=0}^{T-1}1\{x_t=i\}\\ =&\frac{1}{P(O)}\sum_X\sum_{t=0}^{T-1}1\{x_t=i\}P(O,X)\\ =&\frac{1}{P(O)}\sum_{t=0}^{T-1}P(O,x_t=i)\\ =&\frac{1}{P(O)}\sum_{t=0}^{T-1}\sum_{j=1}^NP(O,x_t=i,x_{t+1}=j)\\ =&\frac{1}{P(O)}\sum_{t=0}^{T-1}\sum_{j=1}^N\alpha_t(i)a_{ij}\beta_{t+1}(j)b_{j}(o_{t+1}) \end{array} \end{equation} \begin{equation} \begin{array}{rl} &\sum_XQ(X)\sum_{t=0}^{T}1\{x_t=i,o_t=k\}\\ =&\frac{1}{P(O)}\sum_X\sum_{t=0}^{T}1\{x_t=i,o_t=k\}P(O,X)\\ =&\frac{1}{P(O)}\sum_{t=0}^{T}1\{o_t=k\}P(O,x_t=i)\\ =&\frac{1}{P(O)}\sum_{t=0}^{T}\sum_{j=1}^N1\{o_t=k\}P(O,x_{t}=i,x_{t+1}=j)\\ =&\frac{1}{P(O)}\sum_{t=0}^{T}\sum_{j=1}^N1\{o_t=k\}\alpha_{t}(i)a_{ij}\beta_{t+1}(j)b_{j}(o_{t+1}) \end{array} \end{equation} 將上述三式帶入各參數更新規則中，令\(\gamma_t(i,j)=\alpha_t(i)a_{ij}\beta_{t+1}(j)b_j(o_{t+1})\)，則 \begin{equation} a_{ij}=\frac{\sum_{l=1}^L\sum_{t=0}^{T^l-1}\gamma^l_t(i,j)}{\sum_{l=1}^L\sum_{t=0}^{T^l-1}\sum_{j=1}^N\gamma^l_t(i,j)} \end{equation} \begin{equation} b_{i}(k)=\frac{\sum_{l=1}^L\left(\sum_{t=0}^{T^l-1}\sum_{j=1}^N1\{o_t=k\}\gamma^l_t(i,j)+1\{o_{T^l}=k\}\alpha_{T^l}(i)\right)}{\sum_{l=1}^L\left(\sum_{t=0}^{T^l-1}\sum_{j=1}^N\gamma^l_t(i,j))+\alpha_{T^l}(i)\right)} \end{equation} \begin{equation} \pi_i=\frac{\sum_{l=1}^L\sum_{j=1}^N\gamma^l_0(i,j)}{\sum_{l=1}^L\sum_{i=1}^N\sum_{j=1}^N\gamma^l_0(i,j)} \end{equation}

Scaling

前面的內容已經從理論上解決了HMM相關的問題，可是要實現HMM相關算法卻有不少細節要處理，scaling就是其中之一。在解決前面的任何一個問題時，咱們都要用到衆多範圍在\([0,1]\)之間的機率值的乘積，因此獲得的結果很容易向下溢出。爲了處理這中狀況，在實現HMM算法時有必要引入scaling的技巧，使得計算機處理的中間結果不至於向下溢出，並且最終能夠獲得與原問題相同的解。對於聯合機率\(P(O,X)\)的求解，可採用對數等價轉換 \begin{equation} P(O,X)=\exp\left(\log\pi_{x_0}+\sum_{t=0}^{T-1}\log a_{x_tx_{t+1}}+\sum_{t=0}^T\log b_{x_t}(o_t)\right) \end{equation} Scaling技巧最主要是用在HMM的模型參數估計過程當中，這裏根據文獻[3]中的scaling策略進行簡單分析，最後給出HMM在實際的參數學習過程當中的方法。假設對全部的\(\alpha_t(i)\)和\(\beta_t(i)\)都進行以下放縮 \begin{equation} \begin{array}{cc} &\hat{\alpha}_t(i)=\alpha_t(i)/c_t\\ &\hat{\beta}_t(i)=\beta_t(i)/c_t \end{array} \end{equation} 其中\(c_t=\sum_{i=1}^N\alpha_t(i)\)。前向算法和後向算法中的迭代公式變爲以下形式： \begin{equation} \begin{array}{cc} &\alpha_t(j)=\sum_{i=1}^N\hat{\alpha}_{t-1}(i)a_{ij}b_j(o_t)\\ &\beta_t(i)=\sum_{j=1}^N\hat{\beta}_{t+1}(j)a_{ij}b_j(o_{t+1}) \end{array} \end{equation} 最終很容易獲得 \begin{equation} \begin{array}{cc} &\hat{\alpha}_t(j)=\alpha_{t}(i)/(\prod_{t'=1}^tc_{t'})\\ &\hat{\beta}_{t+1}(j)=\beta_{t+1}(i)/(\prod_{t'=t+1}^Tc_{t'}) \end{array} \end{equation} 因爲各變量放縮後對任意\(t\)都知足\(\sum_{i=1}^N\hat{\alpha}_t(i)=1\)，則有 \begin{equation} \begin{array}{ll} &\sum_{i=1}^N\hat{\alpha}_T(i)=\sum_{i=1}^N\alpha_T(i)/\left(\prod_{t=1}^Tc_t\right)=P(O)/\left(\prod_{t=1}^Tc_t\right)=1\\ &\Longrightarrow \prod_{t=1}^Tc_t=P(O) \end{array} \end{equation} 通過scaling後，參賽更新中的\(\gamma_t(i,j)\)也相應獲得了放縮 \begin{equation} \begin{array}{ll} \hat{\gamma}_t(i,j)&=\hat{\alpha}_t(i)a_{ij}\hat{\beta}_{t+1}(j)b_j(o_{t+1})\\ &=\alpha_t(i)a_{ij}\beta_{t+1}(j)b_j(o_{t+1})/\left(\prod_{t=1}^Tc_t\right)\\ &=\gamma_t(i,j)/\left(\prod_{t=1}^Tc_t\right)=\gamma_t(i,j)/P(O) \end{array} \end{equation} 若是用\(\hat{\gamma}_t(i,j)\)去學習HMM的最優參數，並且要求scaling先後是徹底一致的，則有 \begin{equation} a_{ij}=\frac{\sum_{l=1}^LP(O^l)\sum_{t=0}^{T^l-1}\hat{\gamma}^l_t(i,j)}{\sum_{l=1}^LP(O^l)\sum_{t=0}^{T^l-1}\sum_{j=1}^N\hat{\gamma}^l_t(i,j)} \end{equation} \begin{equation} b_{i}(k)=\frac{\sum_{l=1}^LP(O^l)\left(\sum_{t=0}^{T^l-1}\sum_{j=1}^N1\{o_t=k\}\hat{\gamma}^l_t(i,j)+1\{o_{T^l}=k\}\alpha_{T^l}(i)\right)}{\sum_{l=1}^LP(O^l)\left(\sum_{t=0}^{T^l-1}\sum_{j=1}^N\hat{\gamma}^l_t(i,j))+\alpha_{T^l}(i)\right)} \end{equation} \begin{equation} \pi_i=\frac{\sum_{l=1}^LP(O^l)\sum_{j=1}^N\hat{\gamma}^l_0(i,j)}{\sum_{l=1}^LP(O^l)\sum_{i=1}^N\sum_{j=1}^N\hat{\gamma}^l_0(i,j)} \end{equation}

References

[1] Hidden markov model. http://en.wikipedia.org/wiki/Hidden_Markov_model.

[2] Viterbi algorithm. http://en.wikipedia.org/wiki/Viterbi_algorithm.

[3] Lawrence Rabiner. A tutorial on hidden markov models and selected applications in speech recognition. Proceedings of the IEEE, 77(2):257–286, 1989.

[4] Daniel Ramage. Hidden markov models fundamentals. Lecture Notes. http://cs229. stanford. edu/section/cs229-hmm. pdf, 2007.