序言

…..
本系列對演算法的講解都會從兩篇部分予以呈現：

a. 溼貨部分要淺入淺出，形象生動，讀得明白。
b. 乾貨部分要一文以蔽之，公式羅列，看得通透；

下面是（一）的 b 部分內容

Unigram

Unigram模型認為序列中的每一項都是獨立發生的，所以很自然，假設我們有N個序列，每個序列長度是MnM_n，那麼整個序列的聯合概率分佈就是：

P(X) = \prod_i^N\prod_j^{M_i}p\left(X_i^j\right)

如果xjix_i^j的取值是有限的，屬於集合VV,上式可以轉化為

P(X) = \prod_vp(x=v)^{count(v)}

其中count(v)count(v)表示vv在序列中的出現次數，上式就是Unigram模型的似然函式，並且它的對數似然函式為：

L = lnP(X) = \sum_vcount(v)*ln^{p(x=v)}\\
並且 \sum_vp\left(v\right) = 1

我們將限制條件通過拉格朗日運算元加到對數似然函式裡並分別對引數求導：

\begin{align}
L’ = lnP(X) &= \sum_vcount(v)*ln^{p(v)} \lambda(1-\sum_vp\left(v\right))\\
\frac{\partial L’}{\partial p(v)} &= \frac{count(v)}{p(v)}-\lambda\\
\frac{\partial L’}{\partial \lambda} &= 1-\sum_vp\left(v\right)
\end{align}

令兩個偏導為0，我們先有

p(v)=\frac{count(v)}{\lambda} 並且\\
\sum_vp\left(v\right) = \sum_v\frac{count(v)}{\lambda} =1\\
所以\lambda = count(*)

所以Unigram模型的最大似然估計就是

p(v) = \frac{count(v)}{count(*)}

這是一個非常符合直覺的公式，一個term出現的概率就是它在訓練資料裡的頻率，所以甚至有人會覺得這就是常識，但實際上背後是有數學推導支撐的。

馬爾可夫模型

按照上面類似的邏輯，一階馬爾可夫的最大似然估計就是

p(x^{t}=v_i,x^{t 1}=v_j) = \frac{count(v_i,v_j)}{count(v_i,*)}

這裡不作推導了,有興趣可以自己推導一下。（多說一句，很多書上上面這兩個公式都是直接給出來的，但你自己從頭推導過一次後會感覺概率論這東西還是靠譜的……233）

隱馬爾可夫模型HMM

根據我們a部分的描述，隱馬爾可夫鏈其實是有兩條序列的，我們分別用HH和OO來表示，以一階隱馬爾可夫鏈為例：

那麼對於某個訓練樣本{o,h}\{o,h\}它的聯合概率應該如何表示呢

P(o,h|\theta) = p(h_1|\theta)\prod_{t=2}^Tp(h_t|h_{t-1},\theta)\prod_{t=1}^Tp(o_t|h_t,\theta)

主要是三個部分：1)第一個隱藏層的初始狀態。2)隱藏層的轉移概率。3)從隱藏層到表示層的表現概率，我們分別用特殊的符號來標記對應的引數：

\begin{align}
\pi_{q_i} &= p(h_1=q_i),q \in Q,Q是所有可能的隱藏層狀態\\
A_{q_i,q_j} &= p(h_t=q_j|h_{t-1}=q_i)\\
B_{q_j,s_k} &= p(o_t=s_k|h_t = q_j),s_k \in S,S是所有可能的表現層狀態
\end{align}

所以上面的聯合概率公式可以表示為：

P(o,h|\theta) = \pi_{h_1}\prod_{t=2}^TA_{h_{t-1}h_t}\prod_{t=1}^TB_{h_to_t}

如果我們同時知道{o,h}\{o,h\}，想要求得θ\theta，那其實問題很簡單，無非就是類似Unigram和Bigram統計一下各個term和pair出現的次數，並計算出概率。

但問題的關鍵來了：在實際中我們多半是不知道隱藏層狀態是什麼樣的，否則它們就不叫「隱藏層」了，所以我們只能期望表現序列oo的概率最大：

p(o|\theta) = \sum_hp(o,h|\theta)

記住h可是一個序列，在任何一個位置的取值都有|Q||Q|種，所以對於一個長度為N的序列，我們需要加和的隱藏層候選集大小為|Q|T|Q|^T（T是序列長度），這是一個極其龐大的數字。而且當hh的狀態不同時，從hh表現為oo的概率也不一樣。所以我們想直接去最大化這個聯合分佈（最大似然）是不現實的。

Baum-Welch演算法（EM演算法）

所以在這裡只能通過EM演算法，去一步步迭代計算求得θ\theta，EM演算法首先假設引數是知道的，然後最大化Q函式來更新引數值（Baum-Welch是EM演算法在HMM中的具體實現）
EM演算法作為十大資料探勘演算法之一，後面我們可以單獨開一篇部落格學習，我自己現在也沒有理解得很深入，但在這裡你只需要知道它有兩步，Q&M，Q步驟中計算Q函式，M步驟中最大化Q函式。
在HMM問題中，θold是已知的，θ是我們需要去求解的\theta^{old}是已知的，\theta是我們需要去求解的，第一個公式裡可以認為是Q函式的定義：

\begin{align}
Q(\theta,\theta^{old}) &= \sum_hp(o,h|\theta^{old})ln\ p(o,h|\theta)\\
&= \sum_h{p(h|o,\theta^{old})}{p(o|\theta^{old})}ln\ p(o,h|\theta)\\
p(o|\theta^{old})和\theta無關，所以：\\
&\propto \sum_hp(h|o,\theta^{old})ln\ p(o,h|\theta)\\
帶入前面的p(o,h|\theta)：\\
& = \sum_h p(h|o,\theta^{old})*ln\ \pi_{h_1} \\&\ \ \ \ \sum_h p(h|o,\theta^{old})\sum_{t=2}ln\ A_{h_{t-1},h_t}\\&\ \ \ \ \sum_h p(h|o,\theta^{old})\sum_{t=1}ln\ B_{h_t,o_t}
\end{align}

公式複雜度似乎有些不可控了，我們單把上式的第一部分拎出來：

F = \sum_h p(h|o,\theta^{old})*ln\ \pi_{h_1}

假設我們指定h1=qh_1=q
則F中和h1=qh_1=q 相關的部分就是

\sum_{h \in H’_q} p(h|o,\theta^{old})*ln\ \pi_{q}

H′qH’_q是HH中h1=qh_1=q的所有可能的hh序列構成的集合，此時∑h∈H′p(h|o,θold)\sum_{h \in H’} p(h|o,\theta^{old})發生了微妙的變化：

\sum_{h \in H’} p(h|o,\theta^{old}) = p(h_1=q|o,\theta^{old})

所以原式F：

F = \sum_qp(h_1=q|o,\theta^{old})ln\ \pi_q

類似的方法我們可以吧Q(θ,θold)Q(\theta,\theta^{old})變成下式：

\begin{align}
Q(\theta,\theta^{old}) =& \sum_qp(h_1=q|o,\theta^{old})ln\ \pi_q \\
& \sum_{t=2}\sum_{q_1 \in Q}\sum_{q_2 \in Q}p(h_{t-1}=q_1,h_{t}=q_2|o,\theta^{old})*lnA_{q_1,q_2} \\
& \sum_{t}\sum_qp(h_t=q|o,\theta^{old})ln\ B_{q,o_t}
\end{align}

為了方便後續計算表示，我們用：

\begin{align}
\gamma_t &= p(h_t|o,\theta^{old})\\
\gamma_t(q) &= p(h_t=q|o,\theta^{old})\\
\xi_t &= p(h_{t-1},h_t|o,\theta^{old})\\
\xi_t(q_1,q_2) &= p(h_{t-1}=q_1,h_t=q_2|o,\theta^{old})
\end{align}

所以：

Q(\theta,\theta^{old}) = \sum_q\gamma_1(q)ln\ \pi_q \sum_{t=2}\sum_{q_1 \in Q}\sum_{q_2 \in Q}\xi_t(q_1,q_2)*lnA_{q_1,q_2} \sum_{t}\sum_q\gamma_t(q)ln\ B_{q,o_t}

然後我們最大化Q，分別對每個引數進行求導，將限定（主要是各種概率和為1）以拉格朗日乘子的方式加入，比如對πq\pi_q求導：

\begin{align}
\frac{\partial Q}{\partial \pi_q} &= \frac{\gamma_1(q)}{\pi_q} – \lambda。使它=0，得到\\
\pi_q &= \frac{\gamma_1(q)}{\lambda}\\
\sum_q\pi_q &=\frac{\sum_q\gamma_1(q)}{\lambda}=1\\
所以：\pi_q &= \frac{\gamma_1(q)}{\sum_{q’}\gamma_1(q’)}=\gamma_1(q)，\\（因為\gamma_1本身&也是一個概率分佈，所以其和也為1）
\end{align}

類似的方式可以分別求得：

\begin{align}
A_{q_1,q_2} &= \frac{\sum_{t=2}\xi_t(q_1,q_2)}{\sum_{t=2}\sum_q\xi_t(q,q_2)}=\frac{\sum_{t=2}\xi_t(q_1,q_2)}{\sum_{t=2}\gamma_t(q_2)}\\
B_{q,s} &= \frac{\sum_{t \in T’}\gamma_t(q)}{\sum_{t}\gamma_t(q)},T’是所有o_t=s的集合
\end{align}

上面三個公式就是我們在EM演算法的M步驟中更新引數值的公式。
但在這裡，我們遺留了一個很重要的問題，γ和ξ\gamma和\xi應該如何計算呢，我們只是為了簡化公式引入的這兩個符號，它們背後還是概率公式啊。

前向後向演算法（forward-backward）

所以接下來我們就要解決這個問題，就是如何計算γ\gamma和ξ\xi,也就是在知道了表現序列oo和相關引數θold\theta^{old}後，想要知道和hh相關的一些概率。

上面是一個示意圖，|Q|=3|Q|=3，每個隱藏狀態的取值都有3種情況，所以下面就變成了那樣一個全連線的網路，看起來有點像一個NN了。
有（下面公式省略θold\theta^{old}）：

\begin{align}
\gamma_t &= p(h_t|o).\small(來自上面的定義)\\
&= \frac{p(o|h_t)p(h_t)}{p(o)}\\
&=\frac{p(o_1,o_2…,o_t|h_t)p(o_{t 1}..o_{T}|h_t)p(h_t)}{p(o)}\\
&=\frac{p(o_1,o_2…,o_t,h_t)p(o_{t 1}..o_{T}|h_t)}{p(o)}
\end{align}

上式分母中的部分其實跟θold\theta^{old}沒啥關係，我們把分子的兩部分分別定義為α\alpha和β\beta:

\begin{align}
\alpha(h_t) &= p(o_1,o_2…,o_t,h_t)\\
\beta(h_t) &= p(o_{t 1}..o_{T}|h_t)
\end{align}

其中

\begin{align}
\alpha(h_t)^q &= p(o_1,o_2…,o_t,h_t=q)\\
&= p(o_{1…t}|h_t=q)p(h_t=q)\\
&= p(o_t|h_t=q)p(o_{1…(t-1)}|h_t=q)p(h_t=q)\\
&= p(o_t|h_t=q)p(o_{1…(t-1)},h_t=q)\\
&= p(o_t|h_t=q)\sum_{q’}p(o_{1…(t-1)},h_{t-1}=q’,h_t=q)\\
在這裡h_t只和&h_{t-1}有關係：\\
&=p(o_t|h_t=q)\sum_{q’}p(o_{1…(t-1)},h_{t-1}=q’)p(h_t=q|h_{t-1}=q’)\\ &=p(o_t|h_t=q)\sum_{q’}\alpha(h_{t-1})^{q’}p(h_t=q|h_{t-1}=q’)
\end{align}

上面推導過程用到的性質無非是馬爾可夫假設和條件概率的概念而已，此時我們已經發現α\alpha是可以迭代進行計算的了，整個計算過程是從前往後的，所以是前向演算法，其實也很好理解，比如下圖的a中標紅的節點h3=q2h_3=q_2，其實只能從h2h_2幾個狀態轉移過來，箭頭的引數值就是p(h3=q2|h2=q)p(h_3=q_2|h_{2}=q),也就是A可以指示的。
對於序列的第一個隱藏狀態，它的初始值就是

\alpha(h_1)^q = \pi_q*p(o_1|q)

類似的方法我們可以推導得到β\beta,這是後向演算法的部分,示意圖如下面圖(b)所示：

\begin{align}
\beta(h_t)^q &= p(o_{t 1}..o_{T}|h_t=q)\\
&=\sum_{q’}\beta(h_{t 1})^{q’}p(o_{t 1}|h_{t 1}=q’)p(h_{t 1}=q’|h_t=q)
\end{align}

現在回頭去看就有

\begin{align}
\gamma_t(q) &=\frac{\alpha(h_t)^q\beta(h_t)^q}{p(o)}\\
其中:p(o) &= \sum_q’ \alpha(h_t)^{q’}\beta(h_t)^{q’}.\small(因為\gamma_t是概率分佈）
\end{align}

而

\begin{align}
\xi_t(q_1,q_2) &= p(h_{t-1}=q_1,h_t=q_2|o)
\\
&=\frac{p(o|h_{t-1}=q_1,h_t=q_2)p(h_{t-1}=q_1,h_t=q_2)}{p(o)}\\
&=\frac{p(o_{1…t-1}|h_{t-1}=q_1)p(o_t|h_t=q_2)p(o_{t 1…T}|h_t=q_2)p(h_{t-1}=q_1|h_t=q_2)p(h_t=q_2)}{p(o)}
\\
&\small{第一項和最後一項組合成\alpha,第三項是\beta}
\\
&=\frac{\alpha(h_{t-1})^{q_1}p(o_t|h_t=q_2)\beta(h_t)^{q_2}p(h_{t-1}=q_1|h_t=q_2)}{p(o)}\\
&=\frac{\alpha(h_{t-1})^{q_1}B_{q_2,o_t}\beta(h_t)^{q_2}A_{q_1,q_2}}{p(o)}\\
其中p(o) &= \sum_{q}\sum_{q’}\alpha(h_{t-1})^{q}B_{q’,o_t}\beta(h_t)^{q’}A_{q,q’}.因為\xi也是概率分佈。
\end{align}

可以看到，前向後向演算法我們一開始的意圖已經達到了，而實際上前向後向演算法也可以用來計算出當所有引數給定時，o的概率，這在實際應用中我們判斷一個序列出現的概率時非常有用:

\begin{align}
p(o) &= \sum_q’ \alpha(h_t)^{q’}\beta(h_t)^{q’}\\
\small 或者\\
p(o) &= \sum_{q}\sum_{q’}\alpha(h_{t-1})^{q}B_{q’,o_t}\beta(h_t)^{q’}A_{q,q’}
\end{align}

至此，我們已經得到γ\gamma和ξ\xi了，所以完全可以帶回到EM演算法去迭代地學習到模型中的引數集合θ={A,B,π}\theta=\{A,B,\pi\}了。但在實際的應用中，我們通常希望給出一個觀測序列oo，可以知道對應概率最大的hh是什麼，因為hh有可能是一些有意義的狀態比如單詞詞性、語音詞之類的東西，在溼貨部分的例子裡就是我們希望通過表現外在行為去猜測其最大可能的內在意願，這也是我們人在日常生活中潛意識會幹的事情。

Viterbi演算法

如何解決上一節提到的問題呢，最簡單的方法，我們不是知道γt(q)\gamma_t(q)是t時刻隱狀態為q的概率嗎，對於每一個時刻我們找一個最大的q不就好了。也就是

q_t^* = \mathop{\arg\max}_{q_i}[\gamma_t(q_i)]

這樣的做法被稱為「近似演算法」，但我們知道，所有可能的序列有|Q|T|Q|^T種，近似演算法不能保證求得的解是整體序列概率最大的那一個。說白了，我們希望找到的是下面概率圖（c）（注意此時的連線已經不是轉移概率了）中的一個最優路徑（概率累乘最大）

一提到最優路徑，做過OI/ACM題目的同學可能都會想到動態規劃，沒錯，這裡我們就要通過動態規劃來求解，在序列演算法中被稱為Viterbi演算法。
首先我們定義一個概率計算函式ft(q)f_t(q)表示在ht=qh_t=q的前提下從1到t所有序列中最大概率值：

f_t(q) = \mathop{\max}_{h_{1…t-1}}p(h_t=q,o_{1…t})

在上面那個圖(d)中，我們想求到f3(q1)f_3(q_1),你會發現它只能從前一個隱藏狀態的三種可能性轉移過來，我們只需要看哪個q′q’的ft−1(q′)∗Aq′,qf_{t-1}(q’)*A_{q’,q}更大就好了，至於ft−1(q′)f_{t-1}(q’)怎麼求解我並不關心，交給遞迴就好了，這就是遞迴演算法無後效性的體現。

f_t(q) = \left(\mathop{\max}_q’f_{t-1}(q’)*A_{q’,q}\right)*B_{q,o_t}

對於第一個狀態，初始值為

f_1(q) = \pi_q*B_{q,o_1}

這樣直到最後，我們從fT(q)f_T(q)中找一個概率最大的就是我們求解最優路徑的值了。
等等，我們是想找到最優路徑本身，你算一個最大概率值算怎麼回事啊？其實要得到最優路徑也很簡單，我們只需要在計算ft(q)f_t(q)順帶記錄一下最大的概率是從之前哪個狀態轉移過來的就好了，令：

g_t(q) = \mathop{\arg\max}_{q’}f_{t-1}(q’)*A_{q’,q}

gt(q)g_t(q)表示當t時刻隱藏狀態確認是q的情況下，它的上一時刻狀態應該是什麼，我們最後在fT(q)f_T(q)找到一個概率最大的q，只需要要根據這個q一步步回溯回去就可以找到對應的序列了,如圖e所示，我們在最後一步挑選了q1。

尾巴

隱馬爾可夫模型在語音識別、NLP等領域有著非常廣泛的應用，序列演算法裡未來還有其它跟「隱變數」打交道的東西，隱馬爾可夫應該算經典之一。
之前在不同的書上都看過多次隱馬爾可夫模型，但記憶不深刻，時間稍微一長就忘掉了，所以藉著這個機會從頭到尾梳理一次，按照概率表示)->(EM演算法)->(前向後向演算法)->(Viterbi)這個順序理解相當酣暢，都是有需求才進入下一步的介紹，承上啟下，相信這次要記得久一點了。
在過程中參考了《PRML》和《統計學習方法》，發現書上並不一定按照這個思路來講，不是很好串起來，而且一些概率公式之間跳的太快（估計是我數學比較渣）或者壓根沒有，這次從頭搞一次，我感覺應該是比較豐實了。