[白話解析]用水滸傳爲例學習最大熵馬爾科夫模型

時間 2020-03-08

標籤白話解析水滸傳爲例學習最大熵模型简体版

原文原文鏈接

[白話解析]用水滸傳爲例學習最大熵馬爾科夫模型

0x00 摘要

本文將盡可能使用易懂的方式，儘量不涉及數學公式，而是從總體的思路上來看，運用感性直覺的思考來解釋最大熵馬爾可夫模型。而且從名著中找了個具體應用場景來幫助你們深刻這個概念。html

在機器學習過程當中，會遇到不少晦澀的概念，相關數學公式不少，你們理解起來頗有困難。遇到相似狀況，咱們應該多從直覺角度入手思考，用類別或者舉例來附會，這樣每每會有更好的效果。java

我在講解論述過程當中給本身的要求是：在生活中或者名著中找一個例子，而後用本身的話語闡述出來。python

0x01 背景

在前文[白話解析]以水滸傳爲例學習隱馬爾可夫模型中，咱們用 "梁中書突圍大名府爲例" 講解了隱馬爾可夫模型，可是現實中狀況可能更復雜，好比梁中書突圍遇到了宋江，他再次選擇從宋江處突圍的可能性會變低，由於宋江身邊確定是梁山大部分好漢，突圍難度太大。可是若是遇到史進，危險性就沒有那麼大了。git

這種狀況只用隱馬爾科夫模型就很難解決，須要引入最大熵馬爾可夫模型了。github

最大熵馬爾可夫模型（Maximum Entropy Markov Model，MEMM）的思想把HMM和最大熵結合起來，能夠提取特徵以泛化模型能力，結合上下文依賴，直接判別減小建模負擔。下面就來詳述。算法

1. 隱馬爾科夫模型的缺點

HMM是一種生成模型，定義了聯合機率分佈，其中y和x分別表示觀察序列和相對應的標註序列的隨機變量。爲了可以定義這種聯合機率分佈，生成模型須要枚舉出全部可能的觀察序列，這在實際運算過程當中很困難的，所以咱們須要將觀察序列的元素看做是彼此孤立的個體，即假設每一個元素彼此獨立，任什麼時候刻的觀察結果只依賴於該時刻的狀態。數組

目標函數和預測目標函數不匹配，HMM學到的是狀態和觀察序列的聯合分佈P(Y,X)，而預測問題中，咱們需要的是條件機率P(Y|X)。HMM必須計算出全部的潛在可能路徑的機率值大小（而後再挑機率值最大的那一個做爲最終結果）數據結構

HMM僅僅依賴於每一個狀態和它相應的觀察對象，可是對於某一個觀測值，它可能並不是由一個隱藏狀態決定的，而是兩個以上的隱藏狀態綜合做用而產生的，那麼這時HMM就無能爲力了。app

HMM模型的這個假設前提是在比較小的數據集上是合適的，但實際上在大量真實的語料中觀察序列更多的是一種多重交互特徵形式表現，觀察元素之間普遍存在長程相關性。在命名實體識別的任務中，因爲實體自己結構所具備的複雜性，利用簡單的特徵函數每每沒法涵蓋全部的特性，這時HMM的假設前提使得它沒法使用複雜特徵（它沒法使用多於一個標記的特徵）。框架

2. 最大熵模型的特色

最大熵模型的優勢：首先，最大熵統計模型得到的是全部知足約束條件的模型中信息熵極大的模型；其次，最大熵統計模型能夠靈活地設置約束條件，經過約束條件的多少能夠調節模型對未知數據的擬合程度；再次，他還能天然地解決了統計模型中參數平滑的問題。

最大熵模型的不足：首先，最大熵統計模型中二值化特徵只是記錄特徵的出現是否，而文本分類須要知道特徵的強度，所以，它在分類方法中不是最優的。其次，因爲算法收斂的速度比較慢，因此致使最大熵模型它的計算代價比較大，時空開銷大；再次，數據稀疏問題比較嚴重。

最大熵模型可使用任意的複雜相關特徵，在性能上最大熵分類器超過了bayes分類器。可是，做爲一種分類器模型，這兩種方法有一個共同的缺點：每一個詞都是單獨進行分類的，標記之間的關係沒法獲得充分利用。具備馬爾可夫鏈的HMM模型能夠創建標記之間的馬爾科夫關聯性，這是最大熵模型所沒有的。因此一個很天然的想法就是將二者的優點結合起來，這就獲得了最大熵馬爾可夫模型。

3. MEMM的引入

最大熵馬爾可夫模型（Maximum Entropy Markov Model，MEMM）的思想是利用 HMM 框架預測給定輸入序列的序列標籤，同時結合多項 Logistic 迴歸（又名最大熵，其給出了能夠從輸入序列中提取特徵的類型和數量上的自由度）。這個模型容許狀態轉移機率依賴於序列中彼此之間非獨立的特徵上，從而將上下文信息引入到模型的學習和識別過程當中，提升了識別的精確度，召回率也大大的提升，有實驗證實，這個新的模型在序列標註任務上表現的比 HMM和無狀態的最大熵模型要好得多。

MEMM解決了HMM輸出獨立性假設的問題。由於HMM只限定在了觀測與狀態之間的依賴，而MEMM引入自定義特徵函數，不只能夠表達觀測之間的依賴，還可表示當前觀測與先後多個狀態之間的複雜依賴；（這個複雜依賴是經過特徵函數來作到的，好比該單詞的「前面單詞/後面單詞」分別是什麼）
MEMM考慮到相鄰狀態之間依賴關係。且考慮整個觀察序列，所以MEMM的表達能力更強；
MEMM不考慮P(X)，減輕了建模的負擔。同一時候學到的是目標函數是和預測函數一致；

HMM是生成式模型，參數即爲各類機率分佈元參數，數據量足夠能夠用最大似然估計。

而MEMM是判別式模型。是用函數直接判別，學習邊界，MEMM即經過特徵函數來界定。這裏因爲去掉了獨立性假設，因此不能給出聯合機率分佈，只能求後驗機率，因此是判別模型。但一樣，MEMM也有極大似然估計方法、梯度降低、牛頓迭代、擬牛頓降低、BFGS、L-BFGS等等。

0x02 最大熵模型

由於HMM咱們已經在前文了解，因此本文重點是最大熵模型的介紹。

1. Logistic 迴歸

一般，機器學習分類器經過從全部可能的 y_i 中選擇有最大的 P(y|x) 的那個，來決定將哪一個輸出標籤 y 分配給輸入 x。在分類任務中，logistic 迴歸經過計算給定觀察的屬於每一個可能類別的機率，而後選擇產生最大機率的類別。

Logistic 迴歸是用於分類的一種監督學習算法，它的本質是線性迴歸。Logistic 迴歸用條件極大似然估計進行訓練。這意味着咱們將選擇參數 w，使對給定輸入值 x 在訓練數據中 y 標籤的機率最大化。一般能夠運用隨機梯度降低法、L-BFGS 或共軛梯度法來求此函數的最大值，即找到最優權重。

最大熵模型屬於log-linear model，在給定訓練數據的條件下對模型進行極大似然估計或正則化極大似然估計。

2. 最大熵原理

最大熵原理(Principle of Maximum Entropy)：將已知事實做爲制約條件，求得可以使熵最大化的機率分佈做爲正確的機率分佈。而在機器學習裏，就是若是沒有數聽說明，就不要隨便爲模型加假設。

最大熵原理的實質就是，在已知部分知識的前提下，關於未知分佈最合理的推斷就是：符合已知知識後「最不肯定或最隨機的推斷」是咱們能夠做出的惟一不偏不倚的選擇，任何其它的選擇都意味着咱們增長了其它的約束和假設，這些約束和假設根據咱們掌握的信息沒法做出。

例如，投擲一個骰子，若是問"每一個面朝上的機率分別是多少"，你會說是等機率，即各點出現的機率均爲1/6。由於對這個"一無所知"的色子，什麼都不肯定，而假定它每個朝上機率均等則是最合理的作法。從投資的角度來看，這是風險最小的作法，而從信息論的角度講，就是保留了最大的不肯定性，也就是說讓熵達到最大。

實踐經驗和理論計算都告訴咱們，在徹底無約束狀態下，均勻分佈等價於熵最大（有約束的狀況下，不必定是機率相等的均勻分佈。好比，給定均值和方差，熵最大的分佈就變成了正態分佈）。

3. 最大熵模型

最大熵模型(Maximum Entropy Modeling) : 若是你有數據，請經過你的數據創建先驗信息（在這裏叫作約束條件），剩下的未知部分，請讓他們均勻分佈，也就是讓他們的熵最大。

由於咱們若是要搭建一個最大熵模型來實現分類，那麼咱們定義模型 p(y|x)，這個模型應該是：在 "知足事先已約束" 的條件下的模型中選擇熵最大的模型，即讓不肯定的信息等可能的發生。這樣咱們就獲得了最終的分類模型。

即，給定一個訓練樣本集，咱們但願尋找一個分佈符合以下兩個條件(Given a set of training examples, we wish to find a distribution which)：

知足已知的約束條件（satisfies the input constraints）
最大化其不肯定性（maximizes the uncertainty）

4. 已知的約束條件

上下文 & 特徵函數

咱們作分類問題，看到的數據每每是每一個實例對應一個類別。好比說詞性標註中一個詞對應一個標註。爲了下面討論的方便，將類別稱爲Outcome，將每一個實例的上下文環境叫作Context。實例被稱爲Event，一個實例是Outcome與Context的二元組。

好比，一個多維的向量x，x的每一個維度都是一個特徵，能夠認爲 x 對應了一個 Context（特徵的集合）。而後這條樣本對應了一個label（Outcome）。

問題就是數據並非乖乖地排列好，x的每一個維度都已經取好值等着咱們分類了，因此出現了這個東西：特徵函數。咱們使用特徵函數來表示約束條件，即特徵函數也就是用來指示元素是否屬於某一個子集。

好比記 w 是句子 s 的某個單詞，w是有不少維度的，則用來表示這些維度的特徵函數能夠是：

w在句首
w在句尾
w的前綴是xxx
w的後綴是xxx
w前面的單詞是xxx

能夠看出，這些特徵函數針對這個單詞 w 的判別結果，就構成了 w 的上下文 Context。好比: w不在句首，w在句尾.....

特徵函數

給定一個訓練數據集
\[ T = \{ (x_1, y_1), (x_2, y_2) ... (x_N, y_N) \} \]
爲了表示數據，咱們從數據中抽取了一系列的特徵。通常說的「特徵」都是指輸入的特徵，而最大熵模型中的「特徵」指的是輸入和輸出共同的特徵。每一個特徵與每一個類別都能組成一個特徵，應該用乘法原理計數。也就是說，咱們能夠將X的每一維度特徵下的每一個取值與某個類別配對，並一樣的用一個參數來描繪這個配對的緊密程度。

能夠這麼理解：就是把廣泛意義上的"輸入的特徵"拆分處理變成"輸入和輸出共同的特徵"。

好比原來是：
\[ \{ X = 1 \} ==> \{ Y = Class_i \} \]
如今是
\[ f(x,y) = \begin{cases} 1 & \text{當x,y知足某一事實} \\[2ex] 0 & \text{當不知足該事實} \end{cases} \]
好比每一個特徵函數能夠是從Context到0-1的二值函數。
\[ f(x,y) = \begin{cases} 1 & \text{if $x$ = "is" and y = v} \\[2ex] 0 & \text{otherwise} \end{cases} \]
咱們爲每一個 <feature,label> 對都作一個如上的特徵函數，用來描述數據集數學化。

權重

最大熵模型中的每一個特徵會有一個權重，就是上面說的參數。能夠把它理解成這個特徵所描述的輸入和輸出有多麼傾向於同時出現。若是某類數據中某種輸入和輸出傾向於同時出現，那麼對於這一類數據，表示這一對輸入和輸出的特徵就會有較高的權重。若是認爲這一對徹底不可能成的話即讓這一組的權重爲0。

用途

特徵函數就是用來形式化表示X的相關知識信息：在全部的特徵中，哪些特徵最有表明性，哪些特徵是次要的，須要本身選擇，構造約束集合。一組特徵函數就將Context從上下文空間映射到特徵空間上了。

特徵函數的出現，讓模型獲得了更好的泛化能力。那麼如何讓一個線性模型H(x)也有相似的功能？答案就是特徵函數，讓輸入x先通過一系列特徵函數的處理，變成g(x)，再送給模型分類，好比H(x) = H(g(x))。

特徵函數f(x,y)並不只僅表明計數，它還表明了某種對x和y（尤爲是x）的簡化。一組數據(x,y)通常狀況下並不僅觸發一個特徵。特徵除了「x取特定值，y取特定值」的形式之外，還能夠是「x取某類值，y取某類值」的形式，更通常的是「x和y的取值知足某種關係」的形式。正是這些通常形式的特徵才讓最大熵模型有足夠的泛化能力。當一組數據不僅觸發一個特徵的時候，exp內就不止一個權重求和，就求不出解析解了。

在最大熵模型的視角下，每條輸入的n個「特徵」與k個類別共同組成了nk個特徵，模型中有nk個權重，與特徵一一對應。每一個類別會觸發nk個特徵中的n個。特徵的全體能夠看作是n個特徵函數組成的一個集合。

5. 模型知足已知的約束條件

因爲特徵函數是對創建機率模型有益的特徵，因此應該讓最大熵模型來知足這些約束，經驗分佈與特徵函數結合便能表明機率模型須要知足的約束。

擬合

回到咱們以前給定的訓練數據集
\[ T = \{ (x_1, y_1), (x_2, y_2) ... (x_N, y_N) \} \]
假如數據(x1, y1) 知足了特徵函數 f1(x,y)，那麼若是咱們訓練出的模型 M 也能得出 f1(x1,y1) = 1，則說明 M 能對數據作在(x1, y1)這點上對f1 作了很好的擬合。

機率分佈

對於給定訓練數據集，如今把訓練數據當作由隨機變量 (𝑋,𝑌) 產生。咱們能夠根據訓練數據肯定聯合分佈 P(x,y) 的經驗分佈 𝑃˜(x,y) 和邊緣分佈 P(x)的經驗分佈𝑃˜(x), 即：
\[ 𝑃˜(X=x, Y=y) = \frac{count((X=x, Y=y))}{N} \]

\[ 𝑃˜(X=x) = \frac{count((X=x))}{N} \]

其中，count((X=x, Y=y) 表示訓練數據中樣本 (x,y) 出現的頻數， count((X=x)) 表示訓練數據集中 x 出現的頻數。 N 表示訓練樣本的總容量。

指望

下面介紹兩個指望。

指望一：如今，咱們觀察到一組數據集，經過簡單的統計能夠知道任意一個Context x 和 Outcome y 的組合的聯合機率。有了聯合機率，能夠計算觀察到的某一特徵函數 f 的指望，稱爲 觀察指望/經驗指望。即特徵函數 f(x,y) 關於經驗分佈𝑃˜(x,y) 的指望值。
\[ E_{ref}(f) = \sum_{x,y}p̃(x,y)f(x,y) = \frac{\sum_{x,y}f(x,y)}{N} \]
因爲特徵函數是對創建機率模型有益的特徵，因此應該讓最大熵模型來知足「觀察指望/經驗指望」這一約束，

指望二：假設咱們有一個模型，那麼咱們能夠從這個模型的角度求出這個特徵函數的指望，稱爲 模型指望。即特徵函數 f(x,y) 關於模型 p(y|x) 和經驗分佈 p̃(x) 的指望值。
\[ E_q(f) = \sum_{x,y}p(x,y)f(x,y) \approx\sum_{x,y}p̃(x)p(y|x)f(x,y) \]

E_ref 實際上表明瞭 訓練數據 之中符合特徵函數的數據的佔比。能夠理解爲就是收集到的數據。
E_q 實際上表明的是 模型先驗數據 中符合約束條件的佔比。能夠理解爲是訓練時獲得的數據。

機器學習的目的就是從數據集中學得數據中包含的某種內在信息，所以咱們但願咱們的模型可以很好地反應這些數據中蘊含的現象。那麼從模型角度看到的 f 的指望就應該等於從數據觀察到的 f 的指望。也就是Eq(f) = Eref(f)。或者說，如今先驗數據中的符合約束條件的佔比等於訓練數據中符合特徵函數的數據佔比，實際上就說明了模型已經符合約束條件。

模型集合

假設咱們有n個特徵函數，那麼咱們就有n組等式Eq(fi)=Eref(fi),i∈1,2,…,n。

假設咱們有那麼多的模型，也能夠認爲是機率分佈。他們組成一個空間P，而知足上面一系列特徵函數指望構成的等式的機率分佈構成了 P 的一個子集。
\[ C = \{ p|p∈ P \and E_q(f_i)=E_{ref}(f_i), i∈(1,2,...,n) \} \]
如今要找一個合適的模型描述數據，也就是在 P 中搜索一個p。從知足約束的模型集合C 中找到使得 P(Y|X) 的熵最大的即爲最大熵模型了。

6. 最大化模型不肯定性

對於這個模型p，這個模型須要知足上面一系列特徵函數指望構成的等式（換句話講，是一系列特徵函數指望構成的約束）。從全部特徵函數構成的整個機率分佈出發的話，又須要讓未知事件保持最無序的狀態（等機率分佈的時候就是最無序的）。而衡量無序度的指標就是熵！因此，既有約束條件，又要儘量保持最無序的狀態，儘量將可能性均勻地非配到不肯定的上下文狀況中。那目標狀態就是知足約束條件的狀況下，熵最大的狀態。這裏的熵指的是條件熵（已知x的狀況下y的無序度）

肯定了約束條件，求取最大熵的狀況其實就是：求取有約束條件下的最大熵值。即，最大熵模型就是在知足約束的模型集合中選擇條件機率分佈 p(y|x) 最大的模型，其形式化以下：

根據經驗分佈獲得知足約束集的模型集合 C，即模型應該知足下面兩個約束

\[ E_{ref}(f_i) = E_q(f_i), i =1,2,..,n \\ \sum_y p(y|x) = 1 \]

咱們的訓練目標，即目標函數，也就是在這兩個約束下對下面的式子求最值

\[ H(P) = \{ -\sum_{x,y}p̃(x)p(x|y)logp(y|x) \} \]

上面第一個約束就是 E_ref (訓練數據之中符合特徵函數的數據的佔比) = E_q (模型先驗數據中符合約束條件的佔比)，即模型已經符合約束條件。第二個約束就是機率分佈應該是1。

求最值的式子就是定義在條件機率分佈p(y|x)上的條件熵公式。求最值所得出的解就是最大熵模型學習的模型。

即
\[ max_{𝑃∈𝐶} H(P) = max_{𝑃∈𝐶} \{ -\sum_{x,y}p̃(x)p(x|y)logp(y|x) \} \]

熵越大，可能性就越平均地被分配，於是咱們的最終目標是最大化一個模型的熵。而因爲有前面的約束等式，這個問題變成了一個有約束的最優化問題。

7. 拉格朗日乘子法

MaxEnt 模型最後被形式化爲帶有約束條件的最優化問題。對於"求取有約束條件下的最值"，很容易咱們就想到了拉格朗日乘子法，引入拉格朗日乘子： w0,w1,…,wn，定義拉格朗日函數 𝐿(𝑃,𝑤), 將有約束化爲無約束。

如今問題能夠形式化爲便於拉格朗日對偶處理的極小極大的問題。求解最優的 w∗ 後，便獲得了所要求的MaxEnt 模型。

最終經過一系列極其複雜的運算，獲得了須要極大化的式子(這裏 fi(x,y) 表明特徵函數，wi 表明特徵函數的權值)：
\[ max_{𝑝∈𝐶}∑_{𝑥,𝑦}𝑃˜(𝑥,𝑦)∑_{𝑖=1}^𝑛𝑤_𝑖𝑓_𝑖(𝑥,𝑦)+∑_𝑥𝑃˜(𝑥)𝑙𝑜𝑔𝑍_𝑤(𝑥) \]
這裏 fi(x,y) 表明特徵函數，wi 表明特徵函數的權值， Pw(y|x) 即爲 MaxEnt 模型，將最優解記作 w∗。

8. 極大似然估計

可是上述過程仍是太複雜。有沒有簡單易行的方式呢？答案就是極大似然估計 MLE 了。數學證實能夠得出，最大熵模型學習中的對偶函數極大化等價於最大熵模型的極大似然估計。極大似然估計函數是一個凸函數,這是優化問題中最容易優化的模型，咱們能夠獲得全局最優解。

最大熵模型的似然是使用了(模型學的)真實分佈 p(x,y) 與(來自數據的)經驗分佈p̃(x,y) 的交叉熵來定義。

這裏有訓練數據獲得經驗分佈 𝑃˜(𝑥,𝑦) ，待求解的機率模型 𝑃(𝑌|𝑋) 的似然函數爲：

\[ 𝐿_𝑃˜(𝑃_𝑤)=𝑙𝑜𝑔∏_{𝑥,𝑦}𝑃(𝑦|𝑥)^{𝑃˜(𝑥,𝑦)}=∑_{𝑥,𝑦}𝑃˜(𝑥,𝑦)𝑙𝑜𝑔𝑃(𝑦|𝑥) \]
令 𝑍𝑤(𝑥)表示 𝑒𝑥𝑝(1−𝑤0) ，將Pw(y|x) 帶入公式能夠獲得：
\[ 𝐿_𝑃˜(𝑃_𝑤)=∑_{𝑥,𝑦}𝑃˜(𝑥,𝑦)𝑙𝑜𝑔𝑃(𝑦|𝑥) = \sum_{x,y}𝑃˜(x,y)\sum_{i=1}^nw_if_i(x,y) - \sum_x𝑃˜(x)logZ_w(x) \]
如今只需極大化似然函數便可，順帶優化目標中能夠加入正則項，這是一個凸優化問題，通常的梯度法、牛頓法均可解之，專門的算法有GIS和IIS 算法。

0x03 最大熵模型算法實現

1. 算法總述

如今咱們基本獲得了算法要作的工做，就是從數據中估計出一個特徵函數的權向量λ，最大化"經驗分佈的最大似然估計"。若想讓預測出的結果所有正確的機率最大，根據最大似然估計,就是全部樣本預測正確的機率相乘獲得的P(整體正確)最大,

即, 有n個特徵函數fi(x)，那麼就有n組等式Eq(fi)=Eref(fi),i∈1,2,…,n。這些就是約束條件集合。咱們的訓練目標就是目標函數即熵H(p(y|x))。咱們求「特徵函數f(x,y)和其權重λ的一套組合」，該組合會令熵 H(p(y|x)) 取到最大值。或者說，咱們要的就是根據模型參數來計算指望。最後儘量的優化到使模型指望和先驗指望同樣且熵最大。

算法的流程以下：

初始化λ=0
循環
\[ λ_i^{(t+1)} = λ_i^{(t)} + \frac{1}{C} log \frac {E_{ref}(f_i)}{E_q(f_i)} \]

\[ C = max_{x,y}\sum_{i=1}^nf_i(x,y) \]
重複2到收斂

根據上面算法，在最大熵模型的實現過程當中，咱們須要計算的值包括經驗指望 E_ref(f) 和各輪迭代事後的模型指望E_q(f)。

經驗指望 Eref(fi)=∑x,yp̃ (x,y)fi(x,y)，求Eref(fi)只須要統計訓練數據中符合fi的(x,y)二元組的個數，而後除以訓練實例的個數N。

模型指望 須要首先求p(y|x)。這個條件機率能夠經過簡單地將全部(x,y)符合的fi和對應的參數λi乘起來後相加。歸一化因子是各個Outcome y的p(y|x)的和。在求得p(y|x)後，要求Eq(fi)，只須要枚舉全部符合fi的(x,y)對應的p(y|x)，乘以Context x出現的次數再除以N就能夠。

模型指望有了，經驗指望有了，把他們放到算法裏面去迭代就行了。

如下代碼出自easyME

2. 數據結構和初始化

//咱們作分類問題，看到的數據每每是每一個實例對應一個類別。好比說詞性標註中一個詞對應一個標註。爲了下面討論的方便，將類別稱爲Outcome，將每一個實例的上下文環境叫作Context。 實例被稱爲Event，一個實例是 Outcome y 與Context的二元組
//一個多維的向量x。x的每一個維度都是一個特徵，能夠認爲 x 對應了一個 Context（特徵的集合）。而後這條樣本對應了一個label（Outcome）
struct Event {
    size_t classId; //將類別稱爲Outcome
    std::vector<size_t> context; //將每一個實例的上下文環境叫作Context, 能夠認爲是一系列 feature 的列表, 每一個 feature 被映射成一個數字。
    size_t count; //符合該上下文環境的(x,y)二元組的個數
}

// X的每一維度特徵下的每一個取值與某個類別配對，並一樣的用一個參數（權重）來描繪這個配對的緊密程度
class DataManager {
  typedef std::pair<size_t, double> Pair;
  typedef std::vector<Pair> Param;
  std::vector<Event> mEventSet;
  std::vector<Param> mParamSet; //這個就是最後模型參數（權重）
  //fetId 特徵緯度，classPos 類別
  void DataManager::incLambda(size_t fetId, size_t classPos, double incer) {
      mParamSet[fetId][classPos].second += incer; //更新模型參數（權重）
  }
}

//初始化
bool MaxEntModel::initModel(const char * trainFileName, bool freq, bool isSelectFeature) {
    string line, str;
    vector<size_t> context;
    size_t count = 1;
    // each line is a event, it looks like this: (count) className fetName ... fetName
    while(getline(fin, line)){
        istringstream sin(line);
        // with freq ?
        if(freq) sin >> count;
        sin >> str;
        size_t classId = mClassMap.insertString(str);
        context.clear();
        while(sin >> str){
            size_t fetId = mFetMap.insertString(str);
            context.push_back(fetId); //由此能看出，context是一系列的 feature
        }
        mModelInfo.addEvent(count, classId, context);
    }
    if(!freq) mModelInfo.processEventSet();
    mModelInfo.getAllFeatures();  
}

//在最大熵模型的視角下，每條輸入的n個「特徵」與k個類別共同組成了nk個特徵，模型中有nk個權重，與特徵一一對應。每一個類別會觸發nk個特徵中的n個。特徵的全體能夠看作是n個特徵函數組成的一個集合。
//Event是 實例是Outcome與Context的二元組，count是知足這個組合的數據個數，即，符合該上下文環境的(x,y)二元組的個數
void DataManager::addEvent(size_t count, size_t classId, const std::vector<size_t> & fetVec) {
      mTotEvent += count;
    mEventSet.push_back(Event(count, classId, fetVec));
}

void DataManager::getAllFeatures() {
    size_t eventNum = getEventNum();
    for(size_t i = 0; i < eventNum; ++i){
        int cid = getEventClassId(i);
        for(context_iterator it = getContextBegin(i),
                end = getContextEnd(i);
                it != end; ++it){
            int fid = *it;
            addFeature(cid, fid);
        }
    }
    endAddFeature();
}

void DataManager::addFeature(size_t classId, size_t fetId, double weight = 0.0) {
    if(mParamSet.size() < fetId + 1) mParamSet.resize(fetId + 1);
    if(classId > mMaxCid) mMaxCid = classId;
    if(fetId > mMaxFid) mMaxFid = fetId;
    mParamSet[fetId].push_back(std::make_pair(classId, weight));
}

DataManager::context_iterator DataManager::getContextBegin(size_t eventId) {
    return mEventSet[eventId].context.begin();
}

3. 獲取指望

//Eref(fi)，訓練數據之中符合特徵函數的數據的佔比。能夠理解爲就是收集到的數據。 
//統計訓練數據中符合fi的(x,y)二元組的個數，而後除以訓練實例的個數N
void DataManager::getObserves(vector<vector<double> > & observed) {
    size_t mMaxFid = getFetNum();
    size_t eventNum = getEventNum();
    observed.clear();
    observed.resize(mMaxFid + 1);
    for(size_t i = 1; i <= mMaxFid; ++i){
        DataManager::param_iterator begin = getParamBegin(i);
        DataManager::param_iterator end = getParamEnd(i);
        size_t n = end - begin;
        observed[i].resize(n, 0);
    }
    //每一個event就是一個實例，一個實例是Outcome與Context的二元組
    for(size_t i = 0; i < eventNum; ++i){ 
        size_t classId = getEventClassId(i); //該event實例對應的類 Outcome
        size_t count = getEventCount(i); //知足該event實例的(x,y)的個數
        for(DataManager::context_iterator it = getContextBegin(i),end = getContextEnd(i);
                it != end; ++it){ //遍歷context的全部feature
            int fid = *it;
            int pos = getClassPosition(classId, fid); //找到mParamSet中該類的位置
            if(pos != -1)
                observed[fid][pos] += count; //統計訓練數據中符合該fi的(x,y)二元組的個數，由於每一個event的context包含了多個feature，因此要對特定fid進行累加
        }
    }
}

//E_q(f)，模型先驗數據中符合約束條件的佔比。能夠理解爲是訓練時獲得的數據。
double DataManager::getExpects(vector<vector<double> > & expected) {
    vector<double> probs;
    expected.clear();
    expected.resize(mParamSet.size());
    for(size_t i = 0; i < mParamSet.size(); ++i)
        expected[i].resize(mParamSet[i].size(), 0);
    double logLike = 0;
    for(size_t i = 0, eventNum = getEventNum(); i < eventNum; ++i){ 
        context_iterator ctxtBegin = getContextBegin(i);
        context_iterator ctxtEnd = getContextEnd(i);
        getAllProbs(ctxtBegin, ctxtEnd, probs); //
        size_t count = getEventCount(i); // 知足的個數
        size_t classId = getEventClassId(i); // 類別
        vector<double> newProbs;
        for(size_t i = 0; i < probs.size(); ++i)
            newProbs.push_back(probs[i] * count);
        for(context_iterator it = ctxtBegin; it != ctxtEnd; ++it){
          size_t fid = *it;
          param_iterator pBegin = getParamBegin(fid);
          param_iterator pEnd = getParamEnd(fid);
          for(param_iterator pit = pBegin; pit != pEnd; ++pit){
            int pos = pit - pBegin;
            size_t cid = pit->first;
            expected[fid][pos] += newProbs[cid];
          }
                }
                logLike += log(probs[classId]) * count;
    }
   return logLike;
}

//求p(y|x)。這個條件機率能夠經過簡單地將全部(x,y)符合的fi和對應的參數λi乘起來後相加。 歸一化因子是各個Outcome y 的p(y|x)的和。在求得p(y|x)後，要求Eq(fi)，只須要枚舉全部符合fi的(x,y)對應的p(y|x)，乘以Context x 出現的次數再除以N就能夠。 
size_t DataManager::getAllProbs(
        const context_iterator begin,
        const context_iterator end,
        vector<double> & probs) {
    probs.clear();
    probs.resize(mMaxCid + 1, 0);
    for(context_iterator cit = begin; cit != end; ++cit){
        size_t fid = *cit;
        for(param_iterator it = getParamBegin(fid),
                pend = getParamEnd(fid);
                it != pend; ++it){
            size_t cid = it->first; // Outcome y
            probs[cid] += it->second; //參數λi
        }
    }
    size_t maxK = 0;
    double sum = 0;
    for(size_t i = 1; i <= mMaxCid; i++){
        probs[i] = exp(probs[i]);
        sum += probs[i]; //參數λi乘起來後相加。
        if(probs[i] > probs[maxK])
            maxK = i;
    }
    for(size_t i = 1; i <= mMaxCid; ++i)
        probs[i] /= sum; //除以N
    return maxK;
}

4. 具體訓練過程

void MaxEntTrainer::_initTrainer(void) {
    mEventNum = mModelInfo.getEventNum();
    mMaxFid = mModelInfo.getFetNum();
    mMaxCid = mModelInfo.getClassNum();
    mTotEvent = mModelInfo.getAllEventFreq();
    mModelInfo.getObserves(mObserved);
}

// Calculate the ith GIS parameter updates with Gaussian prior
// using Newton-Raphson method
// the update rule is the solution of the following equation:
//                                   lambda_i + delta_i
// E_ref = E_q * exp (C * delta_i) + ------------------ * N
//                                       sigma_i^2
// note: E_ref and E_q were not divided by N
// this function is copied from Le Zhang's work
double MaxEntTrainer::_newton(double f_q, double f_ref, double lambdaNow, double mEps)
{
    size_t maxiter = 50;
    double x0 = 0.0;
    double x = 0.0;

    for (size_t mIter = 1; mIter <= maxiter; ++mIter) {
        double t = f_q * exp(mSlowF * x0);
        double fval = t + mTotEvent * (lambdaNow + x0) / mSigma2 - f_ref;
        double fpval = t * mSlowF + mTotEvent / mSigma2;
        if (fpval == 0) {
            cerr << "Warning: zero-derivative encountered in newton() method." << endl;
            return x0;
        }
        x = x0 - fval/fpval;
        if (fabs(x-x0) < mEps)
            return x;
        x0 = x;
    }
    cerr << "Failed to converge after 50 iterations in newton() method" << endl;
    exit(-1);
}

bool GisTrainer::train() {
    vector<vector<double> > expected;
    double preLogLike = -1e10;
    for(size_t it = 0; it < mIter; ++it) {
        double newLogLike = mModelInfo.getExpects(expected);
        for(size_t i = 1; i <= mMaxFid; ++i) {
            DataManager::param_iterator begin = mModelInfo.getParamBegin(i);
            DataManager::param_iterator end = mModelInfo.getParamEnd(i);
            for(DataManager::param_iterator it = begin; it != end; ++it) {
                size_t j = it - begin;
                double inc = 0;
                if(mSigma2) {
                    inc = _newton(expected[i][j], mObserved[i][j], it->second);
                } else if(mAlpha) {
                    if(mObserved[i][j] - mAlpha <= 0)
                        continue;
                  inc = (log(mObserved[i][j] - mAlpha) - log(expected[i][j])) / mSlowF;
                    if(it->second + inc <= 0)
                        inc = -it->second;
                } else {
                            inc = (log(mObserved[i][j]) - log(expected[i][j])) / mSlowF;
                        }
                mModelInfo.incLambda(i, j, inc); //更新模型參數（權重）
            }
        }
        if(fabs((newLogLike - preLogLike) / preLogLike) < mEps)
            break;
        preLogLike = newLogLike;
    }
    return true;
}

0x04 用最大熵作文本分類

下面代碼來自最大熵用於文本分類，能夠和上面印證學習。

def get_ctgy(fname):#根據文件名稱獲取類別的數字編號
        index = {'fi':0,'lo':1,'co':2,'ho':3,'ed':4,'te':5,
                 'ca':6,'ta':7,'sp':8,'he':9,'ar':10,'fu':11}
        return index[fname[:2]]

def updateWeight():
        #EP_post是 單詞數*類別 的矩陣
        for i in range(wordNum):
            for j in range(ctgyNum):
                EP_post[i][j] = 0.0 #[[0.0 for x in range(ctgyNum)] for y in range(wordNum)]
        # prob是 文本數*類別 的矩陣，記錄每一個文本屬於每一個類別的機率
        cond_prob_textNum_ctgyNum = [[0.0 for x in range(ctgyNum)] for y in range(textNum)]
        #計算p(類別|文本)

        for i in range(textNum):#對每個文本
                zw = 0.0  #歸一化因子
                for j in range(ctgyNum):#對每個類別
                        tmp = 0.0
                       #texts_list_dict每一個元素對應一個文本，該元素的元素是單詞序號：頻率所組成的字典。
                        for (feature,feature_value) in texts_list_dict[i].items():
                            #v就是特徵f(x,y)，非二值函數，而是實數函數，
                            #k是某文本中的單詞，v是該單詞的次數除以該文本不重複單詞的個數。
                                #feature_weight是 單詞數*類別 的矩陣，與EP_prior相同
                                tmp+=feature_weight[feature][j]*feature_value #feature_weight是最終要求的模型參數，其元素與特徵一一對應，即一個特徵對應一個權值
                        tmp = math.exp(tmp)
                        zw+=tmp #zw關於類別求和
                        cond_prob_textNum_ctgyNum[i][j]=tmp                        
                for j in range(ctgyNum):
                        cond_prob_textNum_ctgyNum[i][j]/=zw
        #上面的部分根據當前的feature_weight矩陣計算獲得prob矩陣（文本數*類別的矩陣，每一個元素表示文本屬於某類別的機率），
        #下面的部分根據prob矩陣更新feature_weight矩陣。

        for x in range(textNum):
                ctgy = category[x] #根據文本序號獲取類別序號
                for (feature,feature_value) in texts_list_dict[x].items(): #該文本中的單詞和對應的頻率
                        EP_post[feature][ctgy] += (cond_prob_textNum_ctgyNum[x][ctgy]*feature_value)#認p(x)的先驗機率相同        
        #更新特徵函數的權重w
        for i in range(wordNum):
                for j in range(ctgyNum):
                        if (EP_post[i][j]<1e-17) |  (EP_prior[i][j]<1e-17) :
                                continue                        
                        feature_weight[i][j] += math.log(EP_prior[i][j]/EP_post[i][j])        

def modelTest():
        testFiles = os.listdir('data\\test\\')
        errorCnt = 0
        totalCnt = 0

        #matrix是類別數*類別數的矩陣，存儲評判結果
        matrix = [[0 for x in range(ctgyNum)] for y in range(ctgyNum)]
        for fname in testFiles: #對每一個文件
                lines = open('data\\test\\'+fname)
                ctgy = get_ctgy(fname) #根據文件名的前兩個字符給出類別的序號
                probEst = [0.0 for x in range(ctgyNum)]         #各種別的後驗機率
                for line in lines: #該文件的每一行是一個單詞和該單詞在該文件中出現的頻率
                        arr = line.split('\t')
                        if not words_dict.has_key(arr[0]) : 
                            continue        #測試集中的單詞若是在訓練集中沒有出現則直接忽略
                        word_id,freq = words_dict[arr[0]],float(arr[1])
                        for index in range(ctgyNum):#對於每一個類別
                            #feature_weight是單詞數*類別墅的矩陣
                            probEst[index] += feature_weight[word_id][index]*freq
                ctgyEst = 0
                maxProb = -1
                for index in range(ctgyNum):
                        if probEst[index]>maxProb:
                            ctgyEst = index
                            maxProb = probEst[index]
                totalCnt+=1
                if ctgyEst!=ctgy: 
                    errorCnt+=1
                matrix[ctgy][ctgyEst]+=1
                lines.close()
        #print "%-5s" % ("類別"),
        #for i in range(ctgyNum):
        #    print "%-5s" % (ctgyName[i]),  
        #print '\n',
        #for i in range(ctgyNum):
        #    print "%-5s" % (ctgyName[i]), 
        #    for j in range(ctgyNum):
        #        print "%-5d" % (matrix[i][j]), 
        #    print '\n',
        print "測試總文本個數:"+str(totalCnt)+"  總錯誤個數:"+str(errorCnt)+"  總錯誤率:"+str(errorCnt/float(totalCnt))

def prepare():
        i = 0
        lines = open('data\\words.txt').readlines()
        #words_dict給出了每個中文詞及其對應的全局統一的序號，是字典類型，示例：{'\xd0\xde\xb5\xc0\xd4\xba': 0}
        for word in lines:
                word = word.strip()
                words_dict[word] = i
                i+=1
        #計算約束函數f的經驗指望EP(f)
        files = os.listdir('data\\train\\') #train下面都是.txt文件
        index = 0
        for fname in files: #對訓練數據集中的每一個文本文件
                file_feature_dict = {}
                lines = open('data\\train\\'+fname)
                ctgy = get_ctgy(fname) #根據文件名的前兩個漢字，也就是中文類別來肯定類別的序號

                category[index] = ctgy #data/train/下每一個文本對應的類別序號
                for line in lines: #每行內容：古蹟  0.00980392156863
                        # line的第一個字符串是中文單詞，第二個字符串是該單詞的頻率
                        arr = line.split('\t')
                        #獲取單詞的序號和頻率
                        word_id,freq= words_dict[arr[0]],float(arr[1])

                        file_feature_dict[word_id] = freq
                        #EP_prior是單詞數*類別的矩陣
                        EP_prior[word_id][ctgy]+=freq
                texts_list_dict[index] = file_feature_dict
                index+=1
                lines.close()
                
def train():
        for loop in range(4):
            print "迭代%d次後的模型效果：" % loop
            updateWeight()
            modelTest()


textNum = 2741  # data/train/下的文件的個數
wordNum = 44120 #data/words.txt的單詞數，也是行數
ctgyNum = 12


#feature_weight是單詞數*類別的矩陣
feature_weight = np.zeros((wordNum,ctgyNum))#[[0 for x in range(ctgyNum)] for y in range(wordNum)]

ctgyName = ['財經','地域','電腦','房產','教育','科技','汽車','人才','體育','衛生','藝術','娛樂']
words_dict = {}

# EP_prior是個12（類）* 44197（全部類的單詞數）的矩陣，存儲對應的頻率
EP_prior = np.zeros((wordNum,ctgyNum))
EP_post = np.zeros((wordNum,ctgyNum))
#print np.shape(EP_prior)
texts_list_dict = [0]*textNum #全部的訓練文本 
category = [0]*textNum        #每一個文本對應的類別

print "初始化:......"
prepare()
print "初始化完畢，進行權重訓練....."
train()

0x05 最大熵馬爾科夫模型 MEMM

1. 建模公式

最大熵馬爾科夫模型能夠看做是HMM與log-linear model結合的一種方式。MEMM利用判別式模型的特色，直接對每個時刻的狀態創建一個分類器，而後將全部的分類器的機率值連乘起來。

與HMM的 o 依賴 i 不同，MEMM的當前狀態依賴於前一狀態與當前觀測；MEMM當前隱藏狀態 i 應該是依賴當前時刻的觀測節點 o 和上一時刻的隱藏節點 i-1。因此MEMM中，給定觀測序列 i1,...in 後，某個狀態序列 in 的條件機率是能夠直接學習的：
\[ P(I|O) = \prod_{t=1}^n P(i_i|i_{i-1},o_i), i = 1, ...,n \]
而且P(i|i',o)這個機率經過最大熵分類器建模（取名MEMM的緣由）:
\[ P(i|i^{''},o) = \frac{1}{Z(o,i^{'})}exp(\sum_a)λ_af_a(o,i) \]
Z(o,i) 這部分是歸一化；f(o,i) 是特徵函數，具體點，這個函數是須要去定義的。MEMM模型就是可以直接容許「定義特徵」。λ是特徵函數的權重，這是個未知參數，須要從訓練階段學習而得。

因此整體上，MEMM的建模公式這樣：
\[ P(I|O) = \prod_{t=1}^n \frac{exp(\sum_a)λ_af_a(o,i)}{Z(o,i^{'}_{i-1})}, i = 1, ...,n \]
公式這部分之因此長成這樣，是由ME模型決定的。

使用維特比算法進行解碼時，MEMM對應的公式爲:
\[ ν_t(j)= max_i^n ν_{t-1}(i)∗P(y_j∣y_i,x_t) \\ 1≤j≤n,1<t<T \]
請務必注意，理解判別模型和定義特徵兩部分含義，這已經涉及到CRF的雛形了。

完整流程:

step1. 先預約義特徵函數 f(o,i)
step2. 在給定的數據上，訓練模型，肯定參數，即肯定了MEMM模型
step3. 用肯定的模型作序列標註問題或者序列求機率問題。

2. MEMM vs HMM

HMM模型是對狀態轉移機率（狀態-狀態）和發射機率（狀態-觀察）直接建模，統計共現機率。
MEMM模型是對狀態轉移機率和發射機率創建聯合機率，統計的是條件機率。但MEMM容易陷入局部最優，是由於MEMM只在局部作歸一化。

舉個例子，對於一個標註任務，「我愛北京天安門「，

標註爲" s s b e b c e"

對於HMM的話，其判斷這個標註成立的機率爲 P= P(s轉移到s)P('我'表現爲s) P(s轉移到b)P('愛'表現爲s) ...*P().訓練時，要統計狀態轉移機率矩陣和表現矩陣。

對於MEMM的話，其判斷這個標註成立的機率爲 P= P(s轉移到s|'我'表現爲s)P('我'表現爲s) P(s轉移到b|'愛'表現爲s)P('愛'表現爲s)..訓練時，要統計條件狀態轉移機率矩陣和表現矩陣。

3. 代碼對比

如下代碼來自 https://github.com/alexkartun/Sequence-Tagging

HMMTag

下面是 HMM 的代碼，計算了Emission / Transaction。

def predict_viterbi(x, trans_c, emiss_c, word_tags_map, interp_weights):
    """
    For each word in vector x predict his tag. Prediction is done by viterbi algorithm. Check all tags options/globally.
    :param x: X vector of words.
    :param trans_c: Transaction counts.
    :param emiss_c: Emission counts.
    :param word_tags_map: Word to tags map.
    :param interp_weights: Interpolation weights.
    :return: Vector of all tags respectively to words in vector x.
    """
    y = []
    v = [{(mle_train.START_SYMBOL, mle_train.START_SYMBOL): 0.0}]
    bp = []
    for ind, word in enumerate(x):
        # Convert word if it was'nt seen in the corpus, to signature word.
        if word not in word_tags_map:
            word = mle_train.subcategorize(word)
        # Pruning of tags to lower amount of possible tags for this word.
        available_tags = word_tags_map[word]

        max_score = {}
        max_tags = {}
        # Calculate for each word best scores/probabilities and best tags for each word.
        for pp_t, p_t in v[ind]:
            for curr_tag in available_tags:
                score = get_score(word, curr_tag, p_t, pp_t, trans_c, emiss_c, interp_weights)
                score += v[ind][(pp_t, p_t)]

                if (p_t, curr_tag) not in max_score or score > max_score[(p_t, curr_tag)]:
                    max_score[(p_t, curr_tag)] = score
                    max_tags[(p_t, curr_tag)] = pp_t

        v.append(max_score)
        bp.append(max_tags)
    # Calculate last 2 best tags.
    max_score = float("-inf")
    prev_last_tag, last_tag = None, None
    for prev_t, curr_t in v[len(x)]:
        score = v[len(x)][(prev_t, curr_t)]
        if score > max_score:
            max_score = score
            last_tag = curr_t
            prev_last_tag = prev_t
    y.append(last_tag)
    if len(x) > 1:
        y.append(prev_last_tag)

    prev_t = last_tag
    prev_prev_t = prev_last_tag
    # By backtracking extract all the path of best tags for each word starting by last 2 tags we calculated above.
    for i in range(len(v) - 2, 1, -1):
        curr_t = bp[i][(prev_prev_t, prev_t)]
        y.append(curr_t)
        prev_t = prev_prev_t
        prev_prev_t = curr_t
    # Reverse the path.
    y = reversed(y)
    return y
  
def get_score(word, curr_tag, p_t, pp_t, trans_c, emiss_c, interp_weights):
    """
    Calculate probability. Prob = e_prob + q_prob.
    :param word: Curr word.
    :param curr_tag: Curr tag.
    :param p_t: Previous tag.
    :param pp_t: Previous of previous tag.
    :param trans_c: Transaction counts.
    :param emiss_c: Emission counts.
    :param interp_weights: Interpolation weights.
    :return:
    """
    e = mle_train.get_e(word, curr_tag, emiss_c, trans_c)
    q = mle_train.get_q(pp_t, p_t, curr_tag, trans_c, interp_weights)
    if e != 0 and q != 0:
        score = math.log(e) + math.log(q)
    else:
        score = float('-inf')
    return score  
  
def get_e(word, tag, emiss_c, trans_c):
    """
    Calculate e probability.
    :param word: Current word to analyze.
    :param tag: Current tag to analyze.
    :param emiss_c: Emission counts.
    :param trans_c: Transaction counts.
    :return: Probability value.
    """
    key = '{} {}'.format(word, tag)
    count_word_tag = float(emiss_c.get(key, 0))
    key = '{}'.format(tag)
    count_tag = float(trans_c.get(key, 0))
    try:
        e_prob = count_word_tag / count_tag
    except ZeroDivisionError:
        e_prob = 0

    return e_prob


def get_q(pp_t, p_t, curr_tag, trans_c, interpolation_weights):
    """
    Calculate q probability.
    :param pp_t: Previous of previous tag to analyze.
    :param p_t: Previous tag to analyze.
    :param curr_tag: Current tag to analyze.
    :param trans_c: Transaction counts.
    :param interpolation_weights: Interpolation weight.
    :return:
    """
    lambda_1 = interpolation_weights[0]
    lambda_2 = interpolation_weights[1]
    lambda_3 = interpolation_weights[2]
    # Calculate trigram prob = count_trigram_abc / count_bigram_ab
    key = '{} {} {}'.format(pp_t, p_t, curr_tag)
    count_trigram_abc = float(trans_c.get(key, 0))
    key = '{} {}'.format(pp_t, p_t)
    count_bigram_ab = float(trans_c.get(key, 0))
    try:
        trigram_prob = count_trigram_abc / count_bigram_ab
    except ZeroDivisionError:
        trigram_prob = 0
    # Calculate bigram prob = count_trigram_bc / count_unigram_b
    key = '{} {}'.format(p_t, curr_tag)
    count_bigram_bc = float(trans_c.get(key, 0))
    key = '{}'.format(p_t)
    count_unigram_b = float(trans_c.get(key, 0))
    try:
        bigram_prob = count_bigram_bc / count_unigram_b
    except ZeroDivisionError:
        bigram_prob = 0
    # Calculate unigram prob = count_unigram_c / num_words
    key = '{}'.format(curr_tag)
    count_unigram_c = float(trans_c.get(key, 0))
    key = '{}'.format(NUM_WORDS_SYMBOL)
    num_words = float(trans_c.get(key, 0))
    try:
        unigram_prob = count_unigram_c / num_words
    except ZeroDivisionError:
        unigram_prob = 0
    # Apply interpolation weight on the probabilities.
    interpolation = lambda_1 * trigram_prob + lambda_2 * bigram_prob + lambda_3 * unigram_prob
    return interpolation

MEMMTag

下面是MEMM的代碼，採用LibLinear作邏輯迴歸(最大熵)。

def predict_viterbi(x, f_map, tags_s, word_t_map, lib_model):
    """
    For each word in vector x predict his tag. Prediction is done by viterbi algorithm. Check all tags options/globally.
    :param x: X vector of all words to be tagged.
    :param f_map: Features map.
    :param tags_s: Tags set.
    :param word_t_map: Word to tags map.
    :param lib_model: Liblinear model.
    :return: Return best predicted list of tags for each word in x respectively.
    """
    y = []
    v = [{(extract.START_SYMBOL, extract.START_SYMBOL): 0.0}]
    bp = []
    for ind, word in enumerate(x):
        # Check if word was seen in the corpus.
        if word not in word_t_map:
            is_rare = True
            available_tags = tags_s
        else:
            is_rare = False
            # Pruning of tags to lower amount of possible tags for this word.
            available_tags = word_t_map[word]

        max_score = {}
        max_tags = {}
        # Calculate for each word best scores/probabilities and best tags for each word.
        for pp_t, p_t in v[ind]:
            for curr_tag in available_tags:
                # 一個word對應的feature列表，可能的tag列表
                word_features = extract.generate_word_features(is_rare, p_t, pp_t, word, ind, x)
                features_vec = features_to_vec(word_features, f_map)
                scores = lib_model.predict(features_vec)
                score = np.amax(scores)
                if (p_t, curr_tag) not in max_score or score > max_score[(p_t, curr_tag)]:
                    max_score[(p_t, curr_tag)] = score
                    max_tags[(p_t, curr_tag)] = pp_t

        v.append(max_score)
        bp.append(max_tags)
    # Calculate last 2 best tags.
    max_score = float("-inf")
    prev_last_tag, last_tag = None, None
    for prev_t, curr_t in v[len(x)]:
        score = v[len(x)][(prev_t, curr_t)]
        if score > max_score:
            max_score = score
            last_tag = curr_t
            prev_last_tag = prev_t

    y.append(last_tag)
    if len(x) > 1:
        y.append(prev_last_tag)

    prev_t = last_tag
    prev_prev_t = prev_last_tag
    # By backtracking extract all the path of best tags for each word starting by last 2 tags we calculated above.
    for i in range(len(v) - 2, 1, -1):
        curr_t = bp[i][(prev_prev_t, prev_t)]
        y.append(curr_t)
        prev_t = prev_prev_t
        prev_prev_t = curr_t
    y = reversed(y)
    return y

提取特徵

def load_model(feature_map_fname):
    """
    Load the model from features map file.
    :param feature_map_fname: Features map file name.
    :return: Tags set, index to tag map, word to tags map and features map.
    """
    features_map = defaultdict(int)
    ind_to_tags_map = defaultdict(str)
    tags_set = set()
    word_to_tags_map = defaultdict(set)

    with open(feature_map_fname, 'r') as f:
        all_data = f.readlines()
    # Booleans to separate the reading process from the file.
    is_features_section = False
    is_words_to_tags_section = False
    is_tags_section = True
    for ind, line in enumerate(all_data):
        line = line.strip()
        if line == '^':
            is_words_to_tags_section = True
            is_tags_section = False
            continue
        if line == '^^':
            is_features_section = True
            is_words_to_tags_section = False
            continue
        if is_tags_section:
            key, value = line.split()
            tags_set.add(key)
            ind_to_tags_map[int(value)] = key
        if is_words_to_tags_section:
            key, value = line.split('^')
            word_to_tags_map[key] = value.split() # 該word可能對應的tag: 動詞，名詞......
        if is_features_section:
            key, value = line.split()
            features_map[key] = int(value) # load事先定義好的特徵，有 string 和 int兩種表示方法
    return features_map, tags_set, ind_to_tags_map, word_to_tags_map

# ExtractFeatures.py, MMEMTag.predict_viterbi會調用
def generate_word_features(is_rare, p_t, pp_t, curr_word, word_ind, words):
    """
    Generate for current word list of features as listed on table one.
    :param is_rare: Boolean that corresponds to type of the word. Rare or not.
    :param p_t: Previous tag.
    :param pp_t: Previous of previous tag.
    :param curr_word: Current word.
    :param word_ind: Current word index.
    :param words: List of all words in the sentence.
    :return: List of all word features.
    """
    word_features = []
    # Check if word is rare.
    if is_rare:
        # Generate the suffixes and prefixes depends on min of (word length or 4).
        for i in range(1, min(5, len(curr_word))):
            word_features.append('prefix' + str(i) + '=' + curr_word[:i]) # 前面的字
        for i in range(1, min(5, len(curr_word))):
            word_features.append('suffix' + str(i) + '=' + curr_word[-i:]) # 後面的字
        # Check with regex if word contains digit.
        if re.search(r'\d', curr_word):
            word_features.append('has_digit=true') # 是否是數字
        # Check with regex if word contains upper case letter.
        if re.search(r'[A-Z]', curr_word):
            word_features.append('has_upper_letter=true') # 大寫開頭
        # Check if word contains hyphen symbol.
        if '-' in curr_word: 
            word_features.append('has_hyphen=true') # 有破折號
    else:
        # Otherwise word is not rare and this word as feature.
        key = 'form={}'.format(curr_word) 
        word_features.append(key)
    # Generate previous tags.
    key = 'pt={}'.format(p_t)
    word_features.append(key)
    key = 'ppt={}^{}'.format(pp_t, p_t)
    word_features.append(key)
    # Generate next words and words that appeared before in the sentence.
    if word_ind > 0:
        key = 'pw={}'.format(words[word_ind - 1])
        word_features.append(key)
    if word_ind > 1:
        key = 'ppw={}'.format(words[word_ind - 2])
        word_features.append(key)
    if word_ind < len(words) - 1:
        key = 'nw={}'.format(words[word_ind + 1])
        word_features.append(key)
    if word_ind < len(words) - 2:
        key = 'nnw={}'.format(words[word_ind + 2])
        word_features.append(key)
    return word_features # word對應的，用string表示的，特徵名字列表

# MEMMTag.py
def features_to_vec(word_features, f_map):
    """
    Convert word features to vector of features indexes whose liblinear model can understand.
    :param word_features: Word features to convert.
    :param f_map: Features map.
    :return: Vector of features.
    """
    features_vec = []
    for feature in word_features:
        if feature in features_map:
            features_vec.append(f_map[feature])

    features_vec = sorted(features_vec)
    return features_vec

邏輯迴歸

LibLinear，這是一個簡單的求解大規模規則化線性分類和迴歸的軟件包。

nr_class:類的個數，若是是迴歸，那麼這個數是2.
nr_feature: 訓練數據的樣本維數（不包括bias項）。
bias: 若是bias >= 0，咱們會在每一個樣本的最後添加一個額外的特徵。
Label: 每一個類的標籤，對迴歸來講，爲空。
w: 一個 nr_w x n權值矩陣。n是nr_feature（特徵維數）或者nr_feature+1（存在bias項時）。若是nr_class=2，而且-s不是4（不是Crammer and Singer的多類SVM），那麼nr_w是1，對於其餘狀況，nr_w等於nr_class。
w 在這裏就是 特徵維數 x 類別數，就對應了前面提到的 "在最大熵模型的視角下，每條輸入的n個特徵與k個類別共同組成了nk個特徵，模型中有nk個權重，與特徵一一對應。每一個類別會觸發nk個特徵中的n個"。

具體代碼：

class LiblinearLogregPredictor(object):
  def __init__(self, model_file_name):
        # dict of feat-num -> numpy array, where array is the per-class values
        self.weights = {}
        with open(model_file_name) as fh:
            solver_type = fh.next().strip()
            nr_classes = fh.next().strip()
            label = fh.next().strip()
            nr_feature = fh.next().strip()
            bias = fh.next().strip()
            w = fh.next().strip()
            assert(w == "w")
            assert(nr_classes.startswith("nr_class"))
            assert(nr_feature.startswith("nr_feature"))
            assert(label.startswith("label"))
            nr_classes = int(nr_classes.split()[-1])
            for i, ws in enumerate(fh, 1):
                ws = [float(x) for x in ws.strip().split()]
                # skip unused features
                if all([x == 0 for x in ws]):
                    continue
                assert len(ws) == nr_classes, "bad weights line %s" % ws
                self.weights[i] = np.array(ws)
            self.bias = float(bias.split()[-1])
            self.labels = label.split()[1:]

    # feature_ids是int類型的特徵列表，weights是權值矩陣
    def predict(self, feature_ids): 
        weights = self.weights
        scores = np.zeros(len(self.labels))
        for f in feature_ids:
            if f in weights:
                scores += weights[f]
        scores = 1/(1+np.exp(-scores))
        scores = scores/np.sum(scores)
        return scores

0x06 水滸傳的繼續應用

嘮叨了這麼半天，總算到實例應用了。仍是以"梁中書突圍大名府爲例"。前文參見[白話解析]以水滸傳爲例學習隱馬爾可夫模型。

隱馬爾可夫模型模型比較簡單，可是現實中狀況可能更復雜，好比梁中書突圍遇到了宋江，他再次選擇從宋江處突圍的可能性會變低，由於宋江身邊確定是梁山大部分好漢，突圍難度太大。可是若是遇到史進，危險性就沒有那麼大了。

如今算法已經有了，其實就看如何構建特徵函數了。從前面的代碼能看出來，獲取word對應的feature的操做以下

word_features = extract.generate_word_features(is_rare, p_t, pp_t, word, ind, x)
features_vec = features_to_vec(word_features, f_map)
scores = lib_model.predict(features_vec)

咱們創建各類相關信息

*********************** 特徵函數
f_1 = 是否是梁山總頭領
f_2 = 是否是打虎英雄
f_3 = 是否是五虎將
f_4 = 是否是八驃騎

*********************** 下一次突圍的門（ 對應詞性標註中的tag ）
t_1 = 逆時針上一門
t_2 = 當前門
t_3 = 順時針下一門
t_3 = 對門
  
*********************** 好漢（ 對應詞性標註中的word ） 
h_1 = 林沖
h_2 = 武松
h_3 = 史進
h_3 = 宋江
  
*********************** 碰見好漢的序列（ 對應詞性標註中的句子 ）  
武松，宋江，史進，林沖

權重矩陣特徵函數維度 x 門類別數目，這個數值屬於隨意炮製的
\[ 權重矩陣 = \left\{ \begin{matrix} \frac12 & \frac14 & \frac18 & \frac18\\ \frac14 & \frac12 & \frac18 & \frac18\\ \frac18 & \frac18 & \frac12 & \frac14\\ \frac18 & \frac18 & \frac14 & \frac12\\ \end{matrix} \right\} \]
能夠帶入算法進行計算。

好吧，其實此次的應用說明比較少 ^_^。