從似然函式到EM演算法(附程式碼實現)
1. 什麼是EM演算法
最大期望演算法(Expectation-maximization algorithm,又譯為期望最大化演算法),是在概率模型中尋找引數最大似然估計或者最大後驗估計的演算法,其中概率模型依賴於無法觀測的隱性變數。
最大期望演算法經過兩個步驟交替進行計算,
第一步是計算期望(E),利用對隱藏變數的現有估計值,計算其最大似然估計值;
第二步是最大化(M),最大化在E步上求得的最大似然值來計算引數的值。M步上找到的引數估計值被用於下一個E步計算中,這個過程不斷交替進行。
極大似然估計用一句話概括就是:知道結果,反推條件θ。
1.1 似然函式
在數理統計學中,似然函式是一種關於統計模型中的引數的函式,表示模型引數中的似然性。“似然性”與“或然性”或“概率”意思相近,都是指某種事件發生的可能性。而極大似然就相當於最大可能的意思。
比如你一位同學和一位獵人一起外出打獵,一隻野兔從前方竄過。只聽一聲槍響,野兔應聲到下,如果要你推測,這一發命中的子彈是誰打的?你就會想,只發一槍便打中,由於獵人命中的概率一般大於你那位同學命中的概率,從而推斷出這一槍應該是獵人射中的。
這個例子所作的推斷就體現了最大似然法的基本思想。
多數情況下我們是根據已知條件來推算結果,而最大似然估計是已經知道了結果,然後尋求使該結果出現的可能性最大的條件,以此作為估計值。
1.3 極大似然函式的求解步驟
假定我們要從10萬個人當中抽取100個人來做身高統計,那麼抽到這100個人的概率就是(概率連乘):
\[L(\theta)=L(x_1,...,x_n|\theta)=\prod_{i=1}^{n}p(x_i|\theta),\theta\in\ominus\]
現在要求的就是這個 \(\theta\) 值,也就是使得 \(L(\theta)\) 的概率最大化,那麼這時的引數 \(\theta\) 就是所求。
為了便於分析,我們可以定義對數似然函式,將其變成連加的形式:
\[H(\theta)=lnL(\theta)=ln\prod_{i=1}^{n}p(x_i|\theta)=\sum_{i=1}^{n}lnp(x_i|\theta)\]
對於求一個函式的極值,通過我們在本科所學的微積分知識,最直接的設想是求導,然後讓導數為0,那麼解這個方程得到的θ就是了(當然,前提是函式L(θ)連續可微)。但,如果θ是包含多個引數的向量那怎麼處理呢?當然是求L(θ)對所有引數的偏導數,也就是梯度了,從而n個未知的引數,就有n個方程,方程組的解就是似然函式的極值點了,最終得到這n個引數的值。
求極大似然函式估計值的一般步驟:
- 寫出似然函式;
- 對似然函式取對數,並整理;
- 求導數,令導數為0,得到似然方程;
- 解似然方程,得到的引數即為所求;
1.4 EM演算法
兩枚硬幣A和B,假定隨機拋擲後正面朝上概率分別為PA,PB。為了估計這兩個硬幣朝上的概率,咱們輪流拋硬幣A和B,每一輪都連續拋5次,總共5輪:
硬幣 | 結果 | 統計 |
---|---|---|
A | 正正反正反 | 3正-2反 |
B | 反反正正反 | 2正-3反 |
A | 正反反反反 | 1正-4反 |
B | 正反反正正 | 3正-2反 |
A | 反正正反反 | 2正-3反 |
硬幣A被拋了15次,在第一輪、第三輪、第五輪分別出現了3次正、1次正、2次正,所以很容易估計出PA,類似的,PB也很容易計算出來(真實值),如下:
PA = (3+1+2)/ 15 = 0.4
PB= (2+3)/10 = 0.5
問題來了,如果我們不知道拋的硬幣是A還是B呢(即硬幣種類是隱變數),然後再輪流拋五輪,得到如下結果:
硬幣 | 結果 | 統計 |
---|---|---|
Unknown | 正正反正反 | 3正-2反 |
Unknown | 反反正正反 | 2正-3反 |
Unknown | 正反反反反 | 1正-4反 |
Unknown | 正反反正正 | 3正-2反 |
Unknown | 反正正反反 | 2正-3反 |
OK,問題變得有意思了。現在我們的目標沒變,還是估計PA和PB,需要怎麼做呢?
顯然,此時我們多了一個硬幣種類的隱變數,設為z,可以把它認為是一個5維的向量(z1,z2,z3,z4,z5),代表每次投擲時所使用的硬幣,比如z1,就代表第一輪投擲時使用的硬幣是A還是B。
- 但是,這個變數z不知道,就無法去估計PA和PB,所以,我們必須先估計出z,然後才能進一步估計PA和PB。
- 可要估計z,我們又得知道PA和PB,這樣我們才能用極大似然概率法則去估計z,這不是雞生蛋和蛋生雞的問題嗎,如何破?
答案就是先隨機初始化一個PA和PB,用它來估計z,然後基於z,還是按照最大似然概率法則去估計新的PA和PB,然後依次迴圈,如果新估計出來的PA和PB和我們真實值差別很大,直到PA和PB收斂到真實值為止。
我們不妨這樣,先隨便給PA和PB賦一個值,比如:
硬幣A正面朝上的概率PA = 0.2
硬幣B正面朝上的概率PB = 0.7
然後,我們看看第一輪拋擲最可能是哪個硬幣。
如果是硬幣A,得出3正2反的概率為 0.20.20.20.80.8 = 0.00512
如果是硬幣B,得出3正2反的概率為0.70.70.70.30.3=0.03087
然後依次求出其他4輪中的相應概率。做成表格如下:
輪數 | 若是硬幣A | 若是硬幣B |
---|---|---|
1 | 0.00512,即0.2 0.2 0.2 0.8 0.8,3正-2反 | 0.03087,3正-2反 |
2 | 0.02048,即0.2 0.2 0.8 0.8 0.8,2正-3反 | 0.01323,2正-3反 |
3 | 0.08192,即0.2 0.8 0.8 0.8 0.8,1正-4反 | 0.00567,1正-4反 |
4 | 0.00512,即0.2 0.2 0.2 0.8 0.8,3正-2反 | 0.03087,3正-2反 |
5 | 0.02048,即0.2 0.2 0.8 0.8 0.8,2正-3反 | 0.01323,2正-3反 |
按照最大似然法則:
第1輪中最有可能的是硬幣B
第2輪中最有可能的是硬幣A
第3輪中最有可能的是硬幣A
第4輪中最有可能的是硬幣B
第5輪中最有可能的是硬幣A
我們就把概率更大,即更可能是A的,即第2輪、第3輪、第5輪出現正的次數2、1、2相加,除以A被拋的總次數15(A拋了三輪,每輪5次),作為z的估計值,B的計算方法類似。然後我們便可以按照最大似然概率法則來估計新的PA和PB。
PA = (2+1+2)/15 = 0.33
PB =(3+3)/10 = 0.6
就這樣,不斷迭代 不斷接近真實值,這就是EM演算法的奇妙之處。
可以期待,我們繼續按照上面的思路,用估計出的PA和PB再來估計z,再用z來估計新的PA和PB,反覆迭代下去,就可以最終得到PA = 0.4,PB=0.5,此時無論怎樣迭代,PA和PB的值都會保持0.4和0.5不變,於是乎,我們就找到了PA和PB的最大似然估計。
總結一下計算步驟:
隨機初始化分佈引數θ
E步,求Q函式,對於每一個i,計算根據上一次迭代的模型引數來計算出隱性變數的後驗概率(其實就是隱性變數的期望),來作為隱藏變數的現估計值:
\[Q_i(z^{(i)})=p(z^{(i)}|x^{(i)};\theta)\]
M步,求使Q函式獲得極大時的引數取值)將似然函式最大化以獲得新的引數值
\[\theta=argmax\sum_{i}\sum_{z^{(i)}}Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\]
然後迴圈重複2、3步直到收斂。
詳細的推導過程請參考文末的參考文獻。
2. 採用 EM 演算法求解的模型有哪些?
用EM演算法求解的模型一般有GMM或者協同過濾,k-means其實也屬於EM。EM演算法一定會收斂,但是可能收斂到區域性最優。由於求和的項數將隨著隱變數的數目指數上升,會給梯度計算帶來麻煩。
3.程式碼實現
高斯混合模型 EM 演算法
4. 參考文獻
如何通俗理解EM演算法
作者:@mantchs
GitHub:https://github.com/NLP-LOVE/ML-NLP
歡迎大家加入討論!共同完善此專案!群號:【541954936】