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.

轉載于:https://www.cnblogs.com/jeromeblog/p/3832751.html