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域名區分的一個變量，有必定的區分型，展示出了必定的規律。