IA9b - K-means分群與EM算法: 實作

這一篇文章延續了上一篇IA9a的內容，但側重在python的實作，所以所有的節數、圖案以及公式編號都會依照IA9a的編號接序下來，所以重複的內容就直接引用，不在這裡重新贅述。

四、實作EM

4.1 高斯混合分佈

讓我們從1.1節的協方差矩陣開始，生成的圖形就如圖4 (腳本內加了seed random方便我們每次追蹤)，另外也看到了執行5次K-means算法後 (再多次也不會變了) 得到的分群結果如圖5，這顯然與我們的預期有極大的落差！
　但我們將按照3.3結給出EM算法的偽代碼，假設資料群是高斯混合分佈的形式，重新執行一遍分群，我們將看到EM算法能很好的分出我們預期的兩群的結果。我們預計會有兩群，所以\(k=1,2\)，且第一步需要初始化參數組\(\boldsymbol{\theta}=\{\mu_k,\Sigma_k,\pi_k\}\)。初始化的過程是給定一個起始點，原則上隨機選擇的方式是合理的，不過為了計算上方便，可以讓初始的協方差矩陣都為單位矩陣，\(\pi_k\)則為每個群均分 (這裡就是每個都是0.5)，如下代碼 
# Initializing parameters
sigma1,sigma2 = np.eye(2),np.eye(2)
pi1 = 0.5
pi2 = 1-pi1
N = len(cluster)
完成初始的參數選擇後，我們用偽代碼的第二步，也就是E-step計算置信度\(\gamma (z_{nk})\)。這裡我將Eq. (10)列出來方便解說 \[ \gamma (z_{nk})=\frac{\pi_k \mathcal{N}(\mathbf{x}|\mu_k,\Sigma_k)}{\sum_{k=1}^2 \pi_k \mathcal{N}(\mathbf{x}|\mu_k,\Sigma_k)}\tag{10}. \] 因為我們是假設兩群都是高斯分佈，所以對第一個群而言的分子就是 
# Starting E-step: Calculating responsibilities
# Numerator for k=1
k1 = pi1*multivariate_normal.pdf(cluster,p1,sigma1)
直接透過python腳本我們很容易理解了分子的意義便是計算每個資料點分別對應到\(\{\mu_1,\Sigma_1\}\)的高斯分佈的機率密度，所以愈靠近\(k=1\)群的點，它分子的值就愈大。接著計算對\(k=2\)以及分母並完整的寫下Eq. (10)的結果 
k2 = pi2*multivariate_normal.pdf(cluster,p2,sigma2)
gamma1 = k1/(k1+k2)
gamma2 = k2/(k1+k2)
不難理解Eq. (10)的分母就是將兩個分子結果相加即可，最後的gamma1和gamma2分別表示\(\gamma(z_{n1})\)和\(\gamma(z_{n2})\)。所以我們可以看下圖7，它是透過初始值所計算出的置信度。顏色代表它們有大於0.5的置信度是屬於同一群。

圖7. 初始條件下所計算的置信度，實線圓圈代表confidence level，內圈是\(1\sigma\)而外圈是\(2\sigma\)

　既然我們已經分出兩群了，我們變透過這兩群的點來計算M-step所需要的值，基本上在M-step裡就是要透過置信度來更新參數組\(\boldsymbol{\theta}\)，計算所需的公式是Eqs. (12-15)，實現的代碼如下 
# Starting M-step: Re-estimate parameters that maximizes expectation
# Re-evaluating N_k, Eq. (13) 
N1,N2 = np.sum(gamma1),np.sum(gamma2)

# Re-evaluating mean mu_k, Eq. (12)
p1 = np.sum(((cluster.T)*gamma1).T,axis=0)/N1
p2 = np.sum(((cluster.T)*gamma2).T,axis=0)/N2

# Re-evaluating pi_k, Eq. (15)
pi1 = N1/N
pi2 = N2/N
上面的部份更新了\(\{\mu_k,\pi_k\}\)，因這兩個計算比較簡單。而協方差矩陣的計算按照Eq. (14)是 \[ \Sigma_k = \frac{1}{N_k}\sum_{n=1}^N \gamma(z_{nk})(x_n-\mu_k)\otimes (x_n-\mu_k).\tag{14} \] 假使我要計算新的\(\Sigma_1\)，那麼它必須逐一計算\(\mathbf{x}\)內每個資料點對\(\mu_1\)的差外積後乘上該點對應的\(\gamma (z_{n1})\)，當計算完\(\mathbf{x}\)內所有的點後，我們會得到\(N\)個\(2\times 2\)的矩陣 (因為我們考慮的是2維的資料)，將這\(N\)個矩陣加總並除以\(N_1\)即得\(\Sigma_1\)，對\(\Sigma_2\)的情形一樣，代碼如下 
# Re-evaluating covariance matrix
sigma1,sigma2 = [],[]
for i in range(N):
    sigma1.append(gamma1[i]*np.outer((cluster[i]-p1),(cluster[i]-p1)))
    sigma2.append(gamma2[i]*np.outer((cluster[i]-p2),(cluster[i]-p2)))
sigma1 = np.sum(sigma1,axis=0)/N1
sigma2 = np.sum(sigma2,axis=0)/N2
OK, very well! 從代碼註解# Starting E-step...開始一路到上面結束我們變完成一次迭代啦，結果如下圖8。

圖8. EM算法經過一次迭代後的結果，為了排版乾淨我們只秀出兩個群對應到\(1\sigma\)所包含的區域

　不難發現我們的\(\mu_k\)已經透過E-step的置信度分群之後更新位置了，但距離正確的分群還有一段距離，所以我們可以加上for迴圈來重複迭帶個幾次，如下 
for iters in range(max_iters):
    ...
其中...就是上面E/M-steps交互計算的代碼部份，最大迭代次數由max_iters決定，當然在本文提供的腳本裏面，你會注意到我還提供了一個停止條件，若當\(\mu_k\)和上一次迭代的距離變化小於0.01時，不管有沒有達到max_iters的次數都會結束計算，這裡不贅述。而這份代碼經過23次迭帶後會碰到停止條件，我們把部份的結果羅列如下圖9。

圖9. 部份迭代後的結果，實線圓圈的區域表\(1\sigma\) CL

　在圖9裡秀出了每一次迭代過程中\(\mu_k\)走過的軌跡，以及該次迭代後修正的\(\Sigma_k\)對應到\(1\sigma\) CL區間，不像K-means (圖5和6) 最終錯誤的分群，EM很正確的分出來我們預期的結果。

4.2 Beyond Gaussians

在上一節中我們用高斯混合為例，發現了EM算法比傳統的K-means更正確的分出了結果，而事實上數學可以證明K-means就是將EM中置信度全取為1的特例，這種我們叫作hard K-means，而EM又叫作soft K-means (McKay 2003)，不過這只是歷史名詞，我們不必太拘泥。
　事實上，分佈\(p(x|\theta)\)不一定要是高斯，同理\(p(z|x)\)也不一定要是Eq. (6)的形式，在這種情況下我們就必須假定一個分佈\(p(x,z)\)並採用第三節的方案，而代碼也不能直接照抄4.1節的結果 (It's Gaussian mixtures ad hoc!)，腳本是死的，方法是活的，了解後便可以理出正確的思路，高斯以外的討論可以見Bishop 2008和Murphy 2012。

五、小結

在IA9裡我們看出EM算法將機率分佈的概念引入分群的方法中，使得它較hard K-means更有彈性，雖然它必須花費較多的計算效能，但卻能得到更接近我們預期的結果，並且這個算法在computer vision裡用來分割圖像或分離前後景也有著重要的應用 (Prince 2012)。但即使這裡引入了機率 (主要是每個資料點對應的pdf) 的計算，我們仍然離真正的隨機性算法八字沒一撇，EM仍是所謂的決定性算法，與IA8用Gibbs採樣去噪本質仍是不一樣的，我會在之後的文章繼續討論這個議題。啊對了，雖然在IA9a的引文裡提到要證明如何從高斯混合EM到K-means，不過就和費馬 (Pierre de Fermat) 說的一樣：空間不夠，而目前我是時間不夠😜，所以只能留待之後再補遺啦！


Source on github: IA9_EM.ipynb

留言

張貼留言

這個網誌中的熱門文章

物理雜記3 - Ising模型與模擬退火

圖片
Ising模型主要是用來描述物質的鐵磁性 (ferromagnetism)，是由德國物理學家Ernst Ising在1920左右提出了處理一維的解析解。該模型包含了每個晶格上原子的自旋 (spin) 參數，兩兩相鄰的原子彼此間可以是偏好處在共同的態上 (都是自旋向上)  或是處在相反的態上 (一上一下)，而外部的磁場會對該晶格上的自旋狀態施加能量，導致其偏向某一個方向。在這篇文章裡，我們不重複Ising當年複雜的數學解析解，而且該解也只能適用在1D的狀態，對於2D以上的系統目前並沒有解析解，所以我們將要另闢蹊徑，採用MCMC模擬的方式來決定一個系統終態每個晶格上自旋的排列情形。 一、Ising模型簡介 考慮一個由圖1所構成的晶格系統，晶格上的點都代表一顆原子，而黑白則表原子的自旋狀態，我們用黑表示向下而白表示向上，所以這個binary系統可以用更簡潔的\(\pm 1\)來表示每個原子可能的自旋 (即這個系統的樣本空間\(\Omega=\{+1,-1\}\))，另外為了方便，我們採用物理學的習慣將自旋表示為希臘字母\(\sigma_{n}\)，下標\(n\)可以想像成在 n 處原子的自旋態。原則上討論的都是一個由費米子 (fermions) 構成的系統，所以自旋應該用\(\pm \frac{1}{2}\)來表述，但在計算機上我們用\(\pm 1\)比較乾脆且簡潔。另外補充一下冷知識自然界已知的物質幾乎由費米子構成，僅有媒介作用力的粒子是玻色子 (bosons)，它們的差別在於費米子有半整數 (即不能約化成整數的分數) 自旋而玻色子則是整數自旋。 圖1. 晶格系統，黑白點分別代表該處原子的自旋狀態，黑表向下，白表向上。紅色的箭頭指出中心原子的4位鄰居\(m\) 　Ising模型假定了每個晶格上的原子的自旋都會受其鄰居 (neighbors \(\mathcal{N}\)) 的狀態影響，例如圖1中間白色的原子雖然他在這一刻是處在\(+1\)的狀態，但根據其週邊四個鄰居 (分別是3白1黑) 的況，它仍有一定的機會會在下一刻翻轉自旋的狀態。 　對一個給定位置 n 的原子，它和它的周圍鄰居構成的局部區域 (local area) 形成的一個系統的能量由Ising模型給出 \[ E_n = -\sum_{m\in \mathcal{N}} J
閱讀完整內容

PS7 - 小專題: 為何高斯分佈擁有最大熵?

圖片
高爾頓板，或稱做quincunx，是由Galton爵士 (1822-1911) 所發明用來驗證中央極限定理 (CLT) 的裝置。它也是第一位討論均值回歸 (regression to the mean) 現象的科學家，但同時亦是優生學 (eugenics) 這樣錯誤觀念的推廣者 在 上篇雜記 的結尾我們提到了對於一個趨於穩態的系統擁有最大的熵，而這也是統計力學的直接應用，用物理的角度來講便是當系統平衡時 Gibbs 自由能不再變化 (\(\Delta G =0\))，而在此時會對應到系統的熵有全域極大值 (global maximum)。故大自然很多現象都已經處於穩態的狀況下，我們對這個現象的母體做抽樣，所得到的觀測樣本應也有最大熵，則描述這些觀測樣本最好的分佈便是符合該條件下擁有最大熵的機率分佈。而指數族的分佈符合最大熵的特點，我們在這個小專題內便嘗試證明高斯與二項式分佈在連續和離散的隨機變量下的條件擁有最大熵的性質。 一、機率熵 在 上篇雜記 裡裡我們探討了夏儂熵，或者稱作資訊量\(I\)的意含，所以一個系統內 (我們這裡就不再用系綜這種彆扭的歷史名詞了) 每個成員出現的機率是\(p_i\)，則對該系統的熵可以用夏儂熵公式來表示 \[ I(p) = -\sum_i p_i \ln p_i. \tag{7.1} \] 所以對某一組 連續 機率分佈函數\(p(x)\)，其中\(x\)是隨機變量，則將(7.1)的求和符號改成積分後有 \[ I(p) = -\int p(x)\ln p(x) dx, \tag{7.2} \] 稱為機率分佈\(p(x)\)的 機率熵 ，而積分的上下限則看隨機變量的範圍，例如高斯分佈就是\(\pm \infty\)。如果今天我們有兩組機率分佈\(p(x)\)和\(q(x)\)，我們便可以定義它們之間的 相對熵 ，或稱作 KL散度 (Kullback-Leibler divergence)，定義為 \[ D_{KL}(q,p)\equiv\int q(x)\ln \left(\frac{q(x)}{p(x)}\right) dx = I(q,p)-I(q)\geq 0, \tag{7.3} \] 是 正定的 (positively defined)，其中 \[ I(q,p)=-\int q(x)\ln p(x) d
閱讀完整內容

IA9a - K-means分群與EM算法: 理論

圖片
市場內分門別類的產品 原則上 K-means分群 是 EM算法 (expectation-maximization algorithm) 的一個特例，相較於EM從似然函數 (likelihood) 的方式出發並導出不同資料點對不同群 (或稱聚類、簇，cluster) 的貢獻度，K-means單純用歐氏距離來分類，雖然這加快了K-means算法的計算效率，但顯然的這麼naïve的方法會產生明顯的缺陷。但仍讓我們會先從簡單的K-means開始，並探討K-means算法的結果在某些情況下和人類肉眼辨識的結論會產生明顯的歧義，接著再推進到通用的EM算法，除了 python 的demo外，我也會證明在數學上如何從EM得到K-means的極限 (僅限高斯混合的情況)，而在這個part內，讓我們先從數學理論的角度來面。 一、K-means分群 今天如果我有一大坨資料共120個點如圖1所示，從資料的散佈趨勢來看，很明顯的分成兩群，那我們是不是可以有一個演算法來自動化將資料點如同我們臆測的一樣分成兩袋？一個最簡單的發法便是使用K-means分群，而這個大寫K表示群的數量，以圖1而言，我們猜測K=2，所以是個2-means分群問題。 圖1. 資料散佈圖 (普魯士藍) 與初始化的紅綠兩點 　所以為了達成分群，我們需要先初始化兩個點，這兩個點代表兩群的中心。事實上，初始化的位置不必要一定剛好在兩群的中心，可以是任意挑選的地方，但為了計算方便，我們會將起始化的兩點開始在靠近資料點的平均值附近，如同圖1中紅綠兩點。我們將這兩個初始化的點叫作\(p_k\)其中\(k=1,2\)表示第\(k\)群。另外將資料點的集合用\(\mathbf{D}=\{d_1,d_2\cdots d_n \cdots d_{120}\}\)表示，所以我們計算每個點\(d_n=(x_{d_n},y_{d_n})\)和\(p_k=(x_{p_k},y_{p_k})\)的歐氏距離 \[ \ell_{n,k}=\sqrt{(x_{d_n}-x_{p_k})^2+(y_{d_n}-y_{p_k})^2}. \tag{1} \] 若該點和\(p_1\)的距離小於和\(p_2\)的距離 (\(\ell_{n,1}<\ell_{n,2}\))，則我們說\(d_n\)比較靠近\(p_1\)且應該
閱讀完整內容