HMM隱馬爾科夫算法(Hidden Markov Algorithm)初探

1. HMM背景

0x1：概率模型 - 用概率分布的方式抽象事物的規律

機器學習最重要的任務，是根據一些已觀察到的證據（例如訓練樣本）來對感興趣的未知變量（例如類別標記）進行估計和推測。

概率模型（probabilistic model）提供了一種描述框架，將學習任務歸結於計算未知變量的概率分布，而不是直接得到一個確定性的結果。

在概率模型中，利用已知變量推測未知變量的分布稱為“推斷（inference）”，其核心是如何基於可觀測變量推測出未知變量的條件分布。

具體來說，假定所關心的變量集合為 Y，可觀測變量集合為 O，其他變量的集合為 R。

生成式模型（generative model）考慮聯合分布 P（Y，R，O）；

判別式模型（discriminative model）考慮條件分布 P（Y，R | O）；

給定一組觀測變量值，推斷就是要由 P（Y，R，O）或 P（Y，R | O）得到條件概率 P（Y | O）。

直接利用概率求和規則消去變量 R 顯然不可行，因為即使每個變量僅有兩種取值的簡單問題，其復雜度已至少是。另一方面，屬性變量之間往往存在復雜的聯系，因此概率模型的學習，即基於訓練樣本來估計變量分布的參數往往相當困難。

為了便於研究高效的推斷和學習算法，需要有一套能夠簡潔緊湊地表達變量間關系的工具。

0x2：概率圖模型（probabilistic graphical model）

概率圖模型是一類用圖來表達變量相關關系的概率模型。它以圖為表示工具，最常見的是用一個結點表示一個或一組隨機變量，結點之間的邊表示變量間的概率相關關系，即“變量關系圖”。

根據邊的性質不同，概率圖模型可大致分為兩類：

1. 使用有向無環圖表示變量間依賴關系：有向圖模型或貝葉斯網（Bayesian network）；
2. 使用無向圖表示變量間的相關關秀：無向圖模型或馬爾科夫網（Markov network）；

0x3：隱馬爾科夫模型

隱馬爾科夫模型（Hiden Markov Model HMM）是結構最簡單的動態貝葉斯網（dynamic Bayesian network），這是一種著名的有向圖模型，在時序數據建模上比較適用。

0x4: 馬爾科夫模型適用的場景 - 處理順序數據

我們知道，在進行貝葉斯估計、SVM等分類的時候，我們常常基於一個大前提假設：數據集里的數據是獨立同分布的。這個假設使得我們將似然函數表示為在每個數據點處計算的概率分布在所有數據點上的乘積。但是在實際應用中，還大量存在另一種場景是獨立同分布不成立的，典型的如順序數據的數據集，這些數據通常產生於沿着時間序列進行的測量（即帶時序特征的樣本），例如某個特定地區的連續若干天的降水量測量，或者每天匯率的值，或者對於語音識別任務在連續的時間框架下的聲學特征，或者一個英語句子中的字符序列。

需要注意的是：時序只是順序數據的一個特例，馬爾科夫模型適用於所有形式的順序數據

0x5：筆者對隱馬爾科夫中“隱序列”的思考

HMM的核心思想認為，任何在現實世界觀測到的現象，其背后一定有一個“本質規律”在支撐這一事實，“本質規律”衍生出被我們觀測到的現象。

其實隱狀態序列只是我們的一個副產品，HMM訓練最重要的事情是推斷概率轉移矩陣以及發射矩陣。

和RNN類似，觀測序列之間不存在依賴，而是在隱狀態序列之間存在依賴，區別在於HMM中這種依賴是有限的（N階），而RNN中這種依賴是一直遞歸到初始狀態的（雖然實際上梯度不可能傳遞到初始狀態）

通過HMM進行解碼問題（例如語音翻譯），即根據觀測序列提取隱狀態序列的這種思想，和將DNN中隱藏層作為輸入向量的抽象特征的思想很類似。

Relevant Link:

《機器學習》周志華

 

2. 從一個虛構的例子引入馬爾科夫模型討論

在開始討論HMM的數學公式以及計算方法之前，我們先來討論一個很有趣的例子。在這里例子中，筆者並不會完全遵從HMM的專業術語，而是采用最符合“直覺”的描述語言。雖然可能定義並不准確，但是希望向讀者朋友傳遞一個推導和理解的思想，我們的HMM為什么會被創造出來，以及HMM具體解決了什么問題。

0x1: 實驗例子說明

擲骰子。假設我手里有三個不同的骰子

1. 第一個骰子是我們平常見的骰子（稱這個骰子為D6），6個面，每個面（1，2，3，4，5，6）出現的概率是1/6
2. 第二個骰子是個四面體（稱這個骰子為D4），每個面（1，2，3，4）出現的概率是1/4
3. 第三個骰子有八個面（稱這個骰子為D8），每個面（1，2，3，4，5，6，7，8）出現的概率是1/8

我們構造的場景就是擲骰子實驗。

1.實驗假設前提

1. 每次從三個骰子里挑一個，挑到每一個骰子的概率未知；
2. 每次選的骰子對下一次選擇是否有影響也是未知；
3. 對每個骰子來說拋到某一個面的概率是相等的（即骰子是沒有動過手腳的）

2. 實驗結果數據采集

我們不斷重復擲骰子，每次得到一個數字，不停的重復該過程，我們會得到一串數字

經過10次實驗，我們得到一串數字：1 6 3 5 2 7 3 5 2 4

3. 實驗目標

得到觀測結果后，我們的目標是進行一些列概率估計：

1. 我們得到了這一串數字，每次拋骰子具體是3個中的哪個骰子是未知的，現在需要估計出整個10次實驗中最可能的骰子類型是什么（即每次實驗對應的骰子）；
2. 我們得到了這一串數字，每次拋骰子具體是3個中的哪個骰子是未知的，現在需要估計出現這個觀測結果（這一串數字）的概率有多大

這2個問題非常具有代表性，它們實際上代表了HMM問題中的解碼和預測兩類核心問題。我們來稍微深入一些討論這個問題。

0x2: 從獨立同分布假設開始 - 極簡化的HMM假設

寫這章的目的是為了能更好的闡述HMM的序列依賴假設，對HMM熟悉的讀者可能會挑戰說，HMM至少是1階的，即至少存在“一步依賴”關系。沒錯！！

筆者在這里作的獨立同部分假設，本質上可以理解為是0階HMM假設，讀者朋友在閱讀完這章后，可以對比下獨立同部分假設和N階序列依賴假設在對時序數據進行模式提取的時候的不同變現，從而在實際項目中選擇是否需要應用HMM去解決你的實際問題。

1. 獨立同分布假設

對上面的問題，我們先建立一個基本假設：

1. 每次實驗選擇什么骰子是完全獨立隨機的（即 1/3）；
2. 每次實驗和上一次實驗也不存在關系，即實驗之間是獨立同分布的；

有了這個基本假設，我們現在可以開始進行最大似然估計（你當然也可以進行極大后驗估計，我們在這里不增加討論的復雜度）

2. 我們得到了這一串數字，每次拋骰子具體是3個中的哪個骰子是未知的，現在需要估計出整個10次實驗中最可能的骰子類型是什么（即每次實驗對應的骰子）

以第一輪實驗說明問題，我們寫出概率表達式：P（S骰子= Sx，Res結果=1 | Model）

P（S骰子= S1，Res結果=1 | Model） = 1 / 6

P（S骰子= S2，Res結果=1 | Model） = 1 / 4

P（S骰子= S3，Res結果=1 | Model） = 1 / 8

這個場景很簡單，求最大似然的解不需要求導即可得是 P（S骰子= S2，Res結果=1 | Model），即第一輪最有可能的骰子類型是2號骰子

同理根據最大似然估計原理可得全部10次實驗的骰子類型序列：

Res：1，S = S2
Res：6，S = S1
Res：3，S = S2
Res：5，S = S1
Res：2，S = S2
Res：7，S = S3
Res：3，S = S2
Res：5，S = S1
Res：2，S = S2
Res：4，S = S2

插一句題外話，再仔細觀察下數據，會發現Res在【1，4】范圍中，最大似然結果都是S2，Res在【5，6】范圍中，最大似然結果都是【S1，S2】范圍中，而Rest在【7，8】，最大似然結果都是【S3】范圍中。這表明我們可以訓練一個決策樹模型去建立一個分類器，並且一定能得到較好的效果。

3. 我們得到了這一串數字，每次拋骰子具體是3個中的哪個骰子是未知的，現在需要估計出現這個觀測結果（這一串數字序列）的概率有多大

還是在獨立同分布假設下，根據概率分布加法原理，我們可得每次實驗出現指定結果的概率：

P（Res結果=1 | Model） = P（S骰子= S1，Res結果=1 | Model） + P（S骰子= S2，Res結果=1 | Model）+ P（S骰子= S3，Res結果=1 | Model） = 1 / 6 + 1 / 4 + 1 / 8 = 13 / 24
P（Res結果=6 | Model） = P（S骰子= S1，Res結果=1 | Model） + P（S骰子= S2，Res結果=6 | Model）+ P（S骰子= S3，Res結果=6 | Model） = 1 / 6 + 0 + 1 / 8 = 7 / 24
P（Res結果=3 | Model） = P（S骰子= S1，Res結果=1 | Model） + P（S骰子= S2，Res結果=3 | Model）+ P（S骰子= S3，Res結果=3 | Model） = 1 / 6 + 1 / 4 + 1 / 8 = 13 / 24
P（Res結果=5 | Model） = P（S骰子= S1，Res結果=1 | Model） + P（S骰子= S2，Res結果=5 | Model）+ P（S骰子= S3，Res結果=5 | Model） = 1 / 6 + 0 + 1 / 8 = 7 / 24
P（Res結果=2 | Model） = P（S骰子= S1，Res結果=1 | Model） + P（S骰子= S2，Res結果=2 | Model）+ P（S骰子= S3，Res結果=2 | Model） = 1 / 6 + 1 / 4 + 1 / 8 = 13 / 24
P（Res結果=7 | Model） = P（S骰子= S1，Res結果=1 | Model） + P（S骰子= S2，Res結果=7 | Model）+ P（S骰子= S3，Res結果=7 | Model） = 0 + 0 + 1 / 8 = 1 / 8
P（Res結果=3 | Model） = P（S骰子= S1，Res結果=1 | Model） + P（S骰子= S2，Res結果=3 | Model）+ P（S骰子= S3，Res結果=3 | Model） = 1 / 6 + 1 / 4 + 1 / 8 = 13 / 24
P（Res結果=5 | Model） = P（S骰子= S1，Res結果=1 | Model） + P（S骰子= S2，Res結果=5 | Model）+ P（S骰子= S3，Res結果=5 | Model） = 1 / 6 + 0 + 1 / 8 = 7 / 24
P（Res結果=2 | Model） = P（S骰子= S1，Res結果=1 | Model） + P（S骰子= S2，Res結果=2 | Model）+ P（S骰子= S3，Res結果=2 | Model） = 1 / 6 + 1 / 4 + 1 / 8 = 13 / 24
P（Res結果=4 | Model） = P（S骰子= S1，Res結果=1 | Model） + P（S骰子= S2，Res結果=4 | Model）+ P（S骰子= S3，Res結果=4 | Model） = 1 / 6 + 1 / 4 + 1 / 8 = 13 / 24

最后，求10次實驗的聯合概率即為出現指定觀測序列的概率：

(13 / 24) * (7/24) *  (13 / 24) * (7/24) * (13 / 24) * (1/8) * (13 / 24) * (7/24) * (13 / 24) * (13 / 24) = 7.83363029759e-05

0x3: 馬爾科夫假設

現在把問題放到馬爾科夫的理論框架下討論。

1. N階序列依賴假設

好，現在我們修改假設前提：

1. 從一次實驗到下一次選擇什么骰子是完全隨機的（即 1/3 ），也即狀態轉移概率為 1/3；
2. 每次實驗之前不再是獨立同分布的，每次實驗都與之前的N次實驗有依賴關聯；

這個假設就是N階序列依賴性假設，每次擲的骰子取決於之前幾次擲的骰子，即擲骰子實驗之間存在轉換概率。

2. 觀測序列

在實驗觀測中，得到的一串數字：1 6 3 5 2 7 3 5 2 4 叫做可見狀態鏈，也即觀測序列。

3. 隱藏狀態序列

在HMM馬爾科夫假設中，每次實驗之間並不是獨立同分布的而是存在時序依賴關系，但是實際問題中我們永遠不知道實際的隱含狀態鏈，因為極大似然估計的也不一定就是真實的情況）。

4. 隱藏狀態轉移矩陣

這導致擲骰子之間的轉換概率（transition probability）分布是未知的（這等價於EM算法中的隱變量是未知的）

5. 輸出概率矩陣

同樣的，盡管可見狀態之間沒有轉換概率，但是隱含狀態向可見狀態之間有一個概率叫做輸出概率（emission probability）。

就我們的例子來說：

六面骰（D6）產生1的輸出概率是1/6。產生2，3，4，5，6的概率也都是1/6；
四面骰（D4）產生1的輸出概率是1/4。產生2，3，4的概率也都是1/4；
八面骰（D8）產生1的輸出概率是1/8。產生2，3，4，5，6，7，8的概率也都是1/8；

6. 我們得到了這一串數字（可見狀態鏈已知），現在需要估計出整個10次實驗中最可能的骰子類型是什么（即求隱狀態序列）

在實驗觀測中，得到的一串數字：1 6 3 5 2 7 3 5 2 4。現在要在馬爾科夫假設的前提下求每次實驗最可能的篩子序列。

這本質上是在根據實現結果（可見狀態鏈），求隱藏狀態序列時序的最大似然估計。很顯然，因為實驗之間不再獨立同分布了，我們無法像上面那樣逐個計算了。

表面上看，其實最簡單而暴力的方法就是窮舉所有可能的骰子序列，然后把每個序列對應的概率算出來（像前面的算式那樣）。然后我們從里面把對應最大概率的序列挑出來就行了。

但是要注意的每次擲骰子實驗之間並不是獨立同分布的，所以所有序列要計算的次數是：1 + 2! + 3! + ... + 序列Length！，如果馬爾可夫鏈不長，當然可行。如果長的話，窮舉的數量太大，就很難完成了。

另外一種很有名的算法叫做Viterbi algorithm. 要理解這個算法，算法的數學公式和推導我們在本章的后面會討論，這里我們用簡單的描述先來嘗試理解它的思想：局部最優遞推思想。

首先，如果我們只擲一次骰子：

實驗觀測結果為1，根據最大似然原理對應的最大概率骰子序列就是D4，因為D4產生1的概率是1/4，高於1/6和1/8，所以P1（D4） = 1 / 3

以P1（D4）作為起點，把這個情況拓展（遞推），我們擲兩次骰子：

實驗結果為1，6，我們計算三個值，分別是第二個骰子是D6，D4，D8的最大概率。

根據最大似然原理，第二個骰子應該是D6，第二個骰子取到D6的最大概率是：
（復用了第一步的結果）
同樣的，我們可以計算第二個骰子是D4或D8時的最大概率。我們發現，第二個骰子取到D6的概率最大。

而使這個概率最大時，第一個骰子為D4。所以最大概率骰子序列就是D4 D6。

繼續拓展，我們擲三次骰子：

同樣，我們計算第三個骰子分別是D6，D4，D8的最大概率。我們再次發現，要取到最大概率，第二個骰子必須為D6（注意，當轉移概率不是1/3等分的時候，當前步的前次依賴不一定是上一次的計算結果，也可能是其他值）。

這時，第三個骰子取到D4的最大概率是
同上，我們可以計算第三個骰子是D6或D8時的最大概率。我們發現，第三個骰子取到D4的概率最大。而使這個概率最大時，第二個骰子為D6，第一個骰子為D4。所以最大概率骰子序列就是D4 D6 D4。

可以看到，擲多少次都可以以此類推。最后得到最有可能的篩子序列。

我們發現，我們要求最大概率骰子序列時要做這么幾件事情

1. 首先，不管序列多長，要從序列長度為1算起，算序列長度為1時取到每個骰子的最大概率；

2. 然后，逐漸增加長度，每增加一次長度，重新算一遍N步長度前到這個長度下最后一個位置取到每個骰子的最大概率（如果是1階HMM只要算上一步即可），
注意，不是說直接把上一步結果拿過來直接用就可以了，在轉移概率不是等比隨機的情況下，上一步不是概率最大的路徑可能在這一步變為最大。
但是好消息是，N步長度前到這個長度長度下的取到每個骰子的最大概率都算過了，計算過程中只要保存下來重新算並不難。
按照這個方法，當我們算到最后一位時，就知道最后一位是哪個骰子的概率最大了。

3. 最后，我們把對應這個最大概率的序列從后往前推出來，就得到最大概率骰子序列

7. 我們得到了這一串數字（可見狀態鏈已知），現在需要估計狀態轉移矩陣以及輸出概率矩陣（即估計模型參數）

換言之，什么樣的模型能最好地描述這個觀測數據，

8. 我們得到了這一串數字（可見狀態鏈已知），每次拋骰子類型已知（轉換概率已知），現在需要估計出現這個觀測結果的概率有多大

其實本質上這相當於隱變量已知，則問題就退化為一個普通的概率求解問題，解決這個問題的算法叫做前向算法（forward algorithm）。

隱狀態鏈以及觀測結果序列如下：

同時假設一個轉換概率：

1. D6的下一個狀態是D4，D6，D8的概率都是1/3
2. D4的下一個狀態是D4，D6，D8的概率都是1/3
3. D8的下一個狀態是D4，D6，D8的概率都是1/3

則出現該觀測結果的概率計算公式為：

P =

P（D6） * P（D6 -> 1） *

P（D6 -> D8）* P（D8 -> 6） *

P（D8 -> D8）* P（D8 -> 3） *

P（D8 -> D6）* P（D6 -> 5） *

P（D6 -> D4）* P（D4 -> 2） *

P（D4 -> D8）* P（D8 -> 7） *

P（D8 -> D6）* P（D6 -> 3） *

P（D6 -> D6）* P（D6 -> 5） *

P（D6 -> D4）* P（D4 -> 2） *

P（D4 -> D8）* P（D8 -> 4） =

1 / 3 * 1 / 6 * 1 / 3 * 1 / 8 * 1 / 3 * 1 / 8 * 1 / 3 * 1 / 6 * 1 / 3 * 1 / 4 * 1 / 3 * 1 / 8* 1 / 3 * 1 / 6 * 1 / 3 * 1 / 6 * 1 / 3 * 1 / 4 * 1 / 3 * 1 / 8 = 1.99389608506e-13

可以看到，HMM衍生於EM，但是又在EM的基礎上增加了一些特性，例如

1. 在HMM中，隱變量是轉換概率。知道了轉換概率，根據實驗觀測結果可以最大似然推測出時序
2. HMM中，隱變量和觀測變量的關系是前者決定后者的關系，即輸出概率

Relevant Link:

http://blog.163.com/zhoulili1987619@126/blog/static/353082012013113191924334/ 
https://www.zhihu.com/question/20962240

 

3. 馬爾科夫模型介紹

0x1：隱馬爾可夫模型圖結構

如上圖所示，隱馬爾科夫模型中的變量可分為兩組：

第一組：狀態變量，其中表示第 i 時刻的系統狀態。通常假定狀態變量是隱藏的、不可被觀測的，因此狀態變量也被稱為隱變量（hidden variable）

第二組：觀測變量，其中表示第 i 時刻的觀測值。

在隱馬爾可夫模型中，系統通常在多個狀態之間轉換，因此狀態變量的取值范圍 Y（狀態空間）通常是有 N 個可能取值的離散空間。

圖中的箭頭表示了變量間的依賴關系，它是一種動態貝葉斯網。

在任一時刻，觀測變量的取值僅依賴於狀態變量，即由確定，與其他狀態變量及觀測變量的取值無關。這是由HMM的概率圖結構決定的。

當然，讀者在研究RNN及其相關擴展的時候，會看到大量其他的變種模型，變量之間的依賴會變得非常的復雜同時也帶來新的能力。但是對於HMM，我們要牢記上面這個概率圖模型，HMM規定了這種概率圖模型，這也限定了HMM能夠具備的能力以及解決問題的范圍。

同時，在 N 階馬爾科夫模型中，t 時刻的隱狀態僅依賴於 t - N 時刻的狀態 Yt-n，與其余狀態無關。這就是所謂的“馬爾柯夫鏈（Markov chain）”。

0x2: 馬爾科夫假設 - 序列依賴

這個小節我們來討論下馬爾科夫假設，它是什么？為什么要這么假設？這種假設的合理性在哪里？

1. 對含有序列特征的數據集建立獨立同分布假設不合理

為了理解馬爾科夫假設，我們先從獨立同分布假設開始（朴素貝葉斯中采用了獨立同分布假設）。

毫無疑問，處理序列數據的最簡單的方式是忽略數據集中的序列性質，而將所有觀測看作是獨立同分布

但是這種方法缺點就是無法利用數據中的序列模式（在很多場景下這種假設都是不成立的）

2. 考慮當前狀態和所有歷史狀態的依賴關系 - 一種極端的假設

我們來對假設進行一些改進，我們假設：每一個當前狀態都與此之前的所有狀態存在依賴關系，即今天是否下雨其實是由今天之前的整個慢慢歷史長河（從宇宙爆發開始）的所有日子是否下雨的結果而決定的（因果論），即：

從某種程度上來說，這種假設沒有問題，理論上它是成立的，但問題在於計算量太過巨大，而且實際上距離太過久遠之前的序列對當前序列的影響實際上是微乎其微的。

我們為了預測明天的天氣，要把地球形成之初到昨天所有的天氣情況都輸入模型進行訓練，這種計算量是不可接受的。

3. 一階馬爾科夫假設 - 只考慮相鄰上一步的序列依賴關系

為了緩解獨立同分布和完整歷史序列依賴的缺點，馬爾科夫假設建立了一個近似逼近，對上面的聯合概率分布，進行一個改進：

我們現在假設右側的每個條件概率分布只與最近的⼀次觀測有關，⽽獨⽴於其他所有之前的觀測，那么我們就得到了⼀階馬爾科夫鏈（first-order Markov chain）

從定義可知，一階隱馬爾可夫鏈作了兩個基本假設：

1）齊次馬爾科夫性假設

假設隱藏的馬爾柯夫鏈（隱狀態序列）在任意時刻 t 的狀態只依賴於前一時刻的狀態，與其他時刻的狀態及觀測無關，也與當前時刻 t 無關：

2）觀測獨立性假設

假設任意時刻的觀測只依賴於該時刻的馬爾柯夫鏈的狀態（觀測與隱狀態一一對應），與其他觀測及狀態無關：

0x3: 隱馬爾可夫模型的定義

隱馬爾可夫模型是關於時序的概率模型。

1. 模型結構信息

一階馬爾科夫模型的結構如下：

2. 狀態轉移概率

模型在各個狀態間轉換的概率，通常記為矩陣 A = ，其中 N 是所有可能的狀態數，

表示在任意時刻 t，若狀態為，則在下一時刻狀態為的概率。

3. 輸出觀測概率

模型根據當前狀態獲得各個觀測值的概率，通常記為矩陣，其中 M 是所有可能的觀測的類型數。

表示在任意時刻 t，若狀態為，則觀測值被獲取的概率。

4. 初始狀態概率

模型在初始時刻各狀態出現的概率，通常記為 = ，其中

表示模型的初始狀態（t = 1）為的概率。

確定上述參數，就能確定一個隱馬爾可夫模型。

隱馬爾可夫模型由初始狀態概率向量、狀態轉移概率矩陣A、和觀測概率矩陣B決定。

和A決定狀態序列；B決定觀測序列。

因此，隱馬爾可夫模型可以用三元符號表示，即：，A，B，稱為隱馬爾可夫模型的三要素。

0x4: 隱馬爾可夫模型的3個基本問題

理解隱馬爾可夫模型的3個基本問題筆者認為非常重要，這3個問題是對現實世界實際問題的數學抽象。

1. 觀測序列概率計算 - 評估模型與觀測序列之前的匹配程度

給定模型和觀測序列，計算在模型下觀測序列 O 出現的概率。

這個問題場景也常使用維特比算法（Viterbi algorithm）、前向后向算法（Forward-Backward algorithm）這類遞推算法來完成概率似然計算。

這個場景是異常檢測中最常見的形式，通常我們根據一組訓練樣本訓練得到一個HMM模型，之后需要根據這個模型預測未來的未知觀測序列，根據概率計算結果評估這些新序列的異常程度。

這類問題還有另一個常用的場景，根據以往的觀測序列來推測當前時刻最有可能的觀測值，顯然可以通過softmax來確定概率最大的下一步序列。

2. 學習問題 - 模型訓練

已知觀測序列，估計模型的參數，使得在該模型下觀測序列概率最大，即用極大似然估計的方法估計參數。

在大多數現實應用中，人工指定模型參數已變得越來越不可行。

對這類問題的求解，需要使用EM算法來完成在存在隱變量條件下的最大似然估計，例如Baum-Welch algorithm。

3. 預測問題/解碼問題 - 反推隱藏狀態序列

已知模型和觀測序列，求對給定觀測序列條件概率最大的隱狀態序列。

即給定觀測序列，求最有可能的對應的隱狀態序列。

這常見於語言翻譯，語言解碼翻譯等問題（從待翻譯語言的空間映射到目標語言的空間中）。

值得注意的是，HMM的隱狀態序列的最大似然估計和傳統的統計概率最大似然估計不同在於HMM中狀態之間還要考慮前序狀態的依賴，即當前的狀態和前面N（N階馬爾科夫模型）個狀態還存在一定的轉移概率關系，這使得HMM的最大似然估計不能簡答的進行計數統計分析，而需要利用維特比算法（Viterbi algorithm）、前向后向算法（Forward-Backward algorithm）這類遞推算法來完成序列最大似然估計。

Relevant Link:

https://www.zhihu.com/question/20962240
http://www.cnblogs.com/skyme/p/4651331.html
http://blog.csdn.net/bi_mang/article/details/52289087
https://yanyiwu.com/work/2014/04/07/hmm-segment-xiangjie.html
https://www.zhihu.com/question/20962240
https://www.zhihu.com/question/52778636
https://www.zhihu.com/question/37391215
https://www.zhihu.com/question/20551866
https://www.zhihu.com/question/28311766
https://www.zhihu.com/question/28311766
https://www.zhihu.com/question/22181185

 

4. 觀測序列概率計算算法 - 馬爾科夫模型中的觀測序列概率計算

概率計算算法的目標是給定模型和觀測序列，計算在模型下觀測序列 O 出現的概率。

0x1: 直接計算法

給定模型和觀測序列，計算觀測序列 O 出現的概率，最直接的方法是按照概率公式直接計算，通過枚舉所有可能的長度為 T 的狀態序列，並求各個狀態序列 I 與觀測序列的聯合概率，然后對所有可能的狀態序列求和，得到

狀態序列的概率是：

對固定的狀態序列，觀測序列的概率是，

O 和 I 同時出現的聯合概率為：

然后，對所有可能的狀態序列 I 求和，得到觀測序列 O 的概率，即

但是，該公式的計算量很大，是階的，這種算法在實際應用中不可行，取而代之的是接下來要討論的前向-后向算法（forward-backward algorithm）。

0x2: 前向-后向算法

1. 前向算法 - 一種遞歸算法

首先定義前向概率，給定隱馬爾可夫模型，定義到時刻 t 部分觀測序列為，且當前時刻狀態為的概率為前向概率，記作：

可以遞推地求得所有時刻的前向概率及觀測序列概率。算法過程如下

（1）初值：，是初始時刻的狀態和觀測的聯合概率

 

（2）遞推：對，有：，乘積就是時刻 t 觀測到並在時刻 t 處於狀態而在時刻 t + 1 到達狀態 的聯合概率。對這個乘積在時刻 t 的所有可能的 N 個狀態求累加和，其結果就是到時刻 t 觀測為並在時刻 t + 1處於狀態的聯合概率

（3）終止：

可以看到，前向算法實際是基於“狀態序列的路徑結構”遞推計算。前向算法高效的關鍵是其局部計算前向概率，然后利用路徑結構將前向概率“遞推”到全局。得到。
具體地，在各個時刻，計算的N個值（i = 1，2，...，N），而且每個的計算都要利用前一時刻N的。

減少計算量的原因在於每一次計算直接引用前一時刻的計算結果，避免重復計算。這樣，利用前向概率計算的計算量是階，而不是直接計算的階

前向算法是一種計算優化算法，本質是把中間計算結果緩存的直接概率計算法。

2. 前向算法的一個具體應用舉例

考慮盒子和球模型

狀態的集合是：

球的顏色對應觀測，觀測的集合是：

狀態序列和觀測序列長度 T = 3

初始概率分布為：

狀態轉移概率矩陣為：

觀測概率矩陣為：

觀測序列集合V為：O = ｛紅，白，紅｝

現在要計算的問題是：出現該觀測序列的最大似然概率是多少？

1）計算初值

在第一輪初始值的計算中，所有狀態都沒有t-1狀態，只有初始狀態自己。

a1（1） = π1 * b1（o1）= 0.2 * 0.5 = 0.1

a1（2） = π2 * b2（o1）= 0.4 * 0.4 = 0.16

a1（3） = π3 * b3（o1）= 0.4 * 0.7 = 0.28

2）遞推計算 t = 2

從t=2開始，計算每個狀態的概率都要累加上一輪所有其他狀態轉移到該狀態的概率和。

= （a1（1） * a11 + a1（2） * a21 + a1（3） * a31） * b1（o2） = （0.1 * 0.5 + 0.16 * 0.3 + 0.28 * 0.2）* 0.5 = 0.077

= （a1（1） * a12 + a1（2） * a22 + a1（3） * a32） * b2（o2） = （0.1 * 0.2 + 0.16 * 0.5 + 0.28 * 0.3）* 0.6 = 0.184 * 0.6 = 0.1104

= （a1（1） * a13 + a1（2） * a23 + a1（3） * a33） * b3（o2） = （0.1 * 0.3 + 0.16 * 0.2 + 0.28 * 0.5）* 0.3 = 0.202 * 0.3 = 0.0606

3）遞推計算 t = 3

= （a2（1） * a11 + a2（2） * a21 + a2（3） * a31） * b1（o3） = （0.077 * 0.5 + 0.1104 * 0.3 + 0.0606* 0.2）* 0.5 = 0.04187

= （a2（1） * a12 + a2（2） * a22 + a2（3） * a32） * b1（o3） = （0.077 * 0.2 + 0.1104 * 0.5 + 0.0606* 0.3）* 0.4 = 0.035512

= （a2（1） * a13 + a2（2） * a23 + a2（3） * a33） * b1（o3） = （0.077 * 0.3 + 0.1104 * 0.2 + 0.0606* 0.5）* 0.7 = 0.052836

可以看到，t = 3的時刻復用了 t = 2的結果。

4）終止

= a3（1） + a3（2）+ a3（3） = 0.04187 + 0.035512 + 0.052836 = 0.130218

最后一步終止只要對對T-1時刻的結果進行累加即可。條件概率可以通過邊緣概率累加計算得到。

3. 后向算法 - 一種遞推算法

先來定義后向概率，給定隱馬爾可夫模型，定義在時刻 t 狀態為的條件下，從 t + 1 到 T 的部分觀測序列為的概率為后向概率，記作：，可以用遞推的方法求得后向概率及觀測序列概率。算法流程如下

（1）初始化后向概率：，對最終時刻的所有狀態規定

（2）后向概率的遞推公式：，為了計算在時刻 t 狀態為條件下時刻 t + 1之后的觀測序列為的后向概率，只需考慮在時刻 t + 1所有可能的N個狀態的轉移概率（即項），以及在此狀態下的觀測的觀測概率（即項），然后考慮狀態之后的觀測序列的后向概率（即項）（遞推下去）

 （3）最終結果：，可以看到后向概率的最終結果表達式就是前向概率的初始化表達式

筆者思考：后向算法就是前向算法的遞歸倒推的改寫，很像C語言中的for循環和遞歸的關系，核心原理都是逐步遞推循環。

前向概率和后向概率可以統一起來，利用前向概率和后向概率的定義可以將觀測序列概率統一寫成：

Relevant Link:

《統計機器學習》李航

 

5. 隱馬爾柯夫模型學習算法 - 模型參數最大似然估計

0x1: 監督學習方法

如果訓練數據包括觀測序列和對應的狀態序列，則可以使用監督學習方法進行隱馬爾可夫模型的訓練學習。

假設訓練數據包含 S 個長度相同的觀測序列和對應的狀態序列，那么可以利用極大似然估計法來估計隱馬爾科夫模型的參數，具體流程如下

1. 轉移概率矩陣的估計

設樣本中時刻 t 處於狀態 i，時刻 t + 1 轉移到狀態 j 的頻數為，那么狀態轉移概率的最大似然估計公式是：

2. 觀測概率矩陣的估計

設樣本中狀態為 j 並觀測為 k 的頻數是，那么狀態為 j 觀測為 k 的概率的最大似然估計是：

3. 初始狀態概率的估計

為 S 個樣本中初始狀態為的頻率

可以看到，最大似然估計本質上就是最簡單的頻數統計法，但是在實際項目中，大多數時候我們是拿不到隱狀態序列。所以接下來要討論的是無監督下的模型參數訓練。

0x2: 無監督學習方法 - Baum-Welch算法（EM算法）

假設給定訓練數據只包含 S 個長度為 T 的觀測序列，而沒有對應的狀態序列（即包含隱變量），目標是學習隱馬爾可夫模型的參數。

套用EM算法的框架，我們將觀測序列數據看作觀測數據O，狀態序列數據看作不可觀測的隱數據 I，那么隱馬爾可夫模型實際上是一個含有隱變量的概率模型：

模型參數的學習可以通過EM算法實現。

1. 確定完全數據的對數似然函數

所有觀測數據寫成，所有隱數據寫成，完整的聯合概率分布是，完全數據的對數似然函數是：

2. EM算法的 E 步：

求 Q 函數 ，其中，是隱馬爾可夫模型參數的本輪次當前估計值，是要極大化的隱馬爾可夫模型參數

又因為，

所以 Q 函數可以寫成：，

式中所有求和都是對所有訓練數據的序列總長度 T 進行的。

3. EM算法的 M 步

極大化本輪次 Q 函數求模型參數A，B，π

由於要極大化的參數在 E 步的式子中單獨出現在3個項中，所有只需對各項分別逐一極大化

（1）第一項：

注意到滿足約束條件（初始概率分布累加和為1），利用拉格朗日乘子法，寫出拉格朗日函數：

對其求偏導數並令其結果為0：，得：

左右兩邊同時對 i 求和得到：，帶入朗格朗日求導為零公式得：。可以看到，初始概率的最大似然估計就是傳統統計法

（2）第二項：

同初始概率分布一樣，轉移概率同樣也滿足累加和為1的約束條件：，通過拉格朗日乘子法可以求出：

（3）第三項：

同樣用拉格朗日乘子法，約束條件是：，求得：

Relevant Link:

 

6. 隱馬爾科夫隱狀態估計算法 - 隱狀態概率矩陣最大似然估計

我們在本章節討論隱馬爾可夫模型隱狀態序列估計的兩種算法：近似算法、維特比算法（viterbi algorithm）

0x1: 近似算法

近似算法的思想是，在每個時刻 t 選擇在該時刻最可能出現的狀態，從而得到一個狀態序列，將它作為預測的結果。給定隱馬爾可夫模型和觀測序列 O，在時刻 t 處於狀態的概率是：，在每一時刻 t 最有可能的狀態是：，從而得到狀態序列

近似算法的優點是計算簡單，其缺點是不能保證預測的狀態序列整體是最優可能的狀態序列（即不能保證全局最優），因為預測的狀態序列可能有實際不發生的部分。同時，根據近似算法得到的狀態序列中有可能存在轉移概率為0的相鄰狀態（對某些i，j，aij = 0時）

0x2: 維特比算法

維特比算法實際是用動態規划解隱馬爾可夫模型預測問題，即用動態規划（dynamic programming）求概率最大路徑（最優路徑）。這時一條路徑對應着一個狀態序列。

根據動態規划原理，最優路徑具有這樣的特征：如果最優路徑在時刻 t 通過，那么這一路徑從結點到終點的部分路徑，對於從到終點的所有可能的部分路徑集合來說，必須是最優的（即每步都要求局部最優）

依據這一原理，我們只需從時刻 t = 1開始，遞推地計算在時刻 t，狀態為 i 的各條部分路徑的最大概率，在每一個 t = T計算最大概率對應的路徑。最終得到最優路徑序列，這就是維特比算法。

為了明確定義維特比算法，我們首先導入兩個變量，（最大概率）和（最大概率對應的狀態）。

定義在時刻 t 狀態為 i 的所有單個單個路徑中概率最大值為：，由此可得變量的遞推公示：

定義在時刻 t 狀態為 i 的所有單個路徑中概率最大的路徑的第 t - 1個結點為：

維特比算法的流程如下：

1. 初始化

2. 遞推：對t = 2，3，...，T

3. 終止

 

4. 最優路徑回溯：對t = T - 1，T - 2，....，1

得到最優路徑：

0x3: 通過一個例子來說明維特比算法

還是盒子球的場景，模型，

，

已知觀測序列，試求最優狀態序列，即最優路徑

我們應用維特比算法來解決這個問題

（1）初始化

在 t = 1時，對每一個狀態 i，i = 1，2，3，求狀態為 i 觀測為紅的概率，記此概率為，則

代入實際數據得：

（2）遞推 t = 2

在時刻 t = 2，對每個狀態i，i = 1，2，3，求在 t = 1時狀態為 j 觀測為紅並在 t = 2時狀態為 i 觀測為白的路徑的最大概率，記此最大概率為，則

同時，對每個狀態i，i = 1，2，3，記錄概率最大路徑的前一個狀態 j：

計算：

考慮上一輪所有可能狀態，本輪的出現狀態1的最大概率： = max｛0.1 * 0.5，0.16 * 0.3，0.28 * 0.2｝* 0.5 = 0.028。

考慮上一輪所有可能狀態，本輪的出現狀態2的最大概率： = max｛0.1 * 0.2，0.16 * 0.5，0.28 * 0.3｝* 0.6 = 0.0504。

考慮上一輪所有可能狀態，本輪的出現狀態3的最大概率： = max｛0.1 * 0.3，0.16 * 0.2，0.28 * 0.5｝* 0.7 = 0.098。

（3）遞推 t = 3

 同樣，在t = 3步，還是逐步引入t = 2步的所有概率最大概率，並結合t = 3步的轉移概率和觀測概率

 

計算：

（4）終止

以表示最優路徑的概率，則，該結果是逐步遞推出來的

最優路徑的終點是：

（5）由最優路徑的終點，逆向逆推最優路徑序列

本輪最大概率的計算公式中，包含了上一輪次最大概率狀態的結果

於是得到最優路徑序列，即最優狀態序列

筆者思考：

無論是計算觀測序列概率的前向-后向算法，還是計算隱狀態序列的維特比算法，本質上都是基於頻數統計的極大似然估計，所依據的理論都是貝葉斯邊緣概率累加公式。

同時我們也看到，每一步的計算需要考慮上一步的所有可能，這包含了一個思想，雖然上一步得到了局部最優但不一定就是全局最優，所以我們還不能完全“信任”上一步的argmax結果
當前本輪遞推要把上一輪的所有可能都參與到本輪的計算中，在本輪重新計算一次似然概率。所以 t+1 步才確定 t 步的最大似然狀態，每一輪的最大似然狀態都要到下一輪才能得到確定。
在更高階的N階馬爾科夫模型中，這個回溯要追溯到更早的 N 步前。

Relevant Link:

 

7. 隱馬爾科夫模型的應用

在運用馬爾科夫模型解決實際的問題的時候，我認為最重要的是正確理解我們要解決問題的本質，進行有效的抽象，將問題的各個維度套用到馬爾科夫模型的三要素上。這相當於數學建模過程，完成建模之后就可以利用馬爾科夫的學習和預測過程來幫助我們解決問題

0x1：序列標注問題

假設有4個盒子，每個盒子里都裝有紅白兩種顏色的球，盒子里的紅白球數由下表給出（等價於已知觀測概率矩陣）

游戲的主辦方按照如下的方法抽球，抽球的過程不讓觀察者看到，產生一個由球的顏色組成的觀測序列：

1. 開始，從4個盒子里隨機選取1個盒子（1/4等概率），從這個盒子里隨機抽出1個球，記錄其顏色后，放回；
2. 然后，從當前盒子隨機轉移到下一個盒子，規則是
　　1）如果當前盒子是盒子1，那么下一個盒子一定是2
　　2）如果當前盒子是2或3，那么分別以概率0.4和0.6轉移到左邊或右邊的盒子
　　3）如果當前盒子是盒子4，那么各以0.5的概率停留在盒子4或者轉移到盒子3
3. 確定轉移的盒子后，再從這個盒子里隨機抽出1個球，記錄其顏色，放回；
4. 如此重復下去，重復進行5次

得到一個球的顏色的觀測序列： O = ｛紅，紅，白，白，紅｝
在這個過程中，觀察者只能觀測到球的顏色的序列，觀測不到球是從哪個盒子取出的，即觀測不到每次實驗抽到盒子的隱狀態序列。

這是一個典型的馬爾科夫模型場景，根據所給條件，可以明確狀態集合、觀測集合、序列長度以及模型的三要素

盒子對應隱狀態，隱狀態的集合是：

球的顏色對應觀測，觀測的集合是：

狀態序列和觀測序列長度 T = 3

初始概率分布為：

狀態轉移概率分布為：

觀測概率分布為：

觀測集合V為：O = ｛紅，紅，白｝

要求解的問題是：在這次實驗中，每次實驗的盒子序列是什么。

0x2：根據觀測數據推測最大似然的隱狀態序列（序列解碼問題）

值得注意的是，在進行實際問題解決的時候，很重要的一點是抽象能力，就是要想清楚哪個數據對應的是觀測數據，哪個數據又是對應的隱狀態序列。

例如下面的根據觀測對象的日常行為，推測當前的天氣的具體問題。

我們對一個觀測對象進行6天的持續觀測，該對象每天的行為有3種：walk：散步；shop：購物；clean：打掃衛生。

同時每天的天氣有2種狀態：Rainy：下雨；Sunny：晴天。

我們已知每天天氣的轉換概率，以及在不同天氣下觀測對象做不同行為的概率分布。

最后我們的目標是根據一系列的對觀測對象在這6天的行為觀測結果：[0, 2, 1, 1, 0, 2]，即[walk，clean，shop，shop，walk，clean]，反推出這6天最可能的每天的天氣。

我們來編碼實現一下

# -*- coding:utf-8 -*-

from __future__ import division
import numpy as np
from hmmlearn import hmm


states = ["Rainy", "Sunny"]   ## 隱藏狀態類型
n_states = len(states)        ## 隱藏狀態數量

observations = ["walk", "shop", "clean"]  ## 可觀察的狀態類型
n_observations = len(observations)        ## 可觀察序列的狀態數量

start_probability = np.array([0.6, 0.4])  ## 初始狀態概率

## 狀態轉移概率矩陣
transition_probability = np.array([
  [0.7, 0.3],   # 晴天傾向保持晴天
  [0.4, 0.6]    # 陰天傾向轉為晴天
])

## 觀測序列概率矩陣
emission_probability = np.array([
  [0.1, 0.4, 0.5],
  [0.6, 0.3, 0.1]
])

# 構建了一個MultinomialHMM模型，這模型包括開始的轉移概率，隱狀態間的轉移矩陣A（transmat），隱狀態到觀測序列的發射矩陣emissionprob，下面是參數初始化
model = hmm.MultinomialHMM(n_components=n_states)
model.startprob_ = start_probability
model.transmat_ = transition_probability
model.emissionprob_ = emission_probability

# predict a sequence of hidden states based on visible states
bob_says = np.array([[0, 2, 1, 1, 0, 2]]).T  ##預測時的可見序列
print bob_says.shape
print bob_says
logprob, alice_hears = model.decode(bob_says, algorithm="viterbi")
print logprob ##該參數反映模型擬合的好壞

## 最后輸出decode解碼結果
print "Bob says(觀測行為序列):", ", ".join(map(lambda x: observations[x[0]], bob_says))
print "Alice hears(天氣序列)(隱狀態序列):", ", ".join(map(lambda x: states[x], alice_hears))

補充一句，如果是在語音識別問題中，觀測數據就是音頻流，需要估計的隱狀態序列就是對應的發音文字

Relevant Link:

http://www.cs.ubc.ca/~murphyk/Bayes/rabiner.pdf
http://www.cnblogs.com/pinard/p/7001397.html
http://blog.csdn.net/u014365862/article/details/50412978
http://blog.csdn.net/yywan1314520/article/details/50533850
http://hmmlearn.readthedocs.io/en/latest/api.html
http://blog.csdn.net/wowotuo/article/details/40783357

0x3：學習問題和概率計算問題搭配使用 - 先根據觀測數據（訓練樣本）train出一個模型（包含模型參數），然后基於這個模型預測新的樣本（觀測變量）的出現概率（即異常判斷）

import numpy as np
from hmmlearn.hmm import GaussianHMM


###############################################################################
# Working with multiple sequences
X1 = [[0.5], [1.0], [-1.0], [0.42], [0.24]]
X2 = [[2.4], [4.2], [0.5], [-0.24]]

# concatenate them into a single array and then compute an array of sequence lengths:
X = np.concatenate([X1, X2])
lengths = [len(X1), len(X2)]

###############################################################################
# Run Gaussian HMM
print "fitting to HMM ..."

# Make an HMM instance and execute fit
# There is 4 type of states
model = GaussianHMM(n_components=4, n_iter=100).fit(X)

# Train the model parameters and hidden state
model_pre = model.fit(X, lengths)

###############################################################################
print "Compute the observation sequence log probability under the model parameters"
print model_pre.score(X1)

print "Compute the log probability under the model and compute posteriors."
print model_pre.score_samples(X1)

我們輸入訓練樣本得到一組模型參數，之后基於這個模型參數預測接下來輸入的觀測序列出現的可能性，這常被用在文本異常檢測的場景中。score方法返回的是LOG的結果，所以是負的，越接近0，表示越符合當前的HMM模型。

模型預測的觀測序列的對數似然概率為：

-2.64250783257

0x4：學習問題和解碼問題搭配使用 - 先根據觀測數據（訓練樣本）train出一個模型（包含模型參數），然后基於這個模型預測在給定觀測序列的情況下，最大似然的隱變量序列（及語音解碼識別、解碼問題）

import numpy as np
from hmmlearn.hmm import GaussianHMM


###############################################################################
# Working with multiple sequences
X1 = [[0.5], [1.0], [-1.0], [0.42], [0.24]]
X2 = [[2.4], [4.2], [0.5], [-0.24]]

# concatenate them into a single array and then compute an array of sequence lengths:
X = np.concatenate([X1, X2])
lengths = [len(X1), len(X2)]

###############################################################################
# Run Gaussian HMM
print "fitting to HMM ..."

# Make an HMM instance and execute fit
# There is 4 type of states
model = GaussianHMM(n_components=4, n_iter=100).fit(X)

# Train the model parameters and hidden state
model_pre = model.fit(X, lengths)


###############################################################################
print "Find most likely state sequence corresponding to X."
print model_pre.decode(X1)

print "pridect the mostly likly hidden state sequence according to the sample"
print model_pre.predict(X1)

print "Compute the posterior probability for each state in the model."
print model_pre.predict_proba(X1)

 

模型預測的隱變量序列概率分布為：

[[  4.28596850e-013   0.00000000e+000   1.00000000e+000   4.82112557e-178]
 [  1.00000000e+000   1.82778069e-311   3.61988083e-021   0.00000000e+000]
 [  1.00000000e+000   0.00000000e+000   7.31100741e-098   0.00000000e+000]
 [  7.17255159e-002   0.00000000e+000   9.28274484e-001   0.00000000e+000]
 [  9.97881894e-001   0.00000000e+000   2.11810616e-003   0.00000000e+000]]

模型預測的隱變量序列的結果為，其實我們根據概率分布取argmax也可以得到和api返回相同的結果

[2 0 0 2 0]

Relevant Link:

http://hmmlearn.readthedocs.io/en/stable/api.html
http://hmmlearn.readthedocs.io/en/stable/tutorial.html 

0x5：基於HMM進行文本異常檢測

1. 以白找黑

如果我們能定義出一個場景，即正常的情況下文本的概率空間是集中在一個相對固定的范圍內的（在一定的字符空間中進行狀態轉移），而異常的樣本的字符空間及其字符間的轉換關系極其不符合

正常的“模式”，這種情況就可以使用HMM。

# -*- coding: utf-8 -*-

import urllib
import urlparse
import re
from hmmlearn import hmm
import numpy as np
from sklearn.externals import joblib


MIN_LEN = 6     # 處理參數值的最小長度
N_Hidden = 10    # 隱狀態個數（從領域經驗來看，隱狀態數量和觀測變量的狀態數量基本上應該持平）
T = -30  # 最大似然概率閾值
# URL中會出現的特殊字符
SEN = [
    '<',
    '>',
    ',',
    ':',
    '\'',
    '/',
    ';',
    '"',
    '{',
    '}',
    '(',
    ')'
]


# 排除中文干擾 只處理127以內的字符
def isascii(istr):
    if re.match(r'^(http)', istr):
        return False
    for i, c in enumerate(istr):
        if ord(c) > 127 or ord(c) < 31:
            return False
        if c in SEN:
            return True
    return True


# 特征離散規范化，減少參數狀態空間
def etl(istr):
    vers = []
    for i, c in enumerate(istr):
        c = c.lower()
        if ord(c) >= ord('a') and ord(c) <= ord('z'):   # ascii字母
            vers.append([ord(c)])
        elif ord(c) >= ord('0') and  ord(c) <= ord('9'):    # 數字
            vers.append([1])
        elif c in SEN:
            vers.append([ord(c)])
        else:
            vers.append([2])
    return np.array(vers)


def extractvec(filename):
    X = [[0]]
    LINES = [['']]
    X_lens = [1]
    with open(filename) as f:
        for line in f:
            query = urllib.unquote(line.strip('\n'))  # url解碼
            if len(line) >= MIN_LEN:
                params = urlparse.parse_qsl(query, True)
                for k, v in params:
                    if isascii(v) and len(v) >= N_Hidden:
                        vers = etl(v)  # 數據處理與特征提取
                        X = np.concatenate([X, vers])  # 每一個參數value作為一個特征向量
                        # # 通過np.concatenate整合成了一個1D長向量，同時需要額外傳入len list來標明每個序列的長度邊界
                        X_lens.append(len(vers))  # 長度
                        LINES.append(v)
    return X, X_lens, LINES


def train(filename):
    X, X_lens, LINES = extractvec(filename)
    remodel = hmm.GaussianHMM(n_components=N_Hidden, covariance_type="full", n_iter=100)
    remodel.fit(X, X_lens)
    joblib.dump(remodel, "../model/xss-train.pkl")
    return remodel


def predict(filename):
    remodel = joblib.load("../model/xss-train.pkl")
    X, X_lens, LINES = extractvec(filename)
    LINES = LINES[1:]
    X = X[1:, :]
    X_lens = X_lens[1:]
    lastindex = 0
    for i in range(len(X_lens)):
        line = np.array(X[lastindex:lastindex+X_lens[i], :])
        lastindex += X_lens[i] 
        pros = remodel.score(line)
        if pros < T:
            print "Find bad sample POR: %d LINE:%s" % (pros, LINES[i])


if __name__ == '__main__':
    # 以白找黑，給HMM輸入純白樣本，讓其記住正常url參數的模型轉化概率
    # remodel = train("../data/web-attack/normal-10000.txt")    # train a GMM model with train sample

    # 輸入待檢測樣本，設定一個閾值(score的分值越小，說明出現的概率越小，也即說明偏離正常樣本的程度)
    predict("../data/xss-20.txt")
    # predict the probability of the new sample according to the model(abnormal detection)

我們將黑白待檢測樣本混合在一起，看HMM的預測效果如何

2. 以黑找黑

# -*- coding:utf-8 -*-

import sys
import urllib
import urlparse
import re
from hmmlearn import hmm
import numpy as np
from sklearn.externals import joblib
import HTMLParser
import nltk

MIN_LEN = 10    # 處理參數值的最小長度
N = 5   # 狀態個數
T = -200    # 最大似然概率閾值
SEN = ['<', '>', ',', ':', '\'', '/', ';', '"', '{', '}', '(', ')']  #其他字符2
index_wordbag = 1     # 詞袋索引
wordbag = {}    # 詞袋


# 數據提取與特征提取，這一步我們不采用單字符的char特征提取，而是根據領域經驗對特定的phrase字符組為基礎單位，進行特征化，這是一種token切分方案

# </script><script>alert(String.fromCharCode(88,83,83))</script>
# <IMG SRC=x onchange="alert(String.fromCharCode(88,83,83))">
# <;IFRAME SRC=http://ha.ckers.org/scriptlet.html <;
# ';alert(String.fromCharCode(88,83,83))//\';alert(String.fromCharCode(88,83,83))//";alert(String.fromCharCode(88,83,83))
# //\";alert(String.fromCharCode(88,83,83))//--></SCRIPT>">'><SCRIPT>alert(String.fromCharCode(88,83,83))</SCRIPT>
tokens_pattern = r'''(?x)
 "[^"]+"
|http://\S+
|</\w+>
|<\w+>
|<\w+
|\w+=
|>
|\w+\([^<]+\) #函數 比如alert(String.fromCharCode(88,83,83))
|\w+
'''


def split_tokens(line, tokens_pattern):
    tokens = nltk.regexp_tokenize(line, tokens_pattern)
    return tokens

def load_wordbag(filename, max=100):
    X = [[0]]
    X_lens = [1]
    tokens_list = []
    global wordbag
    global index_wordbag

    with open(filename) as f:
        for line in f:
            line = line.strip('\n')
            line = urllib.unquote(line)   # url解碼
            #h = HTMLParser.HTMLParser()    # 處理html轉義字符
            #line = h.unescape(line)
            if len(line) >= MIN_LEN:  # 忽略短url value
                line, number = re.subn(r'\d+', "8", line)   # 數字常量替換成8
                line, number = re.subn(r'(http|https)://[a-zA-Z0-9\.@&/#!#\?:=]+', "http://u", line)    # ulr日換成http://u
                line, number = re.subn(r'\/\*.?\*\/', "", line)     # 去除注釋
                tokens = split_tokens(line, tokens_pattern)  # token切分
                tokens_list += tokens

    fredist = nltk.FreqDist(tokens_list)  # 單文件詞頻
    keys = fredist.keys()
    keys = keys[:max]     # 降維，只提取前N個高頻使用的單詞，其余規范化到0
    for localkey in keys:  # 獲取統計后的不重復詞集
        if localkey in wordbag.keys():  # 判斷該詞是否已在詞袋中
            continue
        else:
            wordbag[localkey] = index_wordbag
            index_wordbag += 1
    print "GET wordbag", wordbag
    print "GET wordbag size(%d)" % index_wordbag


def train(filename):
    X = [[-1]]
    X_lens = [1]
    global wordbag
    global index_wordbag

    with open(filename) as f:
        for line in f:
            line = line.strip('\n')
            line = urllib.unquote(line)   # url解碼
            #h = HTMLParser.HTMLParser() # 處理html轉義字符
            #line = h.unescape(line)
            if len(line) >= MIN_LEN:
                line, number = re.subn(r'\d+', "8", line)   # 數字常量替換成8
                line, number = re.subn(r'(http|https)://[a-zA-Z0-9\.@&/#!#\?:]+', "http://u", line)     # ulr日換成http://u
                line, number = re.subn(r'\/\*.?\*\/', "", line)     # 干掉注釋
                words = split_tokens(line, tokens_pattern)
                vers = []
                for word in words:
                    # 根據詞匯編碼表進行index編碼，對不在詞匯表中的token詞不予編碼
                    if word in wordbag.keys():
                        vers.append([wordbag[word]])
                    else:
                        vers.append([-1])

            np_vers = np.array(vers)
            X = np.concatenate([X, np_vers])
            X_lens.append(len(np_vers))

    remodel = hmm.GaussianHMM(n_components=N, covariance_type="full", n_iter=100)
    remodel.fit(X, X_lens)
    joblib.dump(remodel, "../model/find_black_with_black-xss-train.pkl")

    return remodel


def test(remodel, filename):
    with open(filename) as f:
        for line in f:
            line = line.strip('\n')
            line = urllib.unquote(line)     # url解碼
            # 處理html轉義字符
            #h = HTMLParser.HTMLParser()
            #line = h.unescape(line)

            if len(line) >= MIN_LEN:
                line, number = re.subn(r'\d+', "8", line)      # 數字常量替換成8
                line, number = re.subn(r'(http|https)://[a-zA-Z0-9\.@&/#!#\?:]+', "http://u", line)     # ulr日換成http://u
                line, number = re.subn(r'\/\*.?\*\/', "", line)     # 干掉注釋
                words = split_tokens(line)
                vers = []
                for word in words:
                    # test和train保持相同的編碼方式
                    if word in wordbag.keys():
                        vers.append([wordbag[word]])
                    else:
                        vers.append([-1])

                np_vers = np.array(vers)
                pro = remodel.score(np_vers)
                if pro >= T:
                    print "SCORE:(%d) XSS_URL:(%s) " % (pro, line)


if __name__ == '__main__':
    train_filename = "../data/xss-200000.txt"
    # 基於訓練樣本集（純黑，我們這里是以黑找黑），得到詞頻編碼表
    load_wordbag(train_filename, 2000)

    # 基於同樣的訓練樣本集（純黑，我們這里是以黑找黑），訓練HMM模型
    remodel = train(train_filename)

    # 基於訓練得到的HMM模型，對未知樣本的最大似然概率進行預測（即預測未知樣本觀測序列出現的可能性，即找類似的黑樣本）
    test(remodel, "../data/xss-20.txt")

整個過程分為：

泛化 -> 分詞 -> 詞集編碼 -> 詞集模型處理流程

0x6：將HMM的對數似然概率預測值作為特征

我們可以采用類似於DNN中隱層的思想，將HMM的對數似然概率輸出作為一個抽象后的特征值，甚至可以作為新的一份樣本特征輸入到下一層的算法模型中參與訓練。

# -*- coding:utf-8 -*-

import sys
import urllib
import urlparse
import re
from hmmlearn import hmm
import numpy as np
from sklearn.externals import joblib
import HTMLParser
import nltk
import csv
import matplotlib.pyplot as plt

MIN_LEN = 10    # 處理域名的最小長度
N = 8     # 狀態個數
T = -50   # 最大似然概率閾值
FILE_MODEL = "../model/12-4.m"     # 模型文件名


def load_alexa(filename):
    domain_list = []
    csv_reader = csv.reader(open(filename))
    for row in csv_reader:
        domain = row[1]
        if domain >= MIN_LEN:
            domain_list.append(domain)
    return domain_list


def domain2ver(domain):
    ver = []
    for i in range(0, len(domain)):
        ver.append([ord(domain[i])])
    return ver

def train_hmm(domain_list):
    X = [[0]]
    X_lens = [1]
    for domain in domain_list:
        ver = domain2ver(domain)    # 逐字符ascii向量化
        np_ver = np.array(ver)
        X = np.concatenate([X, np_ver])
        X_lens.append(len(np_ver))

    remodel = hmm.GaussianHMM(n_components=N, covariance_type="full", n_iter=100)
    remodel.fit(X, X_lens)
    joblib.dump(remodel, FILE_MODEL)

    return remodel


def load_dga(filename):
    domain_list = []
    # xsxqeadsbgvpdke.co.uk,Domain used by Cryptolocker - Flashback DGA for 13 Apr 2017,2017-04-13,
    # http://osint.bambenekconsulting.com/manual/cl.txt
    with open(filename) as f:
        for line in f:
            domain = line.split(",")[0]
            if domain >= MIN_LEN:
                domain_list.append(domain)
    return domain_list


def test_dga(remodel,filename):
    x = []
    y = []
    dga_cryptolocke_list = load_dga(filename)
    for domain in dga_cryptolocke_list:
        domain_ver = domain2ver(domain)
        np_ver = np.array(domain_ver)
        pro = remodel.score(np_ver)
        x.append(len(domain))
        y.append(pro)
    return x, y


def test_alexa(remodel,filename):
    x=[]
    y=[]
    alexa_list = load_alexa(filename)
    for domain in alexa_list:
        domain_ver = domain2ver(domain)
        np_ver = np.array(domain_ver)
        pro = remodel.score(np_ver)
        x.append(len(domain))
        y.append(pro)
    return x, y


if __name__ == '__main__':
    domain_list = load_alexa("../data/top-1000.csv")
    # remodel = train_hmm(domain_list)    # 以白找黑: 基於alexa域名數據訓練hmm模型
    remodel = joblib.load(FILE_MODEL)
    x_3, y_3 = test_dga(remodel, "../data/dga-post-tovar-goz-1000.txt")
    x_2, y_2 = test_dga(remodel,"../data/dga-cryptolocke-1000.txt")
    fig, ax = plt.subplots()
    ax.set_xlabel('Domain Length')      # 橫軸是域名長度
    ax.set_ylabel('HMM Score')          # 縱軸是hmm對觀測序列計算得到的似然概率
    ax.scatter(x_3, y_3, color='b', label="dga_post-tovar-goz")
    ax.scatter(x_2, y_2, color='g', label="dga_cryptolock") 
    ax.legend(loc='right')
    plt.show()

我們以域名長度為橫軸，以HMM值作為縱軸，從可視化結果可以看到，HMM作為DGA域名區分的一個變量，有一定的區分型，展現出了一定的規律。