EM算法也稱期望最大化（Expectation-Maximum，簡稱EM）算法，它是一個基礎算法，是很多機器學習領域的基礎，比如隱式馬爾科夫算法（HMM），LDA主題模型的變分推斷算法等等。本文對於EM算法，我們主要從以下三個方向學習：

  • 1，最大似然
  • 2，EM算法思想及其推導
  • 3，GMM（高斯混合模型）

1，最大似然概率

我們經常會從樣本觀察數據中，找到樣本的模型參數。最常用的方法就是極大化模型分佈的對數似然函數。怎麼理解呢？下面看我一一道來。

假設我們需要調查我們學習的男生和女生的身高分佈。你會怎麼做啊？你說那麼多人不可能一個一個的去問吧，肯定是抽樣了。假設你在校園隨便捉了100個男生和100個女生。他們共200個人（也就是200個身高的樣本數據，爲了方便表示，下面我們說的“人”的意思就是對應的身高）都在教室了。那麼下一步怎麼辦？你開始喊：“男的左邊，女的右邊”，然後就統計抽樣得到的100個男生的身高。假設他們的身高是服從高斯分佈的。但是這個分佈的均值 μ 和方差 σ 2 我們不知道，這兩個參數就是我們要估計的。記做 θ = [ μ,   δ ] T 。

用數學的語言來說就是：在學校那麼多男生（身高）中，我們獨立地按照概率密度 p(x | θ) ，我們知道了是高斯分佈 N（μ,   δ）的形式，其中的未知參數是 θ = [ μ,   δ ] T 。抽到的樣本集是 X={X1, X2, X3.....Xn}，其中 Xi 表示抽到的第 i 個人的身高，這裏的 N 就是100，表示抽到的樣本個數。

由於每個樣本都是獨立地從 p(x | θ)  中抽取的，換句話說這 100 個男生中的任何一個，都是隨便抽取的，從我們的角度來看這些男生之間是沒有關係的。那麼，我們從學習這麼多男生中爲什麼恰好抽到了這 100個人呢？抽到這 100 個人的概率是多少呢？因爲這些男生（的身高）是服從同一個高斯分佈 p(x | θ)的。那麼我們抽到男生A（的身高）的概率是 p(x A | θ)，抽到男生B的概率是 p(x B | θ)，那因爲他們是獨立的，所以很明顯，我們同時抽到男生A和男生B的概率是 p(x A | θ) * p(x B | θ)，同理，我同時抽到這 100 個男生的概率就是他們各自概率的乘積了。用數學家的口吻說就是從分佈是 p(x | θ) 的總體樣本中抽取到這 100 個樣本的概率，也就是樣本集 X 中各個樣本的聯合概率，用下式表示：

這個概率反映了，在概率密度函數的參數是 θ 時，得到 X 這組樣本的概率。因爲這裏 X 是已知的，也就是說我抽取到這100個人的身高是可以測出來的，也就是已知的。而 θ 是未知的，則上面這個公式只有 θ 是未知的，所以他是關於 θ 的函數。這個函數反映的是在不同的參數 θ 的取值下，取得當前這個樣本集的可能性，因此稱爲參數 θ 相對於樣本集 X 的似然函數（likehood  function），記爲 L(θ)。

這裏出現了一個概念，似然函數。還記得我們的目標嗎？我們需要在已經抽到這一組樣本X的條件下，估計參數  θ 的值。怎麼估計呢？似然函數有什麼用呢？我們下面先了解一下似然的概念。

直接舉個例子：

某位同學與一位獵人一起外出打獵，一隻野兔從前方竄過。只聽一聲槍響，野兔應聲倒下，如果要你推測，這一發命中的子彈是誰打的？你就會想，只發一槍便打中，由於獵人命中的概率一般大於這位同學命中的概率，所以這一槍應該是獵人射中的。

這個例子所作的推斷就體現了極大似然法的基本思想。下面再說一個最大似然估計的例子。

再例如：假設你去賭場了，但是不知道能不能賺錢，你就在門口堵着出來，出來一個人你就問人家賺錢還是賠錢，如果問了五個人，假設五個人都說自己賺錢了，那麼你就會認爲，賺錢的概率肯定是非常大的。

從上面兩個例子，你得到了什麼結論？

總的來說： 極大似然估計就是估計模型參數的統計學方法 。

回到男生身高的例子。在學校那麼多男生中，我一抽就抽到這 100 個男生（表示身高），加入不是其他人，那是不是表示在整個學校中，這 100 個人（的身高）出現的概率最大啊，那麼這個概率怎麼表示呢？就是上面那個似然函數 L(θ)。所以我們就只需要找到一個參數 θ ，其對應的似然函數 L(θ) 最大，也就是說抽到這 100 個男生（的身高）的概率最大，這個叫做 θ 的最大似然估計量，記爲：

有時，可以看到 L(Θ) 是連乘的，所以爲了便於分析，還可以定義對數似然函數，將其變爲連加的：

現在我們知道了，要求出 Θ ，只需要使 Θ 的似然函數 L(Θ) 極大化，然後極大值對應的 Θ 就是我們的估計。這裏就回到了求最值的問題了。怎麼求一個函數的最值？當然是求導，然後令導數爲零，那麼解這個方程得到的 Θ  就是了（前提是函數L(Θ)連續可微）。那如果 Θ 是包含多個參數的向量，那麼如何處理呢？當然是求 L(Θ) 對所有參數的偏導數，也就是梯度了，那麼 n 個未知的參數，就有 n 個方程，方程組的解就是似然函數的極值點了，當然就得到了這 n 個參數了。

最後使用PPT整理一下，更直觀：

所以，其實最大似然估計你可以把它看做是一個反推。多數情況下我們是根據已知條件來推斷結果，而最大似然估計是已經知道了結果，然後尋找使該結果出現的可能性最大的條件，以此作爲估計值。比如，如果其他條件一定的話，抽菸者發生肺癌的危險是不抽菸者的五倍，那麼如果我們知道有個人是肺癌，我想問你這個人是抽菸還是不抽菸。你如何判斷，你可能對這個人一無所知，你所知道的只有一件事，那就是抽菸更容易發生肺癌，那麼你會猜測這個人不抽菸嗎？我相信你更可能會說，這個人抽菸，爲什麼？這就是“最大可能”，我只能說他“最有可能”是抽菸的，“他是抽菸的”這一估計纔是“最有可能”得到“肺癌” 這樣的結果，這就是最大似然估計。

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

而求最大似然估計值的一般步驟

  • 1，寫出似然函數
  • 2，對似然函數取對數，並整理
  • 3，求導數，令導數爲0，得到似然方程
  • 4，解似然方程，得到的參數即爲所求

2，EM算法要解決的問題

上面說到，通過樣本觀察數據中，找到樣本的模型參數。但是在一些情況下，我們得到的觀察數據有未觀察到的隱含數據，此時我們未知的有隱含數據和模型參數，因而無法直接用極大化對數似然函數得到模型分佈的參數。怎麼辦呢？這就是EM算法可以派上用場的地方了。

下面，我們重新回到上面那個身高分佈估計的問題。現在，通過抽取得到的那100個男生的身高和已知的其身高服從高斯分佈，我們通過最大化其似然函數，就可以得到了對應高斯分佈的參數 θ = [ μ,   δ ] T 了。那麼，對於我們學習的女生的身高分佈也可以用同樣的方法得到了。

再回到例子本身，如果沒有“男的左邊，女的右邊”這個步驟，或者說這二百個人已經混到一起了，這時候，你從這二百個人（的身高）裏面隨便給我指一個人（的身高），我都無法確定這個的身高是男生的身高，還是女生的身高。也就是說你不知道抽取的那二百個人裏面的每一個人到底是從男生的那個身高分佈裏面抽取的，還是女生的那個身高分佈抽取的。用數學的語言就是： 抽取的每個樣本都不知道是從那個分佈抽取的 。

這個時候，對於每一個樣本或者你抽取到的人，就有兩個東西需要猜測或者估計的了， 一是這個人是男生還是女的？二是男生和女生對應的身高的高斯分佈的參數是多少 ？

所以問題就難了。

只有當我們知道了那些人屬於同一個高斯分佈的時候，我們才能夠對這個分佈的參數做出靠譜的預測，例如剛開始的最大似然所說的，但現在兩種高斯分佈的人混在一塊了，我們又不知道那些人屬於第一個高斯分佈，那些人屬於第二個，所以就沒法估計這兩個分佈的參數。反過來，只有當我們對這兩個分佈的參數做出了準確的估計的時候，才能知道到底那些人屬於第一個分佈，那些人屬於第二個分佈。

這就成了一個先有雞還是先有蛋的問題了。雞說，沒有我，誰把你生出來的啊。蛋不服，說，沒有我你從哪蹦出來的啊。爲了解決這個你依賴我，我依賴你的循環依賴問題，總得有一方要先打破僵局，說，不管了，我先隨便整一個值出來，看你怎麼變，然後我再根據你的變換調整我的變換，然後如此迭代着不斷互相推導，最終就會收斂到一個解，這就是EM算法的基本思想了。

不知道大家能理解其中的思想嗎，然後再舉個例子。

例如小時候，老媽給一大袋子糖果給你，叫你和你姐姐等分，然後你懶得去點糖果的個數，所以你也就不知道每個人到底該分多少個，咱們一般怎麼做呢？先把一袋糖果目測的分爲兩袋，然後把兩袋糖果拿在左右手，看那個重，如果右手重，那麼很明顯右手這袋糖果多了，然後你再在右手這袋糖果中抓一把放到左手這袋，然後再感受下那個重，然後再從重的那袋抓一小把放入輕的那一袋，繼續下去，直到你感覺兩袋糖果差不多相等了爲止。然後爲了體現公平，你還讓你姐姐先選。

EM算法就是這樣，假設我們想估計A和B兩個參數，在開始狀態下二者都是未知的，但如果知道了A的信息就可以得到B的信息，反過來知道了B也就得到了A。可以考慮首先賦予A某種初值，以此得到B的估計值，然後從B的當前值出發，重新估計A的取值，這個過程一直持續到收斂爲止。

EM 的意思是“Expectation  Maximization”，在我們上面這個問題裏面，我們先隨便猜一下男生（身高）的正態分佈的參數：如均值和方差。例如男生的均值是1米7，方差是 0.1 米（當然了，剛開始肯定沒有那麼準），然後計算出每個人更可能屬於第一個還是第二個正態分佈中的（例如，這個人的身高是1米8，那很明顯，他最大可能屬於男生的那個分佈），這個是屬於Expectation 一步。有了每個人的歸屬，或者說我們已經大概地按上面的方法將這 200 個人分爲男生和女生兩部分，我們就可以根據之前說的最大似然那樣，通過這些被大概分爲男生的 n 個人來重新估計第一個分佈的參數，女生的那個分佈同樣方法重新估計。這個就是 Maximization。然後，當我們更新了這兩個分佈的時候，每一個屬於這兩個分佈的概率又變了，那麼我們就再需要調整 E 步、.... 如此往復，直到參數基本不再發生變化爲止。

這裏把每個人（樣本）的完整描述看做是三元組  yi = {xi,  zi1,  zi2}，其中， xi 是第 i 個樣本的觀測值，也就是對應的這個人的身高，是可以觀測到的值。zi1和zi2表示男生和女生這兩個高斯分佈中哪個被用來產生值 xi，就是說這兩個值標記這個人到底是男生還是女生（的身高分佈產生的）。這兩個值我們是不知道的，是隱含變量。確切的說， zij 在 xi 由第 j個高斯分佈產生時值爲1，否則爲0。例如一個樣本的觀測值爲 1.8，然後他來自男生的那個高斯分佈，那麼我們可以將這個樣本表示爲 {1.8,  1, 0}。如果 zi1 和 zi2的值已知，也就是說每個人我已經標記爲男生或者女生了，那麼我們就可以利用上面說的最大似然算法來估計他們各自高斯分佈的參數。但是他們未知，所以我們只能用 EM算法。

咱們現在不是因爲那個噁心的隱含變量（ 抽取得到的每隔樣本都不知道是那個分佈抽取的 ）使得本來簡單的可以求解的問題變複雜了，求解不了嗎。那怎麼辦呢？人類解決問題的思路都是想能否把複雜的問題簡單化。好，那麼現在把這個複雜的問題逆回來，我假設已經知道這個隱含變量了，哎，那麼求解那個分佈的參數是不是很容易了，直接按上面說的最大似然估計就好了。那你就問我了，這個隱含變量是未知的，你怎麼就來一個假設說已知呢？你這種假設是沒有依據的。所以我們可以先給這個分佈弄一個初始值，然後求這個隱含變量的期望，當成是這個隱含變量的已知值，那麼現在就可以用最大似然求解那個分佈的參數了吧，那假設這個參數比之前的那個隨機的參數要好，它更能表達真實的分佈，那麼我們再通過這個參數確定的分佈去求這個隱含變量的期望，然後再最大化，得到另一個更優的參數，...迭代，就能得到一個皆大歡喜的結果了。

這裏最後總結一下EM算法的思路：

所以EM算法解決問題的思路就是使用啓發式的迭代方法，既然我們無法直接求出模型分佈參數，那麼我們可以先猜想隱含數據（EM算法的E步），接着基於觀察數據和猜測的隱含數據一起來極大化對數似然，求解我們的模型參數（EM算法的M步）。由於我們之前的隱藏數據是猜測的，所以此時得到的模型參數一般還不是我們想要的結果。不過沒關係，我們基於當前得到的模型參數，繼續猜測隱含數據（EM算法的E步），然後繼續極大化對數似然，求解我們的模型參數（EM算法的M步）。以此類推，不斷的迭代下去，直到模型分佈參數基本無變化，算法收斂，找到合適的模型參數。

下面再看看本文中EM算法最後一個例子： 拋硬幣

Nature Biotech在他的一篇EM tutorial文章《Do, C. B., & Batzoglou, S. (2008). What is the expectation maximization algorithm?. Nature biotechnology, 26(8), 897.》中，用了一個投硬幣的例子來講EM算法的思想。

比如兩枚硬幣A和B，如果知道每次拋的是A還是B，那可以直接估計（見下圖a）。

如果不知道拋的是A還是B（這時就刺激了吧，對的，這就是咱們不知情的隱變量），只觀測到5輪循環每輪循環10次，共計50次投幣的結果，這時就沒法直接估計A和B的正面概率。這時，就輪到EM算法出場了（見下圖b）。

可能直接看，這個例子不好懂，下面我們來通俗化這個例子。

還是兩枚硬幣A和B，假定隨機拋擲後正面朝上概率分別爲PA，PB。我們就不拋了，直接使用上面的結果，每一輪都是拋10次，總共5輪：

硬幣A被拋了30次，在第二輪，第三輪，第五輪分別出現了9次，8次，7次正，所以很容易估計出PA，類似的PB 也很容易計算出來，如下：

但是問題來了，如果我們不知道拋的硬幣是A還是B呢（即硬幣種類是隱變量），然後再輪流拋五次，得到如下的結果呢？

下面問題變得有意思了，現在我們的目標還是沒有變，還是估計PA和PB，需要怎麼做呢？

思路如下：

顯然這裏我們多了一個硬幣種類的隱變量，設爲 z，可以把它看做一個10維的向量（z1, z2, ...z10），代表每次投擲時所使用的硬幣，比如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是一個相當靠譜的估計！

這就是說，我們初始化的PA和PB，按照最大似然概率就可以估計出Z，然後基於Z，按照最大似然概率可以反過來估計出P1和P2，當與我們初始化的PA和PB一樣時，說明是P1和P2很可能就是真實的值。這裏包含了兩個交互的最大似然估計。

如果新估計出來的PA和PB和我們初始化的值差別很大，怎麼辦呢？就是繼續用新的P1和P2迭代，直至收斂。

做法如下：

我們不妨，這樣，先隨便給PA和PB賦一個值，比如：

所以我們的初始值 theta 爲：

然後，我們看看第一輪拋擲最可能是哪個硬幣，（五正五負）：

其實上式很簡單就是：PA= 0.6*0.6*0.6*0.6*0.6*0.4*0.4*0.4*0.4*0.4 = 0.0007962624000000002（後面類似）

PB= 0.5**10 = 0.0009765625

然後求出選擇PA的概率和PB的概率值

然後依次求出其他四輪中相應的概率。做成表格如下：

然後我們通過上面的概率值來估計Z，大概是這樣的：

這樣我們通過最大似然估計法則就估算出來PA和PB的值了，如下：

然後我們再次迭代，繼續估算，假設我們估算到第10次，就算出PA和PB的值是正確的了，如下：

下面我們來對比初始化的PA和PB和新估計的PA和PB，和最終的PA和PB：

我們可以發現，我們估計的PA和PB相比於他們的初始值，更接近於他們的真實值，就這樣不斷迭代，不斷接近真實值，這就是EM算法的奇妙之處。

可以期待，我們繼續按照上面的思路，用估計出的PA和PB再來估計Z，再用Z來估計新的PA和PB，反覆迭代下去，就可以最終得到PA和PB的真值了。而這裏我們假設爲 0.8 和 0.52。此時無論怎麼迭代，PA和PB的值都會保持在 0.8和0.52不變，於是，我們就找到了PA和PB的最大似然估計。

下面在推導EM算法之前，我們先學習一下Jensen不等式

3，導數性質和Jensen不等式

在學習Jensen不等式之前，我們先說說二階導數的一些性質。

3.1，二階導數的性質1

如果有一個函數 f(x) 在某個區間 I 上有 f ''(x)（即二階導數） > 0 恆成立，那麼對於區間 I上的任意 x,  y 總有：

如果總有 f''(x) < 0 成立，那麼上式的不等號反向。

幾何的直觀解釋：如果有一個函數 f(x) 在某個區間 I 上有 f ''(x)（即二階導數） > 0 恆成立，那麼在區間 I 上 f(x)的圖像上的任意兩點連出的一條線段，這兩點之間的函數圖像都在該線段的下方，反之在該線段的上方。

3.2，二階導數的性質2：判斷函數極大值以及極小值

結合一階，二階導數可以求函數的極值。當一階導數等於0，而二階導數大於0時，爲極小值點。當一階導數等於0，而二階導數小於0時，爲極大值點；當一階導數和二階導數都等於0，爲 駐點 。

駐點：在微積分，駐點（Stationary Point）又稱爲平穩點、穩定點或臨界點（Critical Point）是函數的一階導數爲零，即在“這一點”，函數的輸出值停止增加或減少。對於一維函數的圖像，駐點的切線平行於x軸。對於二維函數的圖像，駐點的切平面平行於xy平面。值得注意的是，一個函數的駐點不一定是這個函數的極值點（考慮到這一點左右一階導數符號不改變的情況）；反過來，在某設定區域內，一個函數的極值點也不一定是這個函數的駐點（考慮到邊界條件），駐點（紅色）與拐點（藍色），這圖像的駐點都是局部極大值或局部極小值

3.3，二階導數的性質3：函數凹凸性

設 f(x) 在 [a, b] 上連續，在（a, b）內具有一階和二階導數，那麼：

  • （1）若在（a, b）內 f ''(x) > 0，則 f(x) 在[a, b]上的圖形是凹的
  • （2）若在（a, b）內 f ''(x) < 0，則 f(x) 在[a, b]上的圖形是凸的

3.4 期望

在概率論和統計學中，數學期望（mean）（或均值，也簡稱期望）是實驗中每次可能結果的概率乘以其結果的總和，是最基本的數學特徵之一。它反映隨機變量平均取值的大小。

需要注意的是：期望值並不一定等同於常識中的“期望”——“期望值”也許與每一個結果都不相等。期望值是該變量輸出值的平均數。期望值並不一定包含於變量的輸出值集合裏。

大數定理規定：隨着重複次數接近無窮大，數值的算術平均值幾乎肯定地收斂於期望值。

離散型：

如果隨機變量只取得有限個數或無窮能按一定次序一一列出，其值域爲一個或若干個有限或無限區域，這樣的隨機變量稱爲離散型隨機變量。

離散型隨機變量只取得有限個值或無窮能按一定次序一一列出，其值域爲一個或若干個有限或無限區間，這樣的隨機變量稱爲離散型隨機變量。

離散型隨機變量的一切可能的取值 xi 與對應的概率 p(xi) 乘積之和稱爲該離散型隨機變量的數學期望（若該求和絕對收斂），記爲E(x)。它是簡單算術平均的一種推廣，類似加權平均。

公式：

離散型隨機變量X的取值爲X1, X2, X3....Xn， p(X1)，p(X2), p(X3)....P(Xn)爲X對應取值的概率，可以理解爲數據 X1, X2, X3.....Xn出現的頻率 f(Xi)，則：

例子：

某城市有10萬個家庭，沒有孩子的家庭有 1000 個，有一個孩子的家庭有 9 萬個，有兩個孩子的家庭有 6000個，有 3個孩子的家庭有 3000個。

則此城市中任一個家庭中孩子的數目是一個隨機變量，記爲X。它可取值 0, 1， 2， 3 。

其中，X取 0 的概率爲 0.01 ，取 1的概率爲 0.9，取 2的概率爲0.06，取 3的概率爲0.03

則它的數學期望爲：

即此城市一個家庭平均有小孩 1.11 個，當然人不可能用1.11來算，約等於2個。

設 Y 是隨機變量X的函數：

它的分佈律爲：

連續型：

設連續型隨機變量 X 的概率密度函數爲 f(x)，若積分絕對收斂，則稱積分的值爲隨機變量的數學期望，記爲 E(x)，即：

若隨機變量 X的分佈函數 F(x) 可表示成一個非負可積函數 f(x)的積分，則稱 X爲連續型隨機變量，f(x)稱爲 X 的概率密度函數（分佈密度函數）。

數學期望 E(X) 完全由隨機變量 X 的概率分佈所確定。若X服從某一分佈，也稱 E(X) 是這一分佈的數學期望。

定理：

若隨機變量 Y 符合函數 Y = g(x)，且下面函數絕對收斂，則有下下面的期望

該定理的意義在於：我們求 E(Y) 時不需要計算出 Y的分佈律或概率分佈，只要利用X的分佈律或者概率密度即可。

上述定理還可以推廣到兩個或以上隨機變量的函數情況。

設 Z 是隨機變量 X， Y 的函數 Z = g(X, Y)（g是連續函數），Z是一個一維隨機變量，二維隨機變量（X, Y）的概率密度爲 f(x,  y)，則有：

離散型隨機變量與連續型隨機變量的區別

離散型隨機變量與連續型隨機變量都是由隨機變量取值範圍確定。

變量取值只能取離散型的自然數，就是離散型隨機變量。例如，一次擲20個硬幣，k額硬幣正面朝上，k是隨機變量。k的取值只能是自然數0， 1, 2 ，3， 4， 。。。。20。而不能取小數 3.5 ，無理數 √ 20 ，因此k是離散型隨機變量。

如果變量可以在某個區間內取任一實數，即變量的取值可以是連續的，這隨機變量就稱爲連續型變量。例如公共汽車每 15 分鐘一班，某人在站臺等車時間 x 是個隨機變量，x的取值範圍是 [0, 15)，它是一個區間，從理論上說在這個區間內可取任一實數 3.5 ，無理數 ，√ 20 等，因而稱這隨機變量是連續型隨機變量。

3.5  Jensen不等式

如果用圖表示會很清晰：

4，EM算法的推導

從上面的總結，我們知道EM算法就是迭代求解最大值的算法，同時算法在每一次迭代時分爲兩步，E步和M步。一輪輪迭代更新隱含數據和模型分佈參數，直到收斂，即得到我們需要的模型參數。

這時候你就不服了，說你老迭代迭代的，你咋知道新的參數的估計就比原來的好啊？爲什麼這種方法行得通呢？有沒有失效的時候呢？什麼時候失效呢？用到這個方法需要注意什麼問題呢？

爲了說明這一串問題，我們下面學習EM算法的整個推導過程。

首先，我們接着上面的問題聊，就是100名學生的身高問題，總結一下問題，如下：

樣本集 {x(1),   x(2) , .... x(m)}，包含 m 個 獨立 的樣本，其中每個樣本 i 對應的類別 z(i) 是未知的，所以很難用極大似然求解。而我們的目的就是想找到每個樣例隱含的類別z，能使得p(x, z) 最大。

p(x, z) 的最大似然估計如下：

上式中，我們要考慮每個樣本在各個分佈中的情況。本來正常求偏導就可以了，但是因爲有隱含數據的存在，也就是log後面還有求和，這個就沒有辦法直接求出 theta了。

對於每一個樣例 i，讓Qi 表示該樣例隱含變量 z 的某種分佈，Qi 滿足的條件是：

如果 z 是連續性的，那麼Qi是概率密度函數，需要將求和符號換做積分符號。繼續上面的例子，假設隱藏變量 z是身高，那就是連續的高斯分佈；假設隱藏變量是男女，那就是伯努利分佈 。

所以說如果我們確定了隱藏變量z，那求解就容易了，可是怎麼求解呢？這就需要一些特殊的技巧了。

技巧： Jensen不等式應用於凹函數時，不等號方向反向 。

我們知道：

所以對右式分子分母同時乘以 Q(z)：

爲什麼要這麼做呢？說白了就是要湊Jensen不等式（Q(z) 是Z的分佈函數），下面繼續

這裏假設：

則上式可以寫爲：

我們由Jensen不等式可以知道（由於對數函數是凹函數）：

則代入 上式：

所以：

上面這個過程可以看做是對 l(Θ) 求了下界。對於Qi 的選擇，有多種可能，那種更好呢？假設 Θ 已經給定，那麼 l(Θ) 的值就取決於 Q(z)和P(x,  z)了。

4.1 如何優化下界？

因爲下界比較好求，所以我們要優化這個下界來使得似然函數最大。那麼如何優化下界呢？

隨着我們的參數變換，我們的下界是不斷增大的趨勢。我們紅色的線是恒大於零的。簡單來說，就是EM算法可以看做是J 的座標上升法，E步固定 theta，優化 Q，M步固定Q，優化 theta。

迭代到收斂

可能用這個圖更明顯（參考：https://blog.csdn.net/u010834867/article/details/90762296）

我們需要調整上式中Q(z)和P(x,  z) 這兩個概率，使下界不斷上升，以逼近 l(Θ) 的真實值，那麼什麼時候算是調整好了呢？當不等式變成等式時，說明我們調整後的概率能夠等價於  l(Θ) 了。按照這個思路，我們要找到等式成立的條件。根據Jensen不等式，要想讓等式成立，就需要什麼呢？下面繼續說。

4.2  如何能使Jensen不等式等號成立呢？

什麼時候，使得我們的似然函數最大呢？當然是取等號的時候。什麼時候等號成立呢？

由上面Jensen不等式中所說，我們要滿足Jensen不等式的等號，則當且僅當 P(x=E(x)) = 1的時候，也就是說隨機變量X是常量。在這裏，我們則有：

c爲常數，不依賴於 z。對此式子做進一步推導，由於Q(z) 是z的分佈函數，所以滿足：

所有的分子和等於常數C（分母相同），即：

多個等式分子分母相加不變，這個認爲每個樣例的兩個概率比值都是 c，那麼就可以求解 Q(z)了。

4.3  Q(z)的求解

我們知道Q(z) 是z的分佈函數，而且：

由上式就可得C就是 P(xi,  z)對Z求和，我們上面有公式。

從上面兩式，我們可以得到：

Q(z) 代表第 i 個數據是來自 zi 的概率。

至此，我們推出了在固定其他參數 Θ 後，Q(z)的計算公式就是後驗概率，解決了Q(z)如何選擇的問題。這一步就是E步，建立 l(Θ) 的下界。

所以去掉 l(Θ) 式中爲常數的部分，則我們需要極大化的對數似然下界爲：

注意上式也就是我們EM算法的M步，就是在給定 Q(z)後，調整Θ，去極大化 l(Θ)的下界（在固定Q(z)後，下界還可以調整的更大）。

E步我們再強調一次，注意到上式Q(z)是一個分佈，因此 下式可以理解爲 logP(x, z; Θ) 基於條件概率分佈 Q(z)的期望。

至此，我們學習了EM算法那中E步和M步的具體數學含義。

5，EM算法的一般流程

下面我們總結一下EM算法的流程。

輸入：觀察數據 x = {x1,  x2, ... xm}，聯合分佈 p(x, z; θ)，條件分佈 p(z|x; θ)，最大迭代次數 J 。

1）隨機初始化模型參數 θ 的初值 θ 0 。

2）循環重複，即for j from 1 to J  開始EM算法迭代。

E步：計算聯合分佈的條件概率期望

M步：極大化L(θ,  θ j )，得到 θ j+1 ：

3）如果 θ j+1 已經收斂，則算法結束，否則繼續回到步驟E步迭代。

輸出：模型參數 θ。

用PPT總結如下：

6，EM算法的收斂性思考

（這裏直接粘貼劉建平老師的博客了。方便。。）

EM算法的流程並不複雜，但是還有兩個問題需要我們思考：

  • 1） EM算法能保證收斂嗎？
  • 2）EM算法如果收斂，那麼能保證收斂到全局最大值嗎？

首先，我們先看第一個問題，EM算法的收斂性。要證明EM算法收斂，則我們需要證明我們的對數似然函數的值在迭代的過程中一直在增大。即：

由於：

令：

上面兩個式子相減得到：

在上式中分別取 Θ 爲 Θ j 和 Θ j+1 ，並且相減得到：

要證明EM算法的收斂性，我們只需要證明上式的右邊是非負的即可。

由於 Θ j+1  使得L（Θ，Θ j  ）極大，因此有：

上式的意思是說極大似然估計單調遞增，那麼我們會達到最大似然估計的最大值。

而對於第二部分，我們有：

上式中，第二個不等式使用了Jensen不等式，只不過和前面的使用相反而已，最後一個等式使用了概率分佈累積爲1的性質。

至此，我們得到了：

從上面的推導可以看出，EM算法可以保證收斂到一個穩定點，但是卻不能保證收斂到全局的極大值點，因此它是局部最優的算法，當然，如果我們的優化目標 L（Θ，Θ j  ）是凸的，則EM算法可以保證收斂到全局最大值，這點和梯度下降法這樣的迭代算法相同。至此我們也回答了上面的第二個問題。

如果我們從算法思想的角度來思考 EM 算法，我們可以發現我們的算法裏已知的是觀察數據，未知的隱含數據和模型參數。在E步，我們所做的事情是固定模型參數的值，優化隱含數據的分佈，而在M步，我們所做的事情是固定隱含數據分佈，優化模型參數的值。比較下其他的機器學習算法，其實很多算法都有類似的思想。比如SMO算法，座標軸下降法都是由了類似的思想來求解問題。

7，GMM（高斯混合模型）

我們上面已經推導了EM算法，下面再來說一下高斯混合模型。高斯混合模型是什麼呢？顧名思義就是用多個高斯模型來描述數據的分佈，就是說我們的數據可以看做是從多個Gaussian Distribution 中生成出來的。那GMM是由K個Gaussian分佈組成，所以每個Gaussian稱爲一個“Component”。

首先說一句，爲什麼是高斯混合模型呢？而不是其他模型呢？

因爲從中心極限定理知：只要k足夠大，模型足夠複雜，樣本量足夠多，每一塊小區域就可以用高斯分佈描述（簡單說：只要樣本量足夠大時，樣本均值的分佈慢慢變爲正態分佈）。而且高斯函數具有良好的計算性能，所以GMM被廣泛的應用。

有時候單一的高斯分佈不能很好的描述分佈，如下圖，二維高斯密度函數的等概率先爲橢圓（當然後面會做練習）

從上圖我們看到，左圖用單一高斯分佈去描述，顯然沒有右圖用兩個高斯分佈去描述的效果好。

所以我們引入混合高斯模型。

如下圖所示，我們用三個高斯分佈去描述一個二維的數據：

現在我們定義 K 個高斯密度疊加：

對於每個高斯密度函數有自己的均值 μk 和方差 ∑k ， πk作爲混合的比例係數有：

下圖中（a）爲不同混合比例下的高斯概率密度分佈；（b）爲混合狀態下的概率密度分佈；（c）爲概率密度分佈的表面圖：

所以 p(x) 可以改寫爲：

並與公式（1）對比，有：

則後驗概率 p(k|x)，根據貝葉斯理論，可表示爲：

因此GMM由 μ， ∑ ， π 確定，並在有參數 K的存在。

而GMM模型的求解方式和EM算法是一樣的，就是不斷的迭代更新下去。一個最直觀的瞭解算法的思路就是參考K-Means算法。博客地址：

Python機器學習筆記：K-Means算法，DBSCAN算法

這裏再簡單的說一下，在K-Means聚類時，每個聚類簇的質心是隱含數據。我們會假設K個初始化質心，即EM算法的E步；然後計算得到每個樣本最近的質心，並且把樣本聚類到最近的這個質心，即EM算法的M步，重複這個E步和M步，知道質心不再變化爲止，這樣就完成了K-Means聚類。當然K-Means算法是比較簡單的，實際問題並沒有這麼簡單，但是大概可以說明EM算法的思想了。

上面遺留了一個問題，就是單一的高斯分佈不能很好的描述分佈，而兩個高斯分佈描述的好。爲了將這個說的透徹，我們這樣做：我們生成四簇數據，其中兩簇數據離的特別近，然後我們看看K-Means算法是否能將這兩簇離的特別近的數據按照我們的想法區分開，再看看GMM算法呢，如果可以區分，那麼我們將數據拉伸呢？

所以下面我們會做兩個實驗，對於隨機生成的數據，我們會做K-Means算法和GMM算法的效果比對，然後再對數據進行處理，再做K-Means算法和GMM算法的效果比對。

首先生成四簇數據，效果如下：

from sklearn.datasets.samples_generator import make_blobs
from matplotlib import pyplot as plt
from sklearn.cluster import KMeans
from sklearn.mixture import GaussianMixture

# 生成新的數據
X, y_true = make_blobs(n_samples=800, centers=4, random_state=11)
plt.scatter(X[:, 0], X[:, 1])

圖如下：

我們看看K-Means分類效果：

我們再看看GMM算法的效果：

都能分出來，這樣很好，下面我們對數據進行處理： 拉伸 

rng = np.random.RandomState(13)
X_stretched = np.dot(X, rng.randn(2, 2))

我們拉伸後，看看原圖效果：

下面分別使用K-Means算法和GMM算法看效果：

首先是K-Means算法

kmeans = KMeans(n_clusters=4)
kmeans.fit(X)
y_kmeans = kmeans.predict(X)

plt.scatter(X[:, 0], X[:, 1], c=y_kmeans, s=50, cmap='viridis')

效果如下

K-Means算法沒有把拉伸的數據分開，而是按照聚類的思想，將上面的數據聚在一起，將下面的數據聚在一起。

我們再看看GMM算法的效果：

是不是GMM算法區分出來我們想要的效果了。所以在這個問題上GMM算法的優勢就體現出來了。

完整代碼去我的GitHub上拿（地址：https://github.com/LeBron-Jian/MachineLearningNote）。

8，Sklearn實現GMM（高斯混合模型）

Sklearn庫GaussianMixture類是EM算法在混合高斯分佈的實現，現在簡單記錄下其參數說明。

對此源碼中的參數，我們瞭解其意義：

  • 1. n_components:混合高斯模型個數，默認爲1
  • 2. covariance_type:協方差類型，包括{‘full’,‘tied’, ‘diag’, ‘spherical’}四種，分別對應完全協方差矩陣（元素都不爲零），相同的完全協方差矩陣（HMM會用到），對角協方差矩陣（非對角爲零，對角不爲零），球面協方差矩陣（非對角爲零，對角完全相同，球面特性），默認‘full’ 完全協方差矩陣
  • 3. tol：EM迭代停止閾值，默認爲1e-3.
  • 4. reg_covar:協方差對角非負正則化，保證協方差矩陣均爲正，默認爲0
  • 5. max_iter:最大迭代次數，默認100
  • 6. n_init:初始化次數，用於產生最佳初始參數，默認爲1
  • 7. init_params: {‘kmeans’, ‘random’}, defaults to ‘kmeans’.初始化參數實現方式，默認用kmeans實現，也可以選擇隨機產生
  • 8. weights_init:各組成模型的先驗權重，可以自己設，默認按照7產生
  • 9. means_init:初始化均值，同8
  • 10. precisions_init:初始化精確度（模型個數，特徵個數），默認按照7實現
  • 11. random_state :隨機數發生器
  • 12. warm_start :若爲True，則fit（）調用會以上一次fit（）的結果作爲初始化參數，適合相同問題多次fit的情況，能加速收斂，默認爲False。
  • 13. verbose :使能迭代信息顯示，默認爲0，可以爲1或者大於1（顯示的信息不同）
  • 14. verbose_interval :與13掛鉤，若使能迭代信息顯示，設置多少次迭代後顯示信息，默認10次。

下面看看網友對源碼的分析：

（1）首先將模型完全轉換成對數計算，根據高斯密度函數公式分別計算k個組成高斯模型的log值，即logP(x|z)的值
def _estimate_log_gaussian_prob(X, means, precisions_chol, covariance_type):
# 計算精度矩陣的1/2次方log_det（代碼精度矩陣是通過cholesky獲取）
    log_det = _compute_log_det_cholesky(
        precisions_chol, covariance_type, n_features)
# 對應上面四種協方差類型，分別計算精度矩陣與（x-u)相乘那部分log_prob
    if covariance_type == 'full':
        log_prob = np.empty((n_samples, n_components))
        for k, (mu, prec_chol) in enumerate(zip(means, precisions_chol)):
            y = np.dot(X, prec_chol) - np.dot(mu, prec_chol)
            log_prob[:, k] = np.sum(np.square(y), axis=1)

    elif covariance_type == 'tied':
        log_prob = np.empty((n_samples, n_components))
        for k, mu in enumerate(means):
            y = np.dot(X, precisions_chol) - np.dot(mu, precisions_chol)
            log_prob[:, k] = np.sum(np.square(y), axis=1)

    elif covariance_type == 'diag':
        precisions = precisions_chol ** 2
        log_prob = (np.sum((means ** 2 * precisions), 1) -
                    2. * np.dot(X, (means * precisions).T) +
                    np.dot(X ** 2, precisions.T))

    elif covariance_type == 'spherical':
        precisions = precisions_chol ** 2
        log_prob = (np.sum(means ** 2, 1) * precisions -
                    2 * np.dot(X, means.T * precisions) +
                    np.outer(row_norms(X, squared=True), precisions))
# 最後計算出logP(x|z)的值
return -.5 * (n_features * np.log(2 * np.pi) + log_prob) + log_det   

（2）P(x|z)*P(z)計算每個模型的概率分佈P(x,z),求對數則就是相加了
def _estimate_weighted_log_prob(self, X):
 return self._estimate_log_prob(X) + self._estimate_log_weights()

（3）最後開始計算每個模型的後驗概率logP(z|x)，即Q函數
def _estimate_log_prob_resp(self, X):
weighted_log_prob = self._estimate_weighted_log_prob(X)
#計算P(X)
log_prob_norm = logsumexp(weighted_log_prob, axis=1)
 with np.errstate(under='ignore'):
 # 忽略下溢,計算每個高斯模型的後驗概率，即佔比，對數則相減
log_resp = weighted_log_prob - log_prob_norm[:, np.newaxis]
    return log_prob_norm, log_resp
(4)調用以上函數
 #返回所有樣本的概率均值，及每個高斯分佈的Q值
def _e_step(self, X):
log_prob_norm, log_resp = self._estimate_log_prob_resp(X)
    return np.mean(log_prob_norm), log_resp
2.對應M step

def _m_step(self, X, log_resp):
#根據上面獲得的每個高斯模型的Q值(log_resp)。重新估算均值self.means_，協方差self.covariances_，當前符合各高斯模型的樣本數目self.weights_（函數名起的像權重，實際指的是數目）
n_samples, _ = X.shape
self.weights_, self.means_, self.covariances_ = (
            _estimate_gaussian_parameters(X, np.exp(log_resp), self.reg_covar,
                                          self.covariance_type))
#更新當前各高斯模型的先驗概率，即P(Z)
self.weights_ /= n_samples
#根據cholesky分解計算精度矩陣
self.precisions_cholesky_ = _compute_precision_cholesky(
            self.covariances_, self.covariance_type)
然後重複以上過程，就完成了EM算法的實現啦。

下面進入我們的實戰環節。。。

數據（可以去我的GitHub上拿，地址：https://github.com/LeBron-Jian/MachineLearningNote）

數據的背景是這樣的，這裏簡單介紹一下：在某個地區有一個橋，橋東面有一個區域，橋西邊有一個區域。假設我們這兩個地方是存放共享單車的，那給出的數據是每天某個時間點，共享單車被騎走的數據。（此數據是關於時間序列的數據），數據有 42312條，每條數據有兩個特徵（也就是前面提到的兩個地區）。

首先我們展示一下數據，代碼如下：

import pandas as pd

filename = 'Fremont.csv'
data = pd.read_csv(filename, index_col='Date', parse_dates=True)
res = data.head()
print(res)

數據如下：

爲了直觀，我們畫圖，看看原數據長什麼樣子：

這樣是很難看出特徵的吧，原數據是按照小時來統計兩個區域的數據，下面我們對數據進行重採樣，按周進行計算，我們再看看效果：

# 展示原數據
# data.plot()
data.resample('w').sum().plot()

注意這個resample 的意思，是按照時間索引進行合併後求和計算的。就是在給定的時間單位內重取樣，常見的時間頻率爲：

  • A year
  • M month
  • W week
  • D day
  • H hour
  • T minute
  • S second

說明白了意思，下面我們看看按周進行計算的圖：

現在數據直觀上看起來清晰多了吧，當然我們對於時間序列問題，不止可以採取重取樣查看數據，也可以採用滑動窗口，下面嘗試一下采用365天做一個滑動窗口，這裏是窗口的總和，代碼和圖如下：

# 對數據採用滑動窗口，這裏窗口是365天一個
data.resample('D').sum().rolling(365).plot()

圖如下：

下面我們不按照時間來算了，不對這五年的數據進行分析了，我們想看看某一天的數據分佈，如何做呢？

# 查看某一天的數據分佈
data.groupby(data.index.time).mean().plot()
plt.xticks(rotation=45)

圖如下：

我們可以看到在早上七點到九點，West地區騎走的共享單車比較多，而在下午的四點到六點，East騎走的共享單車比較多。

下面我們將這兩個特徵分別命名爲East和West，再看看前五個小時時間變換圖。

# 查看前五個小時時間變換
# pivot table
data.columns = ['West', 'East']
data['Total'] = data['West'] + data['East']
pivoted = data.pivot_table('Total', index=data.index.time, columns=data.index.date)
res = pivoted.iloc[:5, :5]
print(res)

結果如下：

我們展示一下 24個小時的兩地數據圖：

# 畫圖展示一下
pivoted.plot(legend=False, alpha=0.01)
plt.xticks(rotation=45)

圖如下（數據比較多，畫圖比較慢，等等就OK）：

對數據分析完後，我們開始使用模型訓練，這裏爲了能將數據展示在二維圖上，我們做這樣一個處理，首先我們將數據轉置，畢竟將數據分爲很多特徵，而且只有24個數據不好。所以將樣本數分爲1763，特徵爲24個，則效果就比較好了。。。。

而這樣是可以做的，但是方便演示，我們希望將數據轉換爲二維的，那我們需要對數據進行降維，我們把二十四個特徵降維到二維空間，降維代碼如下：

X = pivoted.fillna(0).T.values
print(X.shape)
X2 = PCA(2).fit_transform(X)
print(X2.shape)
plt.scatter(X2[:, 0], X2[:, 1])

轉換爲二維，數據如圖所示：

數據經過PCA降維後，其物理意義很難解釋，我們可以看到圖中都有負值了，所以不要糾結這個問題。

下面看GMM算法的代碼：

gmm = GaussianMixture(2)
gmm.fit(X)
labels_prob = gmm.predict_proba(X)
print(labels_prob)
labels = gmm.predict(X)
print(labels)
plt.scatter(X2[:, 0], X2[:, 1], c=labels, cmap='rainbow')

我們看看預測出來的label概率和label值（GMM算法比較特別，我們可以知道某個類別的概率值）：

使用GMM算法，分類的圖如下：

下面我們將二維數據還原回去，看效果，代碼如下：

fig, ax = plt.subplots(1, 2, figsize=(14, 6))
pivoted.T[labels == 0].T.plot(legend=False, alpha=0.1, ax=ax[0])
pivoted.T[labels == 1].T.plot(legend=False, alpha=0.1, ax=ax[0])
ax[0].set_title('Purple Cluster')
ax[1].set_title('Red Cluster')

圖如下：

完整代碼及其數據，請移步小編的GitHub

傳送門： 請點擊我

如果點擊有誤：https://github.com/LeBron-Jian/MachineLearningNote

PS：這篇文章是我整理自己學習EM算法的筆記，結合老師的PPT，然後找到網友優秀的博客，將自己的理解寫到這裏。

參考文獻：https://www.cnblogs.com/pinard/p/6912636.html

https://www.cnblogs.com/jerrylead/archive/2011/04/06/2006936.html

https://blog.csdn.net/fuqiuai/article/details/79484421

https://blog.csdn.net/u010834867/article/details/90762296

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

https://blog.csdn.net/lihou1987/article/details/70833229

https://blog.csdn.net/v_july_v/article/details/81708386

相關文章