本文寫於2018.12.11，原文地址，技術博客地址

說明：造輪子是為了深入理解演算法原理和參數，實際工作中不用造輪子

01 物以類聚

今天我們學習並實踐KMeans聚類演算法，分成以下幾個部分，跟上節奏燥起來！

1. KMeans演算法理論和代碼實現
2. 改進，BiKMeans演算法理論和代碼實現
3. 實例，上車點規劃
4. 抉擇，如何挑選最佳的聚類簇數？

02 KMeans理論和演算法實現

聚類是一種無監督學習的方法，所謂「無監督」，就是指參與訓練的樣本沒有標籤。

KMeans聚類演算法過程如下： 1. 對於一組數據集，隨機選取k個點作為質心，將數據集中的點歸為離其最近的質心一簇，此時數據集被劃分為k個簇； 2. 對這k個簇，重新計算各簇的質心(均值)； 3. 根據新的質心，按照step1繼續聚類，然後再根據聚類重新計算各簇質心，直到質心不再改變，分類完成。

說白了，就是不斷地聚類、劃分的過程。

通過KMeans原理，可以看到幾個顯而易見的缺點： 1. 簇數量k由用戶指定，無法預先知道最佳k值 >>解法：分為幾簇，最終由輪廓係數S(i)決定，取輪廓係數最大的分類數(05節劇透) 2. 最終質心可能與初始點選擇有關 >> 因此KMeans的結果可能收斂到局部最小值，而不是全局最小值 >> 解法： - BiKMeans（03節） - KMeans++（KMeans++ 演算法在選擇初始質心時並不是隨機選擇，而是選擇盡量相互分離的質心，即，下一個質心點總是離上一個質心點較遠的點）

代碼實現

def loadDataSet(fileName):
dataList=[]
dataMat=[]
fr=open(fileName)
for line in fr.readlines():
curLine=line.strip().split( )
fltLine=list(map(float,curLine))
dataList.append(fltLine)
dataMat=mat(dataList)
return dataMat

def distEclud(vecA,vecB):
return sqrt(sum(power(vecA-vecB,2))) #歐式距離

#為輸入數據集構造k個隨機中心，中心位置在各特徵最大最小值之間
def randCent(dataSet,k):
n=shape(dataSet)[1]
center=mat(zeros((k,n)))
for j in range(n): #對每個特徵
minJ=min(dataSet[:,j])
rangeJ=float(max(dataSet[:,j])-minJ)
center[:,j]=mat(minJ+rangeJ*random.rand(k,1)) #質心第j維坐標在數據集第j維數據之間
return center

def KMeans(dataSet,k,distMeas=distEclud,createCent=randCent):
m=shape(dataSet)[0]
clusterAssment=mat(zeros((m,2))) #用於記錄各樣本當前歸屬於哪個簇以及到該簇質心的歐式距離平方
center=createCent(dataSet,k)
clusterChanged=True

while clusterChanged:
clusterChanged=False
#對每個樣本，計算樣本到各質心的距離，尋找距離最近的質心，將該樣本歸為該質心所在簇
for i in range(m):
minDist=inf;minIndex=-1
for j in range(k): #對每個質心,計算到i樣本的距離
distJI=distMeas(center[j,:],dataSet[i,:])
if distJI<minDist:
minDist=distJI;minIndex=j #i樣本暫屬於j簇，到j簇質心距離為minDist
if clusterAssment[i,0]!=minIndex:
clusterChanged=True #若任一樣本在本次迭代中改變了簇類，則要進行下一次迭代(即，直到任何樣本都不再改變簇類，聚類停止)
clusterAssment[i,:]=minIndex,minDist**2 #記錄樣本i的簇類情況
#print (center)
#更新質心
for cent in range(k):
ptsInClust=dataSet[nonzero(clusterAssment[:,0].A==cent)[0]] #篩選出屬於當前簇類的點
center[cent,:]=mean(ptsInClust,axis=0) #對該簇類各樣本的各列求均值，作為新質心
return center,clusterAssment

用一組數據來測試一下：

dataMat1=loadDataSet(rD:DMpythondataMLiA_SourceCodemachinelearninginactionCh10 estSet.txt)
dataMat2=loadDataSet(rD:DMpythondataMLiA_SourceCodemachinelearninginactionCh10 estSet2.txt)
center_testSet1,clusterAssment_testSet1=KMeans(dataMat1,4)
center_testSet2,clusterAssment_testSet2=KMeans(dataMat2,3)

plt.figure(figsize=(6,6))
plt.scatter(dataMat1[:,0].T.tolist()[0],dataMat1[:,1].T.tolist()[0],c=pink,s=30)
plt.scatter(center_testSet1.T[0].tolist()[0],center_testSet1.T[1].tolist()[0],c=blue,s=50)

結果如下，藍色點為聚類質心

03 BiKMeans理論和演算法實現

剛才我們實現了KMeans演算法，不過也提到，KMeans有個缺點，就是其聚類的最終質心可能與初始點選擇有關，因此KMeans的結果可能收斂到局部最小值，而不是全局最小值。

怎麼解決呢？提示：想想決策樹是如何分支的？

為了得到全局最優解，BiKMeans演算法出現了，它的過程如下（是不是跟決策樹分支有點神似？） 1. 將整個數據集看作一個簇，計算初始質心，即所有數據點各特徵的均值 2. 遍歷各質心，對各質心，將質心所在簇用原始KMeans演算法二分，計算二分後整個數據集的SSE(即平方誤差和，即簇各點到簇質心距離平方和)，找到二分後整體數據集SSE最小的質心，認為此質心是本次劃分的最佳質心，對其進行二分 3. 不斷重複step2，直到質心總數=設置的k

說白了，BiKMeans演算法過程類似於決策樹的分支，通過啟發的方法，每次迭代只分裂當前最佳質心，直到簇數量達到k，這樣的方法可以保證最終劃分得到的質心是全局最優解，而原始KMeans可能會陷入局部最優解。

代碼實現

def biKMeans(dataSet,k,distMeas=distEclud):
m=shape(dataSet)[0]
clusterAssment=mat(zeros((m,2))) #記錄各樣本歸屬和距離平方
center0=mean(dataSet,axis=0).tolist()[0] #初始質心
centerList=[center0] #用於記錄聚類質心坐標
for j in range(m):
clusterAssment[j,1]=distMeas(mat(center0),dataSet[j,:])**2
while len(centerList)<k:
lowestSSE=inf #SSE=Sum of Square Error
#對每個簇，嘗試二分，計算二分後整體數據的SSE，若小於lowestSSE，則將簇二分，如此往複，直到分到k個簇為止
for i in range(len(centerList)):
ptsInCluster=dataSet[nonzero(clusterAssment[:,0].A==i)[0],:] #默認質心索引就是數據對應簇類
centerMat,splitClustAss=KMeans(ptsInCluster,2,distMeas)
sseSplit=sum(splitClustAss[:,1]) #被二分後的簇的平方誤差和
sseNotSplit=sum(clusterAssment[nonzero(clusterAssment[:,0]!=i)[0],1]) #整體數據中，未被二分的簇的平方誤差和
if (sseSplit+sseNotSplit)<lowestSSE:
bestCentToSplit=i
bestNewCents=centerMat
bestClustAss=splitClustAss.copy() #.copy()防止splitClustAss被覆蓋時影響到bestClustAss
lowestSSE=sseSplit+sseNotSplit
#確定好本次迭代被二分的簇後，將被二分的數據對應的簇類更新
bestClustAss[nonzero(bestClustAss[:,0].A!=0)[0],0]=len(centerList) #更新先後順序很重要！
bestClustAss[nonzero(bestClustAss[:,0].A==0)[0],0]=bestCentToSplit
print(The bestCentToSplit is:,bestCentToSplit)
print(The number of samples to split is,len(bestClustAss))
centerList[bestCentToSplit]=bestNewCents[0,:].tolist()[0] #將被二分的原簇質心替換為二分後的質心之一
centerList.append(bestNewCents[1,:].tolist()[0]) #將二分後的另一質心添加在質心列表末尾
#將二分後的簇的數據歸屬更新到總記錄中
clusterAssment[nonzero(clusterAssment[:,0].A==bestCentToSplit)[0],:]=bestClustAss
return mat(centerList),clusterAssment

測試

center_testSet22,clusterAssment_testSet22=biKMeans(dataMat2,3)

plt.figure(figsize=(6,6))
plt.scatter(dataMat2[:,0].T.tolist()[0],dataMat2[:,1].T.tolist()[0],c=pink,s=30)
plt.scatter(center_testSet22.T[0].tolist()[0],center_testSet22.T[1].tolist()[0],c=blue,s=50)

04 實例：上車點規劃

我們在前節實現了KMeans演算法和BiKMeans演算法，現在讓我們來用用這兩個輪子吧。

設想，你邀請了70個朋友參加你的party，他們住在城市的各個地方，你需要在party開始前2小時開車去接他們(假設你的車子可以容納70人，哈哈哈)。

為了時間效率最大化，請問你該如何規劃每個人的上車點呢？總不能你一個一個去家門口接吧。

我們用聚類的方法來規劃！

現在我們獲取了每個朋友住址的經緯度，那麼請開始你的表演。

很簡單，只需要調用我們早好的輪子BiKMeans

#球面距離計算函數：球面餘弦距離(向量夾角*地球半徑)
#求球面上兩向量vecA,vecB的距離
def distSLC(vecA,vecB):
a=sin(vecA[0,1]*pi/180)*sin(vecB[0,1]*pi/180) #由於抓取的經緯度為角度，需要通過 *pi/180來轉換為弧度
b=cos(vecA[0,1]*pi/180)*cos(vecB[0,1]*pi/180)*cos(pi*(vecA[0,0]-vecB[0,0])/180)
return 6371.0*arccos(a+b) #6371為地球半徑，單位為英尺

def clusterClubs(k=5):
dataList=[]
for line in open(rD:DMpythondataMLiA_SourceCodemachinelearninginactionCh10places.txt).readlines():
lineArr=line.strip().split( )
dataList.append([float(lineArr[4]),float(lineArr[3])]) #讀取地址的經緯度
dataMat=mat(dataList)
center,clusterAss=biKMeans(dataMat,k,distMeas=distSLC) #將地址按經緯度聚類
#作圖
fig=plt.figure(figsize=(10,8))
rect=[0.1,0.1,0.8,0.8] #用於設置坐標軸刻度，[]中前兩個值表示左邊起始位置，後兩個值對應坐標長度
scatterMarkers=[^,o,h,8,p,d,v,s,>,<] #用於設置散點圖點的形狀
axprops=dict(xticks=[],yticks=[])
ax0=fig.add_axes(rect,label=ax0,**axprops)

#讀圖，並將圖片顯示在設定好的坐標軸中
imgP=plt.imread(rD:DMpythondataMLiA_SourceCodemachinelearninginactionCh10Portland.png)
ax0.imshow(imgP)
#將地址按聚類結果作散點圖
ax1=fig.add_axes(rect,label=ax1,frameon=False)
for i in range(k):
ptsInCluster=dataMat[nonzero(clusterAss[:,0].A==i)[0],:]
markerStyle=scatterMarkers[i]
ax1.scatter(ptsInCluster[:,0].flatten().A[0],ptsInCluster[:,1].flatten().A[0],
marker=markerStyle,s=90)
#標出質心
ax1.scatter(center[:,0].flatten().A[0],center[:,1].flatten().A[0],marker=+,s=300)
plt.show()

好了，我們寫的這個「上車點規劃器」就可以使用了，只需要輸入聚類簇數量即可 ，下面看看幾個測試結果：

但是到底聚為幾類比較合理呢？這個合理是指即不讓你的朋友跑太遠去坐車，也不用你跑太多地方去接人。

05 如何挑選最佳的聚類簇數

不論是KMeans演算法還是BiKMeans演算法，都仍有一個缺點沒有解決：簇數量k由用戶指定，用戶指定的k不一定是最佳的簇數量。這也是上一節我們沒有解決的問題。

怎麼辦呢？

使用輪廓係數，最佳k值由輪廓係數S(i)決定，取輪廓係數最大的分類數。

什麼是輪廓係數，輪廓係數S(i)用于衡量聚類效果，取值在-1~1之間，越接近1表示聚類效果越好，公式如下，

其中a(i)=i點到同簇各點距離均值；b(i)=min(i點到某個非同簇各點均值)，即i點到其他簇質心距離的最小值，整體數據集的輪廓係數是各點輪廓係數的平均值

基於此，我們可以寫一個計算聚類結果輪廓係數的函數，然後看看輪廓係數與聚類數量k的關係，從而找到最佳的聚類數量k。

輪廓係數計算函數

def outlineOfCluster(filename,maxClusterNum,distMeas=distEclud):
dataList=[]
for line in open(filename).readlines():
lineArr=line.strip().split( )
dataList.append([float(lineArr[4]),float(lineArr[3])]) #讀取地址的經緯度
dataMat=mat(dataList)
#遍歷2~nn個質心，求不同數量的簇的輪廓係數，設置輪廓係數最大的簇數量為k值
sk={}
for k in range(2,maxClusterNum+1):
center,clusterAss=biKMeans(dataMat,k,distMeas)
outline=[] #代表聚類為k個簇時各點的輪廓係數列表
for i in range(k): #遍歷各簇i
ptsInCluster=clusterAss[nonzero(clusterAss[:,0].A==i)[0],:]
for j in range(len(ptsInCluster)): #遍歷i簇的各點j
ptsInClusterNotJ=vstack([ptsInCluster[:j,:],ptsInCluster[j+1:,:]])
ajn=[]
for n in range(len(ptsInClusterNotJ)): #遍歷i簇非j的點
ajn.append(distSLC(ptsInCluster[n],ptsInCluster[j]))
aj=nanmean(ajn)
bjm=[]
centerWithoutI=vstack([center[:i,:],center[i+1:,:]])
for m in range(len(centerWithoutI)): #遍歷非i簇質心
bjm.append(distSLC(centerWithoutI[m],ptsInCluster[j]))
bj=min(bjm)
sj=(bj-aj)/max(bj,aj) #i簇中j點的輪廓係數
outline.append(sj) #將i簇中各點的輪廓係數保存在outline中，outline用於存儲聚類為k個簇時各點的輪廓係數，在遍歷k時需要重置
sk[k]=nanmean(outline) #聚類為k個簇時的輪廓係數是各點輪廓係數的均值
return sk

現在我們可以接近上車點規劃問題的遺留問題了，規劃幾個上車點最合理？

sk=outlineOfCluster(rD:DMpythondataMLiA_SourceCodemachinelearninginactionCh10places.txt,35,distMeas=distSLC)

plt.figure(figsize=(10,5))
plt.plot(list(sk.keys()),list(sk.values()),linewidth=3)
plt.title(聚類簇數-效果,fontsize=18,fontweight=bold)
plt.xlabel(最終聚類簇數量,fontsize=14)
plt.ylabel(輪廓係數,fontsize=14)

可以看到，聚類簇數在達到10個簇之後，輪廓係數的增量就變得非常小了，因此從性價比方面考慮，選擇聚類為10個簇的性價比是最高的，即規劃10個上車點就好。

06 總結

本文實踐了KMeans聚類演算法，分成以下幾個部分，

1. KMeans演算法理論和代碼實現
2. 改進，BiKMeans演算法理論和代碼實現
3. 實例，上車點規劃
4. 抉擇，如何挑選最佳的聚類簇數？

下期我們將實踐Apriori演算法，這是一種挖掘關聯規則得到頻繁項集的演算法，比如可以預測你買了A商品後會買B商品，敬請期待~

07 參考

《機器學習實戰》 Peter Harrington Chapter10

如果覺得本文對你有幫助，點贊讓更多的同學看到吧

推薦閱讀：

相关文章