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

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

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

0x00 摘要

在機器學習過程中，會遇到很多晦澀的概念，相關數學公式很多，大家理解起來很有困難。遇到類似情況，我們應該多從直覺角度入手思考，用類別或者舉例來附會，這樣往往會有更好的效果。

我在講解論述過程中給自己的要求是：在生活中或者名著中找一個例子，然后用自己的話語闡述出來。

0x01 背景

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

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

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

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

HMM是一種生成模型，定義了聯合概率分布，其中y和x分別表示觀察序列和相對應的標注序列的隨機變量。為了能夠定義這種聯合概率分布，生成模型需要枚舉出所有可能的觀察序列，這在實際運算過程中很困難的，因此我們需要將觀察序列的元素看作是彼此孤立的個體，即假設每個元素彼此獨立，任何時刻的觀察結果只依賴于該時刻的狀態。

目標函數和預測目標函數不匹配，HMM學到的是狀態和觀察序列的聯合分布P(Y,X)，而預測問題中，我們須要的是條件概率P(Y|X)。HMM必須計算出所有的潛在可能路徑的概率值大小（然后再挑概率值最大的那一個作為最終結果）

HMM僅僅依賴于每個狀態和它相應的觀察對象，但是對于某一個觀測值，它可能并非由一個隱藏狀態決定的，而是兩個以上的隱藏狀態綜合作用而產生的，那么這時HMM就無能為力了。

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\} \]

可以帶入算法進行計算。

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