基於貝葉斯網（Bayes Netword）圖模型的應用實踐初探

1. 貝葉斯網理論部分

筆者在另一篇文章中對貝葉斯網的理論部分進行了總結，在本文中，我們重點關注其在具體場景里的應用。

 

2. 從概率預測問題說起

0x1：條件概率預測模型之困

我們知道，朴素貝葉斯分類器和Logistic regression模型都是產生概率估計來代替硬性的分類。對於每個類值，它們都是估計某個實例屬於這個類的概率。

實際上，大多數其他機器學習分類器都可以轉化為產生這類信息的模型，例如：

• 通過計算葉子節點上每類的相對頻率，就能從決策樹中得到概率
• 通過檢驗某條規則所覆蓋的實例，就能從決策列表中得到概率

概率估計經常比簡單的預測更為有用，它們可以對所做的預測結果進行排名，使期望成本達到最小化。

本質上說，概率預測屬於生成式模型的一種特殊形式，基於概率的分類模型所估計的是給定其他屬性值時類屬性值的條件概率分布。理想情況下，分類模型是用一種簡潔易懂的形式來表達這個條件分布。

由此看來，朴素貝葉斯分類器、Logistic回歸模型、決策樹等只是用不同的方法表達條件概率的分布。當然，它們的表達能力有所差別：

• 朴素貝葉斯分類器和Logistic回歸模型只能表達簡單的分布
• 決策樹可以近似表達任意分布，但決策樹也有缺陷，它將訓練集分隔成越來越小的數據集，必然造成概率估計可靠性的下降，並且還存在重復子樹問題。

總結來說，傳統的概率預測模型有如下幾個困局與挑戰：

• 欠擬合問題，對復雜場景的表征能力有限
• 過擬合問題，以決策樹為代表的分治預測算法存在”過度編碼問題“，雖然這可以通過剪枝技術得以緩解

0x2：解決條件概率預測的兩個學術分支

解決這一問題的一路研究分支是拋棄概率預測模型，轉而研究和使用預測精度更高的判別式模型，例如：

• SVM（支持向量機模型）：
  • SVM的決策面只和訓練集中的少數實例有關，不會因為某些實例的變化而改變分界面，這從理論上證明了過擬合的發生
  • SVM的決策面考慮了最大化原理類間實例距離，這提供了更強的泛化能力，緩解了欠擬合的發生
• 感知機
• 多層感知機、神經網絡
• ......

另一方面，另一路研究分支是引入更為復雜的概率模型，從統計理論方法上出發，采用圖形方式來表達概率分布，這個結構被稱為貝葉斯網絡（Bayesian network）。畫出的圖形就像是節點的網絡圖，每個節點代表一個屬性，節點間用有向線段連接，但不能形成環，即一個有向無環圖（directed acyclic graph）。

我們接下來討論關於貝葉斯網算法的具體應用方面的知識點。

0x3：什么場景適合貝葉斯網建模

貝葉斯網的主要任務是對隨機變量間的條件概率依賴關系（conditional independence relationships  among the variables）進行建模，它的底層思想是：

• 整體世界觀：”事物之間不是簡單的二元的一對一關系，而是多對多的多元依賴關系“
• 局部世界觀：”在復雜糾纏的關系網內，也存在眾多的局部小關系網，這些關系網內部擁有較高的內聚性，對外保持一定的獨立性。或者說整個關系網是相對稀疏的，某個事件的發生一般只和相關的某些因素有關，不太可能和世界上的所有其他因素都存在關系“

貝葉斯網通過條件概率因子分解的形式來抽象表達這種關系，即：

 

每個局部概率分布（local distribution）都有自己的參數集Θ_Xi，⋃ Θ_X 遠小於 Θ，這源自於貝葉斯網的局部條件獨立性。

如下圖所示，A、S、E、O、R、T 分別代表着6個不同的隨機變量，其對應的貝葉斯網關系可能如下：

貝葉斯網關系圖

 

因子分解式

Relevant Link:     

《數據挖掘 實用機器學習工具與技術》Ian H.Witten Eibe Frank，Mark A.Hall

 

3. 從一個具體的例子建立貝葉斯網預測感性認識

在討論貝葉斯網的建模和訓練之前，筆者希望先通過一個具體的例子，來給大家建立一個對貝葉斯網的感性認識。

一是因為在訓練中涉及到的一些技術和概念和預測是共享的，二是通過對預測的討論，將貝葉斯網這個看起來有些抽象的概念拉回到貝葉斯公式、朴素貝葉斯分類器、馬爾柯夫鏈的理解框架之內，有助於我們建立一些直觀的認知和加速理解。

0x1：一個關於天氣預測的貝葉斯網預測示例

下圖是關於一個天氣數據的貝葉斯網實例，

一個天氣數據的簡單貝葉斯網絡實例

可以看到， 圖中的數據包含了4個屬性：outlook、temperature、humidity、windy，以及類屬性play，分別用結點表示。

• play類屬性：先驗概率分別為0.633和9.367，分別和4個屬性結點存在一條有向邊，實際上，這個貝葉斯網退化成了一個朴素貝葉斯結構
• outlook屬性：在先驗play的條件下，分別有：sunny、overcast、rainy這3個后驗概率屬性
• windy屬性：在先驗play的條件下，分別有：false、true這2個后驗概率屬性
• temperature屬性：在先驗play的條件下，分別有：hot、mild、cool這3個后驗概率屬性
• humidity屬性：在先驗play的條件下，分別有：high、normal這2個后驗概率屬性

從這個例子我們可以看到，描述一個貝葉斯網絡，由兩部分組成：

• 網絡結構：結點間的連接方式，結點間的有向線段代表相關關系，並不一定是因果關系
• 結點屬性的概率值：結點內部屬性的條件概率分布，每行的概率對應於一組父節點屬性值的組合，行中的每個概率代表了這個結點屬性的每個屬性值對應於該屬性值組合的概率。每行中各個概率的綜合始終為1

在繼續討論概率預測之前，我們繼續來看一個稍微復雜一些的例子，上面的例子可能還不能非常明顯地表達出貝葉斯網的概念。

0x2：一個稍微復雜一些的天氣預測例子

我們來看一個和上一小節相同問題的，但是情況更為復雜的網絡例子，

天氣數據的另一個貝葉斯網絡

在這個例子中，“windy”、“temperature”、“humidity”有兩個父節點。在這3個節點中，每個父節點都會在屬性概率表格中產生一列，以聯合條件概率的形式逐行排列，右邊的列數等於節點屬性的屬性值數量。

1. 基於貝葉斯概率鏈式法則計算實例每個類值概率

假設我們現在已經確定了某個數據集對應的貝葉斯網的結構和節點屬性概率，接下來的問題是怎么利用這些表來預測某個實例的每個類值的概率呢？答案很簡單，依然是我們熟悉的貝葉斯鏈式法則。

實例的每個屬性都有確定的屬性值，對網絡中的每個結點，根據父節點屬性值找到相應的行，並查看該行結點屬性值的概率，然后將這些概率相乘，得到的累乘結果就是我們要的實例類值概率。

例如，考慮這樣一個實例：【play=？outlook=rainy，temperature=cool，humidity=high，windy=true】，現在需要計算play=no的概率。

• P（play=no）= 0.367
• P（outlook=rainy | play=no）= 0.385
• P（temperature=cool | play=no，outlook=rainy）= 0.429
• P（humidity=high | play=no，temperature=cool）= 0.250
• P（windy=true | play=no，outlook=rainy）= 0.167

所以有：

• P（play=no，outlook=rainy，temperature=cool，humidity=high，windy=true）= 
  • P（play=no）* 
  • P（windy=true | play=no，outlook=rainy）
  • P（outlook=rainy | play=no）*
  • P（temperature=cool | play=no，outlook=rainy）* 
  • P（humidity=high | play=no，temperature=cool）* 
  • = 0.367 * 0.385 * 0.429 * 0.250 * 0.167 = 0.00253069942125
• 同理有：P（play=yes，outlook=rainy，temperature=cool，humidity=high，windy=true）= 0.0077

顯然，這不是最終的答案，最終的概率之和應該為1，而現在0.0025+0.0077顯然不等於1，那么問題出在哪呢？

實際上，它們是 Pr[play=no，E] 和 Pr[play=yes，E] 的聯合概率，其中 E 是指由這個實例的屬性值所給出的所有證據。聯合概率既度量實例的屬性值出現在 E 中的可能性，也度量相應的類值。只有窮盡了包括類屬性在內的所有可能的屬性值組合空間時，它們的和才為1。顯然這個例子不屬於這種情況。

上面這么說可能還有些抽象，通俗地說就是，Pr[play=no，E] 和 Pr[play=yes，E] 只代表了2個概率鏈條，圖中還有其他的概率鏈條，例如：

【play=？outlook=overcast，temperature=cool，humidity=high，windy=true】

outlook屬性值的不同，使概率鏈條進入了另一個分布空間

2. 貝葉斯網概率計算中鏈式分解法則的解釋

給出了實例類值概率的計算過程和結果后，接下來我們繼續追問一個問題。憑什么可以將這些概率相乘呢？

在馬爾柯夫鏈式公式中我們很好理解，條件概率是以因果鏈條的關系依次展開的，如下圖：

一階馬爾科夫鏈公式

每個節點都只考慮其鄰接父節點對其的影響 

現在到了一個圖結構中，還能依然遵循這種鏈式展開思想嗎？答案是肯定的！！

實際上，貝葉斯網是比馬爾柯夫鏈更泛化的一種通用框架，馬爾柯夫鏈是貝葉斯網在馬爾柯夫假設前提下的一種特殊形式。

以上面天氣的例子為例，按照標准貝葉斯鏈式分解法則，實例類值 P(play=?，E) 的全概率鏈式分解公式如下：

• P(play=?，windy=?，outlook=?，humidity=?，temperature=?) =
  • P(play) 
  • P(windy | play) 
  • P(outlook | windy，play) * 
  • P(temperature | outlook，windy，play) * 
  • P(humidity | temperature，toutlook，windy，play) * 

這是一個標准的因果遞歸鏈分解公式。基於這種分解方式得到的是一個完全的網狀結構。

但是，根據貝葉斯網條件概率分布（CPD: conditional probability distribution）獨立性法則，網絡中存在很多結點間是無相關性關系的，也就是說很多結點間是不存在有向邊的。這種現象稱之為貝葉斯網的緊湊性。

基於這種緊湊性原理，標准的貝葉斯鏈式分解法則可以通過約簡得到一個更緊湊的表示，上面全概率公式可以改寫為：

• P(play=?，windy=?，outlook=?，humidity=?，temperature=?) =
  • P(play) = P(play)
  • P(windy | play) = P(windy | play，outlook)
  • P(outlook | windy，play) = P(outlook | play)
  • P(temperature | toutlook，windy，play) = P(temperature | outlook，play)
  • P(humidity | temperature，toutlook，windy，play) = P(humidity | temperature，play)

和上一小節的實例類值公式進行一些對比會發現，在公式形式上是相同的，

• P（play=no，outlook=rainy，temperature=cool，humidity=high，windy=true）=
  • P（play=no）* 
  • P（windy=true | play=no，outlook=rainy）
  • P（outlook=rainy | play=no）*
  • P（temperature=cool | play=no，outlook=rainy）* 
  • P（humidity=high | play=no，temperature=cool）* 

可以看到，首先，貝葉斯全概率鏈式分解定理是累乘合理性的理論基礎。同時，通過局部條件獨立性原理，繁雜的全概率公式可以得到化簡，得到一個相對精簡的概率累乘式。

這在貝葉斯網理論體系中稱之為因子分解式。

a₁, a₂,.. 代表圖中每個節點，Par₍a_i)就表示除了a_i以外的所有節點。公式中的逗號表示“與”的關系 

乘積步驟的有效性有一個前提假設，就是給定每個父節點的屬性值，知道任何其他非子孫結點的屬性值並不能使該結點的各個可能的屬性值所對應的概率發生變化，也就是所謂的局部獨立性假設（local independencies）。

換句話說，就是非子孫結點並不能提供任何多於父節點所能提供有關該結點屬性值的信息，這點可以表示為下式：

 

給定父節點，每個結點對於它的祖父節點、層祖父節點、以及其他非子孫結點都是條件獨立的，在這種情況下，實例類值的累乘式是有效的。乘積步驟遵循概率論中的鏈式規則。

有一點值得注意的，這個分解表達式對於任何一種屬性排列都是成立的，因為貝葉斯網是一個無環圖，可以將網絡結點進行排列，將上式重寫為：

這就是上一小節應用的乘積規則。

3. 兩種貝葉斯網絡的不同

對比一下這個章節和上個章節的兩個貝葉斯網絡，我們會發現它們是完全不同的。

• 第一個網絡具有更嚴格的獨立假設，它的每個結點的父節點都是實例類屬性結點。實際上，它就是朴素貝葉斯分類器
• 第二個網絡中的條件概率列表含有更多的行，因此網絡需要更多的參數。它可以對潛在領域進行更精確地描述

0x3：警惕因果推斷陷阱

假設貝葉斯網中的有向邊代表的是因果關系，這是一個直覺上的判斷。但是要謹慎！！

在本章的例子中，play屬性值會使outlook的某個具體值所對應的期望值提高，但事實上兩者沒有因果關系，實際上兩者之間存在的是相關關系，相關性不等於因果性。 

為某個領域某件貝葉斯網絡的數據分析師們總是期望有向邊能夠表達因果性，然而，當使用機器學習技術（統計方法）從因果結構未知的數據中推導模型時，所能做的只是根據從數據中觀察到的相關性建立網絡，而從相關性推導因果關系的風險總是很大的。

Relevant Link:     

https://www.cnblogs.com/LittleHann/p/11683607.html#_label3_2_1_1

 

4. 貝葉斯網訓練

通過前面的討論我們知道，決定一個貝葉斯網的兩個核心因素分別是：

• 網絡結點連接結構：一個有向無環圖（DAG），其中每個結點 v_i ∈ V 對應一個隨機變量 X_i
• 結點屬性內的條件概率分布：

其中網絡結點連接結構是最關鍵的因素，因為網絡結點間連接方式決定了貝葉斯網概率鏈式分解公式可以如何約簡。確定了概率鏈式分解公式之后，很容易基於貝葉斯定理反向計算出每個結點的條件概率分布。

整個貝葉斯過程如下圖所示：

我們首先來討論如何基於數據集訓練確定貝葉斯網絡的結構。

0x1：貝葉斯網結點結構學習

在貝葉斯網中，實例的每個屬性（包括類別）都有一個結點。學習網絡結構等於是在可能的有向邊空間中進行搜索，對每組連接方案對應的條件概率表進行評估，並計算結果網絡基於某個數據集的對數似然，以此作為對網絡質量的度量。

要實現網絡結構的學習，需要兩個基本組件：

• “網絡結構質量評估方法”：基於局部條件獨立性進行概率累計運算
• “網絡拓朴空間搜索方法”：如何在指數級的窮舉空間中有效地進行拓朴搜索

1. 基於某個數據集的網絡結構質量評估函數

假設網絡結構（即所有的邊）是已知的，很容易估計條件概率列表中的概率分布，只要計算訓練數據中對應屬性值組合的相對頻率即可。同時為了避免出現零概率，一般使用一個常量來初始化分布的初始先驗概率（概率平滑）。

同時為了避免累乘中出現算數下溢（累乘結果數字變得非常小），以至於不能較好地反映質量，因此我們使用概率對數的總和來代替原先的乘積，最終的結果就是給定數據時，整個網絡的對數似然。有一點需要注意，如果直接基於對數似然率對訓練數據進行了最大化，增加更多的有向邊總是會獲得更好的結果，造成最終網絡過度擬合。有很多方法可以緩解這個問題，包括：

• 采用交叉驗證來估計擬合的良好度
• 增加對網絡參數數量的懲罰項
• 賦予網絡結構一個先驗分布，通過組合先驗概率以及數據與網絡相符的概率，找出可能性最大的網絡。這是網絡評分的“貝葉斯”方法 

有兩個常用的度量方法可用於評估網絡質量：

• Akaike信息准則（Akaike Information Criterion，AIC）：AIC得分 = -LL + K
  • 其中LL表示概率對數似然結果
  • K為參數數量
• 基於最短描述長度原理的MDL度量：MDL得分 = -LL + K/2 * logN
  • 其中N為數據集中的實例數量

這兩個表達式的結果都是正數，因此優化目標就是使它們的得分最小化。

如果使用適當的評分度量，搜尋一個好的網絡結構就成功了一半。我們知道，網絡中某個實例的概率是所有條件概率表中單個概率的乘積，所以數據集的總體概率是把所有的實例概率再累乘起來得到。

由於乘積運算中的項是可以互換的，所以這個乘積可以重寫為把同屬於一個表的各個系數組合在一起的形式。這同樣適合使用對數似然用求和運算代替乘積運算。這意味着似然的優化可以在網絡的每個結點中分別進行（分段優化）。

我們可以通過增加或去除其他結點指向正在進行優化的結點的邊來完成，唯一的限制就只不可引入環。如果使用局部評分度量（如AIC或MDL）方案來代替朴素的對數似然，這個方法也同樣生效，因為懲罰項會被分裂成多個組成部分，每個結點各一個，而且每個結點可以進行獨立的優化。

基於以上關於優化過程的討論，就引出了接下來要講的“有向邊連接空間的搜索策略”，貝葉斯網絡結點支持分段優化的這種特性，使得我們可以應用啟發式的優化算法，對網絡結構進行啟發式地探索。 

2. 網絡結點有向邊連接空間搜索策略

各種貝葉斯網絡學習算法的不同之處在於它們在網絡結構空間的搜索方式不同，本小節我們討論一些經典的算法思想，把握其發展的脈絡。

1）K2算法

K2算法是一個簡單而快速的啟發式算法，簡要描述如下：

• 算法起始於某個給定的屬性（即結點）的排序
• 然后對每個結點依次進行處理，
  • 貪心地增加從先前處理過的結點指向當前結點的邊
  • 每一步過程中都增加那些能使網絡得分達到最高值的邊
• 當不再有改進時，便將注意力轉向下一個結點
• 每個結點的父節點數量可以限制在一個預設的最大值范圍內，這是為防止過度擬合附加的機制

由於只考慮起始於前面已經處理過的節點的邊，並且順序是固定的，所以這個過程不會產生環。

2）使用朴素貝葉斯結構作為K2的起始搜索結構

特別地，朴素貝葉斯分類器是這樣一種網絡，它的邊是由類屬性指向其他每個屬性。當為分類建立網絡時，使用這種網絡作為搜尋起點是很有幫助的。可以使用K2方案來實現，強制把類屬性列為序列中的第一個屬性並合理地設定初始邊。

3）利用馬爾可夫毯特性輔助K2進行結構探索

另一個潛在的有用技巧就是有效的利用結點的馬爾可夫毯（Markov blanket）特性。一個結點的馬爾可夫毯包含該結點的父節點、子結點、子結點的父節點。

可以證明結點與它的馬爾可夫毯內的所有其他結點都是條件獨立的。因此，如果一個結點不包括在類屬性的馬爾可夫毯內，那么該結點所代表的屬性與分類就是毫無關系的。

反過來，如果K2發現一個網絡結構沒有將某個相關屬性包含類屬性的馬爾可夫毯內，或者可以增加一條邊來糾正這個缺點。一個簡單的方法就是增加一條從這個屬性結點指向類結點（或相反）的邊，這要取決於哪種方向可以避免環，同時還取決於網絡結構質量評估函數的判定結果。

4）完全貪心策略的K2搜索算法

不進行結點排序，而是直接采取貪心策略增加或去除任意一對結點間的邊，是一個更為周全卻較慢的K2版本。再進一步就是要同時考慮反轉現有的邊。

需要明白的是，使用任何的貪心算法，結果網絡都只能代表某個評分函數局部最大值，並不一定就是全局最優值。這是啟發式搜索的原生弊端。

通常建議采用不同的隨機初始值多次運行算法，還可以使用如模擬退火、禁忌搜索、或遺傳算法等更為復雜的優化策略。

5）HHMC啟發式搜索算法

一種啟發式最短路徑搜索算法。

Relevant Link:    

https://www.google.com.hk/search?newwindow=1&safe=strict&hl=zh-CN&sxsrf=ACYBGNQ6fvIClIfRwRyXUasIvDRpeZpjEQ%3A1574044570225&source=hp&ei=mgPSXYb0Cs6AtgX18K6QBg&q=mmhc.r+&oq=mmhc.r+&gs_l=psy-ab.3...1738.1738..1925...0.0..0.0.0.......0....2j1..gws-wiz.&ved=0ahUKEwiGqd373PLlAhVOgK0KHXW4C2IQ4dUDCAY&uact=5

0x2：貝葉斯網結點條件概率學習

通過上一章節的討論我們知道，通過啟發式的結構搜索，以及基於使每個實例的聯合概率Pr[a₁，a₂，...，a_n]最大化的評分度量，我們可以得到一個局部最優的網絡結構。

確定了結構后，問題就轉變了一個常規的貝葉斯條件概率估計的運算，這通過簡單的統計概率論知識可以做到。

得到結點概率分布如下圖所示：

然而有一點要注意的是，在分類問題中，真正要最大化的是在給定其他屬性值的情況下類的條件概率，換句話說就是條件似然。

但不幸的是，對貝葉斯網絡的列表中所需要的最大條件似然概率估計沒有解析解（啟發式局部搜索沒有全局解析解）。因此，在實際的情況中，網絡的訓練使用標准最大似然估計，而對某個具體的網絡結構則采用條件似然來評估。

Relevant Link:     

《數據挖掘 實用機器學習工具與技術》Ian H.Witten Eibe Frank，Mark A.Hall

  

5. 貝葉斯網推斷（Bayesnetwork Inference）

貝葉斯網推斷本身也屬於貝葉斯推斷范疇內的一種，通常來說，我們討論和使用較多的是下列兩種形式：

• 貝葉斯網最大后驗概率估計推斷（Bayesnetwork Maximum Posteriori (MAP)  Inference）：從一個已知證據出發，推斷其他未知因素的出現概率
• 貝葉斯網條件概率查詢（Bayesnetwork Conditional Probability Queries）：已知一個實例，推斷其在該貝葉斯網條件下出現的似然概率

0x1：貝葉斯網最大后驗概率估計推斷

我們以一個離散隨機變量組成的貝葉斯網為例，其DAG圖如下：

推斷過程示意圖如下：

 

換句話說，貝葉斯網MAP查詢就是在已有證據的前提條件下，推斷網絡中余下結點最有可能的后驗概率組合，從這個角度來看，貝葉斯網MAP是在進行分類預測任務。

 

0x2：貝葉斯網條件概率查詢

Relevant Link:      

 

6. pomegranate - 一種高效的概率圖建模框架

0x1：基本使用介紹

1. 從已知模型參數中創建概率分布

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

from pomegranate import *
import pygraphviz as pgv
import matplotlib.pyplot as plt


def plot_xy(x, y):
    plt.plot(x, y)
    plt.show()


if __name__ == '__main__':
    mu, sig = 0, 2
    model = NormalDistribution(mu, sig)

    x_y = {}
    for i in range(9999):
        x = float("%.1f" % model.sample())
        if x_y.has_key(x):
            x_y[x] += 1
        else:
            x_y[x] = 1

    x_y = sorted(x_y.items(), key=lambda d: d[0])
    print x_y

    x = [i[0] for i in x_y]
    y = [i[1] for i in x_y]

    print x
    print y

    plot_xy(x, y)

2. 從數據中直接學習概率分布

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

from pomegranate import *
import pygraphviz as pgv
import matplotlib.pyplot as plt
import numpy as np


def plot_xy(x, y):
    plt.plot(x, y)
    plt.show()


def generate_data():
    X = np.random.normal(0, 1, 100)
    model = NormalDistribution.from_samples(X)

    return model


def generate_xy_plot(model):
    x_y = {}
    for i in range(9999):
        x = float("%.1f" % model.sample())
        if x_y.has_key(x):
            x_y[x] += 1
        else:
            x_y[x] = 1

    x_y = sorted(x_y.items(), key=lambda d: d[0])
    print x_y

    x = [i[0] for i in x_y]
    y = [i[1] for i in x_y]

    print x
    print y

    return x,y


if __name__ == '__main__':
    model = generate_data()
    x, y = generate_xy_plot(model)

    plot_xy(x, y)

3. HMM detection

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

from pomegranate import *
import pygraphviz as pgv
import matplotlib.pyplot as plt
import numpy as np


if __name__ == '__main__':
    d1 = DiscreteDistribution({
        'A': 0.25,
        'C': 0.25,
        'G': 0.25,
        'T': 0.25
    })
    d2 = DiscreteDistribution({
        'A': 0.10,
        'C': 0.40,
        'G': 0.40,
        'T': 0.10
    })

    s1 = State(d1, name="background")   # 隱狀態
    s2 = State(d2, name="CG island")    # 顯狀態

    hmm = HiddenMarkovModel("CG-detector")
    hmm.add_states(s1, s2)

    # 狀態概率轉移矩陣
    hmm.add_transition(hmm.start, s1, 0.5)
    hmm.add_transition(hmm.start, s2, 0.5)

    # 隱狀態和顯狀態的概率發射矩陣
    hmm.add_transition(s1, s1, 0.9)
    hmm.add_transition(s1, s2, 0.1)
    hmm.add_transition(s2, s1, 0.1)
    hmm.add_transition(s2, s2, 0.9)
    hmm.bake()

    #
    observed_sequence = "GACTACGACTCGCGCTCGCGCGACGCGCTCGACATCATCGACACGACACTC"

    print "how likely the sequence happen: ", hmm.log_probability(list(observed_sequence))
    print "the corresponding hidden state sequence: ", hmm.predict(list(observed_sequence))

4. 貝葉斯網建模”Monty Hall problem“

關於”Monty Hall problem“的背景描述可以參閱這篇文章， 這里對：

• 參賽選手的選擇
• 獎品所在的門
• 主持人打開的門

這3個隨機變量分別抽象為3個結點，根據題意分析，每個結點內部的條件概率分布都是已知的，所以我們可以直接跳過結構學習和概率學習這2個環節，直接進行建模。

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

from pomegranate import *
import pygraphviz as pgv
import matplotlib.pyplot as plt
import numpy as np


if __name__ == '__main__':
    guest = DiscreteDistribution({'A': 1. / 3, 'B': 1. / 3, 'C': 1. / 3})
    prize = DiscreteDistribution({'A': 1. / 3, 'B': 1. / 3, 'C': 1. / 3})
    monty = ConditionalProbabilityTable(
        [['A', 'A', 'A', 0.0],
         ['A', 'A', 'B', 0.5],
         ['A', 'A', 'C', 0.5],
         ['A', 'B', 'A', 0.0],
         ['A', 'B', 'B', 0.0],
         ['A', 'B', 'C', 1.0],
         ['A', 'C', 'A', 0.0],
         ['A', 'C', 'B', 1.0],
         ['A', 'C', 'C', 0.0],
         ['B', 'A', 'A', 0.0],
         ['B', 'A', 'B', 0.0],
         ['B', 'A', 'C', 1.0],
         ['B', 'B', 'A', 0.5],
         ['B', 'B', 'B', 0.0],
         ['B', 'B', 'C', 0.5],
         ['B', 'C', 'A', 1.0],
         ['B', 'C', 'B', 0.0],
         ['B', 'C', 'C', 0.0],
         ['C', 'A', 'A', 0.0],
         ['C', 'A', 'B', 1.0],
         ['C', 'A', 'C', 0.0],
         ['C', 'B', 'A', 1.0],
         ['C', 'B', 'B', 0.0],
         ['C', 'B', 'C', 0.0],
         ['C', 'C', 'A', 0.5],
         ['C', 'C', 'B', 0.5],
         ['C', 'C', 'C', 0.0]], [guest, prize])

    s1 = Node(guest, name="guest")
    s2 = Node(prize, name="prize")
    s3 = Node(monty, name="monty")

    model = BayesianNetwork("Monty Hall Problem")
    model.add_states(s1, s2, s3)
    model.add_edge(s1, s3)
    model.add_edge(s2, s3)
    model.bake()

    model.plot()

5. 貝葉斯網問題預測

基於訓練得到的貝葉斯網模型進行預測任務，就是沿着DAG中結點的連接拓朴，將所涉及結點的條件概率進行累乘。

以上小節Monty Hall問題為例，累乘的對象就是：

• Guest選擇門編號的條件概率
• 獎品所在門編號的概率
• 主持人Monty在看到Guest作出選擇后，所打開門編號的概率

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

from pomegranate import *
import pygraphviz as pgv
import matplotlib.pyplot as plt
import numpy as np


if __name__ == '__main__':
    guest = DiscreteDistribution({'A': 1. / 3, 'B': 1. / 3, 'C': 1. / 3})
    prize = DiscreteDistribution({'A': 1. / 3, 'B': 1. / 3, 'C': 1. / 3})
    monty = ConditionalProbabilityTable(
        [['A', 'A', 'A', 0.0],
         ['A', 'A', 'B', 0.5],
         ['A', 'A', 'C', 0.5],
         ['A', 'B', 'A', 0.0],
         ['A', 'B', 'B', 0.0],
         ['A', 'B', 'C', 1.0],
         ['A', 'C', 'A', 0.0],
         ['A', 'C', 'B', 1.0],
         ['A', 'C', 'C', 0.0],
         ['B', 'A', 'A', 0.0],
         ['B', 'A', 'B', 0.0],
         ['B', 'A', 'C', 1.0],
         ['B', 'B', 'A', 0.5],
         ['B', 'B', 'B', 0.0],
         ['B', 'B', 'C', 0.5],
         ['B', 'C', 'A', 1.0],
         ['B', 'C', 'B', 0.0],
         ['B', 'C', 'C', 0.0],
         ['C', 'A', 'A', 0.0],
         ['C', 'A', 'B', 1.0],
         ['C', 'A', 'C', 0.0],
         ['C', 'B', 'A', 1.0],
         ['C', 'B', 'B', 0.0],
         ['C', 'B', 'C', 0.0],
         ['C', 'C', 'A', 0.5],
         ['C', 'C', 'B', 0.5],
         ['C', 'C', 'C', 0.0]], [guest, prize])

    s1 = Node(guest, name="guest")
    s2 = Node(prize, name="prize")
    s3 = Node(monty, name="monty")

    model = BayesianNetwork("Monty Hall Problem")
    model.add_states(s1, s2, s3)
    model.add_edge(s1, s3)
    model.add_edge(s2, s3)
    model.bake()

    print model.probability([['A', 'A', 'A'],
                             ['A', 'A', 'B'],
                             ['C', 'C', 'B']])

Relevant Link:     

https://www.ctolib.com/jmschrei-pomegranate.html
https://pomegranate.readthedocs.io/en/latest/whats_new.html
https://blog.laisky.com/p/pygraphviz/
http://www.sohu.com/a/150274125_505915
https://github.com/jmschrei/pomegranate/blob/master/docs/BayesianNetwork.rst
https://pomegranate.readthedocs.io/en/latest/faq.html
https://pomegranate.readthedocs.io/en/latest/BayesianNetwork.html

 

7. bnlearn

Relevant Link:  

http://www.bnlearn.com/examples/
http://www.bnlearn.com/examples/useR19-tutorial/
http://www.bnlearn.com/book-crc/