Statistical algorithm - 讓人著迷的EM演算法

07-09

來自專欄生物信息學-info

如何理解EM演算法？

通常來講，一堆由已知分布得到的數據，如果模型中的變數都是可以觀測到的，為了求其中的參數 $heta$ ，可以直接使用極大似然估計。

然而有這麼一種情況，當模型含有隱含變數的時候，極大似然估計的計算過程就會顯得極為複雜，甚至不可解。這時候就需要用到EM演算法了，那麼，EM演算法的思想究竟是什麼呢？

這裡，推薦一篇寫的很好的博客。

如何感性地理解EM演算法？?

www.jianshu.com

在直觀的理解EM思想後，我下面進行公式推導。從數學上來理解EM演算法。

EM演算法的數學公式推導

假設，有一個含有隱含變數 $Z$ 的模型，其概率密度函數為 $P(y,z| heta)$ 。現在，我們希望得到:

$hat{ heta}=mathop{argmax}_{ heta}logL(Y| heta)$ 【備註：極大似然估計】【1】

由 $P(y_j| heta) = mathop{sum}_{z^{(i)}}P(y_j,z^{(i)}| heta)$ 【備註：邊緣分布】【2】，

則 $logL(Y| heta)$ 有，

$egin{align} logL(Y| heta) & = mathop{sum}_{j}logP(y_j| heta) \ & = mathop{sum}_{j}logmathop{sum}_{z^{(i)}}P(y_j,z^{(i)}| heta) end{align}$ 【3】

到這一步，卡住了，卡死了，公式無從展開了。

我們需要思考。如果極大似然函數如上，求導時，勢必會面臨類似 $log(x_{1}+...+x_{n})$ 這樣的函數，它的導函數很複雜，難以求解。

說到 $x_{1}+...+x_{n}$ 這樣的式子，第一反應便是均值不等式，驗證一下想法，

由 $frac{x_1+x_2+...+x_n}{n} geq sqrt[n]{x_1x_2...x_n}$ ，有，

$log(frac{x_1+x_2+...+x_n}{n} )geq frac{1}{n}log(x_1x_2...x_n)$ ,變形，有,

$log(frac{x_1}{n}+frac{x_2}{n}+...+frac{x_n}{n})geq frac{1}{n}log(x_1)+frac{1}{n}log(x_2)+...+frac{1}{n}log(x_n)$ ,

咦，這是什麼？這個是不是說明， $logE(x) geq E(logx)$ 【備註：Jense不等式】【4】

當然，上述推導，Jense不等式是可以解釋的。但看到log里的連加，想變成log中的連乘，還是均值不等式比較容易聯想到。

也就是說，需要把【3】式中log里的連加變成連乘，需要構造一個期望函數。

對於【3】式，構造期望函數時，很容易想到 $logmathop{sum}_{z^{(i)}}a imes frac{P(y_j,z^{(i)}| heta)}{a}$ 這樣的式子，而 $logmathop{sum}_{z^{(i)}}P(y_j,z^{(i)}| heta)$ 是一個對於每一個z下的關於隨機變數y_j邊緣分布概率的求和公式，因此「a」是一個和z^(i)相關的分布函數，因此有概率密度函數 $Q(z^{(i)})$ ，且， $Q(z^{(i)})ge0,mathop{sum}_{z^{(i)}}Q(z^{(i)})=1$

現在，對於固定下來的每一個z，都能拿到z的概率密度值 $Q(z^{(i)})$ ，和一個形如 $frac{P(y_j,z^{(i)}| heta)}{Q(z^{(i)})}$ 這樣的式子。

根據 $E(g(x))=int g(x)f(x)dx$ 【備註： $Q(z^{(i)})$ 是 $f(x)$ 】和【3】式，有

$egin{align} logL(Y| heta) & = logE(frac{P(y_j,z^{(i)}| heta)}{Q(z^{(i)})}) \&ge mathop{sum}_{j}E(logfrac{P(y_j,z^{(i)}| heta)}{Q(z^{(i)})})\&=mathop{sum}_{j}mathop{sum}_{z^{(i)}}logfrac{P(y_j,z^{(i)}| heta)}{Q(z^{(i)})}*Q(z^{(i)}) end{align}$ 【5】

由於log函數在最後求導中會省略【5】式中分母下的 $Q(z^{(i)})$ ，因此【5】式可簡化為，

$logL(Y| heta) = mathop{sum}_{j}mathop{sum}_{z^{(i)}}Q(z^{(i)})logP(y_j,z^{(i)}| heta)$ 【6】

其中， $Q(z^{(i)})=P(z^{(i)}|y^{(i)}, heta^{(i)})$ 【7】

【6】式為EM演算法中的M步，【7】式為EM演算法中的E步

每一步在做什麼呢，這裡有一個簡單的手繪圖稍作解釋，

通過每一次E步和M步的迭代，可以尋找到一個局部最優點。不同的初始值，可以找到不同的局部最優點，從而找到全局最優點。下面的文章中將會實戰伯努利分布和二項分布來演示EM演算法的使用。