高斯混合模型(Gaussian Mixture Model)是機器學習中一種常用的聚類演算法，本文介紹了其原理，並推導了其參數估計的過程。主要參考Christopher M. Bishop的《Pattern Recognition and Machine Learning》。

以粗體小寫字母表示向量，粗體大寫字母表示矩陣；標量不加粗，大寫表示常數。

1. 高斯分布

高斯分布(Gaussian distribution)，也稱為正態分布(normal distribution)，是一種常用的連續變數分布的模型。若單個隨機變數 x 服從均值為 mu ，方差為 sigma^2 的高斯分布，記為 xsim N left( mu, sigma^2 ight) ，則其概率密度函數為：

p(x|mu,sigma^2) = frac{1}{left(2pisigma^2 ight)^{1/2}} exp{left{ -frac{1}{2sigma^2} (x-mu)^2 ight } }

對於一個 D 維的向量 mathbf{x} ，若其各元素服從均值為向量 mathbf{mu} ，協方差矩陣為 mathbf{Sigma} 的多元高斯分布，記為 mathbf{x} sim N left( mathbf{mu,Sigma} ight) ，則概率密度為：

p(mathbf{x|mu,Sigma}) = frac{1}{left(2pi ight)^{D/2}} frac{1}{mathbf{vert Sigma} vert ^{1/2}} exp{left{ -frac{1}{2} mathbf{ left(x-mu ight)^T Sigma^{-1} left(mathbf{x-mu} ight)} ight } } ag{1}

其中 mathbf{mu} 為 D 維均值向量， mathbf{Sigma} 為 D imes D 的協方差矩陣， vert mathbf{Sigma} vert 表示 mathbf{Sigma} 的行列式。

(1)式中，指數部分的二次型 Delta^2 = mathbf{ left(x-mu ight)^T Sigma^{-1} left(mathbf{x-mu} ight)} 稱為 mathbf{x} 到 mathbf{mu} 的馬哈拉諾比斯距離(馬氏距離，Mahalanobis distance)；當 mathbf{Sigma} 為單位矩陣時退化為歐幾里得距離(Euclidean distance)。多元高斯分布密度函數的等高線即 Delta^2 為常數時 mathbf{x} 的方程，是橢球方程(Ellipsoid - Wikipedia)。

2. 高斯混合模型(Gaussian Mixture Model)

多個高斯分布的線性疊加能擬合非常複雜的密度函數；通過足夠多的高斯分布疊加，並調節它們的均值，協方差矩陣，以及線性組合的係數，可以精確地逼近任意連續密度([1], Section 2.3.9, p111)。

我們考慮 K 個高斯分布的線性疊加，這個高斯混合分布(Gaussian mixture distiburion)的概率密度函數為：

p(mathbf{x}) = sum_{k=1}^{K}{pi_k pleft(mathbf{x | mu_k, Sigma_k} ight)} ag{2}

其中， pleft(mathbf{x | mu_k, Sigma_k} ight) 表示參數為 mathbf{mu_k, Sigma_k} 的高斯分布的概率密度。

我們稱(2)式為一個高斯混合(Mixture of Gaussians, Gaussian Mixture)。其中每個高斯密度函數稱為混合的一個分模型(component)，有自己的均值 mathbf{mu_k} 和協方差矩陣 mathbf{Sigma_k} 。

(2)式中的參數 pi_k 是模型的混合係數(mixing coefficients)。將(2)式左右兩側對 mathbf{x} 積分，得到

sum_{k=1}^{K}{pi_k} = 1

此外，由於 p(mathbf{x})geq 0 ， pleft(mathbf{x | mu_k, Sigma_k} ight)geq 0 ，所以 pi_kgeq 0 。即混合係數應滿足 0leq pi_k leq 1 。

(更一般地，混合模型也可以是其他任意分布的疊加。)

由全概率公式， mathbf{x} 的邊緣分布(marginal distribution)的概率密度為：

p(mathbf x) = sum_{k=1}^{K}{p(k)p(mathbf x|k)}

上式與(2)式等價，其中 p(k) = pi_k 可以看作選擇第 k 個分模型的先驗概率(prior probability)， p(mathbf x|k) 是 mathbf x 對 k 的條件概率密度。

在觀測到 mathbf x 後，其來自第 k 個分模型的後驗概率(posterior probability) p(k |mathbf x) 稱為第 k 個分模型的響應度(responsibility)。

下圖所示為包含兩個一維分模型的高斯混合：

兩個一維分模型的高斯混合

3. 隱變數 & 完全數據

引入一個 K 維的二值型隨機變數 mathbf{z} = left[ z_1, ldots, z_K ight] ，來表示樣本 mathbf{x} 由哪一分模型產生。 mathbf{z} 滿足這樣的條件：z_k in left{ 0, 1 ight} ，且 sum_{k=1}^{K}{z_k} = 1，即其 K 個元素中，有且只有一個為1，其餘為0。 z_k = 1 表示樣本 mathbf{x} 由分模型 k 抽樣得到。可以看出， mathbf{z} 一共有 K 種可能的狀態。

mathbf{z} 的邊緣分布由混合係數給出： pleft( z_k = 1 ight) = pi_k 。也可寫成如下形式：

p(mathbf{z}) = prod_{k=1}^{K} pi_k^{z_k}

考慮由以下方式產生樣本 mathbf{x}：

  1. 先以離散分布 p(mathbf{z}) 抽樣得到變數 mathbf{z} ；
  2. 設根據 mathbf{z} 的取值選擇了第 k 個分模型，則以高斯分布 pleft(mathbf{x | mu_k, Sigma_k} ight) 抽樣得到 mathbf{x} 。

記高斯混合模型的參數為 mathbf{pi} equiv left{ pi_1, ldots, pi_K ight} ， mathbf{mu} equiv left{ mathbf{mu_1, ldots, mu_K} ight} ， mathbf{Sigma} equiv left{ mathbf{Sigma_1, ldots, Sigma_K} ight} ，則這個過程可由如下的graphical model表示[1]：

高斯混合模型的graphical model

變數 mathbf{z} 稱為隱變數(latent variable)，包含 (mathbf{x,z}) 取值的樣本稱為完全數據(complete data)，只含有 mathbf{x} 取值的樣本稱為不完全數據(incomplete data)，

給定 mathbf{z} 後 mathbf{x} 的條件概率密度為：

p(mathbf{x} | z_k=1) = pleft(mathbf{x | mu_k, Sigma_k} ight)

或者寫成如下形式：

p(mathbf{x|z}) = prod_{k=1}^{K} pleft(mathbf{x | mu_k, Sigma_k} ight)^{z_k}

mathbf{x} 的邊緣分布為聯合概率分布 p(mathbf{x,z}) 對 mathbf{z} 的所有可能狀態求和：

egin{align} p(mathbf x) &= sum_{mathbf z} {p(mathbf{x,z})} = sum_{mathbf z} {p(mathbf{z}) p(mathbf{x|z})} \ &= sum_{mathbf z} {left( prod_{k=1}^{K} pi_k^{z_k} prod_{k=1}^{K} pleft(mathbf{x | mu_k, Sigma_k} ight)^{z_k} ight)} \ &= sum_{mathbf z} {left( prod_{k=1}^{K} left[ pi_k pleft(mathbf{x | mu_k, Sigma_k} ight) ight]^{z_k} ight)} \ &= sum_{k=1}^{K}{pi_k pleft(mathbf{x | mu_k, Sigma_k} ight)} end{align}

也可由下面的式子得到：

p(mathbf x) = sum_{j=1}^{K}{p(z_j=1) p(mathbf x | z_j=1)} = sum_{j=1}^{K}{pi_j pleft(mathbf{x | mu_j, Sigma_j} ight)}

這表明 mathbf{x} 的邊緣分布就是高斯混合的形式。

4. 後驗概率 & 響應度

根據貝葉斯定理(Bayes theorem - Wikipedia)，觀測到 mathbf{x} 後，其來自第 k 個分模型的後驗概率(posterior responsibility)為：

egin{align} p(z_k=1|mathbf{x}) &= frac{p(z_k=1)p(mathbf x | z_k=1)}{p(mathbf x)}\ &= frac{pi_k pleft(mathbf{x | mu_k, Sigma_k} ight)} {sum_{j=1}^{K}{pi_j pleft(mathbf{x | mu_j, Sigma_j} ight)} } \ end{align} ag{3}

上式中， p(z_k=1|mathbf{x}) 和 p(z_k=1) 為概率， p(mathbf x | z_k=1) 和 p(mathbf x) 為概率密度。

將上式定義為： gamma(z_k)equiv p(z_k=1|mathbf{x}) ，稱為第 k 個分模型對 mathbf x 的響應度(responsibility)[2]。

對於樣本集 mathbf{X} = left{ mathbf{x_1, ldots, x_N} ight} ，記 mathbf{x_n} 對應的隱變數為 mathbf{z_n} = left[ z_{n1}, ldots, z_{nK} ight] ， n = 1,ldots, N ，則第 k 個分模型對 mathbf{x_n} 的響應度為：

gamma(z_{nk})equiv p(z_{nk}=1|mathbf{x_n}) = frac{pi_k pleft(mathbf{x_n | mu_k, Sigma_k} ight)} {sum_{j=1}^{K}{pi_j pleft(mathbf{x_n | mu_j, Sigma_j} ight)} } ag{4}

5. 對數似然函數 & 最大似然估計

現在有一個樣本集 mathbf{X} = left{ mathbf{x_1, ldots, x_N} ight} ，我們假設 mathbf{X} 是由一個高斯混合模型產生的，要估計這個模型的參數： mathbf{pi} equiv left{ pi_1, ldots, pi_K ight} ， mathbf{mu} equiv left{ mathbf{mu_1, ldots, mu_K} ight} ， mathbf{Sigma} equiv left{ mathbf{Sigma_1, ldots, Sigma_K} ight} 。

樣本集 mathbf{X} (不完全數據)的似然函數(likelihood function)為：

p(mathbf{X|pi,mu,Sigma}) = prod_{n=1}^{N} left[ sum_{k=1}^{K}{pi_k pleft(mathbf{x_n | mu_k, Sigma_k} ight)} ight]

似然函數中的連乘求導比較麻煩，取自然對數將其轉換成對數的和，得到對數似然函數(the log of the likelihood function)：

ln p(mathbf{X|pi,mu,Sigma}) = sum_{n=1}^{N}{lnleft[ sum_{k=1}^{K}{pi_k pleft(mathbf{x_n | mu_k, Sigma_k} ight)} ight]} ag{5}

其中，

p(mathbf {x_n|mu_k,Sigma_k}) = frac{1}{left(2pi ight)^{D/2}} frac{1}{mathbf{vert Sigma_k} vert ^{1/2}} exp{left{ -frac{1}{2} mathbf{ left(x_n-mu_k ight)^T Sigma_k^{-1} left(mathbf{x_n-mu_k} ight)} ight } }

最大似然估計(maximum likelihood estimation)，即通過求似然函數的最大值來估計概率模型的參數。用最大似然估計來計算高斯混合模型的參數，形成如下的優化問題：

max_{mathbf{pi,mu,Sigma}} ln p(mathbf{X|pi,mu,Sigma}) \ s.t. sum_{k=1}^{K}{pi_k}=1

採用拉格朗日乘子法來求極值，拉格朗日函數為：

L(mathbf{pi,mu,Sigma}) = ln p(mathbf{X|pi,mu,Sigma}) + lambda left( sum_{k=1}^{K}{pi_k}-1 ight)

先將 L(mathbf{pi,mu,Sigma}) 對 mathbf {mu_k} 求梯度。因為

frac{partial p(mathbf {x_n|mu_k,Sigma_k})}{partial mathbf{mu_k}} = p(mathbf {x_n|mu_k,Sigma_k}) left[Sigma_k^{-1} left(mathbf{x_n-mu_k} ight) ight]

所以，

egin{align} frac{partial L}{partial mathbf{mu_k}} &= frac{partial ln p(mathbf{X|pi,mu,Sigma}) }{partial mathbf{mu_k}} \ &= frac{partial left(sum_{n=1}^{N}{lnleft[ sum_{k=1}^{K} {pi_k pleft(mathbf{x_n | mu_k, Sigma_k} ight)} ight]} ight)} {partial mathbf{mu_k}} \ &= sum_{n=1}^{N}{ frac{partial lnleft[ sum_{k=1}^{K} {pi_k pleft(mathbf{x_n | mu_k, Sigma_k} ight) } ight]} {partial mathbf{mu_k}}} \ &= sum_{n=1}^{N}{left[ frac{pi_k pleft(mathbf{x_n | mu_k, Sigma_k} ight)}{sum_{j=1}^{K}{pi_j pleft(mathbf{x_n | mu_j, Sigma_j} ight)}} left[mathbf{Sigma_k}^{-1} left(mathbf{x_n-mu_k} ight) ight] ight]} end{align}

注意上式中 frac{pi_k pleft(mathbf{x_n | mu_k, Sigma_k} ight)} {sum_{j=1}^{K}{pi_j pleft(mathbf{x_n | mu_j, Sigma_j} ight)} } 即為響應度 gamma(z_{nk}) ，所以有：

egin{align} frac{partial L}{partial mathbf{mu_k}} &= sum_{n=1}^{N}{left[ gamma(z_{nk})mathbf{Sigma_k}^{-1} left(mathbf{x_n-mu_k} ight) ight]} \ end{align}

令 frac{partial L}{partial mathbf{mu_k}} = mathbf 0 ，則 sum_{n=1}^{N}{left[ gamma(z_{nk})mathbf{Sigma_k}^{-1} left(mathbf{x_n-mu_k} ight) ight]} = mathbf 0 ，

左乘 mathbf{Sigma_k} ，整理得

egin{align} sum_{n=1}^{N}{left[ gamma(z_{nk})mathbf{x_n} ight]} = sum_{n=1}^{N}{left[ gamma(z_{nk})mathbf{mu_k} ight]} = mathbf{mu_k} sum_{n=1}^{N}{gamma(z_{nk})}\ end{align}

定義 N_k = sum_{n=1}^{N}{gamma(z_{nk})} ，則

egin{align} mathbf{mu_k} = frac{1}{N_k} sum_{n=1}^{N}{left[ gamma(z_{nk})mathbf{x_n} ight]} end{align}

N_k 可理解為被分配到第 k 個分模型(聚類)的「有效「的樣本數。

對 mathbf{Sigma_k} 中各元素求偏導：

egin{align} frac{partial L}{partial mathbf{Sigma_k}} &= frac{partial ln p(mathbf{X|pi,mu,Sigma}) }{partial mathbf{Sigma_k}} \ &= frac{partial left(sum_{n=1}^{N}{lnleft[ sum_{k=1}^{K} {pi_k pleft(mathbf{x_n | mu_k, Sigma_k} ight)} ight]} ight)} {partial mathbf{Sigma_k}} \ &= sum_{n=1}^{N}{ frac{partial lnleft[ sum_{k=1}^{K} {pi_k pleft(mathbf{x_n | mu_k, Sigma_k} ight) } ight]} {partial mathbf{Sigma_k}}} \ &= sum_{n=1}^{N}{left[ frac{pi_k }{sum_{j=1}^{K}{pi_j pleft(mathbf{x_n | mu_j, Sigma_j} ight)}} frac{partial pleft(mathbf{x_n | mu_k, Sigma_k} ight)}{partial mathbf{Sigma_k}} ight]} end{align}

接下來涉及矩陣求導(Matrix calculus - Wikipedia)，要複雜一些，這裡不做推導，按參考文獻[1]Chapter 9的公式(9.19)，給出 frac{partial L}{partial mathbf{Sigma_k}} = mathbf{0} 的結果：

mathbf{Sigma_k} = frac{1}{N_k} sum_{n=1}^{N}{gamma(z_{nk})(mathbf{x_n-mu_k})(mathbf{x_n-mu_k})^T}

對 pi_k 求偏導並令其為0：

egin{align} frac{partial L}{partial pi_k} &= frac{partial ln p(mathbf{X|pi,mu,Sigma}) }{partial pi_k} \ &= frac{partial left[sum_{n=1}^{N}{lnleft[ sum_{k=1}^{K} {pi_k pleft(mathbf{x_n | mu_k, Sigma_k} ight)} ight]+lambdaleft(sum_{k=1}^{K}{pi_k-1} ight)} ight]} {partial pi_k} \ &= sum_{n=1}^{N}{frac{pleft(mathbf{x_n | mu_k, Sigma_k} ight)}{sum_{j=1}^{K}{pi_j pleft(mathbf{x_n | mu_j, Sigma_j} ight)}} }+lambda=0 \ end{align}

注意到上式左右兩邊乘以 pi_k 可湊出 gamma(z_{nk}) ，所以有

egin{align} sum_{n=1}^{N}{frac{pi_k pleft(mathbf{x_n | mu_k, Sigma_k} ight)}{sum_{j=1}^{K}{pi_j pleft(mathbf{x_n | mu_j, Sigma_j} ight)}} }+lambda pi_k &= sum_{n=1}^{N}{gamma(z_{nk})}+lambda pi_k \&=N_k + lambda pi_k = 0 end{align}

上式對 k 求和得

sum_{k=1}^{K}{(N_k + lambda pi_k)} = sum_{k=1}^{K}{N_k} + lambda sum_{k=1}^{K}{pi_k} = sum_{k=1}^{K}{N_k} + lambda = 0

又

egin{align} sum_{k=1}^{K}{N_k} &= sum_{k=1}^{K}{sum_{n=1}^{N}{gamma(z_{nk})}} = sum_{n=1}^{N}{sum_{k=1}^{K}{gamma(z_{nk})}} \&= sum_{n=1}^{N}{sum_{k=1}^{K}{frac{pi_k pleft(mathbf{x | mu_k, Sigma_k} ight)} {sum_{j=1}^{K}{pi_j pleft(mathbf{x | mu_j, Sigma_j} ight)} }}} \&= sum_{n=1}^{N}{1} = N end{align}

所以 lambda = -N ，進而 pi_k = -frac{N_k}{lambda} = frac{N_k}{N} 。

綜上，對數似然函數 ln p(mathbf{X|pi,mu,Sigma}) 的極值點滿足的條件為：

left{ egin{align} mathbf{mu_k} &= frac{1}{N_k} sum_{n=1}^{N}{left[ gamma(z_{nk})mathbf{x_n} ight]}\ mathbf{Sigma_k} &= frac{1}{N_k} sum_{n=1}^{N}{gamma(z_{nk})(mathbf{x_n-mu_k})(mathbf{x_n-mu_k})^T} \ pi_k &= frac{N_k}{N} end{align}\ ight. ag{6}

需要注意的是，上式並未給出高斯混合模型的解析解/閉式解(analytical soluiton/closed-form solution)，因為響應度 gamma(z_{nk}) 由式(4)給出，而參數 mathbf{mu_k, Sigma_k} 未知，故 gamma(z_{nk}), N_k 無法計算。

不過，根據(6)式可使用迭代演算法來計算模型參數。

6. EM演算法計算高斯混合模型的參數見後續。

參考文獻

[1] Christopher M. Bishop. Pattern Recognition and Machine Learning, Springer, 2006.

[2] 李航，統計學習方法，清華大學出版社，2012年。


推薦閱讀：
相关文章