1. 程式人生 > >機器學習之最大期望(EM)演算法

機器學習之最大期望(EM)演算法

1.EM演算法簡介

最大期望(Expectation Maximum)演算法是一種迭代優化演算法,其計算方法是每次迭代分為期望(E)步最大(M)步。我們先看下最大期望演算法能夠解決什麼樣的問題。

假如班級裡有50個男生和50個女生,且男生站左,女生站右。我們假定男生和女生的身高分佈分別服從正態分佈。這時我們用極大似然法,分別通過這50個男生和50個女生的樣本來估計這兩個正態分佈的引數,便可知道男女身高分佈的情況。
N1(μ1,σ12);N2(μ2,σ22) N_1(\mu_1,\sigma_1^2)\ ; N_2(\mu_2,\sigma_2^2)

)
但我們面對下類問題如何解決呢?就是這50個男生和50個女生混在一起,我們擁有100個人的身高資料,但卻不知道這100個同學中的每個人是男生還是女生。通常來說,我們只有知道了精確的男女身高的正態分佈引數才能知道每位同學更有可能是男生還是女生,從另一方面來說我們只有知道每個人是男生還是女生才能儘可能準確的估計男女生各自身高的正態分佈引數。但現在兩者都不知道,如何去估計呢?

EM演算法表示我們可以用迭代的方法來解決問題。我們先假設男生身高和女生身高分佈的引數(初始值),然後根據這些引數去判斷每個人是男生還是女生,之後根據標註後的樣本反過來重新估計這些引數。多次重複上述過程,直到穩定。這樣說還是有點抽象,我們先從拋硬幣例項來講解EM演算法流程,然後再講解具體數學推導和原理。

2.EM演算法例項

假如現在我們有兩枚硬幣1和2,隨機拋擲後面朝上概率分別為P1,P2。為了估計兩硬幣概率,我們做如下實驗,每次取一枚硬幣,連擲5次後得到如下結果。

硬幣 結果 統計
1 正正反正反 3正-2反
2 反反正正反 2正-3反
1 正反反反反 1正-4反
2 正反反正正 3正-2反
1 反正正反反 2正-3反

我們可以很方便的估計出硬幣1概率P1=0.4,硬幣2概率P2=0.5。
P1=3+1+215=0.4 P1=\frac{3+1+2}{15}=0.4

P2=2+310=0.5 P2=\frac{2+3}{10}=0.5

=102+3=0.5

下面我們增加問題難度。如果並不知道每次投擲時所使用的硬幣標記,那麼如何估計P1和P2呢?

硬幣 結果 統計
Unknown 正正反正反 3正-2反
Unknown 反反正正反 2正-3反
Unknown 正反反反反 1正-4反
Unknown 正反反正正 3正-2反
Unknown 反正正反反 2正-3反

此時我們加入隱含變數z,可以把它認為是一個5維的向量(z1,z2,z3,z4,z5),代表每次投擲時所使用的硬幣。比如z1就代表第一輪投擲時所使用的是硬幣1還是2,我們必須先估計出z,然後才能進一步估計P1和P2。

我們先隨機初始化一個P1和P2,用它來估計z,然後基於z按照最大似然概率方法去估計新的P1和P2。如果估計出的P1和P2與我們初始化的P1和P2一樣,說明P1和P2很有可能就是真實的值。如果新估計出來的P1和P2和我們初始化時的值差別很大,那麼我們繼續用新的P1和P2迭代,直到收斂。

例如假設P1=0.2和P2=0.7,然後我們看看第一輪投擲的最可能是哪個硬幣。如果是硬幣1,得出3正2反的概率為0.2*0.2*0.2*0.8*0.8=0.00512,如果是硬幣2,得出3正2反的概率為0.7*0.7*0.7*0.3*0.3=0.03087。然後依次求出其他4輪中的相應概率,接下來便可根據最大似然方法得到每輪中最有可能的硬幣。

輪數 若是硬幣1 若是硬幣2 最有可能硬幣
1 0.00512 0.03087 硬幣2
2 0.02048 0.01323 硬幣1
3 0.08192 0.00567 硬幣1
4 0.00512 0.03087 硬幣2
5 0.02048 0.01323 硬幣1

我們把上面的值作為z的估計值(2,1,1,2,1),然後按照最大似然概率方法來估計新的P1和P2。得到
P1=2+1+215=0.33 P1=\frac{2+1+2}{15}=0.33

P2=3+310=0.6 P2=\frac{3+3}{10}=0.6

P1和P2的最大似然估計是0.4和0.5,那麼對比下我們初識化時的P1和P2。

初始化的P1 估計的P1 真實的P1 初始化的P2 估計的P2 真實的P2
0.2 0.33 0.4 0.7 0.6 0.5

可以預估,我們繼續按照上面思路,用估計出的P1和P2再來估計z,再用z來估計新的P1和P2,反覆迭代下去,可以最終得到P1=0.4,P2=0.5。然後無論怎樣迭代,P1和P2的值都會保持0.4和0.5不變,於是我們就找到P1和P2的最大似然估計。

上面我們用最大似然方法估計出z值,然後再用z值按照最大似然概率方法估計新的P1和P2。也就是說,我們使用了最有可能的z值,而不是所有的z值。如果考慮所有可能的z值,對每一個z值都估計出新的P1和P2,將每一個z值概率大小作為權重,將所有新的P1和P2分別加權相加,這樣估計出的P1和P2是否會更優呢?

但所有的z值共有2^5=32種,我們是否進行32次估計呢?當然不是,我們利用期望來簡化運算。

輪數 若是硬幣1 若是硬幣2
1 0.00512 0.03087
2 0.02048 0.01323
3 0.08192 0.00567
4 0.00512 0.03087
5 0.02048 0.01323

利用上面表格,我們可以算出每輪投擲種使用硬幣1或者使用硬幣2的概率。比如第一輪使用硬幣1的概率
z1=0.005120.00512+0.03087=0.14 z_1=\frac{0.00512}{0.00512+0.03087}=0.14
相應的算出其他4輪的概率。

輪數 z_i=硬幣1 z_i=硬幣2
1 0.14 0.86
2 0.61 0.39
3 0.94 0.06
4 0.14 0.86
5 0.61 0.39

上表中表示期望值,例如0.86表示此輪中使用硬幣2的概率是0.86。前面方法我們按照最大似然概率直接將第一輪估計為硬幣2,此時我們更加嚴謹,只說有0.14的概率是硬幣1,有0.86的概率是硬幣2。這樣我們在估計P1或者P2時,就可以用上全部的資料,而不是部分的資料。此步我們實際上是估計出z的概率分佈,稱為E步

按照期望最大似然概率的法則來估計出新的P1和P2。以P1估計為例,第一輪的3正2反相當於有0.14*3=0.42的概率為正,有0.14*2的概率為反。然後依次計算出其他四輪。那麼我們可以得到P1概率,可以看到改變了z的估計方法後,新估計出的P1要更加接近0.4,原因是我們使用了所有拋擲的資料,而不是部分的資料。此步中我們根據E步中求出z的概率分佈,依據最大似然概率法則去估計P1和P2,稱為M步

輪數 正面 反面
1 0.42 0.28
2 1.22 1.83
3 0.94 3.76
4 0.42 0.28
5 1.22 1.93
總計 4.22 7.98

P1=4.224.22+7.98=0.35 P1=\frac{4.22}{4.22+7.98}=0.35
上面我們是通過迭代來得到P1和P2。但是我們同時想要知道,新估計出的P1和P2一定會更接近真實的P1和P2嗎,能夠收斂嗎?迭代一定會收斂到真實的P1和P2嗎?下面我們將從數學推導方法詳解EM演算法。

3.EM演算法推導

對於m個樣本觀察資料x=(x(1),x(2),x(3),,x(m))x=(x^{(1)},x^{(2)},x^{(3)},…,x^{(m)}),找出樣本的模型引數θ\theta,極大化模型分佈的對數似然函式如下所示
θ=argmaxθi=1mlogP(x(i);θ) \theta =\arg \max_{\theta} \sum _{i=1}^{m}logP(x^{(i)};\theta)

如果我們得到的觀察資料有未觀察到的隱含資料z=(z(1),z(2),z(3),,z(m))z=(z^{(1)},z^{(2)},z^{(3)},…,z^{(m)}),此時我們極大化模型分佈的對數似然函式如下
θ=argmaxθi=1mlogP(x(i);θ) \theta =\arg \max_{\theta} \sum _{i=1}^{m}logP(x^{(i)};\theta)

=argmaxθi=1mlogz(i)P(x(i),z(i);θ) =\arg \max_{\theta} \sum _{i=1}^{m}log\sum _{z^{(i)}}P(x^{(i)},z^{(i)};\theta)

上面方程是無法直接求出θ的,因此需要一些特殊技巧,在此我們引入Jensen不等式

設f是定義域為實數的函式,如果對於所有的實數X,f(X)的二次導數大於等於0,那麼f是凸函式。相反,f(X)的二次導數小於0,那麼f是凹函式。

**Jensen不等式定義:**如果f是凸函式,X是隨機變數,那麼E[f(X)]≥f(E(X))。相反,如果f式凹函式,X是隨機變數,那麼E[f(X)]≤f(E[X])。當且僅當X是常量時,上式取等號,其中E[X]表示X的期望。

我們再回到上述推導過程,得到如下方程式。
i=1mlogz(i)P(x(i),z(i);θ) \sum _{i=1}^{m}log\sum _{z^{(i)}}P(x^{(i)},z^{(i)};\theta)

=i=1mlogz(i)Qi(z(i))P(x(i),z(i);θ)Qi(z(i))(1) =\sum _{i=1}^{m}log\sum _{z^{(i)}}Q_i(z^{(i)})\frac{P(x^{(i)},z^{(i)};\theta) }{Q_i(z^{(i)})}\ \ \ \ \ \ \ (1)

i=1mz(i)Qi(z(i))logP(x(i),z(i);θ)Qi(z(i))(2) \ge \sum _{i=1}^{m}\sum _{z^{(i)}}Q_i(z^{(i)})log\frac{P(x^{(i)},z^{(i)};\theta) }{Q_i(z^{(i)})} \ \ \ \ \ \ \ (2)

我們來解釋下怎樣得到的方程式(1)和方程式(2),上面(1)式中引入一個未知的新的分佈Qi(z(i))Q_i(z^{(i)}),第二式用到Jensen不等式。

首先log函式是凹函式,那麼E[f(X)]≤f(E[X]),也就是f(E(X))≥E(f(X))。其中z(i)Qi(z(i))P(x(i),z(i);θ)Qi(z(i))\sum _{z^{(i)}}Q_i(z^{(i)})\frac{P(x^{(i)},z^{(i)};\theta) }{Q_i(z^{(i)})}