這篇部落格整理K均值聚類的内容，包括：
1、K均值聚類的原理；
2、初始類中心的選擇和類别數K的确定；
3、K均值聚類和EM算法、高斯混合模型的關系。
一、K均值聚類的原理
K均值聚類（K-means）是一種基于中心的聚類算法，通過疊代，将樣本分到K個類中，使得每個樣本與其所屬類的中心或均值的距離之和最小。
1、定義損失函數
假設我們有一個資料集｛x1, x2,..., xN｝，每個樣本的特征次元是m維，我們的目标是将資料集劃分為K個類别。假定K的值已經給定，那麼第k個類别的中心定義為μk，k=1,2,..., K，μk是一個m維的特征向量。我們需要找到每個樣本所屬的類别，以及一組向量｛μk｝，使得每個樣本與它所屬的類别的中心μk的距離平方和最小。
首先，這個距離是什麼距離呢？聚類需要根據樣本之間的相似度，對樣本集合進行劃分，将相似度較高的樣本歸為一類。度量樣本之間相似度的方法包括計算樣本之間的歐氏距離、馬氏距離、餘弦距離或相關系數，而K均值聚類是用歐氏距離的平方來度量樣本之間的相似度。歐式距離的平方公式如下：
其中rnk∈{0,1}，n=1,2,...,N，k=1,2,...,K，如果rnk=1，那麼表示樣本xn屬于第k類，且對于j≠k，有rnj=0，也就是樣本xn隻能屬于一個類别。
于是我們需要找到｛rnk｝和｛μk｝的值，使得距離平方之和J最小化。
2、進行疊代
K均值聚類的算法是一種疊代算法，每次疊代涉及到兩個連續的步驟，分别對應rnk的最優化和μk的最優化，也對應着EM算法的E步（求期望）和M步（求極大）兩步。
首先，為μk選擇一些初始值，也就是選擇K個類中心。然後第一步：保持μk固定，選擇rnk來最小化J，也就是把樣本指派到與其最近的中心所屬的類中，得到一個聚類結果；第二步：保持rnk固定，計算μk來最小化J，也就是更新每個類别的中心。不斷重複這兩個步驟直到收斂。
具體來說：
E步：在類中心μk已經确定的情況下，最優化rnk
這一步比較簡單，我們可以對每個樣本xn獨立地進行最優化。将某個樣本配置設定到第k個類别，如果這個樣本和第k個類别的距離最小，那麼令rnk=1。對N個樣本都這樣進行配置設定，自然就得到了使所有樣本與類中心的距離平方和最小的｛rnk｝，進而得到了一個聚類結果。
M步：确定了資料集的一種劃分，也就是｛rnk｝确定後，最優化μk
目标函數J是μk的一個二次函數，令它關于μk的導數為零，就可以使目标函數達到最小值，即
這個μk就是類别k中所有樣本的均值，是以把這個算法稱為K均值聚類。
重新為樣本配置設定類别，再重新計算每個類别的均值，不斷重複這兩個步驟，直到聚類的結果不再改變。
需要注意以下幾點：
①K均值聚類算法可能收斂到目标函數J的局部極小值，不能保證收斂到全局最小值；
②在聚類之前，需要對資料集進行标準化，使得每個樣本的均值為0，标準差為1；
③初始類中心的選擇會直接影響聚類結果，選擇不同的初始類中心，可能會得到不同的聚類結果。
④K均值聚類算法的複雜度是O(mnk)，m是樣本的特征次元，n是樣本個數，k是類别個數。
二、初始類中心的選擇和類别數K的确定
K均值聚類算法的思想比較簡單，不涉及到什麼數學知識，關鍵點在于初始類中心的選擇和類别數K的确定，這對聚類的結果有比較大的影響。
（一）初始類中心的選擇
1、第一種方法
用層次聚類算法進行初始聚類，然後用這些類别的中心作為K均值聚類的初始類中心。層次聚類的複雜度為O(mn3)，m是樣本的特征次元，n是樣本個數，複雜度也是蠻高的，那為什麼用層次聚類的結果作為初始類呢？我想是因為層次聚類的結果完全是由算法确定的，完全沒有人工的幹預，是一個客觀的結果，這樣就把K均值聚類的初始類選擇問題，由主觀确定變成了客觀決定。
2、第二種方法
首先随機選擇一個點作為第一個初始類中心點，然後計算該點與其他所有樣本點的距離，選擇距離最遠的點作為第二個初始類的中心點，以此類推，直到選出K個初始類中心點。
（二）類别數K的确定
1、輪廓系數
輪廓系數（Silhouette Coefficient）可以用來判定聚類結果的好壞，也可以用來确定類别數K。好的聚類要保證類别内部樣本之間的距離盡可能小（密集度），而類與類之間樣本的距離盡可能大（離散度），輪廓系數就是一個用來度量類的密集度和離散度的綜合名額。
輪廓系數的計算過程和使用如下：
①計算樣本xi到同類Ck其他樣本的平均距離ai，将ai稱為樣本xi的簇内不相似度，ai越小，說明樣本xi越應該被配置設定到該類。
②計算樣本xi到其他類Cj所有樣本的平均距離bij，j=1, 2 ,..., K，j≠k，稱為樣本xi與類别Cj的不相似度。定義樣本xi的簇間不相似度：bi =min{bi1, bi2, ...,bij,..., biK}，j≠k，bi越大，說明樣本xi越不屬于其他簇。
③根據樣本xi的簇内不相似度ai和簇間不相似度bi，定義樣本xi的輪廓系數si，作為樣本xi分類結果的合理性的度量。
④輪廓系數範圍在[-1,1]之間，該值越大，聚類結果越好。si接近1，則樣本xi被配置設定到類别Ck的結果比較合理；si接近0，說明樣本xi在兩個類的邊界上；si接近-1，說明樣本xi更應該被配置設定到其他類别。
⑤計算所有樣本的輪廓系數si的均值，得到聚類結果的輪廓系數S，作為聚類結果合理性的度量。輪廓系數越大，聚類結果越好。
⑥使用不同的K值進行K均值聚類，計算各自聚類結果的輪廓系數S，選擇較大的輪廓系數所對應的K值。
2、肘部法則
K均值聚類的損失函數是所有樣本到類别中心的距離平方和J，也就是誤差平方和SSE：
從類别數K=1開始，一開始随着類别數的增大，且類别數K小于真實的類别數時，SSE會快速地下降。而當類别數到達真實類别數的臨界點後，SSE開始緩慢下降，也就是說SSE和K的關系是一個肘部形狀，這個肘部所對應的K值可以認為是合适的類别數。
那麼我們選擇不同的K值，用K均值聚類訓練多個模型，然後計算模型的SSE，選擇SSE開始緩慢下降時的K值作為聚類的類别數。如下圖，可以選擇4或5作為聚類的類别數。
三、K均值聚類與高斯混合模型的關系
關于K均值聚類與EM算法、高斯混合模型的關系，主要有以下三點：
1、K均值聚類是一種非機率的聚類算法，屬于硬聚類方法，也就是一個樣本隻能屬于一個類（類與類之間的交集為空）。相比之下，高斯混合模型（GMM）是一種基于機率的聚類算法，屬于軟聚類方法，每個樣本按照一個機率分布，屬于多個類。
2、K均值聚類在一次疊代中的兩個步驟，可以看做是EM算法的E步和M步，而且K均值聚類可以看做是用EM算法對⾼斯混合模型進行參數估計的⼀個特例，也就是高斯混合模型中分模型的方差σ2相等，為常數，且σ2→0時的極限情況。
3、K均值聚類和基于EM算法的高斯混合模型，對參數的初始化值比較敏感，由于K均值聚類的計算量遠小于基于EM算法的高斯混合模型，是以通常運⾏K均值算法找到⾼斯混合模型的⼀個初始化值，再使用EM算法進行調節。具體而言，用K均值聚類劃分的K個類别中，各類别中樣本所占的比例，來初始化K個分模型的權重；用各類别中樣本的均值來初始化K個高斯分布的期望；用各類别中樣本的方差來初始化K個高斯分布的方差。
（一）從圖形來了解
為了了解以上這幾點，尤其是第2點，我們可以先從圖形來看。假設高斯混合模型由4個高斯分布混合而成，高斯分布的密度函數如下。這裡和《聚類之高斯混合模型與EM算法》的符号表示一緻，y為樣本。
（二）從公式來了解
用EM算法對高斯混合模型進行極大似然估計，在E步，我們需要基于第i輪疊代的參數θ(i)=(αk, μk, σk)來計算γjk，γjk是第j個樣本yj來自于第k個高斯分布分模型的機率，k=1,2,...,K。在高斯混合模型中，γj是一個K維的向量，也就是第j個樣本屬于K個類的機率。假設分模型的方差σk2都相等，且是一個常數，不需要再估計，那麼在EM算法的E步我們計算γjk
也就是樣本yj屬于第k類的機率近似為1，屬于其他類别的機率近似為0，也就成為了一種硬聚類，也就是K均值聚類。
其實在σ2→0時的極限情況下，最大化高斯混合模型的完全資料的對數似然函數的期望，等價于最小化K均值聚類的目标函數J。
比如有4個高斯分布，樣本yj的γj為[0.55, 0.15, 0.2, 0.1]，那麼屬于第1類的機率γj1最大。而當分模型的方差σ2→0時，樣本yj的γj可能為[0.98, 0.01, 0.005, 0.005]，也就是該樣本直接被配置設定到了第1類，成為了硬聚類。
附：4個高斯分布的機率密度函數的圖形代碼
import matplotlib.pyplot as plt
import math
import numpy as np
# 均值
u1 = 0
u2 = 2
u3 = 4
u4 = 6
# 标準差
sig1 = math.sqrt(1)
sig2 = math.sqrt(2)
sig3 = math.sqrt(3)
sig4 = math.sqrt(1)
def x(u,sig):
return np.linspace(u - 5*sig, u + 5*sig, 100)
x1 = x(u1,sig1)
x2 = x(u2,sig2)
x3 = x(u3,sig3)
x4 = x(u4,sig4)
# 機率密度
def y(x,u,sig):
return np.exp(-(x - u) ** 2 /(2* sig**2))/(math.sqrt(2*math.pi)*sig)
y1 = y(x1,u1,sig1)
y2 = y(x2,u2,sig2)
y3 = y(x3,u3,sig3)
y4 = y(x4,u4,sig4)
plt.plot(x1,y1, "r-")
plt.plot(x2, y2, "g-")
plt.plot(x3, y3, "b-")
plt.plot(x4, y4, "m-")
plt.xticks(range(-6,16,2))
plt.show() 參考資料：
1、李航：《統計學習方法》（第二版）
2、《Pattern Recognition and Machine Learning》
3、https://www.jianshu.com/p/335b376174d4