Python機(jī)器學(xué)習(xí)算法實(shí)現(xiàn)

Author：louwill

Machine Learning Lab

? ? ?

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

極大似然估計(jì)

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

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

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

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

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

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

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

EM算法引入

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

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

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

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

EM算法

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

• 初始化模型參數(shù)，進(jìn)行迭代：

• E步：記為第次迭代參數(shù)的估計(jì)值，在第次迭代的E步，計(jì)算函數(shù)：
  其中為給定觀測數(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算法的一個經(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é)果，至于這個結(jié)果是由硬幣A拋出來的還是由硬幣B拋出來的，我們無從知曉。所以這個過程中依概率選擇哪一個硬幣拋擲就是一個隱變量。因此我們需要使用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ù)值：

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

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

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

import numpy as npdef em(data, thetas, max_iter=50, eps=1e-3):'''data：觀測數(shù)據(jù)thetas：估計(jì)參數(shù)max_iter：最大迭代次數(shù)eps：收斂閾值'''# 初始化似然函數(shù)值ll_old = -np.inftyfor i in range(max_iter):### E步：求隱變量分布# 對數(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:breakll_old = ll_newreturn thetas

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

# 觀測數(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.5 thetas = np.array([[0.6, 0.4], [0.5, 0.5]]) eps = 0.01 max_iter = 50 thetas = em(observed_data, thetas, max_iter=50, eps=1e-3)

迭代過程如下：

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

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

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

參考資料：

李航 統(tǒng)計(jì)學(xué)習(xí)方法 第二版

https://zhuanlan.zhihu.com/p/36331115

往期精彩：

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法26：隨機(jī)森林

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法25：CatBoost

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法24：LightGBM

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法23：kmeans聚類

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法22：最大熵模型

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法21：馬爾科夫鏈蒙特卡洛

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法20：LDA線性判別分析

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法19：PCA降維

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法18：奇異值分解SVD

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法17：XGBoost

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法16：Adaboost

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法15：GBDT

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法14：Ridge嶺回歸

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法13：Lasso回歸

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法12：貝葉斯網(wǎng)絡(luò)

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法11：樸素貝葉斯

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法10：線性不可分支持向量機(jī)

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法8-9：線性可分支持向量機(jī)和線性支持向量機(jī)

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法7：神經(jīng)網(wǎng)絡(luò)

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法6：感知機(jī)

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法5：決策樹之CART算法

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法4：決策樹之ID3算法

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法3：k近鄰

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法2：邏輯回歸

數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法1：線性回歸

往期精彩回顧適合初學(xué)者入門人工智能的路線及資料下載機(jī)器學(xué)習(xí)及深度學(xué)習(xí)筆記等資料打印機(jī)器學(xué)習(xí)在線手冊深度學(xué)習(xí)筆記專輯《統(tǒng)計(jì)學(xué)習(xí)方法》的代碼復(fù)現(xiàn)專輯 AI基礎(chǔ)下載機(jī)器學(xué)習(xí)的數(shù)學(xué)基礎(chǔ)專輯獲取一折本站知識星球優(yōu)惠券，復(fù)制鏈接直接打開：https://t.zsxq.com/yFQV7am本站qq群1003271085。加入微信群請掃碼進(jìn)群：

總結(jié)

以上是生活随笔為你收集整理的【机器学习基础】数学推导+纯Python实现机器学习算法27：EM算法的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。