EM算法解析以及EM應用於GMM

目錄

• 參考blog and 視頻
• EM算法的定義
  • 一、極大似然
    • 1.1 似然函數
    • 1.2 似然函數舉例：已知樣本X，求參數θ
    • 1.3 極大似然即最大可能
  • 二、EM算法的理解
    • 2.1 極大似然估計的復雜情況
    • 2.2 EM算法中的隱變量
    • 2.3 EM算法的例子：拋硬幣
  • 三、EM算法的公式推導
    • 3.1 EM算法的目標函數
    • 3.2 Jensen不等式
    • 3.3 EM算法的流程及證明
    • 3.4 整體推導過程
  • 四、EM算法應用到高斯混合模型中
    • 4.1 高斯混合模型簡單介紹
    • 4.2 GMM的EM算法
    • 4.3 Python實現

參考blog and 視頻

高斯混合模型 EM 算法的 Python 實現
如何通俗理解EM算法
機器學習-白板推導系列(十一)-高斯混合模型GMM

EM算法的定義

---

最大期望算法（Expectation-maximization algorithm，又譯為期望最大化算法），是在概率模型中尋找參數最大似然估計或者最大后驗估計的算法，其中概率模型依賴於無法觀測的隱性變量。

最大期望算法經過兩個步驟交替進行計算，

第一步是計算期望（E），利用對隱藏變量的現有估計值，計算其最大似然估計值；
第二步是最大化（M），最大化在E步上求得的最大似然值來計算參數的值。M步上找到的參數估計值被用於下一個E步計算中，這個過程不斷交替進行。

---

一、極大似然

1.1 似然函數

---

在數理統計學中，似然函數是一種關於統計模型中的參數的函數，表示模型參數中的似然性。“似然性”與“或然性”或“概率”意思相近，都是指某種事件發生的可能性。
而極大似然估計就相當於最大可能的意思
比如我的一位同學和一位獵人一起外出打獵，一只野兔從前方竄過，只聽一聲槍響，野兔應聲倒下，如果要我推測，這一發命中的子彈是誰打的？我就會想，只發一槍便打中，由於獵人命中的概率一般大於我的同學命中的概率，從而推斷出這一槍應該是獵人射中的。
這個例子所做出的推斷就體現了最大似然法的基本思想。

---

多數情況下，我們根據已知條件來推算結果，而最大似然估計是已經知道了結果，然后通過這個已經發生的記過來尋求出現這個結果的可能性最大的條件，以此來作為估計值；即通過 現有結果 推出 條件的一個算法。（類似與后驗概率）
概率是根據條件推測結果，而似然則是根據結果反推條件；在這種意義上，似然函數可以理解為條件概率的逆反。
假定已知某個參數B時，推測事件A發生的概率寫作：

通過貝葉斯公式，可以推導出

現在，我們反過來：事件A已經發生了，通過似然函數L(B|A),估計參數B的可能性。

1.2 似然函數舉例：已知樣本X，求參數θ

  自從Google的圍棋機器人AlphaGo通過4：1戰勝人類世界冠軍李世石之后，人工智能的大潮便一發不可收拾，在無人駕駛、人臉識別、安防監控、醫療影像等各領域大行其道。而專注AI教育的七月在線也已於2017年年底積累了10萬AI學員。
  假定我們需要統計7月在線10萬學員中男生女生的身高分布，怎么統計呢？考慮到10萬的數量巨大，所以不可能一個一個的取統計；所以，通過隨機抽樣，從10萬學院中隨機抽取100個男生，100個女生，然后依次統計他們各自的身高。
  對於這個問題，我們可以通過數學建模抽象整理：

1. 首先我們從10萬學院中抽取到100個男生/女生的身高，組成樣本集X,X={x1, x2, ..., xN}，其中xi表示抽到的第i個人的身高，N=100，表示抽到的樣本個數。
2. 假設男生的身高服從正態分布N1，女生的身高則服從另一個正態分布：N2.
3. 但是這兩個分布的均值和方差都不知道，設為未知參數θ=[u, ∂]T
4. 現在需要使用極大似然(MLE)，通過着100個男生/女生的身高結果，即通過樣本集X來估計兩個正態分布的位置參數θ；簡而言之，通過已知結果來推測θ，即求p(θ|x)中概率最大的θ是啥。
     因為這些男生的升高服從同一個高斯分布，那么抽到男生A的概率是p(xA|θ)，抽到男生B的概率是p(xB|θ)，所以同時抽到男生A和男生B的概率是p(xA|θ)*p(xB|θ)。
     同理，從分布是p(x|θ)的總體樣本中同時抽到這100個男生樣本的概率，也就是樣本集X中100個樣本的聯合概率(即他們各自概率的乘積)，可用下面的式進行表示：
   

  插一句，有的文章會用p(x|θ)，有的文章則會用p(x;θ)。不過，不管使用哪種表示方法，本質都是一樣的。如果涉及到Bayes公式的話，用前者表示p(x|θ)會好一些。
  在七月在線的那么多男學員中，我一抽就抽到了着100個男生，而不是其他人，那么說明在整個學校中這100個人的身高出現的概率是最大的，這個概率就是上面這個似然函數L(θ)。那么現在我想要確定的θ，就是讓L(θ)這個函數值最大的時候θ的取值。

1.3 極大似然即最大可能

  假定我們找到一個參數
，能夠讓似然函數L(θ)最大(也就是說確定一個θ使得抽到這100個男生的身高概率最大)，則
應該是"最可能"的參數值，相當於θ的極大似然估計量，記作：

這里的L(θ)是連乘的，為了便於分析，我們可以定義對數似然函數，將其變成連加的：

  現在需要使似然函數L(θ)極大化，相當於極大化H(θ),最后極大值對應的θ就是我們估計的參數。
  對於求一個函數的極值，通過我們在本科所學習的微積分知識，最直接的設想就是求導，然后讓導數為0，最后求解出這個方程的0。但，如果0是包含多個參數的向量應該怎么處理呢？當然是求L(0)對所有參數的偏導數，也就是梯度了，從而n個未知的參數，就有n個方程，方程組的解就是似然函數的極值點了，最終得到這個n個參數的值。
求極大似然函數估計值的一般步驟

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

二、EM算法的理解

2.1 極大似然估計的復雜情況

我們已經知道，極大似然估計用一句話來概括就是：知道結果，反推條件0。
例如，在上文中，為了統計七月在線的10萬學院中男生女生的身高分布

---

• 首先我們從10萬學院中抽取到100個男生/女生的身高，組成樣本集X,X={x1,x2,...,xN}，其中xi表示抽到的第i個人的身高，N=100，表示抽到的樣本個數。
• 假定男生的身高服從正態分布N1,女生的身高服從正態分布N2；但是這兩個分布的均值和方差都不知道，設為未知參數θ=[u, ∂]T
• 現在需要使用極大似然估計法(MLE)，通過這100個男生/女生的身高結果來估計這兩個正態分布的位置參數θ，問題定義為已知X，求θ；換而言之就是求p(θ|x)最大的θ的數值。

---

極大似然估計的目標是求解實現結果的最佳θ，但極大似然估計需要面臨的概率分布只有一個或者知道結果是通過哪個概率分布實現的，只不過我們不知道這個概率分布的參數。
但現在我們讓情況復雜一些，比如這100個男生和100個女生混在一起了。我們擁有着200個人的身高數據，卻不知道着200個人每一個是男生還是女生，此時的男女性別就像一個隱變量。
這時候情況就有點兒尷尬，因為通常來說，我們只有知道了精確的男女身高的正態分布參數，我們才能知道每一個人更有可能是男生還是女生。反過來，我們只有知道了每個人是男生還是女生才能盡可能准確的估計男女各自身高的正態分布的參數。
而EM算法就是為了解決"極大似然估計"這種更復雜而存在的。

2.2 EM算法中的隱變量

理解EM算法中的隱變量很關鍵，就如理解KMP算法中的理解Next數組的意義一樣。
一般用Y表示觀測的隨機變量的數據，Z表示隱隨機變量的數據（因為我們觀測不到結果是從哪個概率分布中得出的，所以將這個叫做隱變量）。於是Y和Z連在一起被稱為完全數據，僅Y一個稱為不完全數據。
這時有沒有發現EM算法面臨的問題主要就是：有一個隱變量的數據Z；而如果Z已知的話，那問題就可用極大似然估計求解了。於是乎，怎么把Z變成已知的呢？

2.3 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還是(這就是我們不知道的隱變量，label：A,B)，只觀測5輪循環，每輪循環10此，一共50次拋硬幣的結果，這時候就沒有辦法直接估計A和B正面的概率；此時，就輪到EM算法出場了(如圖b)

可能咋一看，沒看懂；沒關系，我們來通俗化這個例子。
還是兩枚硬幣A和B，假定隨機拋擲后正面朝上分別為PA,PA.為了估計這2個硬幣朝上的概率，咱們輪流拋硬幣A和B，每一輪都連續拋5次，總共5輪：

硬幣	結果	統計
A	正正反正反	3正-2反
B	反反正正反	2正-3反
A	正反反反反	1正-4反
B	正反反正正	3正-2反
A	反正正反反	2正-3反
硬幣A被拋了15次，在第一輪、第三輪和第五輪分別出現了3次正、1次正、2次正，所以很容易估計出PA，類似的PB也很容易計算出來，如下：
PA = (3 + 1 + 2) / 15 = 0.4
PB = (2 + 3) / 10 - 0.5
問題來了，如果我們不知道拋的硬幣是A還是B(即硬幣種類是隱變量)，然后再輪流五輪，得到如下結果
硬幣	結果	統計
----	----	----
Unknown	正正反正反	3正-2反
Unknown	反反正正反	2正-3反
Unknown	正反反反反	1正-4反
Unknown	正反反正正	3正-2反
Unknown	反正正反反	2正-3反
OK，問題現在變得更加有意思了；現在我們的目標沒有變化，依然是取估計PA和PB，需要如何取做呢？

---

顯然，此時我們多了一個硬幣種類的隱變量，設為z，可以把它認為是一個5維的向量
(z1, z2, z3, z4, z5)，代表，每次投擲時所使用的硬幣，比如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賦初值，比如：
硬幣A正面朝上的概率PA = 0.2
硬幣B正面朝上的概率PB = 0.7

然后，我們看看第一輪拋擲最可能是哪個硬幣。
如果是硬幣A，得出3正2反的概率為 0.20.20.20.80.8 = 0.00512
如果是硬幣B，得出三正2反的概率為 0.70.70.70.30.3 = 0.03087
然后一次求出其他四輪中的相應概率。做成表格如下（標粗表示其概率更大）：

輪數	若是硬幣A	若是硬幣B
1	0.00512	0.03087
2	0.02048	0.01323
3	0.08192	0.00567
4	0.00512	0.03087
5	0.02048	0.01323
按照最大似然法則：
第1輪中最有可能的是硬幣B
第2輪中最有可能的是硬幣A
第3輪中最有可能的是硬幣A
第4輪中最有可能的是硬幣B
第5輪中最有可能的是硬幣A
我們就把概率更大，即更可能是A的，就第2，3，5輪出現正的次數2，1，2相加除以A被拋的總次數15次（A拋了三輪，每輪5次），B同理 。
PA = (2 + 1 + 2) / 15 = 0.33
PB = (3 + 3) / 10 = 0.6

不妨設真實的PA和PB分別是0.4和0.5。那么對比下面我們初始化的PA和PB以及新估計出來的PA和PB。

初始化的PA	估計的PA	真實的PA	初始化的PB	估計的PB	真實的PB
0.2	0.33	0.4	0.7	0.6	0.5

我們估計的PA和PB相比於它們的初是值，更接近它們的真實值了！就這樣，不斷迭代 不斷接近真實值，這就是EM算法的奇妙支出。
可以期待，我們繼續按照上面的思路，用估計出的PA和PB再來估計z，再用z來估計新的PA和PB，反復迭代下去，就可以最終得到PA = 0.4, PB = 0.5, 此時無論怎樣迭代，PA和PB的數值都會保持0.4和0.5不變，於是乎，我們就找到了PA和PB的最大似然估計。

三、EM算法的公式推導

3.1 EM算法的目標函數

還記得1.2節開頭說的把？

---

從分布是p(x|θ)的總體樣本中抽取到着100個樣本的概率，也就是樣本集X中各個樣本的聯合概率，用下式表示：

---

  假設我們有一個樣本集{x(1),...,x(m)}，包含m個獨立的樣本，現在我們想找到每個樣本隱含的類別z，能使得p(x,z)最大。p(x,z)的極大似然估計如下：

  第一步是對極大似然取對數，第二步是對每個樣例的每個可能類別z求聯合分布概率和。但是直接求一般比較困難，因為有隱變量z存在，但是一般確定了z后，求解就容易了。
  對於參數估計，我們本質上還是想獲得一個使似然函數最大化的那個參數θ，現在與極大似然不同的知識似然函數式中多了一個位置的變量z，見下式(1)。也就是說我們的目標是尋找到合適的θ和z，從而容L(θ)最大。那我們也許會想，你就是多了一個未知變量而已，我們也可以分別對位置的θ和z分別求偏導，再令其等於0，求解出來不也是一樣的嗎？

  本質上我們需要最大化(1)式，也就是似然函數，但是可以看到里面有"和的對數"，求導后形式會非常復雜(可以想象一下log(f1(x)+f2(x)+...)復合函數的求導)，所以很難求解得到未知參數z和θ。
  我們把分子和分母同時乘上一個相等的函數(即隱變量Z的概率分布Qi(Z(i)))，其概率之和等於1，即

，即得到上圖中的(2)式，但(2)式還是有"和的對數"，還是求解不了，咋辦呢？接下來，就是見證奇跡的時刻了，我們通過Jensen不等式可得到(3)式，此時(3)式變成了"對數的和"，如此一來求導就容易了。
  從(2)式到(3)式，神奇的地方有兩點：

1. 當我們在求(2)式的極大值時，(2)式不太好計算，我們便想(2)式大於某個方便計算的下界(3)式，而我們盡可能讓下界(3)式最大化，直到(3)式的計算結果等於(2)式，便也就間接求得了(2)式的極大值；
2. Jensen不等式，促進神奇發生的Jensen不等式到底是什么來歷呢？

3.2 Jensen不等式

f是定義域為實數的函數

1. 如果對於所有的實數x,f(x)的二次導數f''(x)>=0,那么f是凸函數
2. 當x是向量時，如果其hessian矩陣H是半定的(H>=0),那么f是凸函數
3. 如果f''(x)>0或者H>0，那么f是嚴格凸函數。

Jensen不等式表述如下：

---

如果f是凸函數，X是隨機變量，那么：E[f(X)]>=f(E[X])，通俗的說法是函數的期望大於期望的函數。
特別地，如果f是嚴格凸函數，當且僅當P(X=EX)=1，即X是常量時，上式取等號，即E(f(X))=f(EX)。

---

如下圖所示：

  圖中，實線f是凸函數(因為函數上方的區域是一個凸集)，X是隨機變量，有0.5的概率是a，有0.5的概率是b(就像拋硬幣一樣)。X的期望值就是a和b的中值，可以明顯的看出，E[f(X)] >= f(E[X])。
當然，當f是(嚴格)凹函數時，不等式方向，也就是E[f(X)]<=f(E(X))
回到公式(2)，因為f(x)=log(x)為凹函數(其二次導數為1/x2<0)

(2)式中的

就是

的期望，為什么，可以回想期望公式中的Lazy Statistician規則，如下：

---

設Y是隨機變量X的函數，Y=g(X)(g是連續函數)，那么

1. X是離散型隨機變量，它的分布律為P(X=xk)=pk,k=1,2,...。若
   

則有：

2. X是連續型隨機變量，它的概率密度為f(x)，若
   

絕對收斂，則有：

---

考慮到E(X)=∑xp(x),F(x)是X的函數，則E[f(X)] = ∑xp(x)，又

所以就可以得到公式(3)的不等式了：

  到這里，現在式(3)就比較容易求導了，但是式(2)和式(3)是不等號，而不是等號，導致式(2)的最大值不是式(3)的最大值，而我們想得到式(2)的最大值，應該怎么辦呢？
  上面的式(2)和式(3)不等式可以寫成：似然函數L(θ)>=J(z,Q)，如3.1節最后所說，我們可以通過不斷地最大化這個下屆J，來使地L(θ)不斷提高，最后達到L(θ)的最大值。

  見上圖，我們固定θ，調整Q(z)使下屆J(z,Q)達到最大值(θt到θt+1)，然后再固定θ，調整Q(z).....直到收斂到似然函數L(θ)的最大值處的θ*。
  這里又兩個問題：

1. 什么時候下屆J(z,Q)與L(θ)在此點θ處相等？
2. 為什么一定會收斂？
   首先第一個問題，在Jensen不等式中說到，當自變量X是常數的時候，等式成立。換言之，為了讓(3)式取等號，我們需要：
   

其中c為常熟，不依賴與z(i)。對該式做個變換，並對所有的z求和，得到

因為前面提到的，隱變量Z的概率分布，其概率之和等於1，所以可以推導出：

通過

可以求出來Q(zi)的值，即：

又因為：

所以Q(Zi)為：

  瞬間豁然開朗！至此，我們推出了在固定參數θ后，使用下界拉升的Q(z)的計算公式就是條件概率，解決了Q(z)如何選擇的問題。這一步就是E步，建立L(θ)的下界。
  接下來的M步，就是在給定Q(z)后，調整θ，去極大化L(θ)的下界J(z,Q)，畢竟在固定Q(z)后，下屆還可以調整的更大。
  這就是EM算法的步驟。

3.3 EM算法的流程及證明

  我們來總結一下，期望最大化EM算法是一種從不完全數據或又數據丟失的數據集(存在隱含變量)中求解概率模型參數的最大似然估計方法。
  換言之，當我們不知道隱變量z的具體取值時(比如是硬幣A還是硬幣B)，也就不好判斷硬幣A或硬幣B正面朝上的概率θ，既然這樣，那就：

---

1. 隨機初始化分布參數θ
2. 然后循環重復直到收斂{

   1. (E步，求Q函數)對於每一個i，計算根據上一次迭代的模型參數來計算隱性變量的后驗概率(其實就是隱性變量的期望)，來作為隱藏變量的現估計值：
      

   2. (M步，求使用Q函數獲得極大時的參數取值)將似然函數最大化以獲得新的參數值
      

}

---

3.4 整體推導過程

四、EM算法應用到高斯混合模型中

4.1 高斯混合模型簡單介紹

首先簡單介紹一些，高斯混合模型(GMM, Gaussian Mixture Model)是多個高斯模型的線性疊加，高斯混合模型的概率分布可以表示如下：

其中，K 表示模型的個數，αk 是第 k 個模型的系數，表示出現該模型的概率，ϕ(x;μk,Σk) 是第 k 個高斯模型的概率分布。

這里討論的是多個隨機變量的情況，即多元高斯分布，因此高斯分布中的參數不再是方差 σk，而是協方差矩陣 Σk 。

我們的目標是給定一堆沒有標簽的樣本和模型個數K，以此求得混合模型的參數，然后就可以用這個模型來對樣本進行聚類。

4.2 GMM的EM算法

我們直到，EM算法是通過不斷迭代來求得最佳參數的。在執行該算法之前，需要先給出一個初始化的模型參數。

我們讓每個模型的 μ 為隨機值，Σ 為單位矩陣，α 為 1K，即每個模型初始時都是等概率出現的。

EM算法可以分為E步和M步。

E步
直觀理解就是我們已經知道了樣本xi，那么它是由哪個模型產生的呢？我們這里求的就是：樣本xi來自於第k個模型的概率，我們把這個概率稱為模型k對樣本xi的"責任"，也叫"響應度"，記作γik，計算公式如下：

M步
根據樣本和當前 γ 矩陣重新估計參數，注意這里 x 為列向量：

4.3 Python實現

在給出代碼前，先作一些說明。

• 在對樣本應用高斯混合模型的EM算法前，需要先進行數據預處理，即把所有樣值都縮放到0和1之間。
• 初始化模型參數時，要確保任意兩個模型之間參數沒有完全相同，否則迭代到最后，兩個模型的參數也將完全相同，相當於一個模型。
• 模型的個數必須大於1.當K等於1時相當於將樣本聚成一類，沒有任何意義。

generate_data.py

import numpy as np
import matplotlib.pyplot as plt

cov1 = np.mat("0.3 0;0 0.1")
cov2 = np.mat("0.2 0;0 0.3")
mu1 = np.array([0, 1])
mu2 = np.array([2, 1])

sample = np.zeros((100, 2))
sample[:30, :] = np.random.multivariate_normal(mean=mu1, cov=cov1, size=30)
sample[30:, :] = np.random.multivariate_normal(mean=mu2, cov=cov2, size=70)
np.savetxt("./data/sample.data", sample)

plt.plot(sample[:30, 0], sample[:30, 1], "bo")
plt.plot(sample[30:, 0], sample[30:, 1], "rs")
plt.show()

代碼解析

---

1. np.mat()：生成矩陣
   np.array()：創建數組
   mat可以從字符串或列表中生成；array只能從列表中生成

2. 返回來一個給定形狀和類型的用0填充的數組；

3. np.random.multivariate_normal方法用於根據實際情況生成一個多元正態分布矩陣
   其中mean和cov為必要的傳參而size，check_valid以及tol為可選參數。
   mean：mean是多維分布的均值維度為1；
   cov：協方差矩陣（協方差基本概念戳這里），注意：協方差矩陣必須是對稱的且需為半正定矩陣；
   size：指定生成的正態分布矩陣的維度（例：若size=(1, 1, 2)，則輸出的矩陣的shape即形狀為 1X1X2XN（N為mean的長度））。
   check_valid：這個參數用於決定當cov即協方差矩陣不是半正定矩陣時程序的處理方式，它一共有三個值：warn，raise以及ignore。當使用warn作為傳入的參數時，如果cov不是半正定的程序會輸出警告但仍舊會得到結果；當使用raise作為傳入的參數時，如果cov不是半正定的程序會報錯且不會計算出結果；當使用ignore時忽略這個問題即無論cov是否為半正定的都會計算出結果。3種情況的console打印結果如下：

4. generate_data.py：制定了2個維度為2的正態分布；前面30條數據是N1,后邊70條數據是N2。

---

gmm.py

# -*- coding: utf-8 -*-
# ----------------------------------------------------
# Copyright (c) 2017, Wray Zheng. All Rights Reserved.
# Distributed under the BSD License.
# ----------------------------------------------------

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

DEBUG = True

######################################################
# 調試輸出函數
# 由全局變量 DEBUG 控制輸出
######################################################
def debug(*args, **kwargs):
    global DEBUG
    if DEBUG:
        print(*args, **kwargs)


'''
# 第 k 個模型的高斯分布密度函數
# 每 i 行表示第 i 個樣本在各模型中的出現概率
# 返回一維列表；代表着當前這個高斯分布下，所有樣本的概率是多少
'''
def phi(Y, mu_k, cov_k):
    norm = multivariate_normal(mean=mu_k, cov=cov_k)
    return norm.pdf(Y)


'''
# E 步：計算每個模型對樣本的響應度
# Y 為樣本矩陣，每個樣本一行，只有一個特征時為列向量（100，2）
# mu 為均值多維數組，每行表示一個樣本各個特征的均值；2個期望矩陣，維度為(2, 2)
# cov 為協方差矩陣的數組，alpha 為模型響應度數組；含有K個數值的數組，表示每個模型的概率，初始化為1/K (1, K)
'''
def getExpectation(Y, mu, cov, alpha):
    # 樣本數
    N = Y.shape[0]
    # 模型數
    K = alpha.shape[0]

    # 為避免使用單個高斯模型或樣本，導致返回結果的類型不一致
    # 因此要求樣本數和模型個數必須大於1
    assert N > 1, "There must be more than one sample!"
    assert K > 1, "There must be more than one gaussian model!"

    # 響應度矩陣，行對應樣本，列對應響應度
    gamma = np.mat(np.zeros((N, K)))

    # 計算各模型中所有樣本出現的概率，行對應樣本，列對應模型
    prob = np.zeros((N, K))
    for k in range(K):     #對同一列，的每一行進行操作；給同一個K的不同樣本進行賦值，表示再當前這個高斯分布下，這些樣本的概率值
        prob[:, k] = phi(Y, mu[k], cov[k])
    prob = np.mat(prob)   #(100, 2)

    # 計算每個模型對每個樣本的響應度
    for k in range(K):
        gamma[:, k] = alpha[k] * prob[:, k]   # 對於每個樣本而言，每個模型出現的概率是alpha(k)
    for i in range(N):
        gamma[i, :] /= np.sum(gamma[i, :])    #歸一化
    return gamma


'''
# M 步：迭代模型參數
# Y 為樣本矩陣(樣本數,特征數)(100,2)
# gamma 為響應度矩陣(樣本數,模型數)(100,2)
'''
def maximize(Y, gamma):
    # 樣本數和特征數
    N, D = Y.shape
    # 模型數
    K = gamma.shape[1]

    #初始化參數值
    mu = np.zeros((K, D))
    cov = []
    alpha = np.zeros(K)

    # 更新每個模型的參數
    for k in range(K):
        # 第 k 個模型對所有樣本的響應度之和
        Nk = np.sum(gamma[:, k])
        # 更新 mu
        # 對每個特征求均值
        for d in range(D):
            mu[k, d] = np.sum(np.multiply(gamma[:, k], Y[:, d])) / Nk
        # 更新 cov
        cov_k = np.mat(np.zeros((D, D)))
        for i in range(N):
            cov_k += gamma[i, k] * (Y[i] - mu[k]).T * (Y[i] - mu[k]) / Nk
        cov.append(cov_k)
        # 更新 alpha
        alpha[k] = Nk / N
    cov = np.array(cov)
    return mu, cov, alpha
'''
# 傳入一個矩陣
# 數據預處理
# 將所有數據都縮放到 0 和 1 之間
'''
def scale_data(Y):
    # 對每一維特征分別進行縮放
    for i in range(Y.shape[1]):
        max_ = Y[:, i].max()
        print(i)
        min_ = Y[:, i].min()
        Y[:, i] = (Y[:, i] - min_) / (max_ - min_)
    debug("Data scaled.")
    return Y

'''
# 初始化模型參數
# shape 是表示樣本規模的二元組，(樣本數, 特征數)
# K 表示模型個數
'''
def init_params(shape, K):
    N, D = shape   # N:樣本數，D：特征數
    mu = np.random.rand(K, D)   # (模型數，特征數)，生成K個期望矩陣，維度為(2, 2)
    cov = np.array([np.eye(D)] * K)    # 生成K個協方差矩陣，即2個單位矩陣
    alpha = np.array([1.0 / K] * K)    # 生成含有K個數值的數組，表示每個模型的概率，初始化為1/K
    debug("Parameters initialized.")
    debug("mu:", mu, "cov:", cov, "alpha:", alpha, sep="\n")
    return mu, cov, alpha


'''
# 高斯混合模型 EM 算法
# 給定樣本矩陣 Y，計算模型參數（100，2）
# K 為模型個數（可能的模型數量，在這里表示高斯模型數量
# times 為迭代次數
'''
def GMM_EM(Y, K, times):
    Y = scale_data(Y)     # 將Y中的數據縮放到 0-1 之間
    mu, cov, alpha = init_params(Y.shape, K)   #Y.shape表示數據的維度，K表示模型的數量
    for i in range(times):
        gamma = getExpectation(Y, mu, cov, alpha)   #得到一個(100,2),gamma(i,j)：第i個樣本是第j個模型的概率
        mu, cov, alpha = maximize(Y, gamma)
    debug("{sep} Result {sep}".format(sep="-" * 20))
    debug("mu:", mu, "cov:", cov, "alpha:", alpha, sep="\n")
    return mu, cov, alpha

main.py

# -*- coding: utf-8 -*-
# ----------------------------------------------------
# Copyright (c) 2017, Wray Zheng. All Rights Reserved.
# Distributed under the BSD License.
# ----------------------------------------------------

import matplotlib.pyplot as plt
from gmm import *

# 設置調試模式
DEBUG = True

# 載入數據
Y = np.loadtxt("./data/sample.data")
matY = np.matrix(Y, copy=True)

# 模型個數，即聚類的類別個數
K = 2

# 計算 GMM 模型參數
'''
1. 這里傳入的三個參數是：matY：所有觀測的數據(x1,x2)組成的矩陣，維度為：(100,2)
2. K：聚類的數量，表示這里由兩種可能的高斯模型
3. 100表示數據的數量，一共有100條數據
4. GMM_EM：這個函數需要返回三個參數：mu,cov,alpha
    1) mu：表示期望
    2) cov：表示協方差矩陣
    3)alpha：表示每個樣本屬於哪個高斯分布的概率 alpha = (alpha1, ..., alphak)表示當前樣本屬於第i個高斯分布的概率值.
'''
mu, cov, alpha = GMM_EM(matY, K, 100)

# 根據 GMM 模型，對樣本數據進行聚類，一個模型對應一個類別
N = Y.shape[0]
# 求當前模型參數下，各模型對樣本的響應度矩陣
gamma = getExpectation(matY, mu, cov, alpha)
# 對每個樣本，求響應度最大的模型下標，作為其類別標識
category = gamma.argmax(axis=1).flatten().tolist()[0]
# 將每個樣本放入對應類別的列表中
class1 = np.array([Y[i] for i in range(N) if category[i] == 0])
class2 = np.array([Y[i] for i in range(N) if category[i] == 1])

# 繪制聚類結果
plt.plot(class1[:, 0], class1[:, 1], 'rs', label="class1")
plt.plot(class2[:, 0], class2[:, 1], 'bo', label="class2")
plt.legend(loc="best")
plt.title("GMM Clustering By EM Algorithm")
plt.show()

---

代碼解析

1. argmax()函數

import numpy as np
a = np.array([3, 1, 2, 4, 6, 1])
#print(np.argmax(a))
#4
#argmax返回的是最大數的索引.argmax有一個參數axis,默認是0,表示第幾維的最大值.

2. flatten()函數

#a.flatten()：a是個數組，a.flatten()就是把a降到一維，默認是按行的方向降
from numpy import *
a=array([[1,2],[3,4],[5,6]])
print(a)
'''
array([[1, 2],
    [3, 4],
    [5, 6]])'''
a.flatten() #默認按行的方向降維
#array([1, 2, 3, 4, 5, 6])
a.flatten('F') #按列降維
#array([1, 3, 5, 2, 4, 6])
a.flatten('A') #按行降維
#array([1, 2, 3, 4, 5, 6])

3. tolist()函數

#tolist()作用：將矩陣（matrix）和數組（array）轉化為列表。
from numpy import *
a1 = [[1,2,3],[4,5,6]] #列表
a2 = array(a1) #數組
print(a2)
# array([[1, 2, 3],
#        [4, 5, 6]])
a3 = mat(a1) #矩陣
print(a3)
# matrix([[1, 2, 3],
#         [4, 5, 6]])
a4 = a2.tolist()
print(a4)
# [[1, 2, 3], [4, 5, 6]]
a5 = a3.tolist()
print(a5)
# [[1, 2, 3], [4, 5, 6]]
print(a4 == a5)
# True

---

聚類結果：