一文說懂EM演算法及其在HMM和GMM中的應用

作者：由一隻懶羊發表于收藏時間：2020-07-15

一、EM演算法要解決的問題

EM演算法就是最大期望演算法，用於解決無法觀測隱性變數的機率模型求引數的問題。這句話是什麼意思呢？舉個例子，如果一個學校只有男生，假設男生身高符合正態分佈，此時需要根據統計得到的男生身高計算出正態分佈模型中的均值和方差，那麼我們可以直接計算。但是如果學校既有男生也有女生，而且因為統計時的疏漏，無法區分統計的身高是男生還是女生的身高，此時要計算模型引數，就需要EM演算法了。

EM演算法依舊透過最大似然估計法來計算機率模型的引數。如果我們知道在狀態z下導致事件y發生的機率，那麼我們在計算事件的最大機率時，可以把表示式寫為：

然後令其求導=0的方式求出

$\theta$

，也就是求事件發生的機率最大時的引數。但是如果有多個狀態都可能導致事件y，而我們沒法確定到底每一件狀態導致了事件y，那麼我們該如何計算這個引數使得事情發生的機率最大呢？我們可以將表示式寫成：

$l(\theta)=log P(y|\theta)=log\sum_{z}^{}{P(y|z,\theta)}=log\sum_{z}^{}{P(z|\theta)P(y|z,\theta)}$

之所以是求和符號，就是因為希望透過求Y的邊緣函式表示其機率。

$l(\theta)=log P(Y|\theta)=\prod_{i=1}^{N}log P(y_{i}|\theta)=\prod_{i=1}^{N}log\sum_{j=1}^{k}{P(y_{i}|z_{ij},\theta)}=\prod_{i=1}^{N}log\sum_{j=1}^{k}{P(z_{ij}|\theta)P(y_{i}|z_{ij},\theta)}$

二、EM演算法的E步是怎麼得出的？

我們在使用EM演算法計算時，常規步驟通常為E步設立公式，然後M步求導，得到引數值再代入E步的公式，如此反覆計算。可是我們E步中的公式是怎麼來的呢？下面就進行解釋

因為不知道

$P(z|\theta)$

，我們沒有辦法直接求得

$\theta$

。但是我們發現

$l(\theta)$

是單調遞增的，而且有上限（因為機率在0到1之間），那麼就可以得知

$l(\theta)$

會收斂到某一個值，那麼就意味雖然不能求出最優解，但是我們可以先隨意初始化一個

$\theta$

，然後透過不斷最佳化

$\theta$

使得

$l(\theta)$

不斷的變大，最終達到收斂狀態，得到其下界以及這個下界對應的引數

$\theta$

。既然如此，那麼就應該滿足

$l(\theta)-l(\theta_{i})>=0$

。我們來看它們的差值：

$l(\theta)-l(\theta_{i})= log\sum_{z}^{}{P(z|\theta)P(y|z,\theta)}-log P(y|\theta_{i})\\ = log\sum_{z}^{}{P(z|y,\theta_{i})\frac{P(y|z,\theta)P(z|\theta)}{P(z|y,\theta_{i})}}-log P(y|\theta_{i})\\$

由於log中帶有求和函式，使得式子難以處理，所以我們利用琴生不等式來對前半部分縮放。所謂琴生不等式就是：

經過不等式變換可得：

$\begin{align} &l(\theta)-l(\theta_{i})>=\sum_{z}^{}{P(z|y,\theta_{i})log\frac{P(y|z,\theta)P(z|\theta)}{P(z|y,\theta_{i})}}-log P(y|\theta_{i})\\ & l(\theta)-l(\theta_{i})>= \sum_{z}^{}{P(z|y,\theta_{i})log\frac{P(y|z,\theta)P(z|\theta)}{P(z|y,\theta_{i})}}-\sum_{z}^{}P(z|y,\theta_{i})log P(y|\theta_{i})\\ & l(\theta)-l(\theta_{i})>=\sum_{z}^{}{P(z|y,\theta_{i})log\frac{P(y|z,\theta)P(z|\theta)}{P(z|y,\theta_{i})P(y|\theta_{i})}} \end{align}$

注意後半部分之所以可以加上

$\sum_{z}^{}{P(z|y,\theta_{i})}$

，是因為它的和為1，而且

$logP(y|\theta_{i})$

是一個常數。我們將公式稍作變形：

$l(\theta)>=l(\theta_{i})+\sum_{z}^{}{P(z|y,\theta_{i})log\frac{P(y|z,\theta)P(z|\theta)}{P(z|y,\theta_{i})P(y|\theta_{i})}}=B(\theta,\theta_{i})$

雖然我們無法直接求

$l(\theta)$

，但是可以不斷的增大

$B(\theta,\theta_{i})$

來逼近

$l(\theta)$

。如果我們將

$\theta$

取值為上一次迭代的值

$\theta_{i}$

，那麼：

$B(\theta_{i},\theta_{i})=l(\theta_{i})，因為\frac{P(y|z,\theta)P(z|\theta)}{P(z|y,\theta_{i})P(y|\theta_{i})}=1，log1=0$

逼近的效果如下圖所示：

也就是說，我們要求一個引數，使得B最大，最接近原來的機率，即

$\begin{align}\\ &\theta_{i+1}=arg\max_{\theta}B(\theta,\theta_{i})\\ &=arg\max_{\theta}[l(\theta_{i})+\sum_{z}^{}{P(z|y,\theta_{i})log\frac{P(y|z,\theta)P(z|\theta)}{P(z|y,\theta_{i})P(y|\theta_{i})}}]\\ &\simeq arg\max_{\theta}[\sum_{z}^{}{P(z|y,\theta_{i})log{P(y|z,\theta)P(z|\theta)}}]\\ &=arg\max_{\theta}[\sum_{z}^{}{P(z|y,\theta_{i})log{P(y,z|\theta)}}]\\ &=arg\max_{\theta}E_{Z\sim P(z|y,\theta_{i})}log{P(y,z|\theta)}\\ &=arg\max_{\theta}Q(\theta,\theta_{i}) \end{align}\\$

在這個過程中，我們去掉常量，發現要想求得最大機率，最終只和在上一個引數

$\theta_{i}$

背景下的機率分佈（即已知事件y發生，判斷z是哪種隱狀態的機率）和現有的事件發生的機率有關。

這就是為什麼我們使用EM演算法，一開始計算期望時就用

$\sum_{z}^{}{P(z|y,\theta_{i})log{P(y,z|\theta)}}$

來計算的原因！而我們回頭看，整個公式最為關鍵的地方就是為前半部分式子添加了

$P(z|y,\theta_{i})$

，為下面的合併兩項奠定了基礎。

三、EM演算法在求解HMM和GMM中的求解

3.1 EM演算法在拋硬幣模型中的引數求解過程

問題背景是如果有三枚硬幣abc，他們出現正面的機率分別為π，p，q，我們先拋擲a硬幣，如果是正面，則用b拋擲，如果是反面，則用c拋擲，然後記錄下b或者c丟擲來的結果。現在得到觀測結果為1、1、0、1、0、0、1、0、1、1。問如何估計三個硬幣正面出現的機率，即三硬幣模型的引數。

注意每一次的拋擲結果是獨立的，所以事件的總機率是各個事件的連乘。由於我們要求的完全資料對數似然函式明確，因此我們直接進入到E步，得到事情最大機率的表示式：

$\begin{align} Q(\theta,\theta_{i})&=\sum_{z}^{} P(z_{j}|Y,\theta_{i}) log{P(Y,z_{j}|\theta)}\\ &=\sum_{z}^{} [\prod_{j=1}^{N} {P(z_{j}|y_{j},\theta_{i})][log\prod_{j=1}^{N} {P(y_{j},z_{j}|\theta)}} ]\\ &=\sum_{z}^{} [\prod_{j=1}^{N} P(z_{j}|y_{j},\theta_{i})] [\sum_{j=1}^{N}P(y_{j},z_{j}|\theta)]\\ &=\sum_{z}^{} A B \end{align}$

我們僅看B這一項，將其分解成兩部分，第一部分是第一項，第二部分是2到N項，計算可得：

$\begin{align} Q(\theta,\theta_{i})&=\sum_{z}^{} [\prod_{j=1}^{N} P(z_{j}|y_{j},\theta_{i}) logP(y_{1},z_{1}|\theta)]+ \sum_{z}^{} [\prod_{j=1}^{N} P(z_{j}|y_{j},\theta_{i}) \sum_{j=2}^{N}logP(y_{j},z_{j}|\theta)]\\ &=C+D \end{align}$

我們發現D中的項也可以分解為N-1項形式如C的式子，因此我們只考慮C即可。Z代表的是選擇哪一枚硬幣，我們考慮其中一枚硬幣，其中每一個硬幣又都有兩種狀態，就是正面和反面。

$\begin{align} C&=\sum_{z_{1……n}}^{} [\prod_{j=1}^{N} P(z_{j}|y_{j},\theta_{i}) logP(y_{1},z_{1}|\theta)] \quad (1)\\ &=\sum_{z_{1……n}}^{} [\prod_{j=2}^{N} P(z_{j}|y_{j},\theta_{i}) P(z_{1}|y_{1},\theta_{i}) logP(y_{1},z_{1}|\theta)]=\sum_{z_{1……n}}^{} [\prod_{j=2}^{N}EFG] \quad (2)\\ &=\sum_{z_{2……n}}^{} [\prod_{j=2}^{N} P(z_{j}|y_{j},\theta_{i}) P(z_{1}=0|y_{1},\theta_{i}) logP(y_{1},z_{1}=0|\theta)]+\sum_{z_{1……n}}^{} [\prod_{j=2}^{N} P(z_{j}|y_{j},\theta_{i}) P(z_{1}=1|y_{1},\theta_{i}) logP(y_{1},z_{1}=1|\theta)] \quad (3)\\ &=[P(z_{1}=0|y_{1},\theta_{i}) logP(y_{1},z_{1}=0|\theta)+P(z_{1}=1|y_{1},\theta_{i}) logP(y_{1},z_{1}=1|\theta)]\sum_{z_{2……n}}^{} [\prod_{j=2}^{N}P(z_{j}|y_{j},\theta_{i}) \quad (4)\\ &=\sum_{z_{1}}^{}{P(z_{1}|y_{1},\theta_{i}) logP(y_{1},z_{1}|\theta)} \sum_{z_{2……n}}^{} [\prod_{j=2}^{N}P(z_{j}|y_{j},\theta_{i}) \quad (5)\\ &=\sum_{z_{1}}^{}{P(z_{1}|y_{1},\theta_{i}) logP(y_{1},z_{1}|\theta)} \quad (6) \end{align}$

其實第二步最為關鍵，因為將累乘符號分解以後，（2）中的F和G有關，而E則和G無關，又因為E是指硬幣z在不同取值下的機率，其實等於1，這樣n-1個1累乘法依然是1，所以最終得到（6）這一表達式。

我們知道

$Q(\theta,\theta_{i})$

可以分解成C+D，而每一個C都會得到類似（6）中的表示式，那麼

$Q(\theta,\theta_{i})$

的表示式可以表示成：

$\begin{align} Q(\theta,\theta_{i}) &=\sum_{z_{1}}^{}{P(z_{1}|y_{1},\theta_{i}) logP(y_{1},z_{1}|\theta)} +\sum_{z_{2}}^{}{P(z_{2}|y_{2},\theta_{i}) logP(y_{2},z_{2}|\theta)} +……+\sum_{z_{n}}^{}{P(z_{n}|y_{n},\theta_{i}) logP(y_{n},z_{n}|\theta)} \\ &=\sum_{j=1}^{N} [{\sum_{z_{j}}^{}{P(z_{j}|y_{j},\theta_{i}) logP(y_{j},z_{j}|\theta)}] } \\ &=\sum_{j=1}^{N} [P(z_{j}=1|y_{j},\theta_{i})lnP(y_{j},z_{j}=1|\theta) + P(z_{j}=0|y_{j},\theta_{i})lnP(y_{j},z_{j}=0|\theta)] \end{align}$