Markov Chain

馬爾科夫鏈(Markov chain)是一個具有馬氏性的隨機過程，其時間和狀態(tài)參數(shù)都是離散的。馬爾科夫鏈可用于描述系統(tǒng)在狀態(tài)空間中的各種狀態(tài)之間的轉移情況，其中下一個狀態(tài)僅依賴于當前狀態(tài)。因為系統(tǒng)是隨機變化的，所以不可能百分百預測出未來某個時刻的系統(tǒng)狀態(tài)，但是我們可以預測出未來時刻系統(tǒng)處在某個狀態(tài)的概率。 下面我們從實際生活中的天氣預測問題入手解析馬爾科夫鏈。現(xiàn)將天氣的狀態(tài)粗分為三種：1-雨雪天氣、2-多云、3-天晴。假設明天的天氣情況僅和今天的天氣有關，根據(jù)大量的氣象數(shù)據(jù)，我們統(tǒng)計出了今明兩天之間的天氣變化規(guī)律，如圖1中的表格所示。天氣狀態(tài)變化對應的概率圖模型為圖1(b)，其中狀態(tài)值之間的邊的權值表示轉移概率，且每個狀態(tài)的所有指出去的邊的權值之和為1。如果今天為雨天，那么明天下雨的概率為0.4；如果今天為多云，明天出現(xiàn)多云天氣的概率為0.6；如果今天為晴天，明天最有可能出現(xiàn)的天氣依然是晴天，其概率為0.8。基于今天的天氣狀況和明天可能出現(xiàn)的天氣狀況，我們可以估算一個星期后的各種天氣狀況對應的概率。根據(jù)該實例，我們可知運用馬爾科夫模型預測未來天氣狀況需要幾個要素：1)當前的天氣(初始狀態(tài))；2)連續(xù)兩天的天氣之間轉變模式(轉移概率)。

我們進一步用數(shù)學語言對馬爾科夫鏈相應的知識點進行抽象和概括。馬爾科夫鏈是由隨機變量組成的序列，每個隨機變量可能的值組成該馬爾科夫鏈的狀態(tài)空間\(\mathcal{S}\)，且該序列滿足馬氏性，馬氏性的數(shù)學表述如下： \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\)處于狀態(tài)\(i\)的前提下，下一步轉移到狀態(tài)\(j\)的概率。當\(a_{ij}(t)\)與時刻\(n\)無關時，馬爾科夫鏈具有平穩(wěn)轉移概率，下面我們只討論這種情況。設\(A\)為轉移概率組成的矩陣，且狀態(tài)空間\(\mathcal{S}=\{1,2,\cdots,N\}\)，則狀態(tài)轉移矩陣\(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\)的馬爾科夫鏈狀態(tài)序列對應的概率為： \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} 其中\(zhòng)(\pi_i\)表示初始狀態(tài)為\(i\)的概率。如果今天是晴天，那么后面一個星期的天氣狀況觀察序列\(zhòng)(O\)是"sun-sun-rain-rain-sun-cloudy-sun"的概率是多大呢？在前面給出的天氣變化模型中，天氣狀態(tài)有三種\(s1\)(雨雪)、\(s2\)(多云)和\(s3\)(晴)，根據(jù)這個公式，我們可以計算出對應的概率值： \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\)處于狀態(tài)\(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} 結合馬爾科夫鏈狀態(tài)序列概率的求解規(guī)則，有 \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\)各有\(zhòng)(N\)種可能的取值，計算每個狀態(tài)序列的時間復雜度為\(O(t)\)，因此上式總的時間復雜度為\(O(tN^t)\)。
2. 第一種解法復雜度太高，因為針對所有的公共子路徑都會重復計算，如果在圖2中的每個結點都存儲概率值\(P(X_i=x_i)\)，就可以消除重復計算，算法的復雜度也會降低。基于這個想法，我們給出如下的計算規(guī)則： \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有\(zhòng)(N(t-1)+1\)個這樣的結點需要計算，因此總的時間復雜度為\(O(tN^2)\)。這種方法無論是用遞歸還是非遞歸都很容易實現(xiàn)，只不過用非遞歸方式實現(xiàn)時，需要按時間遞增的順序計算\(P(X_i=x_i)\)且需要存儲每個時刻的\(N\)種可能的狀態(tài)值對應的概率。

上述問題涉及到馬爾科夫鏈的\(n\)步轉移概率，表示隨機移動的質點在兩個前后相差\(n\)個時刻的狀態(tài)值分別為\(i\)和\(j\)的概率，數(shù)學描述形式如下： \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\)步狀態(tài)轉移矩陣\(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} 其中\(zhòng)(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} 根據(jù)全概率公式及馬氏性，可證明第三條性質： \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\)步轉移概率矩陣后，求某個狀態(tài)在特定時刻出現(xià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)\)。 假設現(xiàn)有\(zhòng)(L\)條馬爾可夫序列\(zhòng)(\{X^1,X^2,\cdots X^L\}\)，如何基于這些訓練數(shù)據(jù)學習一個最合理的Markov模型呢？我們采用最大似然估計的方法估計Markov模型的最優(yōu)參數(shù)。假設訓練數(shù)據(jù)之間相互獨立，則整個訓練集上的似然函數(shù)的對數(shù)形式如下： \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\)條訓練數(shù)據(jù)\(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乘子法求得使上述似然函數(shù)最大且滿足約束條件的參數(shù)，構造的Lagrange函數(shù)如下： \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函數(shù)分別對參數(shù)\(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函數(shù)分別對參數(shù)\(\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} 根據(jù)上面的參數(shù)學習規(guī)則，很容易看出來這些參數(shù)的學習都是基于統(tǒng)計的，\(a_{ij}\)等于所有從狀態(tài)\(i\)跳轉到狀態(tài)\(j\)的次數(shù)與所有狀態(tài)為\(i\)的比值，而\(\pi_i\)則為\(L\)條訓練數(shù)據(jù)中以狀態(tài)\(i\)開始的比率。

Hidden Markov Model

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

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

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

接下來，我們對HMM進行形式化描述。HMM的模型參數(shù)為 \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} 參數(shù)的具體說明如下：

1. 隱含狀態(tài)的數(shù)目\(N\)。盡管HMM中的狀態(tài)是不可見的，但在很多實際應用中，狀態(tài)是具有一定物理意義的。例如，在瓷缸問題中，每個狀態(tài)代表一個特定的瓷缸。一般而言，狀態(tài)間都是相互連通的，但也有例外。狀態(tài)集合表示為\(\mathcal{S}=\{s_1,s_2,\cdots,s_N\}\)，\(t\)時刻的狀態(tài)標號用\(x_t\)表示，其中\(zhòng)(x_t\in\{1,2,\cdots,N\}\)。
2. 每個狀態(tài)可能對應的觀測值\(M\)。觀測值對應模型的輸出，在瓷缸問題中就是球的顏色。觀察集合表示為\(\mathcal{Y}=\{y_1,y_2,\cdots,y_M\}\)，\(t\)時刻的觀測值標號用\(o_t\)表示，其中\(zhòng)(o_t\in\{1,2,\cdots,M\}\)。 \item 狀態(tài)轉移概率矩陣\(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. 觀測值關于狀態(tài)的輸出概率矩陣\(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. 初始狀態(tài)的概率分布\(\pi=\{\pi_i\}\)，其中 \begin{equation} \pi_i=P(x_0=i),1\leq i\leq N \end{equation}

Inference in HMM

圍繞著HMM有三個基本問題，這也是HMM在現(xiàn)實應用中發(fā)揮作用的關鍵性問題[3]。

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

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

Evaluation

給定模型參數(shù)\(\lambda\)和長度為\(T\)的觀測序列\(zhòng)(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} 其中\(zhòng)(X=x_0x_1\cdots x_{T-1}\)為觀測序列對應的未知的狀態(tài)序列，且有 \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ài)規(guī)劃算法，這里討論前向和后向兩種動態(tài)規(guī)劃的策略。 假設\(t\)時刻的狀態(tài)\(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} 根據(jù)前向算法的推理規(guī)則不斷計算，最后合并\(T-1\)時刻的所有結果即可得到觀測序列的概率： \begin{equation} P(O|\lambda)=\sum_{i=1}^N\alpha_{T-1}(i) \end{equation} 對應的算法描述參見1，算法時間復雜讀為\(O(TN^2)\)，空間復雜度為\(O(N)\)。

假設在\(t\)時刻的狀態(tài)\(x_t=j\)的前提下，\(t+1\)至\(T-1\)的時間段內(nèi)對應的觀察序列為\(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\)。根據(jù)后向算法的推理規(guī)則不斷計算，最后合并\(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

給定模型參數(shù)\(\lambda\)和長度為\(T\)的觀測序列\(zhòng)(O=o_0o_1\cdots o_{T-1}\)，求最優(yōu)的狀態(tài)序列\(zhòng)(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} 基于動態(tài)規(guī)劃的Viterbi算法[2]可用于尋找一條最優(yōu)的隱含狀態(tài)序列。假設已知的前\(t\)個觀測值有多條以狀態(tài)值\(i\)收尾的狀態(tài)序列，這些狀態(tài)序列與觀測序列最大的聯(lián)合概率為 \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\)個觀察值對應的狀態(tài)序列以\(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} 根據(jù)上式不斷向前計算\(\delta_{t}(j)\)直至最后時刻\(T-1\)，此時可知最優(yōu)狀態(tài)序列\(zhòng)(X^\star\)對應的概率為 \begin{equation} P(O,X^\star|\lambda)=\underset{1\leq i\leq N}{\max}\delta_{T-1}(i) \end{equation} 那么在\(T-1\)時刻對應的最優(yōu)狀態(tài)\(x^\star_{T-1}\)為 \begin{equation} x^\star_{T-1}=\underset{1\leq i\leq N}{arg\max}\delta_{T-1}(i) \end{equation} 假設我們已經(jīng)知道最優(yōu)狀態(tài)序列在\(t+1\)至\(T-1\)時時間段內(nèi)對應局部狀態(tài)序列，且\(t+1\)時刻的狀態(tài)值為\(x_{t+1}=j\),那么\(t\)時刻的狀態(tài)\(x_t\)為何值時為最優(yōu)呢？由于\(x_{t+1}\cdots x_{T-1}\)均已確定，意味著最優(yōu)狀態(tài)序列從\(t+1\)時刻開始的狀態(tài)轉移概率和輸出概率對剩下的各狀態(tài)變量而言是常數(shù)了，又\(P(X,O)\)可通過初始狀態(tài)概率及一系列的狀態(tài)轉移概率和輸出概率相乘得到，則 \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\)及以后所有時刻的局部最優(yōu)狀態(tài)序列后，\(t\)時刻的最優(yōu)狀態(tài)值為 \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} 根據(jù)上述分析，在遍歷完整個觀測序列后，\(\Psi\)記錄下的狀態(tài)間的連接關系如圖6所示。在Viterbi算法計算過程中，關鍵是要記錄所有時間段內(nèi)的狀態(tài)鏈連接關系\(\{\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的模型參數(shù)\(\lambda=\{A,B,\pi\}\)為何值時才能是觀測序列集合\(\mathcal{O}\)出現(xiàn)的概率最大？接下來，我們基于前面的前向算法、后向算法，用Lagrange乘子處理參數(shù)的約束條件，最后用EM(Expection Maximization)算法的思想完成模型的訓練問題[4]。 根據(jù)前面知識，觀察序列和隱含狀態(tài)序列的聯(lián)合概率計算形式如下： \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不等式為對數(shù)似然函數(shù)構造一個下界函數(shù)，E步驟中在固定模型參數(shù)的情況下通過使隱含狀態(tài)序列滿足特定的概率分布使得下界函數(shù)與真實目標函數(shù)相等，M步驟則在固定隱含狀態(tài)序列概率分布情況下優(yōu)化模型參數(shù)使得下界函數(shù)取得局部最優(yōu)解；如此不斷重復前面兩個步驟，直至收斂為止。雖然EM算法只能找到局部最優(yōu)解，幸運的是大多時候局部最優(yōu)解已經(jīng)可以滿足我們的要求。假設訓練集中的觀測序列相互獨立，則EM算法中的下界函數(shù)\(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)\)最大的參數(shù)，而且參數(shù)必須滿足下列約束 \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函數(shù) \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} 下面用最大似然的參數(shù)估計方法計算最優(yōu)參數(shù)。Lagrange函數(shù)分別對狀態(tài)轉移概率\(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}\)的參數(shù)求解如下 \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函數(shù)分別對輸出概率\(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)\)的計算規(guī)則為 \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函數(shù)分別對初始狀態(tài)概率\(\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} 結合上述兩式，可得初始狀態(tài)概率\(\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} 是否通過上述的參數(shù)計算規(guī)則，我們就能得到了最優(yōu)參數(shù)呢？因為隱含狀態(tài)序列\(zhòng)(X^l\)未知，所有以無法計算\(Q(X^l)\)。我們先來看兩個基本問題，在已知模型參數(shù)\(\lambda\)和觀測序列\(zhòng)(O\)的前提下，求任意單個時刻或連續(xù)兩個時刻的隱含狀態(tài)為特定值的概率

• 任意連續(xù)兩個時刻的隱含狀態(tài)為特定值的概率 \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}
• 任意單個時刻的隱含狀態(tài)為特定值的概率 \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步驟中，我們有\(zhòng)(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} 將上述三式帶入各參數(shù)更新規(guī)則中，令\(\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

前面的內(nèi)容已經(jīng)從理論上解決了HMM相關的問題，但是要實現(xiàn)HMM相關算法卻有很多細節(jié)要處理，scaling就是其中之一。在解決前面的任何一個問題時，我們都要用到眾多范圍在\([0,1]\)之間的概率值的乘積，所以得到的結果很容易向下溢出。為了處理這中情況，在實現(xiàn)HMM算法時有必要引入scaling的技巧，使得計算機處理的中間結果不至于向下溢出，而且最終可以得到與原問題相同的解。對于聯(lián)合概率\(P(O,X)\)的求解，可采用對數(shù)等價轉換 \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的模型參數(shù)估計過程中，這里根據(jù)文獻[3]中的scaling策略進行簡單分析，最后給出HMM在實際的參數(shù)學習過程中的方法。假設對所有的\(\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} 其中\(zhòng)(c_t=\sum_{i=1}^N\alpha_t(i)\)。前向算法和后向算法中的迭代公式變?yōu)槿缦滦问剑?\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} 經(jīng)過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的最優(yōu)參數(shù)，而且要求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.

Hidden Markov Model