【機(jī)器學(xué)習(xí)基礎(chǔ)】數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法27：EM算法

? ? 從本篇開始，整個(gè)機(jī)器學(xué)習(xí)系列還剩下最后三篇涉及導(dǎo)概率模型的文章，分別是EM算法、CRF條件隨機(jī)場(chǎng)和HMM隱馬爾科夫模型。本文主要講解一下EM（Expection maximization），即期望最大化算法。EM算法是一種用于包含隱變量概率模型參數(shù)的極大似然估計(jì)方法，所以本文從極大似然方法說起，然后推廣到EM算法。

極大似然估計(jì)

? ? 統(tǒng)計(jì)學(xué)專業(yè)的朋友對(duì)于極大似然估計(jì)一定是很熟悉了。極大似然是一種統(tǒng)計(jì)參數(shù)估計(jì)方法，對(duì)于某個(gè)隨機(jī)樣本滿足某種概率分布，但其中的統(tǒng)計(jì)參數(shù)未知，我們通過若干次試驗(yàn)結(jié)果來估計(jì)參數(shù)的值的方法。

? ? 舉個(gè)例子來說。比如說我們想了解一下某校學(xué)生的身高分布。我們先假設(shè)該校學(xué)生身高服從一個(gè)正態(tài)分布，其中的分布參數(shù)和未知。全校大幾萬的學(xué)生，我們要一個(gè)個(gè)去實(shí)測(cè)肯定不現(xiàn)實(shí)。所以我們決定用統(tǒng)計(jì)抽樣的方法，隨機(jī)選取100名學(xué)生來看一下身高。

? ? 要通過這100人的身高來估算全校學(xué)生的身高，我們需要明確下面幾個(gè)問題。第一個(gè)就是抽到這100人的概率是多少，因?yàn)槊總€(gè)人的選取都是獨(dú)立的，所以選到這100人的概率可以表示為單個(gè)概率的乘積：

? ? 上式即為似然函數(shù)。通常為了計(jì)算方便，我們會(huì)對(duì)似然函數(shù)取對(duì)數(shù)：

? ? 第二個(gè)問題在于我們要解釋一下為什么能夠剛好抽到這100人。所以按照極大似然估計(jì)的理論，在學(xué)校這么多人中，我們恰好抽到這100人而不是另外的100人，正是因?yàn)檫@100人出現(xiàn)的概率極大，即其對(duì)應(yīng)的似然函數(shù)極大：

? ? 第三個(gè)問題在于如何求解。這個(gè)好辦，直接對(duì)求其參數(shù)的偏導(dǎo)數(shù)并令為0即可。

? ? 所以極大似然估計(jì)法可以看作由抽樣結(jié)果對(duì)條件的反推，即已知某個(gè)參數(shù)能使這些樣本出現(xiàn)的概率極大，我們就直接把該參數(shù)作為參數(shù)估計(jì)的真實(shí)值。

EM算法引入

? ? 上述基于全校學(xué)生身高服從一個(gè)分布的假設(shè)過于籠統(tǒng)，實(shí)際上該校男女生的身高分布是不一樣的。其中男生的身高分布為，女生的身高分布為。現(xiàn)在我們估計(jì)該校學(xué)生的分布，就不能簡(jiǎn)單的一刀切了。

? ? 你可以說我們分別抽選50個(gè)男生和50個(gè)女生，對(duì)其分開進(jìn)行估計(jì)。但大多數(shù)情況下，我們并不知道抽樣得到的這個(gè)樣本是來自于男生還是女生。如果說學(xué)生的身高的觀測(cè)變量（Observable Variable），那么樣本性別就是一種隱變量（Hidden Variable）。

? ? 在這種情況下，我們需要估計(jì)的問題包括兩個(gè)：一個(gè)是這個(gè)樣本是男的還是女的，二是男生和女生對(duì)應(yīng)身高的正態(tài)分布參數(shù)分別是多少。這種情況下常規(guī)的極大似然估計(jì)就不太好使了，要估計(jì)男女身高分布，那必須先估計(jì)該學(xué)生是男還是女，反過來要估計(jì)該學(xué)生是男還是女，又得從身高來判斷（男生身高相對(duì)較高，女生身高相對(duì)較矮）。但二者相互依賴，直接用極大似然估計(jì)沒法算。

? ? 針對(duì)這種包含隱變量的參數(shù)估計(jì)問題，一般使用EM算法來進(jìn)行求解。針對(duì)上述身高估計(jì)問題，EM算法的求解思想認(rèn)為：既然兩個(gè)問題相互依賴，肯定是一個(gè)動(dòng)態(tài)求解過程。不如我們直接給定男生女生身高的分布初始值，根據(jù)初始值估計(jì)每個(gè)樣本是男還是女的概率（E步），然后據(jù)此使用極大似然估計(jì)男女生的身高分布參數(shù)（M步），之后動(dòng)態(tài)迭代調(diào)整直到滿足終止條件為止。

EM算法

? ? 所以EM算法的應(yīng)用場(chǎng)景就是要解決包含隱變量的概率模型參數(shù)估計(jì)問題。給定觀測(cè)變量數(shù)據(jù)，隱變量數(shù)據(jù)，聯(lián)合概率分布以及關(guān)于隱變量的條件分布，使用EM算法對(duì)模型參數(shù)進(jìn)行估計(jì)流程如下：

• 初始化模型參數(shù)，進(jìn)行迭代：
• E步：記為第次迭代參數(shù)的估計(jì)值，在第次迭代的E步，計(jì)算函數(shù)：
  其中為給定觀測(cè)數(shù)據(jù)和當(dāng)前參數(shù)估計(jì)下隱變量數(shù)據(jù)的條件概率分布；
• M步：求使函數(shù)最大化的參數(shù)，確定第次迭代的參數(shù)估計(jì)值：
• 重復(fù)迭代E步和M步直至收斂。

? ? 由EM算法過程我們可以看到，其關(guān)鍵在于E步要確定函數(shù)，E步在固定模型參數(shù)的情況下來估計(jì)隱變量分布，而M步則是固定隱變量來估計(jì)模型參數(shù)。二者交互進(jìn)行，直至滿足算法收斂條件。

三硬幣模型

? ? EM算法的一個(gè)經(jīng)典例子是三硬幣模型。假設(shè)有A，B，C三枚硬幣，其出現(xiàn)正面的概率分別為，和。使用三枚硬幣進(jìn)行如下試驗(yàn)：先拋擲硬幣A，根據(jù)其結(jié)果來選擇硬幣B或者C，假設(shè)正面選B，反面選C，然后記錄硬幣結(jié)果，正面記為1，反面記為0，獨(dú)立重復(fù)5次試驗(yàn)，每次試驗(yàn)重復(fù)拋擲B或者C10次。問如何估計(jì)三枚硬幣分別出現(xiàn)正面的概率。

三硬幣模型

? ? 由于我們只能觀察到最后的拋擲結(jié)果，至于這個(gè)結(jié)果是由硬幣A拋出來的還是由硬幣B拋出來的，我們無從知曉。所以這個(gè)過程中依概率選擇哪一個(gè)硬幣拋擲就是一個(gè)隱變量。因此我們需要使用EM算法來進(jìn)行求解。

? ? E步：初始化B和C出現(xiàn)正面的概率為和，估計(jì)每次試驗(yàn)中選擇B還是C的概率（即硬幣A是正面還是反面的概率），例如選擇B的概率為：

? ? 相應(yīng)的選擇C的概率為。計(jì)算出每次試驗(yàn)選擇B和C的概率，然后根據(jù)試驗(yàn)數(shù)據(jù)進(jìn)行加權(quán)求和。

? ? M步：更新模型參數(shù)的新估計(jì)值。根據(jù)函數(shù)求導(dǎo)來確定參數(shù)值：

? ? 對(duì)上式求導(dǎo)并令為0可得第一次迭代后的參數(shù)值：，然后重復(fù)進(jìn)行第二輪、第三輪...直至模型收斂。

EM算法簡(jiǎn)易實(shí)現(xiàn)

? ? 下面我們用numpy來實(shí)現(xiàn)一個(gè)簡(jiǎn)單的EM算法過程來求解三硬幣問題。完整代碼如下：

import numpy as np
def em(data, thetas, max_iter=50, eps=1e-3):    '''    data：觀測(cè)數(shù)據(jù)    thetas：估計(jì)參數(shù)    max_iter：最大迭代次數(shù)    eps：收斂閾值    '''    # 初始化似然函數(shù)值    ll_old = -np.infty    for i in range(max_iter):        ### E步：求隱變量分布        # 對(duì)數(shù)似然        log_like = np.array([np.sum(data * np.log(theta), axis=1) for theta in thetas])        # 似然        like = np.exp(log_like)        # 求隱變量分布        ws = like/like.sum(0)        # 概率加權(quán)        vs = np.array([w[:, None] * data for w in ws])        ### M步：更新參數(shù)值        thetas = np.array([v.sum(0)/v.sum() for v in vs])        # 更新似然函數(shù)        ll_new = np.sum([w*l for w, l in zip(ws, log_like)])        print("Iteration: %d" % (i+1))        print("theta_B = %.2f, theta_C = %.2f, ll = %.2f"               % (thetas[0,0], thetas[1,0], ll_new))        # 滿足迭代條件即退出迭代        if np.abs(ll_new - ll_old) < eps:            break        ll_old = ll_new    return thetas

基于我們編寫的EM算法函數(shù)來對(duì)三硬幣問題進(jìn)行求解：

# 觀測(cè)數(shù)據(jù)，5次獨(dú)立試驗(yàn)，每次試驗(yàn)10次拋擲的正反次數(shù)# 比如第一次試驗(yàn)為5次正面5次反面observed_data = np.array([(5,5), (9,1), (8,2), (4,6), (7,3)])# 初始化參數(shù)值，即硬幣B的正面概率為0.6，硬幣C的正面概率為0.5thetas = np.array([[0.6, 0.4], [0.5, 0.5]])eps = 0.01max_iter = 50thetas = em(observed_data, thetas, max_iter=50, eps=1e-3)

迭代過程如下：

Iteration: 1theta_B = 0.71, theta_C = 0.58, ll = -32.69Iteration: 2theta_B = 0.75, theta_C = 0.57, ll = -31.26Iteration: 3theta_B = 0.77, theta_C = 0.55, ll = -30.76Iteration: 4theta_B = 0.78, theta_C = 0.53, ll = -30.33Iteration: 5theta_B = 0.79, theta_C = 0.53, ll = -30.07Iteration: 6theta_B = 0.79, theta_C = 0.52, ll = -29.95Iteration: 7theta_B = 0.80, theta_C = 0.52, ll = -29.90Iteration: 8theta_B = 0.80, theta_C = 0.52, ll = -29.88Iteration: 9theta_B = 0.80, theta_C = 0.52, ll = -29.87Iteration: 10theta_B = 0.80, theta_C = 0.52, ll = -29.87Iteration: 11theta_B = 0.80, theta_C = 0.52, ll = -29.87Iteration: 12theta_B = 0.80, theta_C = 0.52, ll = -29.87

? ? 可以看到，算法在第七次時(shí)達(dá)到收斂，最后硬幣B和C的正面概率分別為0.8和0.52。

? ? 關(guān)于EM算法，本文沒有太多深入的解釋，關(guān)于似然函數(shù)下界的推導(dǎo)，EM算法的多種解釋等，感興趣的朋友可以自行查找資料學(xué)習(xí)。