python 数学期望_python机器学习笔记：EM算法
完整代碼及其數據，請移步小編的GitHub
傳送門：請點擊我
如果點擊有誤：https://github.com/LeBron-Jian/MachineLearningNote
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(xA | θ)，抽到男生B的概率是 p(xB | θ)，那因為他們是獨立的，所以很明顯，我們同時抽到男生A和男生B的概率是 p(xA | θ) * p(xB | θ)，同理，我同時抽到這 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算法。博客地址：
這里再簡單的說一下，在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')
圖如下：
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
總結
以上是生活随笔為你收集整理的python 数学期望_python机器学习笔记：EM算法的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： Introduction to ML
- 下一篇： 2020级C语言大作业 - 三国杀