統計學習筆記九----EM算法


前言

  EM算法是一種迭代算法,1977年由Dempster等人總結提出,用於含有隱變量(hidden variable)的概率模型參數的極大似然估計,或極大后驗概率估計。EM算法的每次迭代由兩步組成:E步,求期望(expection);M步,求極大值(maximization),所以這一算法稱為期望極大算法(exception maximization algorithm),簡稱EM算法。

極大似然估計

  極大似然估計,只是一種概率論在統計學中的應用,它是參數估計的方法之一。說的是已知某個隨機樣本滿足某種概率分布,但是其中具體的參數不清楚,參數估計就是通過若干次實驗,觀察其結果,利用結果推出參數的大概值。最大似然估計是建立在這樣的思想上:已知某個參數能使這個樣本出現的概率值最大,我們當然不會再去選擇其他小概率的樣本,所以干脆就把這個參數作為估計的真實值。

最大似然估計你可以把它看作是一個反推。多數情況下我們是根據已知條件來推算結果,而最大似然估計是已經知道了結果,然后尋求使該結果出現的可能性最大的條件,以此作為估計值。比如,如果其他條件一定的話,抽煙者發生肺癌的危險時不抽煙者的5倍,那么如果現在我已經知道有個人是肺癌,我想問你這個人抽煙還是不抽煙。你怎么判斷?你可能對這個人一無所知,你所知道的只有一件事,那就是抽煙更容易發生肺癌,那么你會猜測這個人不抽煙嗎?我相信你更有可能會說,這個人抽煙。為什么?這就是“最大可能”,我只能說他“最有可能”是抽煙的,“他是抽煙的”這一估計值才是“最有可能”得到“肺癌”這樣的結果。這就是最大似然估計。

求最大似然函數估計值的一般步驟:

(1)、寫出似然函數;
(2)、對數似然函數取對數,並整理;
(3)、求導數,令導數為0,得到似然方程;
(4)、解似然方程,得到的參數即為所求。

最大(極大)似然估計也是統計學習中經驗風險最小化(RRM)的例子。如果模型為條件概率分布,損失函數定義為對數損失函數,經驗風險最小化就等價於最大似然估計。

綜之,最大似然估計就是給定模型(含有未知參數)和樣本集的情況下,用來估計模型參數的方法。其基本思想是找到最佳的模型參數,使得模型實現對樣本的最大程度的擬合,也就是使樣本集出現的可能性最大。這一方法是基於這樣的思想:我們所估計的模型參數,要使得產生這個給定樣本的可能性最大。在最大釋然估計中,我們試圖在給定模型的情況下,找到最佳的參數,使得這組樣本出現的可能性最大。舉個極端的反面例子,如果我們得到一個中國人口的樣本,男女比例為3:2,現在讓你估計全國人口的真實比例,你肯定不會估計為男:女=1:0。因為如果是1:0,不可能得到3:2的樣本。我們大多很容易也估計為3:2,為什么?樣本估計總體,其實它背后的思想就是最大似然。

如果你還是不是很理解最大似然,可以參考以下實例鏈接:
參考一
參考二
參考三

EM算法的引入

我們先來看看一下書上的引入實例:

這里寫圖片描述

我們把中間的推導公式寫一下:




由以上三個等式可以得到:

注意:這里的隨機變量y是觀測變量,表示一次實驗觀測的結果是1或者0;隨機變量z是隱變量,表示的實未觀測到的擲硬幣A的結果;是模型參數。這一模型是以上數據的生成模型。注意,隨機變量y的數據可以觀測,隨機變量z的數據不可觀測到。

如果我們將這n次獨立實驗,將觀測數據表示為向量,將未觀測到的數據表示為向量,則觀測數據的似然函數為:


由於觀測到的n個向量是獨立同分布的,所以概率是連乘的形式:
這里寫圖片描述

考慮求模型參數的極大似然估計,也就是說求得一個最優參數使得以上條件概率最大(或者對數條件概率最大)。即:

這個問題沒有解析解,只有通過迭代的方法求解。EM算法就是可以用於求解這個問題的一種迭代算法。求解過程如下所示,其實,理解一下求解過程,你會發現EM算法的思想與坐標上升算法或者SMO算法思想是一致的,其實他們思想的本質是相同的。

這里寫圖片描述

這里寫圖片描述

一般的,用Y表示觀測隨機變量的數據,Z表示隱藏隨機變量的數據。Y和Z連在一起稱為完全數據(complete-data),觀測數據Y又稱為不完全數據(imcomplete-data)。假設給定觀測數據Y,其概率分布是,其中是需要估計的模型參數,那么不完全數據Y的似然函數是,對數似然函數是;假設Y和Z的聯合概率分布是,那么完全數據的對數似然函數是.

EM算法通過迭代求的極大似然估計,每次迭代包含兩步:E步,求期望,M步,求極大化。下面來介紹EM算法。

我們還是直接上圖吧,(那些公式我是真的不想編輯了….)

這里寫圖片描述

這里寫圖片描述

EM算法的推導

假設我們有一個樣本集{x(1),…,x(m)},包含m個獨立的樣本。但每個樣本i對應的類別z(i)是未知的(相當於聚類),也即隱含變量。故我們需要估計概率模型p(x,z)的參數θ,但是由於里面包含隱含變量z,所以很難用最大似然求解,但如果z知道了,那我們就很容易求解了。

對於參數估計,我們本質上還是想獲得一個使得似然函數最大化的那個參數θ,現在與最大化似然不同的只是似然函數式中多了一個未知的變量z。也就是說我們的目標是:找到合適的θ和z讓L(θ)最大。那我們也許會想,你就是多了一個未知的變量而已啊,我們也可以分別對未知的變量θ和z分別求偏導,然后再令其為0,求解出來不也一樣嗎?

本質上我們是需要最大化(1)式(對(1)式,我們回憶下聯合概率密度下某個變量的邊緣概率密度函數的求解,注意這里z也是隨機變量。對每一個樣本i的所有可能類別z求等式右邊的聯合概率密度函數和,也就得到等式左邊為隨機變量x的邊緣概率密度),也就是似然函數,但是可以看到里面有“和的對數”,求導后形式會非常復雜(自己可以想象下log(f1(x)+ f2(x)+ f3(x)+…)復合函數的求導),所以很難求解得到未知參數z和θ。那OK,我們可否對(1)式做一些改變呢?我們看(2)式,(2)式只是分子分母同乘以一個相等的函數,還是有“和的對數”啊,還是求解不了,那為什么要這么做呢?咱們先不管,看(3)式,發現(3)式變成了“對數的和”,那這樣求導就容易了。我們注意點,還發現等號變成了不等號,為什么能這么變呢?這就是Jensen不等式的大顯神威的地方。

 jensen不等式

  設f是定義域為實數的函數,如果對於所有的實數x,f(x)的二次導數大於等於0,那么f是凸函數。當x是向量時,如果其hessian矩陣H是半正定的,那么f是凸函數。如果只大於0,不等於0,那么稱f是嚴格凸函數。

Jensen不等式表述如下:
如果f是凸函數,X是隨機變量,那么:E[f(X)]>=f(E[X])
特別地,如果f是嚴格凸函數,當且僅當X是常量時,上式取等號。

如果用圖表示會很清晰:

圖中,實線f是凸函數,X是隨機變量,有0.5的概率是a,有0.5的概率是b。(就像擲硬幣一樣)。X的期望值就是a和b的中值了,圖中可以看到E[f(X)]>=f(E[X])成立。

當f是(嚴格)凹函數當且僅當-f是(嚴格)凸函數。

Jensen不等式應用於凹函數時,不等號方向反向。
回到公式(2),因為f(x)=log x為凹函數(其二次導數為-1/x^2<0)。

(2)式中的期望,(考慮到是X的函數,則,又,所以就可以得到公式(3)的不等式了:

OK,到這里,現在式(3)就容易地求導了,但是式(2)和式(3)是不等號啊,式(2)的最大值不是式(3)的最大值啊,而我們想得到式(2)的最大值,那怎么辦呢?

現在我們就需要一點想象力了,上面的式(2)和式(3)不等式可以寫成:似然函數L(θ)>=J(z,Q),那么我們可以通過不斷的最大化這個下界J,來使得L(θ)不斷提高,最終達到它的最大值。

見上圖,我們固定θ,調整Q(z)使下界J(z,Q)上升至與L(θ)在此點θ處相等(綠色曲線到藍色曲線),然后固定Q(z),調整θ使下界J(z,Q)達到最大值(),然后再固定θ,調整Q(z)……直到收斂到似然函數L(θ)的最大值處的θ*。這里有兩個問題:什么時候下界J(z,Q)與L(θ)在此點θ處相等?為什么一定會收斂?

首先第一個問題,在Jensen不等式中說到,當自變量X是常數的時候,等式成立。而在這里,即:

再推導下,由於(因為Q是隨機變量z(i)的概率密度函數),則可以得到:分子的和等於c(分子分母都對所有z(i)求和:多個等式分子分母相加比值不變,這個認為每個樣例的兩個概率比值都是c),則:

至此,我們推出了在固定參數θ后,使下界拉升的Q(z)的計算公式就是后驗概率,解決了Q(z)如何選擇的問題。這一步就是E步,建立L(θ)的下界。接下來的M步,就是在給定Q(z)后,調整θ,去極大化L(θ)的下界J(在固定Q(z)后,下界還可以調整的更大)。那么一般的EM算法的步驟如下:

EM的算法流程:

初始化分布參數θ;

重復以下步驟直到收斂:

E步驟:根據參數初始值或上一次迭代的模型參數來計算出隱性變量的后驗概率,其實就是隱性變量的期望。作為隱藏變量的現估計值:

M步驟:將似然函數最大化以獲得新的參數值:

這個不斷的迭代,就可以得到使似然函數L(θ)最大化的參數θ了。那就得回答剛才的第二個問題了,它會收斂嗎?

感性的說,因為下界不斷提高,所以極大似然估計單調增加,那么最終我們會到達最大似然估計的最大值。理性分析的話,就會得到下面的東西:

具體如何證明的,看推導過程參考

接下來我貼幾張圖片(最后的參考文獻部分也給出了鏈接),是我關注的一個人寫的,有助於大家的理解:


利用坐標上升法理解EM算法

圖中的直線式迭代優化的路徑,可以看到每一步都會向最優值前進一步,而且前進路線是平行於坐標軸的,因為每一步只優化一個變量。

這猶如在x-y坐標系中找一個曲線的極值,然而曲線函數不能直接求導,因此什么梯度下降方法就不適用了。但固定一個變量后,另外一個可以通過求導得到,因此可以使用坐標上升法,一次固定一個變量,對另外的求極值,最后逐步逼近極值。對應到EM上,E步:固定θ,優化Q;M步:固定Q,優化θ;交替將極值推向最大。

助於理解EM算法的例子:

實例一

這是一個拋硬幣的例子,H表示正面向上,T表示反面向上,參數θ表示正面朝上的概率。硬幣有兩個,A和B,硬幣是有偏的。本次實驗總共做了5組,每組隨機選一個硬幣,連續拋10次。如果知道每次拋的是哪個硬幣,那么計算參數θ就非常簡單了,如下圖所示。

這里寫圖片描述

如果不知道每次拋的是哪個硬幣呢?那么,我們就需要用EM算法,基本步驟為:1、給θA和θB一個初始值;2、(E-step)估計每組實驗是硬幣A的概率(本組實驗是硬幣B的概率=1-本組實驗是硬幣A的概率)。分別計算每組實驗中,選擇A硬幣且正面朝上次數的期望值,選擇B硬幣且正面朝上次數的期望值;3、(M-step)利用第三步求得的期望值重新計算θA和θB;4、當迭代到一定次數,或者算法收斂到一定精度,結束算法,否則,回到第2步。

稍微解釋一下上圖的計算過程。初始值θA=0.6,θB=0.5。

圖中的0.45是怎么得來的呢?由兩個硬幣的初始值0.6和0.5,容易得出投擲出5正5反的概率是, 0.45就是0.449近似而來的,表示第一組實驗選擇的硬幣是A的概率為0.45。圖中的2.2H,2.2T是怎么得來的呢? ,表示第一組實驗選擇A硬幣且正面朝上次數的期望值是2.2。其他的值依次類推。

Python 實現

現在我們用Python實現上述例子中的第一次迭代過程,也就是對應上圖的(1)、(2)、(3)步,如下所示:

# -*- encoding:utf-8 -*-

import numpy
from scipy.stats import binom


#模擬實現實例中的一次迭代過程:擲硬幣
def em_single(priors,observations):
'''
EM算法的一次迭代
:param priors: 初始化參數theta_A,theta_B [theta_A,theta_B]
:param observations: 觀測矩陣:M*N:M代表實驗組數,N:代表每組實驗的次數
:return:[new_theta_A,new_theta_B]
'''


counts={'A':{'H':0,'T':0},'B':{'H':0,'T':0}}
theta_A=priors[0]
theta_B=priors[1]


# E step:

for observation in observations:

len_observation = len(observation)
num_heads = observation.sum()#正面
num_tails = len_observation-num_heads#反面

#二項分布求解公式,拋硬幣是一個二項分布
contribution_A = binom.pmf(num_heads,len_observation,theta_A)
contribution_B = binom.pmf(num_heads,len_observation,theta_B)

#將兩個概率正規化,得到數據來自硬幣A,B的概率。
weight_A = contribution_A / (contribution_A + contribution_B)
weight_B = contribution_B / (contribution_A + contribution_B)

#更新在當前參數下A,B硬幣產生的正反面次數
counts['A']['H'] += weight_A * num_heads
counts['A']['T'] += weight_A * num_tails
counts['B']['H'] += weight_B * num_heads
counts['B']['T'] += weight_B * num_tails

# M step:
new_theta_A = counts['A']['H'] / (counts['A']['H'] + counts['A']['T'])
new_theta_B = counts['B']['H'] / (counts['B']['H'] + counts['B']['T'])
return [new_theta_A,new_theta_B]

if __name__ == '__main__':

#采集數據集,硬幣投擲結果,1:表示H 正面;0:表示T 反面
observations = numpy.array([[1,0,0,0,1,1,0,1,0,1],
[1,1,1,1,0,1,1,1,1,1],
[1,0,1,1,1,1,1,0,1,1],
[1,0,1,0,0,0,1,1,0,0],
[0,1,1,1,0,1,1,1,0,1]])
priors=[0.6,0.5]

new_thetas=em_single(priors,observations)
print new_thetas

實驗結果:[0.71301223540051617, 0.58133930831366265]和途中的結果是一致的。

下面我們再給出EM算法的主循環:也就是在收斂前迭代的調用單次EM方法。

'''
#EM算法的主循環
給定循環的兩個終止條件:模型參數變化小於閾值;
循環達到最大次數
'''

def em(observations,prior,threshold = 1e-6,iterations=10000):
"""
EM算法
:param observations :觀測數據
:param prior:模型初值
:param tol:迭代結束閾值
:param iterations:最大迭代次數
:return:局部最優的模型參數
"""

iteration = 0
while iteration < iterations:
new_prior = em_single(prior,observations)
delta_change = numpy.abs(prior[0]-new_prior[0])
if delta_change < threshold:
break
else:
prior = new_prior
iteration +=1
#返回最終的權重,以及迭代次數
return (new_prior,iteration)

實驗結果:

最終得到硬幣A正面朝上的概率參數為: 0.796788759383
最終得到硬幣B正面朝上的概率參數為: 0.519583935675
最終的迭代次數為: 14

參考鏈接:
http://blog.csdn.net/u011300443/article/details/46763743

http://blog.csdn.net/abcjennifer/article/details/8170378

http://blog.csdn.net/zouxy09/article/details/8537620/

https://www.zhihu.com/question/27976634/answer/39132183


注意!

本站转载的文章为个人学习借鉴使用,本站对版权不负任何法律责任。如果侵犯了您的隐私权益,请联系我们删除。



 
粤ICP备14056181号  © 2014-2021 ITdaan.com