在知乎上看到一篇有價值的文章，分享給學習機器學習、深度學習、智能科學的朋友：
高斯混合模型（Gaussian Mixed Model）指的是多個高斯分布函數的線性組合，理論上GMM可以擬合出任意類型的分布，通常用於解決同一集合下的數據包含多個不同的分布的情況（或者是同一類分布但參數不一樣，或者是不同類型的分布，比如正態分布和伯努利分布）。
如圖1，圖中的點在我們看來明顯分成兩個聚類。這兩個聚類中的點分別通過兩個不同的正態分布隨機生成而來。但是如果沒有GMM，那麼只能用一個的二維高斯分布來描述圖1中的數據。圖1中的橢圓即為二倍標準差的正態分布橢圓。這顯然不太合理，畢竟肉眼一看就覺得應該把它們分成兩類。
圖 1這時候就可以使用GMM了！如圖2，數據在平面上的空間分布和圖1一樣，這時使用兩個二維高斯分布來描述圖2中的數據，分別記為 和 . 圖中的兩個橢圓分別是這兩個高斯分布的二倍標準差橢圓。可以看到使用兩個二維高斯分布來描述圖中的數據顯然更合理。實際上圖中的兩個聚類的中的點是通過兩個不同的正態分布隨機生成而來。如果將兩個二維高斯分布 和 合成一個二維的分布，那麼就可以用合成後的分布來描述圖2中的所有點。最直觀的方法就是對這兩個二維高斯分布做線性組合，用線性組合後的分布來描述整個集合中的數據。這就是高斯混合模型（GMM）。
圖 22. 高斯混合模型
設有隨機變量 ，則混合高斯模型可以用下式表示：
其中 稱為混合模型中的第 個分量（component）。如前面圖2中的例子，有兩個聚類，可以用兩個二維高斯分布來表示，那麼分量數 . 是混合係數（mixture coefficient），且滿足：
可以看到 相當於每個分量 的權重。
例如，圖3中將兩個高斯混合分布進行線性組合，
圖 3GMM常用於聚類。如果要從 GMM 的分布中隨機地取一個點的話，實際上可以分為兩步：首先隨機地在這個 Component 之中選一個，每個 Component 被選中的概率實際上就是它的係數 ，選中 Component 之後，再單獨地考慮從這個 Component 的分布中選取一個點就可以了──這裡已經回到了普通的 Gaussian 分布，轉化為已知的問題。
將GMM用於聚類時，假設數據服從混合高斯分布（Mixture Gaussian Distribution），那麼只要根據數據推出 GMM 的概率分布來就可以了；然後 GMM 的 個 Component 實際上對應 個 cluster 。根據數據來推算概率密度通常被稱作 density estimation 。特別地，當我已知（或假定）概率密度函數的形式，而要估計其中的參數的過程被稱作參數估計。
例如圖2的例子，很明顯有兩個聚類，可以定義 . 那麼對應的GMM形式如下：
上式中未知的參數有六個： . 之前提到GMM聚類時分為兩步，第一步是隨機地在這 個分量中選一個，每個分量被選中的概率即為混合係數 . 可以設定 ，表示每個分量被選中的概率是0.5，即從中抽出一個點，這個點屬於第一類的概率和第二類的概率各佔一半。但實際應用中事先指定 的值是很笨的做法，當問題一般化後，會出現一個問題：當從圖2中的集合隨機選取一個點，怎麼知道這個點是來自 還是 呢？換言之怎麼根據數據自動確定 和 的值？這就是GMM參數估計的問題。要解決這個問題，可以使用EM算法。通過EM算法，我們可以迭代計算出GMM中的參數：
3. EM算法（Expectation Maximization Algorithm）
EM算法是一種迭代優化策略，由於它的計算方法中每一次迭代都分兩步，其中一個為期望步（E步），另一個為極大步（M步），所以算法被稱為EM算法（Expectation Maximization Algorithm）。EM算法受到缺失思想影響，最初是為了解決數據缺失情況下的參數估計問題，其算法基礎和收斂有效性等問題在Dempster，Laird和Rubin三人於1977年所做的文章Maximum likelihood from incomplete data via the EM algorithm中給出了詳細的闡述。其基本思想是：首先根據己經給出的觀測數據，估計出模型參數的值；然後再依據上一步估計出的參數值估計缺失數據的值，再根據估計出的缺失數據加上之前己經觀測到的數據重新再對參數值進行估計，然後反覆迭代，直至最後收斂，迭代結束。
令 是觀察到的樣本數據集合， 為丟失的數據集合，完整的樣本集合 。有
由於 未知，在給定參數 時，似然函數可以看作 的函數，似然函數
由於 未知，因此我們需要尋找到一個在 的所有可能情況下，平均意義下的似然函數最大值，即似然函數對 的期望的最大值。EM算法分為E步和M步：
EM算法偽代碼如下：
EM算法性質：EM算具有收斂性；EM算法只能保證收斂於似然函數的局部最大值點（極值點），而不能保證收斂於全局最優點。
4. GMM的貝葉斯理解
在介紹GMM參數估計之前，我們先改寫GMM的形式，改寫之後的GMM模型可以方便地使用EM估計參數。GMM的原始形式如下：
前面提到的 可以看成是第 類被選中的概率。我們引入一個新的 維隨機變量 . 只能取0或1兩個值； 表示第 類被選中的概率，即 ；如果 表示第 類沒有被選中的概率。更數學化一點， 要滿足以下兩個條件：
例如圖2中的例子，有兩類，則 的維數是2. 如果從第一類中取出一個點，則 ，如果從第二類中取出一個點，則 的概率就是 ，假設 之間是獨立同分布的，我們可以寫出 的聯合概率分布形式：
因為 只能取0或1，且 中只能有一個 為1而其它 全為0，所以上式是成立的。
圖2中的數據可以分為兩類，顯然，每一類中的數據都是服從正態分布的。這個敘述可以用條件概率來表示：
即第 類中的數據服從正態分布。進而上式有可以寫成如下形式：
上面分別給出了 和 的形式，根據條件概率公式，可以求出 的形式：
可以看到GMM模型的(1)式與(10)式有一樣的形式，且(10)式中引入了一個新的變量 ，通常稱為隱含變量（latent variable）。對於圖2中的數據，『隱含』的意義是：我們知道數據可以分成兩類，但是隨機抽取一個數據點，我們不知道這個數據點屬於第一類還是第二類，它的歸屬我們觀察不到，因此引入一個隱含變量 來描述這個現象。
注意到在貝葉斯的思想下， 是先驗概率， 是似然概率，很自然我們會想到求出後驗概率 ：
上式中我們定義符號 來表示第 個分量的後驗概率。在貝葉斯的觀點下， 可視為 的先驗概率。
上述內容改寫了GMM的形式，並引入了隱含變量 和已知 後的的後驗概率 ，這樣做是為了方便使用EM算法來估計GMM的參數。
5. EM算法估計GMM參數
EM算法（Expectation-Maximization algorithm）分兩步，第一步先求出要估計參數的粗略值，第二步使用第一步的值最大化似然函數。因此要先求出GMM的似然函數。
假設 ，對於圖2， 是圖中所有點（每個點有在二維平面上有兩個坐標，是二維向量，因此 等都用粗體表示）。GMM的概率模型如(1)式所示。GMM模型中有三個參數需要估計，分別是 ， 和 . 將(1)式稍微改寫一下：
為了估計這三個參數，需要分別求解出這三個參數的最大似然函數。先求解 的最大似然函數。對(12)式取對數後再對 求導並令導數為0即得到最大似然函數。
注意到上式中分數的一項的形式正好是(11)式後驗概率的形式。兩邊同乘 ，重新整理可以得到：
其中：
(12)式和(15)式中， 表示點的數量。 表示點 屬於聚類 的後驗概率。則 可以表示屬於第 個聚類的點的數量。那麼 表示所有點的加權平均，每個點的權值是
，跟第 個聚類有關。
同理求 的最大似然函數，可以得到：
最後剩下 的最大似然函數。注意到 有限制條件 ，因此我們需要加入拉格朗日算子： 。
求上式關於 的最大似然函數，得到：
上式兩邊同乘 ，可以得到 ，進而可以得到 更簡潔的表達式：
EM算法估計GMM參數即最大化(14)，(16)和(18)。需要用到(11)，(14)，(16)和(18)四個公式。我們先指定 ， 和 的初始值，帶入(14)中計算出 ，然後再將 帶入(14)，(16)和(18)，求得 ， 和 ；接著用求得的 ， 和 再帶入(11)得到新的 ，再將更新後的 帶入(14)，(16)和(18)，如此往復，直到算法收斂。
6. EM算法（偽代碼）
定義分量數目 ，對每個分量 設置 ， 和 的初始值，然後計算(1)式 的對數似然函數。
E step
根據當前的 、 、 計算後驗概率 .
M step
根據E step中計算的 再計算新的 、 、 .
其中：
計算(1)式 的對數似然函數
檢查參數是否收斂或對數似然函數是否收斂，若不收斂，則返回第2步。
來源：https://zhuanlan.zhihu.com/p/45793456
下面是補充的EM的python實現代碼：
from numpy import *
import numpy as np
import matplotlib.pyplot as plt
import random
def create_sample_data(m, n):
mat_y = mat(zeros((m, n)))
for i in range(m):
for j in range(n):
#通過產生隨機數，每一行表示一次實驗結果
mat_y[i, j] = random.randint(0, 1)
return mat_t
#EM算法
def em(arr_y, theta, tol, iterator_num):
PI = 0
P = 0
Q = 0
m,n = shape(arr_y)
mat_y = arr_y.getA()
for i in range(iterator_num):
miu = []
PI = copy(theta[0])
P = copy(theta[1])
Q = copy(theta[2])
for j in range(m):
miu_value = (PI * (P**mat_y[j]) * ((1 - P)**(1 - mat_y[j]))) / \
(PI * (P**mat_y[j]) * ((1 - P)**(1 - mat_y[j])) + (1 - PI) * (Q**mat_y[j]) * ((1 - Q)**(1 - mat_y[j])))
miu.append(miu_value)
sum1 = 0.0
for j in range(m):
sum1 += miu[j]
theta[0] = sum1 / m
sum1 = 0.0
sum2 = 0.0
for j in range(m):
sum1 += miu[j] * mat_y[j]
sum2 += miu[j]
theta[1] = sum1 / sum2
sum1 = 0.0
sum2 = 0.0
for j in range(m):
sum1 += (1 - miu[j]) * mat_y[j]
sum2 += (1 - miu[j])
theta[2] = sum1 / sum2
print("----")
print(theta)
if(abs(theta[0] - PI) <= tol and abs(theta[1] - P) <= tol \
and abs(theta[2] - Q) <= tol):
print("break")
break
return PI,P,Q
def main():
#mat_y = create_sample_data(100, 1)
mat_y = mat(zeros((10, 1)))
mat_y[0,0] = 1
mat_y[1,0] = 1
mat_y[2,0] = 0
mat_y[3,0] = 1
mat_y[4,0] = 0
mat_y[5,0] = 0
mat_y[6,0] = 1
mat_y[7,0] = 0
mat_y[8,0] = 1
mat_y[9,0] = 1
theta = [0.5, 0.6, 0.5]
print(mat_y)
PI,P,Q = em(mat_y, theta, 0.001, 100)
print(PI, P, Q)
main()
輸出：
[[1.]
[1.]
[0.]
[1.]
[0.]
[0.]
[1.]
[0.]
[1.]
[1.]]
----
[array([0.50505051]), array([0.648]), array([0.55102041])]
----
[array([0.50505051]), array([0.648]), array([0.55102041])]
break
[0.50505051] [0.648] [0.55102041]
來源：https://www.jianshu.com/p/154ee3354b59