[Scikit-learn] Dynamic Bayesian Network - HMM

Warning The sklearn.hmm module has now been deprecated due to it no longer matching the scope and the API of the project. It is scheduled for removal in the 0.17 release of the project.

From: http://scikit-learn.sourceforge.net/stable/modules/hmm.html

The sklearn.hmm module has now been deprecated due to it no longer matching the scope and the API of the project. It is scheduled for removal in the 0.17 release of the project.

 

hmmlean

This module has been moved to a seperate repository: https://github.com/hmmlearn/hmmlearn

hmmlearn doc: http://hmmlearn.readthedocs.io/en/latest/

 

其餘參考連接：

隱馬爾科夫模型HMM的前向算法和後向算法

HMM的Baum-Welch算法和Viterbi算法公式推導細節

涉及的三個問題都是什麼，以及解決每一個問題的主流算法：

1. 機率計算問題即模型評價問題——前向算法和後向算法
2. 學習問題即參數估計問題——Baum-Welch算法
3. 預測問題即解碼問題——Viterbi算法

---

 

Ref：https://www.zhihu.com/question/20962240/answer/33438846

擲骰子

仍是用最經典的例子，擲骰子。假設我手裏有三個不一樣的骰子。

 

 

可見/隱含狀態鏈

假設咱們開始擲骰子，咱們先從三個骰子裏挑一個，挑到每個骰子的機率都是1/3。

而後咱們擲骰子，獲得一個數字，1，2，3，4，5，6，7，8中的一個。

 

不停的重複上述過程，咱們會獲得一串數字，每一個數字都是1，2，3，4，5，6，7，8中的一個。

例如咱們可能獲得這麼一串數字（擲骰子10次）：1 6 3 5 2 7 3 5 2 4

這串數字叫作 可見狀態鏈。

 

 

可是在隱馬爾可夫模型中，咱們不只僅有這麼一串可見狀態鏈，還有一串 隱含狀態鏈。在這個例子裏，這串隱含狀態鏈就是你用的骰子的序列。

好比，隱含狀態鏈有多是： D6 -> D8 -> D8 -> D6 -> D4 -> D8 -> D6 -> D6 -> D4 -> D8

 

轉換機率/輸出機率

隱含狀態（骰子）之間存在 轉換機率（ Transition Probability）。

 

在咱們這個例子裏，D6的下一個狀態是D4，D6，D8的機率都是1/3。

D4，D8的下一個狀態是D4，D6，D8的轉換機率也都同樣是1/3。

這樣設定是爲了最開始容易說清楚，可是咱們實際上是能夠隨意設定轉換機率的。好比，咱們能夠這樣定義，D6後面不能接D4，D6後面是D6的機率是0.9，是D8的機率是0.1。這樣就是一個新的HMM。

一樣的， 儘管可見狀態之間沒有轉換機率，

 

可是隱含狀態和可見狀態之間有一個機率叫做  輸出機率（emission probability）。

就咱們的例子來講，六面骰（D6）產生1的輸出機率是1/6。產生2，3，4，5，6的機率也都是1/6。咱們一樣能夠對輸出機率進行其餘定義。好比，我有一個被賭場動過手腳的六面骰子，擲出來是1的機率更大，是1/2，擲出來是2，3，4，5，6的機率是1/10。

 

提出問題

 

其實對於HMM來講，若是提早知道全部隱含狀態之間的轉換機率和全部隱含狀態到全部可見狀態之間的輸出機率，作模擬是至關容易的。

可是應用HMM模型時候呢，每每是缺失了一部分信息的，

• • 有時候你知道骰子有幾種，每種骰子是什麼，可是不知道擲出來的骰子序列；
  • 有時候你只是看到了不少次擲骰子的結果，剩下的什麼都不知道。

若是應用算法去估計這些缺失的信息，就成了一個很重要的問題。這些算法我會在下面詳細講。


回到正題，和HMM模型相關的算法主要分爲三類，分別解決三種問題：

1）知道骰子有幾種（隱含狀態數量），每種骰子是什麼（轉換機率），根據擲骰子擲出的結果（可見狀態鏈），我想知道每次擲出來的都是哪一種骰子（隱含狀態鏈）。
這個問題呢，在語音識別領域呢，叫作 解碼問題。

這個問題其實有兩種解法，會給出兩個不一樣的答案。每一個答案都對，只不過這些答案的意義不同。

• 第一種解法求最大似然狀態路徑，說通俗點呢，就是我求一串骰子序列，這串骰子序列產生觀測結果的機率最大。（這裏只講這個）
• 第二種解法呢，就不是求一組骰子序列了，而是求每次擲出的骰子分別是某種骰子的機率。好比說我看到結果後，我能夠求得第一次擲骰子是D4的機率是0.5，D6的機率是0.3，D8的機率是0.2。

2）仍是知道骰子有幾種 （隱含狀態數量） ，每種骰子是什麼 （轉換機率） ，根據擲骰子擲出的結果 （可見狀態鏈） ，我想知道擲出這個結果的機率。
看似這個問題意義不大，由於你擲出來的結果不少時候都對應了一個比較大的機率。問這個問題的目的呢，實際上是檢測觀察到的結果和已知的模型是否吻合。

若是不少次結果都對應了比較小的機率，那麼就說明咱們已知的模型頗有多是錯的，有人偷偷把咱們的骰子給換了。

3）知道骰子有幾種 （隱含狀態數量） ，不知道每種骰子是什麼 （轉換機率） ，觀測到不少次擲骰子的結果 （可見狀態鏈） ，我想反推出每種骰子是什麼 （轉換機率） 。
這個問題很重要，由於這是最多見的狀況。不少時候咱們只有可見結果，不知道HMM模型裏的參數，咱們須要從可見結果估計出這些參數，這是建模的一個必要步驟。

 

問題闡述完了，下面就開始說解法。（0號問題在上面沒有提，只是做爲解決上述問題的一個輔助）

---

0. 一個簡單問題

「知道骰子有幾種，每種骰子是什麼，每次擲的都是什麼骰子，根據擲骰子擲出的結果，求產生這個結果的機率。」

 

解法無非就是機率相乘：

 

【已知信息完整，求可見狀態鏈的機率】

 

 

1. 看見不可見的，破解骰子序列

 

這裏我說的是第一種解法，解 最大似然路徑問題。(brute force固然不行)

舉例來講，我知道我有三個骰子，六面骰，四面骰，八面骰。我也知道我擲了十次的結果（1 6 3 5 2 7 3 5 2 4），我不知道每次用了那種骰子，我想知道最有可能的骰子序列。
其實最簡單而暴力的方法就是窮舉全部可能的骰子序列，而後 依照第零個問題的解法把每一個序列對應的機率算出來。而後咱們從裏面把對應最大機率的序列挑出來就好了。若是馬爾可夫鏈不長，固然可行。若是長的話，窮舉的數量太大，就很難完成了。

另一種頗有名的算法叫作 Viterbi algorithm. 要理解這個算法，咱們先看幾個簡單的列子。

（1）首先，若是咱們只擲一次骰子：
　　
　　由於D4產生1的機率是1/4，高於1/6和1/8，so，對應的最大機率骰子序列就是D4，

（2）把這個狀況拓展，咱們擲兩次骰子：
　　
　　結果爲1，6.這時問題變得複雜起來，咱們要計算三個值，分別是第二個骰子是D6，D4，D8的最大機率。

顯然，要取到最大機率，第一個骰子必須爲D4。第二個骰子取D6：

 

一樣的，咱們能夠計算第二個骰子是D4或D8時的最大機率。咱們發現，第二個骰子取到D6的機率最大。

 

（3）繼續拓展，咱們擲三次骰子：
　　  

　　一樣，咱們計算第三個骰子分別是D6，D4，D8的最大機率。咱們再次發現，要取到最大機率，第二個骰子必須爲D6。這時，第三個骰子取到D4的最大機率是：

 

 

同上，咱們能夠計算第三個骰子是D6或D8時的最大機率。咱們發現，第三個骰子取到D4的機率最大。而使這個機率最大時，第二個骰子爲D6，第一個骰子爲D4。因此最大機率骰子序列就是D4 D6 D4。

寫到這裏，你們應該看出點規律了。既然擲骰子一二三次能夠算，擲多少次均可以以此類推。咱們發現，咱們要求最大機率骰子序列時要作這麼幾件事情。首先，無論序列多長，要從序列長度爲1算起，算序列長度爲1時取到每一個骰子的最大機率。而後，逐漸增長長度，每增長一次長度，從新算一遍在這個長度下最後一個位置取到每一個骰子的最大機率。由於上一個長度下的取到每一個骰子的最大機率都算過了，從新計算的話其實不難。當咱們算到最後一位時，就知道最後一位是哪一個骰子的機率最大了。而後，咱們要把對應 這個最大機率的序列從後往前推出來。

 

【但其實，狀態轉移時是消耗能量的，因此機率相對會小些 compared with 狀態不變化】

 

 

2. 誰動了個人骰子？

好比說你懷疑本身的六面骰被賭場動過手腳了，有可能被換成另外一種六面骰： 【異常色子的信息也已知】

這種六面骰擲出來是1的機率更大，是1/2，擲出來是2，3，4，5，6的機率是1/10。你怎麼辦麼？答案很簡單：

• 算一算正常的三個骰子擲出一段序列的機率，
• 再算一算不正常的六面骰和另外兩個正常骰子擲出這段序列的機率。

若是前者比後者小，你就要當心了。

好比說擲骰子的結果是：

 

要算用正常的三個骰子擲出這個結果的機率，其實就是將全部可能狀況的機率進行加和計算。一樣，簡單而暴力的方法就是把窮舉全部的骰子序列，仍是計算每一個骰子序列對應的機率，可是這回，咱們不挑最大值了，而是把全部算出來的機率相加，獲得的總機率就是咱們要求的結果。這個方法依然不能應用於太長的骰子序列（馬爾可夫鏈）。

咱們會應用一個和前一個問題相似的解法，只不過 前一個問題關心的是 機率最大值，這個問題關心的是 機率之和。解決這個問題的算法叫作 前向算法（forward algorithm）。

首先，若是咱們只擲一次骰子：

看到結果爲1。產生這個結果的總機率能夠按照以下計算，總機率爲0.18：

把這個狀況拓展，咱們擲兩次骰子：

看到結果爲1，6.產生這個結果的總機率能夠按照以下計算，總機率爲0.05：


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

 

看到結果爲1，6，3. 產生這個結果的總機率能夠按照以下計算，總機率爲0.03：

一樣的，咱們一步一步的算，有多長算多長，再長的馬爾可夫鏈總能算出來的。

用一樣的方法，也能夠算出不正常的六面骰和另外兩個正常骰子擲出這段序列的機率，

而後咱們 比較一下這兩個機率大小，就能知道你的骰子是否是被人換了。

 

【有點動態規劃的意思，不斷的利用子問題的解】

---

以上是一些基礎問題。

模型這東西，不結合實際features搞一搞，不利於深刻理解。

 

HMM的理解 in NLP

From: https://zhuanlan.zhihu.com/p/26811689

應用到詞性標註中，v表明詞語，是能夠觀察到的。q表明咱們要預測的詞性（一個詞可能對應多個詞性）是隱含狀態。

應用到分詞中，v表明詞語，是能夠觀察的。q表明咱們的標籤（B,E這些標籤，表明一個詞語的開始，或者中間等等）

應用到命名實體識別中，v表明詞語，是能夠觀察的。q表明咱們的標籤（標籤表明着地點詞，時間詞這些）

 

From: NLP - 05. Tagging Problems, and HMM

HMM 算是生成模型中的一種。

這裏用到的是Trigram Hidden Markov Models: 受最近的兩個單詞的影響，也就是分析三個單詞之間的關係。

 

 

寫成一個概括的形式，以下：

x_i: Term.

y_i: Named Entity Recognition.

 

HMM的三大問題

我愛天然語言處理 - HMM

Link: http://www.52nlp.cn/category/hidden-markov-model/page/4

 

Link 三大問題

1. Evaluation problem - 計算問題

   • Evaluation problem answers the question: what is the probability that a particular sequence of symbols is produced by a particular model?
   • For evaluation we use two algorithms: the forward algorithm or the backwards algorithm (DO NOT confuse them with the forward-backward algorithm).

2. Decoding problem - 解碼問題

   • Decoding problem answers the question: Given a sequence of symbols (your observations) and a model, what is the most likely sequence of states that produced the sequence.
   • For decoding we use the Viterbi algorithm.

3. Training problem - 學習問題

   • Training problem answers the question: Given a model structure and a set of sequences, find the model that best fits the data.
   • For this problem we can use the following 3 algorithms:
     1. MLE (maximum likelihood estimation)
     2. Viterbi training (DO NOT confuse with Viterbi decoding)
     3. Baum Welch = forward-backward algorithm

  

考慮一個村莊，全部村民都健康或發燒，只有村民醫生才能肯定每一個人是否發燒。醫生經過詢問患者的感覺來診斷髮燒。村民只能回答說他們以爲正常，頭暈或感冒。

醫生認爲，他的患者的健康情況做爲離散的馬可夫鏈。 「健康」和「發燒」有兩個狀態，但醫生不能直接觀察他們;健康與發燒的狀態是隱藏的。天天都有機會根據患者的健康情況，病人會告訴醫生他/她是「正常」，「感冒」仍是「頭昏眼花」。（正常，感冒，仍是暈眩是咱們前面說的觀測序列）

觀察（正常，感冒，暈眩）以及隱藏的狀態（健康，發燒）造成隱馬爾可夫模型（HMM），並能夠用Python編程語言表示以下：

obs    = ('normal', 'cold', 'dizzy') states = ('Healthy', 'Fever')
 start_p = {'Healthy': 0.6, 'Fever': 0.4}　　#醫生對患者首次訪問時HMM所處的狀態的信念（他知道患者每每是健康的） trans_p = { 'Healthy' : {'Healthy': 0.7, 'Fever': 0.3}, 'Fever' : {'Healthy': 0.4, 'Fever': 0.6} } emit_p = { 'Healthy' : {'normal': 0.5, 'cold': 0.4, 'dizzy': 0.1}, 'Fever' : {'normal': 0.1, 'cold': 0.3, 'dizzy': 0.6} }

那麼用圖來表示上面的例子能夠以下表示：

 

第一個問題 - 計算問題

求，給定模型的狀況下，求某種觀測序列出現的機率。

好比，給定的HMM模型參數已知，求出三天觀察是 (Dizzy,Cold,Normal) 的機率是多少？

對應的HMM模型參數已知的意思，就是說的A,B,pi矩陣是已經知道的。

 

第二個問題 - 學習問題

已知觀測序列是（Dizzy,Cold,Normal）,須要求出HMM的參數問題。

也就是咱們上面說的 A, B, PI 三個矩陣參數

重點，採用EM逼近。

Note, since the EM algorithm is a gradient-based optimization method, it will generally get stuck in local optima. You should in general try to run fit with various initializations and select the highest scored model.

You can use the monitor_ attribute to diagnose convergence.

 

第三個問題 - 解碼問題

咱們知道了觀測序列是（Dizzy,Cold,Normal）,也知道了HMM的參數，讓咱們求出形成這個觀測序列最有可能對應的狀態序列。

好比說是（Healthy，Healthy，Fever）仍是(Healthy,Healthy,Healthy)。

Jeff: 直觀的brute force方法，求出全部的組合，而後求機率，機率最大的序列就是解。

 

實踐出真知

""" Sampling from HMM ----------------- This script shows how to sample points from a Hiden Markov Model (HMM): we use a 4-components with specified mean and covariance. The plot show the sequence of observations generated with the transitions between them. We can see that, as specified by our transition matrix, there are no transition between component 1 and 3. """
print(__doc__) import numpy as np import matplotlib.pyplot as plt from hmmlearn import hmm ############################################################## # Prepare parameters for a 4-components HMM # Initial population probability
startprob = np.array([0.6, 0.3, 0.1, 0.0]) # The transition matrix, note that there are no transitions possible # between component 1 and 3
transmat = np.array([[0.7, 0.2, 0.0, 0.1], [0.3, 0.5, 0.2, 0.0], [0.0, 0.3, 0.5, 0.2], [0.2, 0.0, 0.2, 0.6]]) # The means of each component
means = np.array([[0.0,  0.0], [0.0, 11.0], [9.0, 10.0], [11.0, -1.0]]) # The covariance of each component
covars = .5 * np.tile(np.identity(2), (4, 1, 1)) # Build an HMM instance and set parameters
model = hmm.GaussianHMM(n_components=4, covariance_type="full") # Instead of fitting it from the data, we directly set the estimated # parameters, the means and covariance of the components
model.startprob_ = startprob model.transmat_ = transmat model.means_ = means model.covars_ = covars ###############################################################

# Generate samples
X, Z = model.sample(500) # Plot the sampled data
plt.plot(X[:, 0], X[:, 1], ".-", label="observations", ms=6, mfc="orange", alpha=0.7) # Indicate the component numbers
for i, m in enumerate(means): plt.text(m[0], m[1], 'Component %i' % (i + 1), size=17, horizontalalignment='center', bbox=dict(alpha=.7, facecolor='w')) plt.legend(loc='best') plt.show()

Result:  

 

Finance module的庫有問題： link

The finance module has been split out into its own package https://github.com/matplotlib/mpl_finance and deprecated in matplotlib. Any further development will take place in https://github.com/matplotlib/mpl_finance can you report the issue there? The finance module has been deprecated because non of the matplotlib developers are interested in maintaining it and mpl_finance is currently almost un maintained. It might also be worth checking out https://pypi.python.org/pypi/yahoo-finance as an alternative solution

""" Gaussian HMM of stock data -------------------------- This script shows how to use Gaussian HMM on stock price data from Yahoo! finance. For more information on how to visualize stock prices with matplotlib, please refer to ``date_demo1.py`` of matplotlib. """

from __future__ import print_function import datetime import numpy as np from matplotlib import cm, pyplot as plt from matplotlib.dates import YearLocator, MonthLocator try: from matplotlib.finance import quotes_historical_yahoo_ochl except ImportError: # For Matplotlib prior to 1.5.
    from matplotlib.finance import ( quotes_historical_yahoo as quotes_historical_yahoo_ochl ) from hmmlearn.hmm import GaussianHMM print(__doc__) ############################################################################### # Get quotes from Yahoo! finance
quotes = quotes_historical_yahoo_ochl("^GSPC", datetime.date(2012, 1, 1), datetime.date(2015, 1, 6)) sp = quotes_historical_yahoo_ochl("^GSPC", datetime.date(2012, 1, 1), datetime.date(2015, 1, 6), asobject=True, adjusted=True) # Unpack quotes
dates   = np.array([q[0] for q in quotes], dtype=int) close_v = np.array([q[2] for q in quotes]) volume = np.array([q[5] for q in quotes])[1:] # Take diff of close value. Note that this makes # ``len(diff) = len(close_t) - 1``, therefore, other quantities also # need to be shifted by 1.
diff    = np.diff(close_v) dates = dates[1:] close_v = close_v[1:] # Pack diff and volume for training.
X = np.column_stack([diff, volume]) ############################################################################### # Run Gaussian HMM
print("fitting to HMM and decoding ...", end="") # Make an HMM instance and execute fit
model = GaussianHMM(n_components=4, covariance_type="diag", n_iter=1000).fit(X) # Predict the optimal sequence of internal hidden state
hidden_states = model.predict(X) print("done") ############################################################################### # Print trained parameters and plot
print("Transition matrix") print(model.transmat_) print() print("Means and vars of each hidden state") for i in range(model.n_components): print("{0}th hidden state".format(i)) print("mean = ", model.means_[i]) print("var = ", np.diag(model.covars_[i])) print() fig, axs = plt.subplots(model.n_components, sharex=True, sharey=True) colours = cm.rainbow(np.linspace(0, 1, model.n_components)) for i, (ax, colour) in enumerate(zip(axs, colours)): # Use fancy indexing to plot data in each state.
    mask = hidden_states == i ax.plot_date(dates[mask], close_v[mask], ".-", c=colour) ax.set_title("{0}th hidden state".format(i)) # Format the ticks.
 ax.xaxis.set_major_locator(YearLocator()) ax.xaxis.set_minor_locator(MonthLocator()) ax.grid(True) plt.show()

 

Out:

Transition matrix
[[  9.79220773e-01   2.57382344e-15   2.72061945e-03   1.80586073e-02]
 [  1.12216188e-12   7.73561269e-01   1.85019044e-01   4.14196869e-02]
 [  3.25313504e-03   1.12692615e-01   8.83368021e-01   6.86228435e-04]
 [  1.18741799e-01   4.20310643e-01   1.18670597e-18   4.60947557e-01]]

Means and vars of each hidden state
0th hidden state
mean =  [  2.33331888e-02   4.97389989e+07]
var =  [  6.97748259e-01   2.49466578e+14]

1th hidden state
mean =  [  2.12401671e-02   8.81882861e+07]
var =  [  1.18665023e-01   5.64418451e+14]

2th hidden state
mean =  [  7.69658065e-03   5.43135922e+07]
var =  [  5.02315562e-02   1.54569357e+14]

3th hidden state
mean =  [ -3.53210673e-01   1.53080943e+08]
var =  [  2.55544137e+00   5.88210257e+15]

 

 

 

---

「紙上得來終覺淺，絕知此事要躬行」，在繼續翻譯《HMM學習最佳範例》以前，這裏先補充幾個不一樣程序語言實現的HMM版本，主要參考了維基百科。讀者有興趣的話能夠研究一下代碼，這樣對於HMM的學習會深入不少！

 

C語言版：
一、 HTK（Hidden Markov Model Toolkit）：
　　HTK是英國劍橋大學開發的一套基於C語言的隱馬爾科夫模型工具箱，主要應用於語音識別、語音合成的研究，也被用在其餘領域，如字符識別和DNA排序等。HTK是重量級的HMM版本。
　　HTK主頁：http://htk.eng.cam.ac.uk/
二、 GHMM Library：
　　The General Hidden Markov Model library (GHMM) is a freely available LGPL-ed C library implementing efficient data structures and algorithms for basic and extended HMMs.
　　GHMM主頁：http://www.ghmm.org/
三、 UMDHMM（Hidden Markov Model Toolkit）：
　　Hidden Markov Model (HMM) Software: Implementation of Forward-Backward, Viterbi, and Baum-Welch algorithms.
　　這款屬於輕量級的HMM版本。
　　UMDHMM主頁：http://www.kanungo.com/software/software.html

Java版：
四、 Jahmm Java Library (general-purpose Java library)：
　　Jahmm (pronounced "jam"), is a Java implementation of Hidden Markov Model (HMM) related algorithms. It's been designed to be easy to use (e.g. simple things are simple to program) and general purpose.
　　Jahmm主頁：http://code.google.com/p/jahmm/

Malab版：
五、 Hidden Markov Model (HMM) Toolbox for Matlab：
　　This toolbox supports inference and learning for HMMs with discrete outputs (dhmm's), Gaussian outputs (ghmm's), or mixtures of Gaussians output (mhmm's).
　　Matlab-HMM主頁：http://www.cs.ubc.ca/~murphyk/Software/HMM/hmm.html

Common Lisp版：
六、CL-HMM Library (HMM Library for Common Lisp)：
　　Simple Hidden Markov Model library for ANSI Common Lisp. Main structures and basic algorithms implemented. Performance speed comparable to C code. It's licensed under LGPL.
　　CL-HMM主頁：http://www.ashrentum.net/jmcejuela/programs/cl-hmm/

Haskell版：
七、The hmm package （A Haskell library for working with Hidden Markov Models）：
　　A simple library for working with Hidden Markov Models. Should be usable even by people who are not familiar with HMMs. Includes implementations of Viterbi's algorithm and the forward algorithm.
　　Haskell-HMM主頁：http://hackage.haskell.org/cgi-bin/hackage-scripts/package/hmm
　　注：Haskell是一種純函數式編程語言，它的命名源自美國數學家Haskell Brooks Curry，他在數學邏輯方面上的工做使得函數式編程語言有了普遍的基礎。

　　是否還有C++版、Perl版或者Python版呢？歡迎補充！

注：原創文章，轉載請註明出處「我愛天然語言處理」：www.52nlp.cn

本文連接地址：http://www.52nlp.cn/several-different-programming-language-version-of-hmm