1.1 極大似然估計
一般來說，很多互相獨立的事件都發生了的機率是非常低的，而如果這種情況發生了，我們有理由認為這個組合事件發生的機率在所有可能的組合事件中是最大的。假設我們有一個随機信源 S {S} S，其分布類型已知，但參數 θ {\bf{\theta }} θ 未知。那麼，如果我們對信源 S {S} S 進行多次獨立同分布采樣，并且得到了一系列的随機樣本，其樣本值為 X = { x 1 , x 2 , … , x N } {\bf{X}} = \left\{ {{{\bf{x}}_1},{\rm{ }}{{\bf{x}}_2},{\rm{ }} \ldots ,{\rm{ }}{{\bf{x}}_N}} \right\} X={x1,x2,…,xN}，我們就有理由認為 X {\bf{X}} X 出現的機率是最大的。[注： 一般機率論書籍裡以大寫字母為随機變量， 小寫字母為随機變量的值， 這裡為了簡便， 以大寫字母作為采樣值的集合。 采樣值一般為向量， 有時為了友善推導可能隻用到标量]。是以，參數 θ {\bf{\theta }} θ 的優化就可表示為式（1-1），即
這就是極大似然估計（Maximum Likelihood Estimation, MLE）算法的基本思路。因為連乘不友善求導，而且機率的連乘容易導緻小數位溢出，我們令
因為 log {\log} log 函數是嚴格單調的，是以兩者的極值點也是一樣的。那麼求解式（1-3）即可得到最優的參數，
當然這通常隻能保證是局部最優的。
以一進制單高斯分布為例，
根據式（1-3）可得
可以看到，對于單高斯這種分布函數，通過 MLE 來求解最優參數的計算過程還是比較簡單的。
1.2 高斯混合模型
上面我們讨論了 MLE 求解單高斯分布參數的方法，然而，考慮這麼一種情況，我們先從一個有限整數随機數發生器中随機抽取一個數 y {y} y，例如 1 , 2 , . . . , M {1, 2, ..., M} 1,2,...,M，然後根據抽取到的結果決定使用具有不同參數的分布，例如高斯分布，來生成最終的輸出 x {\bf{x}} x，其中
那麼輸出 x {\bf{x}} x 的機率或機率密度可表示為
這就是所謂的多模态（Multimodal）機率模型。如果所用到的分布為高斯分布，則稱為高斯混合模型（Gaussian Mixture Model, GMM），如圖 1-1 所示。可以看到，如果隻用一個高斯分布來描述則會造成比較大的誤差。
圖1-1 高斯混合模型
這時，如果我們想通過極大似然估計來優化各個分布函數的參數，根據式（1-2）可得，
因為 log {\log} log 函數裡面包含了加法，一般來說我們是很難得到式（1-3）的解析解的，是以得借助一些非線性優化方法，例如梯度下降（上升）法等等，但這些方法還是需要用到式（1-8）的梯度，牛頓法甚至還需要用到二階偏導數，除非是計算機自動微分，依靠人工推導還是十分麻煩的。而且，這些方法的收斂精度和收斂速度十分依賴于泰勒公式展開的近似程度以及學習速率的設定，是以具有一定的不穩定性。
然而，高斯分布有着其獨特的性質。對于單高斯模型，根據 MLE 我們可以得到其兩個最優估計的參數剛好為采樣資料的均值和方差，即如式（1-5）所示。那麼對于 GMM 來說，如果我們知道每個樣本的類别，那是否也可以根據式（1-5）對每個高斯分布的參數進行估計？當然，我們是不可能知道這些樣本的類别的，但是，如果每個高斯分布的參數已知，通過貝葉斯定理我們可以得到這些樣本屬于某個類别的後驗機率，即
那麼，我們就可以對各個高斯分布的參數進行軟估計了，即
慢着！這繞了一圈下來不就成了雞生蛋、蛋生雞的問題了？在求某個樣本屬于各個類别的後驗機率時我們需要先知道每個高斯分布的參數，但這些參數又是基于前面求得的後驗機率的。是以，很自然地，我們想到了一種疊代優化的方法，即先随機選擇各個高斯分布的參數，然後通過式（1-9）和（1-10）得到新的參數，依此重複，最後或許我們能夠得到足夠好的估計結果。注意，是“或許”！我們怎麼能夠保證這種疊代方法一定能夠使得參數往極大似然的方向前進呢？我們這裡甚至連梯度都沒有提到！為了解決這個問題，我們需要引入期望最大化（Expectation Maximization, EM）算法。通過 EM 算法，我們将證明，前面所述的疊代方法是可以保證達到一個局部最優點的。這也是為什麼 GMM 應用能夠這麼廣泛的一個原因，除了這種分布在自然中廣泛存在以外，其參數優化的過程也是比較簡單的。
1.3 期望最大化算法
從直覺上，我們認為可通過式（1-9）和（1-10）的疊代來得到高斯混合模型參數的一個最優估計，但是我們目前沒辦法從理論上去證明這種疊代的有效性，而 EM 算法就是解決這個問題的法寶。對于大多數多模态系統來說，我們隻能擷取最終的輸出 x {\bf{x}} x，而無法知道該輸出所屬的類别 y {y} y。是以，我們一般将 y {y} y 稱為隐含變量， x {\bf{x}} x 則稱為不完整資料（Incomplete Data），而 z = ( x , y ) {{\bf{z}} =({\bf{x}}, y)} z=(x,y) 才稱為完整資料（Complete Data）。EM 算法就是要解決這種包含不完整資料的極大似然估計問題。實際上，隐含變量 y {y} y 不一定是一個多模态分布中的類别資訊，其甚至可以隻是某個樣本值 x {\bf{x}} x 中的某個次元或特征缺失的值，但這裡我們隻讨論 y {y} y 屬于類别資訊的情況。
1.3.1 疊代政策
EM 算法首先定義了一種能保證收斂的疊代過程，如圖 1-2 所示。給定一個需要尋找極大值點的函數 f ( x ) {f(x)} f(x)，找到一個輔助函數 A ( x , x t ) {A(x, x^t)} A(x,xt)，其中 x t {x^t} xt 為 t {t} t 時刻所找到的最優參數值，也就是輔助函數是與目前最優解相關的，并且滿足
那麼，找到 A ( x , x t ) {A(x, x^t)} A(x,xt) 的極大值點，則有
這樣，我們就可以保證整個疊代過程是不會惡化的，最終疊代停滞的時候就說明 f ( x ) {f(x)} f(x)到達了一個極大值點。要證明這一點，令
那麼 x t {x^t} xt 一定是 h ( x ) {h(x)} h(x) 的一個極小值點，是以 f ( x ) {f(x)} f(x) 和 A ( x , x t ) {A(x, x^t)} A(x,xt) 在 x t {x^t} xt 處的導數是相等的，也就是兩函數在 x t {x^t} xt 處相切。是以，疊代停滞則意味着
因為一般的優化方法隻能保證 x t + 1 {x^{t+1}} xt+1 是一個局部極大值，是以 EM 算法的結果也隻能最多保證是一個局部最優解。為了擷取更好的結果，通常會選擇不同的初始點進行疊代，加大了找到全局最優解的幾率。
圖1-2 EM 算法的疊代過程
然而怎樣找到一個合适的輔助函數 A ( x , x t ) {A(x, x^t)} A(x,xt)？這似乎毫無頭緒。而且我們不能讓 A ( x , x t ) {A(x, x^t)} A(x,xt) 的形式過于複雜，例如在 log {\log} log 函數中包含了加法項，因為我們還是得求解 A ( x , x t ) {A(x, x^t)} A(x,xt) 的極值點。當然，即便沒有找到 A ( x , x t ) {A(x, x^t)} A(x,xt) 的極值點，隻要滿足 A ( x t + 1 , x t ) > A ( x t , x t ) {A(x^{t+1}, x^{t})>A(x^t, x^t)} A(xt+1,xt)>A(xt,xt)，我們至少還是能夠往好的方向前進的。這樣問題最終還是落在了如何找到合适的 A ( x , x t ) {A(x, x^t)} A(x,xt) 上，我們希望其能夠像單模态分布的似然函數那樣 log {\log} log 函數裡面不包含加法項，這樣求解極值點也會簡便很多。
1.3.2 Jensen不等式
在推導輔助函數 A ( x , x t ) {A(x, x^t)} A(x,xt) 的具體形式之前，我們有必要先介紹一下 Jensen 不等式。首先，一個函數 f ( x ) {f(x)} f(x) 在區間 [ a , b ] {[a, b]} [a,b] 内稱為凸（Convex）函數，如果
如圖 1-3 所示。注意函數的凹凸從直覺上是以函數以上部分的空間形狀來定義的，是以數學上的凸函數看上去實際上凹形的。
圖1-3 凸函數
那麼，對于一個凸函數，根據 Jensen 不等式有
我們可以使用歸納法來證明。 n = 1 {n=1} n=1 時很明顯成立，對于 n = 2 {n=2} n=2 以上的，有
得證。因為 − log ( x ) {-\log(x)} −log(x) 在 ( 0 , ∞ ) {(0, \infty)} (0,∞)是嚴格的凸函數，根據 Jensen 不等式有
利用這個不等式，我們将找到我們所需要的那個輔助函數 A ( x , x t ) {A(x, x^t)} A(x,xt)。
1.3.3 EM算法推導
回到多模态分布的參數優化，有了 Jensen 不等式和式（1-18），我們很自然地可能想改寫式（1-7），即
然而這樣我們無法建立起參數的疊代過程，是以還得另辟蹊徑。為了構造輔助函數，我們引入另外一個機率分布 q ( y ) {q(y)} q(y)，其中 y = 1 , 2 , . . . , M {y=1, 2, ..., M} y=1,2,...,M。那麼有
這似乎把問題複雜化了，而且看上去我們同樣也沒法進行參數的疊代更新。但是，因為 q ( y ) {q(y)} q(y) 是可以自由選擇的，我們或許可以通過 q ( y ) {q(y)} q(y) 建立一個從舊參數到新參數的橋梁，進而使得疊代過程能夠運作起來。
回顧式（1-11），輔助函數除了需要滿足不大于原函數的條件以外，還必須在舊參數處等于原函數。是以，我們首先看看式（1-20）兩側的差距，有
也就是說，式（1-20）中的不等式兩側的差距剛好為 q ( y ) {q(y)} q(y) 和 p ( y ∣ x , Θ ) {p(y|\bf{x, \Theta})} p(y∣x,Θ) 的 KL 散度。我們同樣可以利用 Jensen 不等式證明 KL 散度是非負的，并且當且僅當 q ( y ) {q(y)} q(y) 和 p ( y ∣ x , Θ ) {p(y|\bf{x, \Theta})} p(y∣x,Θ) 相同的時候為 0。因為我們隻要求原函數和輔助函數在舊參數處的值相等，是以我們令
于是我們得到
并且有
是以，下一步的重點是求解輔助函數的極值點，我們希望輔助函數能夠讓這個計算過程變得更加簡單。由式（1-23）可得
那麼，我們實際隻需要求解 Q ( Θ , Θ o l d ) {Q({\bf{\Theta, \Theta}}^{old})} Q(Θ,Θold) 的極值點即可。而且這時 log {\log} log 函數裡面已經沒有加法項了，對于高斯分布來說這能夠極大地簡化極值點的求解。
另外，細心的你可能已經發現了，單讨論一個樣本， Q ( Θ , Θ o l d ) {Q({\bf{\Theta, \Theta}}^{old})} Q(Θ,Θold) 實際上就是 z i {{\bf{z}}_i} zi 也就是 完整資料的 log {\log} log 似然的條件期望值，英文表述為 Conditional expectation of the complete data log lokelihood，這就是 EM 算法的核心意義，也就是說将不完整資料的極大似然估計轉換為完整資料的條件期望對數似然的最優化。因為完整資料的機率一般是條件機率的連乘，而不是不完整資料的在各種可能事件下的機率疊加，進而很好地避免了 log {\log} log 函數内包含加法的情況。 Q ( Θ , Θ o l d ) {Q({\bf{\Theta, \Theta}}^{old})} Q(Θ,Θold) 優化的具體過程我們先留到後面與 GMM 的參數估計一起詳細讨論，現在我們可以總結 EM 算法的基本流程如下：
- (1) 選擇一個初始的參數集 Θ o l d {{\bf{\Theta}}^{old}} Θold，但為了避免遇到較差的局部最優解，一般會選擇多組的初始參數集獨立優化。
- (2) E Step：根據式（1-9）計算後驗機率 p ( y ∣ x i , Θ o l d ) {p(y|{\bf{x}}_i, {\bf{\Theta}}^{old})} p(y∣xi,Θold)，構造 Q ( Θ , Θ o l d ) {Q({\bf{\Theta, \Theta}}^{old})} Q(Θ,Θold)。
- (3) M Step：求解 Θ n e w = a r g m a x Q ( Θ , Θ o l d ) {{\bf{\Theta}}^{new} = argmax Q({\bf{\Theta, \Theta}}^{old})} Θnew=argmaxQ(Θ,Θold)。
- (4) 檢查參數值的變化，如 ∣ ∣ Θ o l d − Θ n e w ∣ ∣ {||{\bf{\Theta}}^{old}-{\bf{\Theta}}^{new}||} ∣∣Θold−Θnew∣∣，若參數值收斂，算法結束；否則更新 Θ o l d ← Θ n e w {{\bf{\Theta}}^{old}\gets{\bf{\Theta}}^{new}} Θold←Θnew，傳回第（2）步。
1.3.4 GMM的參數優化
了解了 EM 算法的基本原理後，我們現在就可以去推導 GMM 的參數優化流程了，進而證明式（1-9）和（1-10）的疊代方式是否有效。為了簡便，這裡隻讨論一進制的高斯分布，但同理也可以推廣到多元的情況。
首先，一進制 GMM 可表示為
這看上去似乎挺複雜的。根據式（1-25），我們構造需要優化的函數 Q ( Θ , Θ o l d ) {Q({\bf{\Theta, \Theta}}^{old})} Q(Θ,Θold)，其中後驗機率 p ( y ∣ x i , Θ o l d ) {p(y|x_i, {\bf{\Theta}}^{old})} p(y∣xi,Θold) 可根據式（1-9）計算，可得
前面我們一直假設 y {y} y 的機率分布是已知的，但我們通常隻能得到最終的輸出 x {x} x，是以 α y {\alpha_y} αy 實際上也是需要進行估計的參數，不過将其加到 Θ {\bf{\Theta}} Θ 中同樣是滿足 EM 算法的要求的，為了友善這裡還是對其單獨讨論，因為從式（1-27）也可以看出 α y {\alpha_y} αy 和 Θ {\bf{\Theta}} Θ 極值點的計算是互相獨立的。接下來，我們先對 y {y} y 的機率分布進行估計。因為 y {y} y 所有取值的機率之和為 1，是以這是一個等式限制的最優化問題。我們引入拉格朗日算子 λ {\lambda} λ，那麼其最優解滿足
注意， p ( y ∣ x i , Θ o l d ) {p(y|x_i, {\bf{\Theta}}^{old})} p(y∣xi,Θold) 是基于舊參數計算的，是以在對新參數求導時相當于常數。那麼求解式（1-28）可得
這樣我們就可以知道，對于 Q ( Θ , Θ o l d ) {Q({\bf{\Theta, \Theta}}^{old})} Q(Θ,Θold) 來說， y {y} y 等于某個值的機率的最優估計實際上就是各個樣本在目前參數下屬于該類的後驗機率的平均值。同樣，我們可以計算各個高斯分布的參數對于 Q ( Θ , Θ o l d ) {Q({\bf{\Theta, \Theta}}^{old})} Q(Θ,Θold) 的最優估計。因為
根據式（1-27）可得
同理可得
這樣，我們就比較輕松地得到了對于 Q ( Θ , Θ o l d ) {Q({\bf{\Theta, \Theta}}^{old})} Q(Θ,Θold) 的各個高斯分布的參數的最優估計，而這剛好與我們最初猜想的式（1-10）是一緻的，即首先根據現有參數計算出每個樣本屬于某個類的後驗機率，然後以這個機率作為一種權重對該類的單高斯分布的參數進行估計，得到新的參數集合，并利用這個新的參數集進行疊代，進而在EM算法的理論支援下最終收斂到一個局部最優解，使得采樣到的不完整資料的出現機率達到最大。
綜上所述，使用 EM 算法求解 GMM 的參數的流程如圖 1-4 所示。
圖1-4 EM 算法優化 GMM 參數流程
腳本示例
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
def get_gmm_data(n, means, covs, alpha):
"""
生成 GMM 随機樣本，這裡采用了比較直覺的方法，即首先根據 alpha 機率随機選擇樣本所屬的高斯分布，
然後再根據該高斯分布的參數生成随機生成最終的樣本。這種方法相對比較低效，需要 n 次循環。
Parameters
----------
n : int
生成随機樣本的個數。
means : list of 1D array
各個高斯分布的均值，必須使用 list 來存放，或者使用 m x k 的二維矩陣，m 為高斯分布的個數。
covs : list of 2D array
各個高斯分布的協方差矩陣，必須使用 list 來存放，或者使用使用 m x k x k 的三維矩陣，k 為樣本的維數。
alpha : list or 1D array
樣本屬于各個高斯分布的機率，總和必須為 1。
Returns
-------
k x n 的二維數組，其中 k 為樣本維數，可以從 means 推導。
注意，為了友善進行矩陣的運算，腳本中所有樣本都以列向量的形式存在，即每一列代表一個樣本。
"""
models = len(alpha)
xs = []
for i in range(n):
j = np.random.choice(range(models), p=alpha)
xs.append(np.random.multivariate_normal(means[j].reshape(-1), covs[j]))
return np.array(xs).T
def gauss_prob(x, mean, cov):
"""
根據給定 mean, cov 的高斯分布的機率密度函數計算樣本 x 的機率密度。注意隻适用于單高斯分布。
Parameters
----------
x : scalar or array
可以為單個樣本，也可以為多個樣本。樣本既可以是标量，也可以是列向量（單個樣本也應該是 k x 1 的形式）。
mean : scalar or array
scalar 代表單變量的随機分布，array 代表多變量的随機分布。
可以以行向量或列向量形式輸入，list 形式也可以，實際上包括 scalar 都會首先被轉換為 k x 1 的列向量形式。
cov : scalar or 2D array
如果 mean 是 scalar, cov 也必須為 scalar，否則應該為 k x k 的二維矩陣。
Returns
-------
如果 x 是單個 scalar 樣本，那麼也傳回一個 scalar，否則都傳回一個 shape 為 (n,) 的 1D array。
"""
return_scalar = True if np.shape(x) is () else False
mean = np.reshape(mean, [-1, 1])
dim = mean.shape[0]
x = np.reshape(x, [dim, -1])
cov = np.reshape(cov, [dim, dim])
# 參考多元的高斯機率密度函數公式，注意這裡可以一次性求得所有樣本的機率密度
tmp = (x - mean).T
tmp = np.sum(tmp.dot(np.linalg.inv(cov)) * tmp, axis=1)
prob = ((2*np.pi)**dim * np.linalg.det(cov))**(-0.5) * np.exp(-0.5 * tmp)
if return_scalar:
prob = prob[0]
return prob
def gmm_em(xs, models):
"""
基于 EM 算法估計 GMM 參數的具體實作。
Parameters
----------
xs : 2D array
已有的樣本，應該為 k x n 的形式，k 為樣本維數，n 為樣本數量。
注意即便 k=1 樣本也應該使用列向量的形式，因為代碼大部分采用了矩陣運算。
models : int
高斯分布個數。
Returns
-------
alpha : 1D array
樣本屬于各個高斯分布的機率
means : lsit of 1D array
各個高斯分布的均值，以 list 形式存儲
covs : list of 2D array
各個高斯分布的協方差矩陣，以 list 形式存儲
"""
# 根據樣本的總體均值和協方差矩陣随機初始化各個高斯分布的參數，這種方式僅供參考
dims, n = xs.shape[:2]
mean = np.mean(xs, axis=1, keepdims=True)
cov = (xs - mean).dot((xs - mean).T) / n
alpha = np.random.rand(models) + 0.1
alpha = alpha / np.sum(alpha)
means, covs = [], []
rnd_range = (np.max(xs, axis=1) - np.min(xs, axis=1)) * 0.8 # 樣本資料的範圍乘以一個縮放系數
for i in range(models):
means.append(mean.ravel() + (np.random.rand(dims) - 0.5) * rnd_range) # 總體均值加上一個随機數
covs.append(np.copy(cov)) # 直接使用總體協方差矩陣
means, covs = np.array(means), np.array(covs)
""" EM 算法流程 """
for epoch in range(100):
# 根據各個高斯分布的參數求得每個樣本屬于各個高斯分布的後驗機率，即 E step
probs = np.array([gauss_prob(xs, _mean, _cov) for _mean, _cov in zip(means, covs)])
post_probs = np.diag(alpha).dot(probs)
post_probs = post_probs / np.sum(post_probs, axis=0)
# 根據求得的後驗機率更新各個高斯分布的參數，即 M step
new_alpha = np.sum(post_probs, axis=1)
new_means = xs.dot(post_probs.T) / new_alpha
new_covs = []
for i in range(models):
tmp = xs - new_means[:, [i]]
new_covs.append(tmp.dot(np.diag(post_probs[i])).dot(tmp.T) / new_alpha[i])
new_alpha = new_alpha / n
new_means = new_means.T
new_covs = np.array(new_covs)
# 判斷參數的更新是否停滞，這裡使用了無窮階範數，即 max(abs(x))，也可以改為 2nd 範數
diff_alpha = np.linalg.norm(np.ravel(alpha - new_alpha), np.inf)
diff_means = np.linalg.norm(np.ravel(means - new_means), np.inf)
diff_covs = np.linalg.norm(np.ravel(covs - new_covs), np.inf)
alpha = new_alpha
means = new_means
covs = new_covs
if diff_alpha <= 1e-3 and diff_means <= 1e-3 and diff_covs <= 1e-3:
print('Solved in {} epoches.'.format(epoch))
break
return alpha, list(means), list(covs)
if __name__ == '__main__':
means = [np.array([5, 6]), np.array([1, 2])]
covs = [np.array([[1, 0.1], [0.1, 0.5]]), np.array([[0.8, 1.2], [1.2, 3.5]])]
alpha = np.array([0.3, 0.7])
n = 2000
xs = get_gmm_data(n, means, covs, alpha)
#plt.scatter(xs[0], xs[1])
#plt.show()
_alpha, _means, _covs = gmm_em(xs, 2)
print(_alpha)
print(_means)
print(_covs)
延伸閱讀
http://www.seanborman.com/publications/EM_algorithm.pdf
https://statlect.com/fundamentals-of-probability/Jensen-inequality
https://statlect.com/fundamentals-of-probability/Kullback-Leibler-divergence