高斯混合模型（GMM）及其EM算法的了解

一個例子

高斯混合模型（Gaussian Mixed Model）指的是多個高斯分布函數的線性組合，理論上GMM可以拟合出任意類型的分布，通常用于解決同一集合下的資料包含多個不同的分布的情況（或者是同一類分布但參數不一樣，或者是不同類型的分布，比如正态分布和伯努利分布）。

如圖1，圖中的點在我們看來明顯分成兩個聚類。這兩個聚類中的點分别通過兩個不同的正态分布随機生成而來。但是如果沒有GMM，那麼隻能用一個的二維高斯分布來描述圖1中的資料。圖1中的橢圓即為二倍标準差的正态分布橢圓。這顯然不太合理，畢竟肉眼一看就覺得應該把它們分成兩類。

圖1

這時候就可以使用GMM了！如圖2，資料在平面上的空間分布和圖1一樣，這時使用兩個二維高斯分布來描述圖2中的資料，分别記為 N ( μ 1 , Σ 1 ) \mathcal{N}(\boldsymbol{\mu}_1, \boldsymbol{\Sigma}_1) N(μ1​,Σ1​)和 N ( μ 2 , Σ 2 ) \mathcal{N}(\boldsymbol{\mu}_2, \boldsymbol{\Sigma}_2) N(μ2​,Σ2​). 圖中的兩個橢圓分别是這兩個高斯分布的二倍标準差橢圓。可以看到使用兩個二維高斯分布來描述圖中的資料顯然更合理。實際上圖中的兩個聚類的中的點是通過兩個不同的正态分布随機生成而來。如果将兩個二維高斯分布 N ( μ 1 , Σ 1 ) \mathcal{N}(\boldsymbol{\mu_1}, \boldsymbol{\Sigma}_1) N(μ1​,Σ1​)和 N ( μ 2 , Σ 2 ) \mathcal{N}(\boldsymbol{\mu}_2, \boldsymbol{\Sigma}_2) N(μ2​,Σ2​)合成一個二維的分布，那麼就可以用合成後的分布來描述圖2中的所有點。最直覺的方法就是對這兩個二維高斯分布做線性組合，用線性組合後的分布來描述整個集合中的資料。這就是高斯混合模型（GMM）。

圖2

高斯混合模型（GMM）

設有随機變量 X \boldsymbol{X} X，則混合高斯模型可以用下式表示：

p ( x ) = ∑ k = 1 K π k N ( x ∣ μ k , Σ k ) p(\boldsymbol{x}) = \sum_{k=1}^K\pi_k \mathcal{N}(\boldsymbol{x}|\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) p(x)=k=1∑K​πk​N(x∣μk​,Σk​)

其中 N ( x ∣ μ k , Σ k ) \mathcal{N}(\boldsymbol{x}|\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) N(x∣μk​,Σk​)稱為混合模型中的第 k k k個分量（component）。如前面圖2中的例子，有兩個聚類，可以用兩個二維高斯分布來表示，那麼分量數 K = 2 K = 2 K=2. π k \pi_k πk​是混合系數（mixture coefficient），且滿足：

∑ k = 1 K π k = 1 \sum_{k=1}^K\pi_k = 1 k=1∑K​πk​=1

0 ≤ π k ≤ 1 0 \leq \pi_k \leq 1 0≤πk​≤1

實際上，可以認為 π k \pi_k πk​就是每個分量 N ( x ∣ μ k , Σ k ) \mathcal{N}(\boldsymbol{x}|\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) N(x∣μk​,Σk​)的權重。

GMM的應用

GMM常用于聚類。如果要從 GMM 的分布中随機地取一個點的話，實際上可以分為兩步：首先随機地在這 K 個 Component 之中選一個，每個 Component 被選中的機率實際上就是它的系數 π k \pi_k πk​ ，選中 Component 之後，再單獨地考慮從這個 Component 的分布中選取一個點就可以了──這裡已經回到了普通的 Gaussian 分布，轉化為已知的問題。

将GMM用于聚類時，假設資料服從混合高斯分布（Mixture Gaussian Distribution），那麼隻要根據資料推出 GMM 的機率分布來就可以了；然後 GMM 的 K 個 Component 實際上對應 K K K個 cluster 。根據資料來推算機率密度通常被稱作 density estimation 。特别地，當我已知（或假定）機率密度函數的形式，而要估計其中的參數的過程被稱作『參數估計』。

例如圖2的例子，很明顯有兩個聚類，可以定義 K = 2 K=2 K=2. 那麼對應的GMM形式如下：

p ( x ) = π 1 N ( x ∣ μ 1 , Σ 1 ) + π 2 N ( x ∣ μ 2 , Σ 2 ) p(\boldsymbol{x}) =\pi_1 \mathcal{N}(\boldsymbol{x}|\boldsymbol{\mu}_1, \boldsymbol{\Sigma}_1) + \pi_2 \mathcal{N}(\boldsymbol{x}|\boldsymbol{\mu}_2, \boldsymbol{\Sigma}_2) p(x)=π1​N(x∣μ1​,Σ1​)+π2​N(x∣μ2​,Σ2​)

上式中未知的參數有六個： ( π 1 , μ 1 , Σ 1 ; π 2 , μ 2 , Σ 2 ) (\pi_1, \boldsymbol{\mu}_1, \boldsymbol{\Sigma}_1; \pi_2, \boldsymbol{\mu}_2, \boldsymbol{\Sigma}_2) (π1​,μ1​,Σ1​;π2​,μ2​,Σ2​). 之前提到GMM聚類時分為兩步，第一步是随機地在這 K K K個分量中選一個，每個分量被選中的機率即為混合系數 π k \pi_k πk​. 可以設定 π 1 = π 2 = 0.5 \pi_1 = \pi_2 = 0.5 π1​=π2​=0.5，表示每個分量被選中的機率是0.5，即從中抽出一個點，這個點屬于第一類的機率和第二類的機率各占一半。但實際應用中事先指定 π k \pi_k πk​的值是很笨的做法，當問題一般化後，會出現一個問題：當從圖2中的集合随機選取一個點，怎麼知道這個點是來自 N ( x ∣ μ 1 , Σ 1 ) N(\boldsymbol{x}|\boldsymbol{\mu}_1, \boldsymbol{\Sigma}_1) N(x∣μ1​,Σ1​)還是 N ( x ∣ μ 2 , Σ 2 ) N(\boldsymbol{x}|\boldsymbol{\mu}_2, \boldsymbol{\Sigma}_2) N(x∣μ2​,Σ2​)呢？換言之怎麼根據資料自動确定 π 1 \pi_1 π1​和 π 2 \pi_2 π2​的值？這就是GMM參數估計的問題。要解決這個問題，可以使用EM算法。通過EM算法，我們可以疊代計算出GMM中的參數： ( π k , x k , Σ k ) (\pi_k, \boldsymbol{x}_k, \boldsymbol{\Sigma}_k) (πk​,xk​,Σk​).

GMM參數估計過程

GMM的貝葉斯了解

在介紹GMM參數估計之前，先改寫GMM的形式，改寫之後的GMM模型可以友善地使用EM估計參數。GMM的原始形式如下：

p ( x ) = ∑ k = 1 K π k N ( x ∣ μ k , Σ k ) (1) p(\boldsymbol{x}) = \sum_{k=1}^K\pi_k \mathcal{N}(\boldsymbol{x}|\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \tag{1} p(x)=k=1∑K​πk​N(x∣μk​,Σk​)(1)

前面提到 π k \pi_k πk​可以看成是第 k k k類被選中的機率。我們引入一個新的 K K K維随機變量 z \boldsymbol{z} z. z k ( 1 ≤ k ≤ K ) z_k (1 \leq k \leq K) zk​(1≤k≤K)隻能取0或1兩個值； z k = 1 z_k = 1 zk​=1表示第 k k k類被選中的機率，即： p ( z k = 1 ) = π k p(z_k = 1) = \pi_k p(zk​=1)=πk​；如果 z k = 0 z_k = 0 zk​=0表示第 k k k類沒有被選中的機率。更數學化一點， z k z_k zk​要滿足以下兩個條件：

z k ∈ { 0 , 1 } z_k \in \{0,1\} zk​∈{0,1}

∑ K z k = 1 \sum_K z_k = 1 K∑​zk​=1

例如圖2中的例子，有兩類，則 z \boldsymbol{z} z的維數是2. 如果從第一類中取出一個點，則 z = ( 1 , 0 ) \boldsymbol{z} = (1, 0) z=(1,0)；，如果從第二類中取出一個點，則 z = ( 0 , 1 ) \boldsymbol{z} = (0, 1) z=(0,1).

z k = 1 z_k=1 zk​=1的機率就是 π k \pi_k πk​，假設 z k z_k zk​之間是獨立同分布的（iid），我們可以寫出 z \boldsymbol{z} z的聯合機率分布形式，就是連乘：

p ( z ) = p ( z 1 ) p ( z 2 ) . . . p ( z K ) = ∏ k = 1 K π k z k (2) p(\boldsymbol{z}) = p(z_1) p(z_2)...p(z_K) = \prod_{k=1}^K \pi_k^{z_k} \tag{2} p(z)=p(z1​)p(z2​)...p(zK​)=k=1∏K​πkzk​​(2)

因為 z k z_k zk​隻能取0或1，且 z \boldsymbol{z} z中隻能有一個 z k z_k zk​為1而其它 z j ( j ≠ k ) z_j (j\neq k) zj​(j​=k)全為0，是以上式是成立的。

圖2中的資料可以分為兩類，顯然，每一類中的資料都是服從正态分布的。這個叙述可以用條件機率來表示：

p ( x ∣ z k = 1 ) = N ( x ∣ μ k , Σ k ) p(\boldsymbol{x}|z_k = 1) = \mathcal{N}(\boldsymbol{x}|\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) p(x∣zk​=1)=N(x∣μk​,Σk​)

即第 k k k類中的資料服從正态分布。進而上式有可以寫成如下形式：

p ( x ∣ z ) = ∏ k = 1 K N ( x ∣ μ k , Σ k ) z k (3) p(\boldsymbol{x}| \boldsymbol{z}) = \prod_{k=1}^K \mathcal{N}(\boldsymbol{x}|\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)^{z_k} \tag{3} p(x∣z)=k=1∏K​N(x∣μk​,Σk​)zk​(3)

上面分别給出了 p ( z ) p(\boldsymbol{z}) p(z)和 p ( x ∣ z ) p(\boldsymbol{x}| \boldsymbol{z}) p(x∣z)的形式，根據條件機率公式，可以求出 p ( x ) p(\boldsymbol{x}) p(x)的形式：

p ( x ) = ∑ z p ( z ) p ( x ∣ z ) = ∑ z ( ∏ k = 1 K π k z k N ( x ∣ μ k , Σ k ) z k ) = ∑ k = 1 K π k N ( x ∣ μ k , Σ k ) (4) \begin{aligned} p(\boldsymbol{x}) &= \sum_{\boldsymbol{z}} p(\boldsymbol{z}) p(\boldsymbol{x}| \boldsymbol{z}) \\ &= \sum_{\boldsymbol{z}} \left(\prod_{k=1}^K \pi_k^{z_k} \mathcal{N}(\boldsymbol{x}|\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)^{z_k} \right ) \\ &= \sum_{k=1}^K \pi_k \mathcal{N}(\boldsymbol{x}|\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \end{aligned} \tag{4} p(x)​=z∑​p(z)p(x∣z)=z∑​(k=1∏K​πkzk​​N(x∣μk​,Σk​)zk​)=k=1∑K​πk​N(x∣μk​,Σk​)​(4)

（注：上式第二個等号，對 z \boldsymbol{z} z求和，實際上就是 ∑ k = 1 K \sum_{k=1}^K ∑k=1K​。又因為對某個 k k k，隻要 i ≠ k i \ne k i​=k，則有 z i = 0 z_i = 0 zi​=0，是以 z k = 0 z_k=0 zk​=0的項為1，可省略，最終得到第三個等号）

可以看到GMM模型的(1)式與(4)式有一樣的形式，且(4)式中引入了一個新的變量 z \boldsymbol{z} z，通常稱為隐含變量（latent variable）。對于圖2中的資料，『隐含』的意義是：我們知道資料可以分成兩類，但是随機抽取一個資料點，我們不知道這個資料點屬于第一類還是第二類，它的歸屬我們觀察不到，是以引入一個隐含變量 z \boldsymbol{z} z來描述這個現象。

注意到在貝葉斯的思想下， p ( z ) p(\boldsymbol{z}) p(z)是先驗機率， p ( x ∣ z ) p(\boldsymbol{x}| \boldsymbol{z}) p(x∣z)是似然機率，很自然我們會想到求出後驗機率 p ( z ∣ x ) p(\boldsymbol{z}| \boldsymbol{x}) p(z∣x)：

γ ( z k ) = p ( z k = 1 ∣ x ) = p ( z k = 1 ) p ( x ∣ z k = 1 ) p ( x , z k = 1 ) = p ( z k = 1 ) p ( x ∣ z k = 1 ) ∑ j = 1 K p ( z j = 1 ) p ( x ∣ z j = 1 )  (全機率公式) = π k N ( x ∣ μ k , Σ k ) ∑ j = 1 K π j N ( x ∣ μ j , Σ j )  (結合(3)、(4)) (5) \begin{aligned} \gamma(z_k) &= p(z_k = 1| \boldsymbol{x}) \\ &= \frac{p(z_k = 1) p(\boldsymbol{x}|z_k = 1)}{p(\boldsymbol{x}, z_k = 1)} \\ &= \frac{p(z_k = 1) p(\boldsymbol{x}|z_k = 1)}{\sum_{j=1}^K p(z_j = 1) p(\boldsymbol{x}|z_j = 1)} \text{ (全機率公式)}\\ &= \frac{\pi_k \mathcal{N}(\boldsymbol{x|\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k})}{\sum_{j=1}^K \pi_j \mathcal{N}(\boldsymbol{x|\boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j})} \text{ (結合(3)、(4))} \end{aligned} \tag{5} γ(zk​)​=p(zk​=1∣x)=p(x,zk​=1)p(zk​=1)p(x∣zk​=1)​=∑j=1K​p(zj​=1)p(x∣zj​=1)p(zk​=1)p(x∣zk​=1)​ (全機率公式)=∑j=1K​πj​N(x∣μj​,Σj​)πk​N(x∣μk​,Σk​)​ (結合(3)、(4))​(5)

（第2行，貝葉斯定理。關于這一行的分母，很多人有疑問，應該是 p ( x , z k = 1 ) p(\boldsymbol{x}, z_k = 1) p(x,zk​=1)還是 p ( x ) p(\boldsymbol{x}) p(x)，按照正常寫法，應該是 p ( x ) p(\boldsymbol{x}) p(x)。但是為了強調 z k z_k zk​的取值，有的書會寫成 p ( x , z k = 1 ) p(\boldsymbol{x}, z_k = 1) p(x,zk​=1)，比如李航的《統計學習方法》，這裡就約定 p ( x ) p(\boldsymbol{x}) p(x)與 p ( x , z k = 1 ) p(\boldsymbol{x}, z_k = 1) p(x,zk​=1)是等同的）

上式中我們定義符号 γ ( z k ) \gamma(z_k) γ(zk​)來表示來表示第 k k k個分量的後驗機率。在貝葉斯的觀點下， π k \pi_k πk​可視為 z k = 1 z_k=1 zk​=1的先驗機率。

上述内容改寫了GMM的形式，并引入了隐含變量 z \boldsymbol{z} z和已知 x {\boldsymbol{x}} x後的的後驗機率 γ ( z k ) \gamma(z_k) γ(zk​)，這樣做是為了友善使用EM算法來估計GMM的參數。

EM算法估計GMM參數

EM算法（Expectation-Maximization algorithm）分兩步，第一步先求出要估計參數的粗略值，第二步使用第一步的值最大化似然函數。是以要先求出GMM的似然函數。

假設 x = { x 1 , x 2 , . . . , x N } \boldsymbol{x} = \{\boldsymbol{x}_1, \boldsymbol{x}_2, ..., \boldsymbol{x}_N\} x={x1​,x2​,...,xN​}，對于圖2， x \boldsymbol{x} x是圖中所有點（每個點有在二維平面上有兩個坐标，是二維向量，是以 x 1 , x 2 \boldsymbol{x}_1, \boldsymbol{x}_2 x1​,x2​等都用粗體表示）。GMM的機率模型如(1)式所示。GMM模型中有三個參數需要估計，分别是 π \boldsymbol{\pi} π， μ \boldsymbol{\mu} μ和 Σ \boldsymbol{\Sigma} Σ. 将(1)式稍微改寫一下：

p ( x ∣ π , μ , Σ ) = ∑ k = 1 K π k N ( x ∣ μ k , Σ k ) (6) p(\boldsymbol{x}|\boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma}) = \sum_{k=1}^K\pi_k \mathcal{N}(\boldsymbol{x}|\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \tag{6} p(x∣π,μ,Σ)=k=1∑K​πk​N(x∣μk​,Σk​)(6)

為了估計這三個參數，需要分别求解出這三個參數的最大似然函數。先求解 μ k \mu_k μk​的最大似然函數。樣本符合iid，(6)式所有樣本連乘得到最大似然函數，對(6)式取對數得到對數似然函數，然後再對 μ k \boldsymbol{\mu_k} μk​求導并令導數為0即得到最大似然函數。

0 = − ∑ n = 1 N π k N ( x n ∣ μ k , Σ k ) ∑ j π j N ( x n ∣ μ j , Σ j ) Σ k − 1 ( x n − μ k ) (7) 0 = -\sum_{n=1}^N \frac{\pi_k \mathcal{N}(\boldsymbol{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)}{\sum_j \pi_j \mathcal{N}(\boldsymbol{x}_n|\boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)} \boldsymbol{\Sigma}_k^{-1} (\boldsymbol{x}_n - \boldsymbol{\mu}_k) \tag{7} 0=−n=1∑N​∑j​πj​N(xn​∣μj​,Σj​)πk​N(xn​∣μk​,Σk​)​Σk−1​(xn​−μk​)(7)

注意到上式中分數的一項的形式正好是(5)式後驗機率的形式。兩邊同乘 Σ k \boldsymbol{\Sigma}_k Σk​，重新整理可以得到：

μ k = 1 N k ∑ n = 1 N γ ( z n k ) x n (8) \boldsymbol{\mu}_k = \frac{1}{N_k}\sum_{n=1}^N \gamma(z_{nk}) \boldsymbol{x}_n \tag{8} μk​=Nk​1​n=1∑N​γ(znk​)xn​(8)

其中：

N k = ∑ n = 1 N γ ( z n k ) (9) N_k = \sum_{n=1}^N \gamma(z_{nk}) \tag{9} Nk​=n=1∑N​γ(znk​)(9)

(8)式和(9)式中， N N N表示點的數量。 γ ( z n k ) \gamma(z_{nk}) γ(znk​)表示點 n n n（ x n \boldsymbol{x}_n xn​）屬于聚類 k k k的後驗機率。則 N k N_k Nk​可以表示屬于第 k k k個聚類的點的數量。那麼 μ k \boldsymbol{\mu}_k μk​表示所有點的權重平均，每個點的權值是 ∑ n = 1 N γ ( z n k ) \sum_{n=1}^N \gamma(z_{nk}) ∑n=1N​γ(znk​)，跟第 k k k個聚類有關。

同理求 Σ k \boldsymbol{\Sigma}_k Σk​的最大似然函數，可以得到：

Σ k = 1 N k ∑ n = 1 N γ ( z n k ) ( x n − μ k ) ( x n − μ k ) T (10) \boldsymbol{\Sigma}_k = \frac{1}{N_k} \sum_{n=1}^N \gamma(z_{nk}) (x_n - \boldsymbol{\mu}_k)(x_n - \boldsymbol{\mu}_k)^T \tag{10} Σk​=Nk​1​n=1∑N​γ(znk​)(xn​−μk​)(xn​−μk​)T(10)

最後剩下 π k \pi_k πk​的最大似然函數。注意到 π k \pi_k πk​有限制條件 ∑ k = 1 K π k = 1 \sum_{k=1}^K\pi_k = 1 ∑k=1K​πk​=1，是以我們需要加入拉格朗日算子：

ln ⁡ p ( x ∣ π , μ , Σ ) + λ ( ∑ k = 1 K π k − 1 ) \ln p(\boldsymbol{x}|\boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma}) + \lambda(\sum_{k=1}^K \pi_k -1) lnp(x∣π,μ,Σ)+λ(k=1∑K​πk​−1)

求上式關于 π k \pi_k πk​的最大似然函數，得到：

0 = ∑ n = 1 N N ( x n ∣ μ k , Σ k ) ∑ j π j N ( x n ∣ μ j , Σ j ) + λ (11) 0 = \sum_{n=1}^N \frac{\mathcal{N}(\boldsymbol{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)}{\sum_j \pi_j \mathcal{N}(\boldsymbol{x}_n|\boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)} + \lambda \tag{11} 0=n=1∑N​∑j​πj​N(xn​∣μj​,Σj​)N(xn​∣μk​,Σk​)​+λ(11)

上式兩邊同乘 π k \pi_k πk​，我們可以做如下推導：

0 = ∑ n = 1 N π k N ( x n ∣ μ k , Σ k ) ∑ j π j N ( x n ∣ μ j , Σ j ) + λ π k (11.1) 0 = \sum_{n=1}^N \frac{\pi_k \mathcal{N}(\boldsymbol{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)}{\sum_j \pi_j \mathcal{N}(\boldsymbol{x}_n|\boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)} + \lambda \pi_k \tag{11.1} 0=n=1∑N​∑j​πj​N(xn​∣μj​,Σj​)πk​N(xn​∣μk​,Σk​)​+λπk​(11.1)

結合公式(5)、(9)，可以将上式改寫成：

0 = N k + λ π k (11.2) 0 = N_k + \lambda \pi_k \tag{11.2} 0=Nk​+λπk​(11.2)

注意到 ∑ k = 1 K π k = 1 \sum_{k=1}^K \pi_k = 1 ∑k=1K​πk​=1，上式兩邊同時對 k k k求和。此外 N k N_k Nk​表示屬于第個聚類的點的數量（公式(9)）。對 N k ， N_k， Nk​，從 k = 1 k=1 k=1到 k = K k=K k=K求和後，就是所有點的數量 N N N:

0 = ∑ k = 1 K N k + λ ∑ k = 1 K π k (11.3) 0 = \sum_{k=1}^K N_k + \lambda \sum_{k=1}^K \pi_k \tag{11.3} 0=k=1∑K​Nk​+λk=1∑K​πk​(11.3)

0 = N + λ (11.4) 0 = N + \lambda \tag{11.4} 0=N+λ(11.4)

進而可得到 λ = − N \lambda = -N λ=−N，帶入(11.2)，進而可以得到 π k \pi_k πk​更簡潔的表達式：

π k = N k N (12) \pi_k = \frac{N_k}{N} \tag{12} πk​=NNk​​(12)

EM算法估計GMM參數即最大化(8)，(10)和(12)。需要用到(5)，(8)，(10)和(12)四個公式。我們先指定 π \boldsymbol{\pi} π， μ \boldsymbol{\mu} μ和 Σ \boldsymbol{\Sigma} Σ的初始值，帶入(5)中計算出 γ ( z n k ) \gamma(z_{nk}) γ(znk​)，然後再将 γ ( z n k ) \gamma(z_{nk}) γ(znk​)帶入(8)，(10)和(12)，求得 π k \pi_k πk​， μ k \boldsymbol{\mu}_k μk​和 Σ k \boldsymbol{\Sigma}_k Σk​；接着用求得的 π k \pi_k πk​， μ k \boldsymbol{\mu}_k μk​和 Σ k \boldsymbol{\Sigma}_k Σk​再帶入(5)得到新的 γ ( z n k ) \gamma(z_{nk}) γ(znk​)，再将更新後的 γ ( z n k ) \gamma(z_{nk}) γ(znk​)帶入(8)，(10)和(12)，如此往複，直到算法收斂。

EM算法

1. 定義分量數目 K K K，對每個分量 k k k設定 π k \pi_k πk​， μ k \boldsymbol{\mu}_k μk​和 Σ k \boldsymbol{\Sigma}_k Σk​的初始值，然後計算(6)式的對數似然函數。

2. E step

   根據目前的 π k \pi_k πk​、 μ k \boldsymbol{\mu}_k μk​、 Σ k \boldsymbol{\Sigma}_k Σk​計算後驗機率 γ ( z n k ) \gamma(z_{nk}) γ(znk​)

   γ ( z n k ) = π k N ( x n ∣ μ n , Σ n ) ∑ j = 1 K π j N ( x n ∣ μ j , Σ j ) \gamma(z_{nk}) = \frac{\pi_k\mathcal{N}(\boldsymbol{x}_n| \boldsymbol{\mu}_n, \boldsymbol{\Sigma}_n)}{\sum_{j=1}^K \pi_j \mathcal{N}(\boldsymbol{x}_n| \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)} γ(znk​)=∑j=1K​πj​N(xn​∣μj​,Σj​)πk​N(xn​∣μn​,Σn​)​

3. M step

   根據E step中計算的 γ ( z n k ) \gamma(z_{nk}) γ(znk​)再計算新的 π k \pi_k πk​、 μ k \boldsymbol{\mu}_k μk​、 Σ k \boldsymbol{\Sigma}_k Σk​

   μ k n e w = 1 N k ∑ n = 1 N γ ( z n k ) x n Σ k n e w = 1 N k ∑ n = 1 N γ ( z n k ) ( x n − μ k n e w ) ( x n − μ k n e w ) T π k n e w = N k N \begin{aligned} \boldsymbol{\mu}_k^{new} &= \frac{1}{N_k} \sum_{n=1}^N \gamma(z_{nk}) \boldsymbol{x}_n \\ \boldsymbol{\Sigma}_k^{new} &= \frac{1}{N_k} \sum_{n=1}^N \gamma(z_{nk}) (\boldsymbol{x}_n - \boldsymbol{\mu}_k^{new}) (\boldsymbol{x}_n - \boldsymbol{\mu}_k^{new})^T \\ \pi_k^{new} &= \frac{N_k}{N} \end{aligned} μknew​Σknew​πknew​​=Nk​1​n=1∑N​γ(znk​)xn​=Nk​1​n=1∑N​γ(znk​)(xn​−μknew​)(xn​−μknew​)T=NNk​​​

   其中：

   N k = ∑ n = 1 N γ ( z n k ) N_k = \sum_{n=1}^N \gamma(z_{nk}) Nk​=n=1∑N​γ(znk​)

4. 計算(6)式的對數似然函數

   ln ⁡ p ( x ∣ π , μ , Σ ) = ∑ n = 1 N ln ⁡ { ∑ k = 1 K π k N ( x k ∣ μ k , Σ k ) } \ln p(\boldsymbol{x}|\boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma}) = \sum_{n=1}^N \ln \left\{\sum_{k=1}^K \pi_k \mathcal{N}(\boldsymbol{x}_k| \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)\right\} lnp(x∣π,μ,Σ)=n=1∑N​ln{k=1∑K​πk​N(xk​∣μk​,Σk​)}

5. 檢查參數是否收斂或對數似然函數是否收斂，若不收斂，則傳回第2步。

GMM聚類的可分性評價

使用GMM得到聚類結果後如何定量評價兩個類别的可分性呢？可以通過計算兩個或多個類别分布的重疊度來評價模型的可分性。這裡介紹一種高斯混合模型的重疊度計算方法：計算高斯混合模型的可分性和重疊度（Overlap Rate, OLR）。

Reference

1. 漫談 Clustering (3): Gaussian Mixture Model
2. Draw Gaussian distribution ellipse
3. Pang-Ning Tan 等, 資料挖掘導論（英文版）, 機械工業出版社, 2010
4. Christopher M. Bishop etc., Pattern Recognition and Machine Learning, Springer, 2006
5. 李航, 統計學習方法, 清華大學出版社, 2012