寫在前面

老習慣，正文之前瞎扯一通。HMM學了很久，最初是在《統計學自然語言處理》裡面就學到了相關內容，並且知道HMM CRF一直都是NLP比較底層比較基礎且較為有效的演算法模型（雖然感覺還是挺難的），之前僅僅侷限在瞭解前向演算法和維特比演算法上。也沒有去寫程式碼，只知道個大概思路。最近從52nlpHMM系列講解再次入手，結合多篇部落格、github專案以及李航的《統計學習方法》比較全面的對HMM做了一次學習，要求對自己強制輸出，所以在整體公式推導沒有什麼大問題之後，昨天花了一天完善了程式碼，今天來做一個全面的講解，為人為己。
本文還是堅持自己的風格，講解和公式穿插進行，數學公式永遠是最精煉的語言~^_^~

本文目標

• Why - 什麼場景下需要HMM模型
• What - HMM模型的相關概念定義
• How - HMM模型中的3個經典問題
• Code - python實現一個HMM基礎框架以及簡單應用
• 總結與後續問題

Why - 什麼場景下需要HMM模型

X是一個時間和狀態都是離散隨機過程，Xn是n時刻下的狀態，如果Xn + 1對於過去狀態的條件概率分佈僅是Xn的一個函式，即我們通常說的狀態僅與上一個狀態有關，則該過程是一個一階馬爾科夫過程。公式如下：
P(Xn+1=x∣X0,X1,X2,…，Xn)=P(Xn+1=x∣Xn)

上面回顧了一下馬爾科夫模型，現在來看一個實際的例子：
現在有某地一週的天氣預報：
1 Sun 2 Sun 3 Rain 4 Sun 5 Rain 6 Rain 7 Sun
現在我們假設知道某地以往所有的天氣，現在我們假設這是在沙漠地區，那麼我們有一個2階的狀態轉移矩陣：

         sun rain
sun      0.9  0.1
rain     0.8  0.2

那麼我們現在可以很輕鬆的判斷這一週出現這種天氣的概率是多少。但是突然有一天，你被關進了一個小黑屋，你不能直接知道外面的天氣，你只能看看小屋子裡面苔蘚的乾溼情況，並且你知道天氣能影響苔蘚的乾溼情況。那麼現在如果你想了解天氣情況以及相關的資訊，你該怎麼辦？

這裡先來小小總結一下，我們現實生活中有很多這種類似的例子，上述例子中天氣是我們不能直接觀測的（被關小黑屋的情況）稱為隱藏狀態，苔蘚的乾溼被稱為觀察狀態，觀察到的狀態序列與隱藏過程有一定的概率關係，這時我們使用隱馬爾科夫模型對這樣的過程建模，這個模型包含了一個底層隱藏的隨時間改變的馬爾科夫過程，以及一個與隱藏狀態某種程度相關的可觀察到的狀態集合。我們生活中有很多這種無法被直接觀測的例子，比如語音裝置中的某些內部產出，或者是投骰子中某些骰子的狀態選擇（4,6,8面骰子的案例，網上很多這裡不贅述），或者是詞性標註中的詞語詞性等等。整個HMM模型都是圍繞觀測狀態和隱藏狀態建模和處理的。我們先留個印象，讓我們繼續往後看。

What - HMM模型的相關概念定義

隱馬爾可夫模型（Hidden Markov Model，HMM）是統計模型，它用來描述一個含有隱含未知引數的馬爾可夫過程。以上描述來自百度百科。
HMM模型大概長像如下，偷的52NLP的圖，這個模型包含了一個底層隱藏的隨時間改變的馬爾科夫過程，通常是一階的，對應下圖就是天氣之間的狀態轉移概率關係。以及一個與隱藏狀態某種程度相關的可觀察到的狀態集合，對應就是天氣如何影響苔蘚，即我們的術語叫做發射概率或者混淆概率。

HMM模型的5元組

每個模型都有自己相關的概念，弄清演算法之前我們先來看看這個模型的基本概念，5元組。
（S,K,π,A,B）：以上分別對應了HMM中的5個重要概念
S：隱藏狀態的集合（Sun Cloud Rain），N為隱狀態個數 N=3
K : 輸出狀態或者說觀測狀態的集合（Soggy Damp Dryish Dry），M為觀測狀態個數M=4
π : 對應隱藏狀態的初始化概率 （sun : 0.5 cloud : 0.3 rain : 0.2）
A : 隱藏狀態的狀態轉移概率，是一個N*N的概率矩陣
B : 隱藏狀態到觀測狀態的混淆矩陣，是一個N*M的發射概率矩陣
另外符號體系中會有對於每個序列的狀態序列和觀測序列，如一週的觀測和狀態，定義下
X : 某特定觀測序列
O : 某特定隱藏狀態序列
仍然是盜的圖，方便大家看一下，和我上面寫的基本一致。來自《統計學自然語言處理》

HMM中的3個經典問題

關於HMM這個模型建立起來之後，其實有很多問題可以去討論，但是一般講學中我們討論3個問題，即評估，預測，學習，下面我們來一起看一波吧~
這裡為了解釋方便，請大家無條件的接受以下兩個基礎條件：
1.已知對應問題的觀測狀態集合K
2.已知對應問題的隱藏狀態集合S
如果沒有以上兩個問題，可能就不太屬於今天要討論的範疇了，可能屬於聚類後者其他研究，總之這裡我們不予考慮。
3個經典問題如下：
·觀察序列概率的計算 評估（前向演算法）
·隱藏狀態序列最大概率解 預測（維特比演算法）
·馬爾科夫引數求解（π,A,B） 學習（EM演算法 前後向演算法）

How - HMM模型中的3個經典問題

評估

評估描述

給定觀測序列O（o1,o2,…,oT）和模型u = (π,A,B),求出P（O | u）,即給定模型下觀測序列的概率是多少？對應於之前的例子，就是給定天氣的轉移矩陣，天氣和苔蘚的發射矩陣，以及天氣的初始化列表（這些都是已知的，以前統計好的，具體方法這裡不用糾結）。然後求出給定一週苔蘚的狀態，你判斷這個狀態存在的概率有多大（這個評估這裡只是介紹方法，想想感覺這個案例這裡沒有什麼特別大的實際意義）。

評估理論推導

解決該問題有個很直觀的想法就是把所有的隱藏狀態都走一遍，然後看對應觀測狀態概率有多大，一起加起來就是這個狀態的可能性。我們用數學式子表示如下：
自己寫公式還是很費力的，第二個公式中bxtxt+1ot這種寫法是因為有些HMM模型的發射概率是在發射弧上面，即和該狀態與下狀態有關，所以寫成這種樣子，有時候如果只與當前狀態有關可以寫成btot的形式。
第一個公式：利用全概率公式求解所有可能隱狀態序列下的概率之和。
第二個公式：已知狀態下序列的概率。
第三個公式：任意隱藏序列的概率。
第四個公式：利用每個概率表示公式1，這裡bxtot表示了發射只與當前狀態有關，與2略不同，只是多了一個假設條件便於表示。另外此處說明下由於序列長度不同，該公式可能與其他某些書中公式有點差異，但基本思想一致，只不過具體表現上針對不同情況略有不同。最後一個bXTOT是在連乘之後的，不在求積符號裡面！

好了，公式也推到了，寫寫程式碼就能跑了（完結撒花~）！
其實並沒有進度條君命還長呢，我們仔細看看上面公式，計算一下時間複雜度。一共N^T次方的可能序列，好了打住，不用往後看了，這已經指數級別了。我們可以看到計算公式裡面實際是由大量冗餘乘法計算的，現在給大家介紹動態規劃的前向演算法來巧妙的解決實際計算問題。

評估實際演算法：前向計算

我們定義前向變數運算元αt（i）=P(O1,O2,…,Ot,Xt = Si | u),前向變量表示t時刻，si狀態結束是之前路徑上的總概率。可知，對αT（i）求和便能得到評估結果。而此時時間複雜度因為有了動態規劃，乘法次數在T*N^2的級別，進步不少哦，只需要額外的少量空間（T*N）即可。

python前向演算法程式碼

OK，這一部分差不多到這裡，已經可以評估出這個觀測序列存在的概率了，下面附上一點點python程式碼，以下符號體系和上述相同，應該比較好理解，重點在於熟悉numpy對矩陣的操作，這使得python的程式碼看起來十分的簡潔！

    #計算公式中的alpha二維陣列
    def _forward(self,observationsSeq):
        T = len(observationsSeq)
        N = len(self.pi)
        alpha = np.zeros((T,N),dtype=float)
        alpha[0,:] = self.pi * self.B[:,observationsSeq[0]]  #numpy可以簡化迴圈
        for t in range(1,T):
            for n in range(0,N):
                alpha[t,n] = np.dot(alpha[t-1,:],self.A[:,n]) * self.B[n,observationsSeq[t]] #使用內積簡化程式碼
        return alpha

預測

預測描述

給定觀測序列O（o1,o2,…,oT）和模型u = (π,A,B),求出最大概率的隱藏序列X（X1,X2…,XT），那麼解法思路和上述一樣，只不過將求和變成求最大值，並且相同的思路我們可以利用動態規劃來解決這個問題，這裡該方法有一個較為有名的演算法“維特比演算法”。下面引用了數學之美的一段話。

維特比演算法是一個特殊但應用最廣的動態規劃演算法，利用動態規劃，可以解決任何一個圖中的最短路徑問題。而維特比演算法是針對一個特殊的圖——籬笆網路的有向圖（Lattice )的最短路徑問題而提出的。 它之所以重要，是因為凡是使用隱含馬爾可夫模型描述的問題都可以用它來解碼，包括今天的數字通訊、語音識別、機器翻譯、拼音轉漢字、分詞等。——《數學之美》

維特比演算法

argmax(P（X|O,u）),我們想求出最有可能的X序列（X1，X2,….,XT），我們依葫蘆畫瓢寫出對應的維特比中間變數θt（j） = max(P(X1,X2,..Xt=Sj,O1,O2,…,Ot|u)) （這個公式確實有點難寫出來，不過是想還是很容易理解）。這裡我們根據上一時刻的概率，根據轉移概率求出此刻最可能的情況，以此遞迴，最終能找到最優解。下面是具體推導過程，該過程中引入了一個變數來儲存該節點的前一個路過節點，程式碼中用argmax(θ（j-1）)表示。
推導公式如下，解釋下：
第一個公式：狀態1時刻概率都是初始概率，這裡注意有人用1有人用0，完全看是不是寫程式碼方便了
第二個公式：狀態t時刻為t-1時刻θ變數*轉移概率的最大值
第三個公式：狀態t時刻回溯上一個最大概率路徑。

θ1(j) = πj 
θt(j) = max( θt-1(i) * aij * bj ot)   i∈[1,N]
φt(j) = argmax( θt-1(i) * aij )   i∈[1,N]

給張示意圖吧，有關更多維特比的內容，可以去網上搜索更多資料。

python 維特比演算法程式碼

下面仍然給一段python程式碼：

 #維特比演算法進行預測，即解碼，返回最大路徑與該路徑概率
    def viterbi(self,observationsSeq):
        T = len(observationsSeq)
        N = len(self.pi)
        prePath = np.zeros((T,N),dtype=int)
        dpMatrix = np.zeros((T,N),dtype=float)
        dpMatrix[0,:] = self.pi * self.B[:,observationsSeq[0]]

        for t in range(1,T):
            for n in range(N):
                probs = dpMatrix[t-1,:] * self.A[:,n] * self.B[n,observationsSeq[t]]
                prePath[t,n] = np.argmax(probs)
                dpMatrix[t,n] = np.max(probs)

        maxProb = np.max(dpMatrix[T-1,:])
        maxIndex = np.argmax(dpMatrix[T-1,:])
        path = [maxIndex]

        for t in reversed(range(1,T)):
            path.append(prePath[t,path[-1]])

        path.reverse()
        return maxProb,path

學習

說了半天，終於說到以前沒有說的問題了，之前的兩個問題其實都比較容易理解，學習是HMM經典問題裡面比較難的問題，要求是給定一個觀察序列，求出模型的引數（π,A,B）。想想就很頭疼，這怎麼求，有一個辦法就是拿已經標記好的去統計，這是一種方法。另外就是我們今天要講的前後向演算法了。這個演算法又叫baum-welch演算法，是在EM演算法提出之前就已經存在了，之後被證明為一種特殊的EM演算法。一般的部落格會從後向變數開始往後講，但是我覺得從理解的角度出發，我們還是先來看看這個EM演算法。

EM演算法例項理解

首先我們想一下學習引數我們需要做些什麼，有哪些難點。首先我們不知道隱藏狀態序列是什麼樣子的，再一個我們要求轉移矩陣之類。我們想從概率最大入手或者說極大似然，但是隻有知道u=(π,A,B)我們才能求出最可能的序列，有了最可能的序列我們才有可能估算最可能的引數集。這彷彿是一個迴圈的“雞生蛋，蛋生雞”的故事。此問題就是EM演算法的應用場景，那麼我們舉個簡單的例子，看一下EM演算法的思想。
此處會借用http://blog.csdn.net/zouxy09/article/details/8537620/博主的例子，特此說明，為了保證博文的完整性，這裡簡要說明。更詳細說明請看上面連結。

假設現在有N1個男生，N2個女生。並且男女身高集X1,X2均服從高斯分佈X1~N（μ1，б1^2），X2~N（μ2，б2^2），其中均值和標準差μi,бi(i∈1，2)均是未知引數。我們現在隨機抽取了n個人，假設抽取人身高序列為Oi(i∈1,..,n)。那麼我們希望在u=(π,A,B)引數下使得該序列出現可能性最大。用數學公式表示如下：

現在就是想辦法對求出上式u,最直接的思路就是對（π,A,B）分別求導，在分佈已知（隱藏狀態確認）的情況下還是比較好求，前提是對似然函式L（u）求lg,這樣可以將上面的式子程式設計∑，方便計算。下面給出極大似然的一般步驟：

求最大似然函式估計值的一般步驟：

（1）寫出似然函式；

（2）對似然函式取對數，並整理；

（3）求導數，令導數為0，得到似然方程；

（4）解似然方程，得到的引數即為所求；

那麼現在的問題就是，我們不知道這個人是男的還是女的，我們就不知道怎麼去算這個概率了。那麼既然是一個迴圈，我們隨便找一個點切入好了。假設u0=(π0,A0,B0),那麼我們可以計算出對應一個身高資料是男生還是女生的概率大小，然後我們根據這個比例去算該情況下最可能的引數u1。解釋下下面的公式，其中X為隱藏狀態序列，O為觀測序列，在觀測和之前引數條件下求現在條件期望的最大值對應的新引數。這裡可能有點繞口，大家可以仔細想下，就是對應不同的X序列我們求出一個狀態概率，這個概率來自於上一次的引數。然後根據這個序列我們求出新引數下的概率，然後對其相乘求和，這就是一個期望的概念。這個公式我儘量的想去形象化的理解，但總感覺說不太通，這裡我就不多解釋了，總之，在數學推導上，這個式子是完全可以用EM演算法得到的。雖然這個存在於EM之前，我們這裡就不糾結了，大家預設就好。

到這裡我們實際上就可以利用約束最優化的問題去求解HMM中的引數學習問題了。現在雖然是數學公式，而且看得有點頭暈，但是最少我們可以求出一個新引數了，此處並沒有證明這個演算法是全域性最優或者是單調的，這個不知道當時那個作者有沒有考慮這個問題，總之後來人們證明了baum-welch演算法是EM演算法的一個特例。並且證明了單調性，這個之後再講，我們先看看baum-welch演算法是如何操作的。

baum-welch演算法的思路

注意：因為實在不想自己抄一遍公式，所以後文中符號體系與上述內容略有不符，O表示觀測序列，I表示隱藏狀態序列，λ表示引數即上述u。
有上面這個公式之後，我這裡直接借用李航書中的計算，因為涉及到最優化問題，暫時還沒有系統的學習，所以直接上截圖吧。這裡特別說明（10.33）書中是有特別註釋的，其省略了常數因子1/P(O|λ)，所以寫成了聯合概率的形式。


同理之後得出下面結論：



好了，看著上面的程式碼，都是在以上期望公式+約束最大化產生的結果，我們已經可以快完成整個演算法了，現在的最後一步就是如何去計算公式中的幾個概率。

P(O,i1=i|λ)
p(O,it=i,it+1=j|λ)  

終於到了講前後向演算法的時候了，個人認為前後向演算法的引入就是為了計算上面概率的，而不是一開始就引入了後向變數，這個因果關係上我認為這樣說比較符合當時的推導思路吧。

後向演算法與前向演算法相同，這裡給出後向變數的定義:
βt(i) = P(Ot+1,OT|it=i,λ) 該公式表示在it時刻狀態為i的情況下，後面到T時刻觀測序列的概率情況。
由於最後一個時刻沒有OT+1,所以直接硬性規定β（T） = 1,遞迴公式如下：
β1(i) = 1 i∈[1,N]
βt(i) = ∑j  aij*bjot+1βt+1(j)  i∈[1,N] t∈[1,T-1]

推導P(O,i1=i|λ)
αt(i) = P(O1,O2...Ot,it=i|λ)
βt(i) = P(Ot+1，...,OT|it=i,λ)
P(O,it=i|λ) = P(O1,O2...Ot,it=i,Ot+1,..OT |λ)  
                 = P(O1,O2...Ot,it=i|λ) * P(Ot+1，...,OT|it=i,O1,O2...Ot,λ)   
                 = P(O1,O2...Ot,it=i|λ) * P(Ot+1，...,OT|it=i,λ)
                 = αt(i) * βt(i) 
這裡推導的公式，有一個獨立性的假設，這裡並不太清楚，感覺可能觀測序列就是獨立的吧？不然這個等式也就變不過來了。

好勒，我們再來看兩個變數，這兩個變數分別表示，從某t時刻i狀態的概率，和某t時候i狀態到下個狀態j的概率。這是Baum-welch演算法裡面描述的兩個期望。


根據以上兩個期望，我們帶入重估函式裡面，得到如下結果。

上面式子雖然是數學推導而來，但是同樣具有特別強的概率含義。

這次真的完結撒花了！小結一下，我們從EM演算法的思想入手，考慮了我們的迭代函式，對就是這個期望
然後我們直接對他求導，差分開來，利用拉格朗日乘子法求出了各個重估函式，並且他們具有很明顯的概率意義。之後我們引入了前後向演算法來表示P（O,it=i|λ）概率，這一步完全是方便計算。最後我們引入兩個中間期望變數γ和ε。好了下面給出程式碼啦^_^~

python程式碼baum-welch演算法

 #計算公式中的beita二維陣列
    def _backward(self,observationsSeq):
        T = len(observationsSeq)
        N = len(self.pi)
        beta = np.zeros((T,N),dtype=float)
        beta[T-1,:] = 1
        for t in reversed(range(T-1)):
            for n in range(N):
                beta[t,n] = np.sum(self.A[n,:] * self.B[:,observationsSeq[t+1]] * beta[t+1,:])
        return beta

    #前後向演算法學習引數
    def baumWelch(self,observationsSeq,criterion=0.001):
        T = len(observationsSeq)
        N = len(self.pi)

        while True:
            # alpha_t(i) = P(O_1 O_2 ... O_t, q_t = S_i | hmm)
            # Initialize alpha
            alpha = self._forward(observationsSeq)

            # beta_t(i) = P(O_t+1 O_t+2 ... O_T | q_t = S_i , hmm)
            # Initialize beta
            beta = self._backward(observationsSeq)

            #根據公式求解XIt(i,j) = P(qt=Si,qt+1=Sj | O,λ)
            xi = np.zeros((T-1,N,N),dtype=float)
            for t in range(T-1):
                denominator = np.sum( np.dot(alpha[t,:],self.A) * self.B[:,observationsSeq[t+1]] * beta[t+1,:])
                for i in range(N):
                    molecular = alpha[t,i] * self.A[i,:] * self.B[:,observationsSeq[t+1]] * beta[t+1,:]
                    xi[t,i,:] = molecular / denominator

            #根據xi就可以求出gamma，注意最後缺了一項要單獨補上來
            gamma = np.sum(xi,axis=2)
            prod = (alpha[T-1,:] * beta[T-1,:])
            gamma = np.vstack((gamma, prod /np.sum(prod)))

            newpi = gamma[0,:]
            newA = np.sum(xi,axis=0) / np.sum(gamma[:-1,:],axis=0).reshape(-1,1)
            newB = np.zeros(self.B.shape,dtype=float)

            for k in range(self.B.shape[1]):
                mask = observationsSeq == k
                newB[:,k] = np.sum(gamma[mask,:],axis=0) / np.sum(gamma,axis=0)

            if np.max(abs(self.pi - newpi)) < criterion and \
                                    np.max(abs(self.A - newA)) < criterion and \
                                    np.max(abs(self.B - newB)) < criterion:
                break

            self.A,self.B,self.pi = newA,newB,newpi

Baum-welch演算法可行性

Baum-welch演算法的公式實際就是EM的Q函式，下面大致看下怎麼推到的吧，其中用到了Jensen不等式,還是截圖吧，只能理解，還不能把推倒理由講徹底，這就是數學的魅力吧！

Jensen公式很重要的一點就是把log給變了位置，這樣後續計算也是十分有利的，方便引入約束最優化問題。但是Baum-welch過程中貌似直接用期望公式的變化達到了這一點，本人暫時還沒想特別明白。之後有人證明了EM演算法的區域性最優性，這是一個收斂的演算法，並且區域性最優。大家知道這些就差不多了，更為詳細的還是去看數學書吧。

詞性標註例項

（此時，這部落格已經寫了8個小時了。）
上面講完了所有的程式碼和推導，現在舉個簡單的小例子來看看HMM的效果。
先坑，有空補

一些問題

多樣本標註：
現在我們可以根據一個觀測值去計算引數，但是對於多個序列，我們如何去操作呢，這是一個很實際的問題，能不能像神經網路裡面一下，每個樣本修改一點引數，到達較好的效果，網上有一些說法，可以取均值或者其他，這個如果有大神瞭解，歡迎指導。

Gaussian-HMM:
hmmlearn裡面碰到了Gaussian-HMM的東西，一般我們假設觀測和隱狀態是離散的，不過可能並不都是這樣。這是考慮觀測連續且服從高斯分佈的情況

github上面有比較完整的程式碼，想看的可以看下下，暫時沒有提供完整的應用程式碼，只是簡單的演算法實現，以及基礎工具。
https://github.com/Continue7777/HMM/