01 物以類聚
經(jīng)過(guò)半年的不懈努力，我們已經(jīng)學(xué)習(xí)并實(shí)踐了經(jīng)典的分類算法和經(jīng)典的回歸算法耕赘，下面我們開始學(xué)習(xí)經(jīng)典的聚類算法(興奮~~~)
目前打算對(duì)三種聚類算法進(jìn)行學(xué)習(xí)和代碼實(shí)操(俗稱“造輪子”)：
- KMeans
- Apriori
- FP-Growth
今天我們學(xué)習(xí)并實(shí)踐KMeans聚類算法楞黄，分成以下幾個(gè)部分又谋，跟上節(jié)奏燥起來(lái)橘荠！
- KMeans算法理論和代碼實(shí)現(xiàn)
- 改進(jìn)礁遣，BiKMeans算法理論和代碼實(shí)現(xiàn)
- 實(shí)例雪侥，上車點(diǎn)規(guī)劃
- 抉擇碗殷，如何挑選最佳的聚類簇?cái)?shù)？
其中一二三部分參考了Peter Harrington的《機(jī)器學(xué)習(xí)實(shí)戰(zhàn)》速缨，第四部分是純手工打造的輪子锌妻，還請(qǐng)多多指教。
02 KMeans理論和算法實(shí)現(xiàn)
聚類是一種無(wú)監(jiān)督學(xué)習(xí)的方法旬牲，所謂“無(wú)監(jiān)督”仿粹，就是指參與訓(xùn)練的樣本沒有標(biāo)簽。
KMeans聚類算法過(guò)程如下：
1. 對(duì)于一組數(shù)據(jù)集原茅，隨機(jī)選取k個(gè)點(diǎn)作為質(zhì)心吭历，將數(shù)據(jù)集中的點(diǎn)歸為離其最近的質(zhì)心一簇，此時(shí)數(shù)據(jù)集被劃分為k個(gè)簇擂橘；
2. 對(duì)這k個(gè)簇晌区，重新計(jì)算各簇的質(zhì)心(均值)；
3. 根據(jù)新的質(zhì)心通贞，按照step1繼續(xù)聚類朗若，然后再根據(jù)聚類重新計(jì)算各簇質(zhì)心，直到質(zhì)心不再改變昌罩，分類完成哭懈。
說(shuō)白了，就是不斷地聚類茎用、劃分的過(guò)程遣总。
通過(guò)KMeans原理，可以看到幾個(gè)顯而易見的缺點(diǎn)：
1. 簇?cái)?shù)量k由用戶指定绘搞，無(wú)法預(yù)先知道最佳k值 >>解法：分為幾簇彤避，最終由輪廓系數(shù)S(i)決定，取輪廓系數(shù)最大的分類數(shù)(05節(jié)劇透)
2. 最終質(zhì)心可能與初始點(diǎn)選擇有關(guān) >> 因此KMeans的結(jié)果可能收斂到局部最小值夯辖，而不是全局最小值 >> 解法：
- BiKMeans（03節(jié)）
- KMeans++（KMeans++ 算法在選擇初始質(zhì)心時(shí)并不是隨機(jī)選擇琉预，而是選擇盡量相互分離的質(zhì)心，即蒿褂，下一個(gè)質(zhì)心點(diǎn)總是離上一個(gè)質(zhì)心點(diǎn)較遠(yuǎn)的點(diǎn)）
代碼實(shí)現(xiàn)
def loadDataSet(fileName):
dataList=[]
dataMat=[]
fr=open(fileName)
for line in fr.readlines():
curLine=line.strip().split('\t')
fltLine=list(map(float,curLine))
dataList.append(fltLine)
dataMat=mat(dataList)
return dataMat
def distEclud(vecA,vecB):
return sqrt(sum(power(vecA-vecB,2))) #歐式距離
#為輸入數(shù)據(jù)集構(gòu)造k個(gè)隨機(jī)中心圆米，中心位置在各特征最大最小值之間
def randCent(dataSet,k):
n=shape(dataSet)[1]
center=mat(zeros((k,n)))
for j in range(n): #對(duì)每個(gè)特征
minJ=min(dataSet[:,j])
rangeJ=float(max(dataSet[:,j])-minJ)
center[:,j]=mat(minJ+rangeJ*random.rand(k,1)) #質(zhì)心第j維坐標(biāo)在數(shù)據(jù)集第j維數(shù)據(jù)之間
return center
def KMeans(dataSet,k,distMeas=distEclud,createCent=randCent):
m=shape(dataSet)[0]
clusterAssment=mat(zeros((m,2))) #用于記錄各樣本當(dāng)前歸屬于哪個(gè)簇以及到該簇質(zhì)心的歐式距離平方
center=createCent(dataSet,k)
clusterChanged=True
while clusterChanged:
clusterChanged=False
#對(duì)每個(gè)樣本卒暂，計(jì)算樣本到各質(zhì)心的距離，尋找距離最近的質(zhì)心娄帖，將該樣本歸為該質(zhì)心所在簇
for i in range(m):
minDist=inf;minIndex=-1
for j in range(k): #對(duì)每個(gè)質(zhì)心,計(jì)算到i樣本的距離
distJI=distMeas(center[j,:],dataSet[i,:])
if distJI<minDist:
minDist=distJI;minIndex=j #i樣本暫屬于j簇也祠，到j(luò)簇質(zhì)心距離為minDist
if clusterAssment[i,0]!=minIndex:
clusterChanged=True #若任一樣本在本次迭代中改變了簇類，則要進(jìn)行下一次迭代(即近速，直到任何樣本都不再改變簇類诈嘿，聚類停止)
clusterAssment[i,:]=minIndex,minDist**2 #記錄樣本i的簇類情況
#print (center)
#更新質(zhì)心
for cent in range(k):
ptsInClust=dataSet[nonzero(clusterAssment[:,0].A==cent)[0]] #篩選出屬于當(dāng)前簇類的點(diǎn)
center[cent,:]=mean(ptsInClust,axis=0) #對(duì)該簇類各樣本的各列求均值，作為新質(zhì)心
return center,clusterAssment
用一組數(shù)據(jù)來(lái)測(cè)試一下：
dataMat1=loadDataSet(r'D:\DM\python\data\MLiA_SourceCode\machinelearninginaction\Ch10\testSet.txt')
dataMat2=loadDataSet(r'D:\DM\python\data\MLiA_SourceCode\machinelearninginaction\Ch10\testSet2.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)
結(jié)果如下削葱，藍(lán)色點(diǎn)為聚類質(zhì)心
03 BiKMeans理論和算法實(shí)現(xiàn)
剛才我們實(shí)現(xiàn)了KMeans算法奖亚，不過(guò)也提到，KMeans有個(gè)缺點(diǎn)析砸，就是其聚類的最終質(zhì)心可能與初始點(diǎn)選擇有關(guān)昔字，因此KMeans的結(jié)果可能收斂到局部最小值，而不是全局最小值首繁。
怎么解決呢作郭？提示：想想決策樹是如何分支的？
為了得到全局最優(yōu)解弦疮，BiKMeans算法出現(xiàn)了夹攒，它的過(guò)程如下（是不是跟決策樹分支有點(diǎn)神似？）
1. 將整個(gè)數(shù)據(jù)集看作一個(gè)簇挂捅，計(jì)算初始質(zhì)心芹助，即所有數(shù)據(jù)點(diǎn)各特征的均值
2. 遍歷各質(zhì)心，對(duì)各質(zhì)心闲先，將質(zhì)心所在簇用原始KMeans算法二分状土，計(jì)算二分后整個(gè)數(shù)據(jù)集的SSE(即平方誤差和，即簇各點(diǎn)到簇質(zhì)心距離平方和)伺糠，找到二分后整體數(shù)據(jù)集SSE最小的質(zhì)心蒙谓，認(rèn)為此質(zhì)心是本次劃分的最佳質(zhì)心，對(duì)其進(jìn)行二分
3. 不斷重復(fù)step2训桶，直到質(zhì)心總數(shù)=設(shè)置的k
說(shuō)白了累驮，BiKMeans算法過(guò)程類似于決策樹的分支，通過(guò)啟發(fā)的方法舵揭，每次迭代只分裂當(dāng)前最佳質(zhì)心谤专，直到簇?cái)?shù)量達(dá)到k，這樣的方法可以保證最終劃分得到的質(zhì)心是全局最優(yōu)解午绳，而原始KMeans可能會(huì)陷入局部最優(yōu)解置侍。
代碼實(shí)現(xiàn)
def biKMeans(dataSet,k,distMeas=distEclud):
m=shape(dataSet)[0]
clusterAssment=mat(zeros((m,2))) #記錄各樣本歸屬和距離平方
center0=mean(dataSet,axis=0).tolist()[0] #初始質(zhì)心
centerList=[center0] #用于記錄聚類質(zhì)心坐標(biāo)
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
#對(duì)每個(gè)簇，嘗試二分，計(jì)算二分后整體數(shù)據(jù)的SSE蜡坊，若小于lowestSSE杠输，則將簇二分，如此往復(fù)秕衙，直到分到k個(gè)簇為止
for i in range(len(centerList)):
ptsInCluster=dataSet[nonzero(clusterAssment[:,0].A==i)[0],:] #默認(rèn)質(zhì)心索引就是數(shù)據(jù)對(duì)應(yīng)簇類
centerMat,splitClustAss=KMeans(ptsInCluster,2,distMeas)
sseSplit=sum(splitClustAss[:,1]) #被二分后的簇的平方誤差和
sseNotSplit=sum(clusterAssment[nonzero(clusterAssment[:,0]!=i)[0],1]) #整體數(shù)據(jù)中蠢甲，未被二分的簇的平方誤差和
if (sseSplit+sseNotSplit)<lowestSSE:
bestCentToSplit=i
bestNewCents=centerMat
bestClustAss=splitClustAss.copy() #.copy()防止splitClustAss被覆蓋時(shí)影響到bestClustAss
lowestSSE=sseSplit+sseNotSplit
#確定好本次迭代被二分的簇后，將被二分的數(shù)據(jù)對(duì)應(yīng)的簇類更新
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] #將被二分的原簇質(zhì)心替換為二分后的質(zhì)心之一
centerList.append(bestNewCents[1,:].tolist()[0]) #將二分后的另一質(zhì)心添加在質(zhì)心列表末尾
#將二分后的簇的數(shù)據(jù)歸屬更新到總記錄中
clusterAssment[nonzero(clusterAssment[:,0].A==bestCentToSplit)[0],:]=bestClustAss
return mat(centerList),clusterAssment
測(cè)試
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 實(shí)例：上車點(diǎn)規(guī)劃
我們?cè)谇肮?jié)實(shí)現(xiàn)了KMeans算法和BiKMeans算法鹦牛，現(xiàn)在讓我們來(lái)用用這兩個(gè)輪子吧。
設(shè)想若河，你邀請(qǐng)了70個(gè)朋友參加你的party勺择，他們住在城市的各個(gè)地方冒窍，你需要在party開始前2小時(shí)開車去接他們(假設(shè)你的車子可以容納70人揍拆，哈哈哈)溶锭。
為了時(shí)間效率最大化秦躯，請(qǐng)問(wèn)你該如何規(guī)劃每個(gè)人的上車點(diǎn)呢策精？總不能你一個(gè)一個(gè)去家門口接吧冀墨。
我們用聚類的方法來(lái)規(guī)劃鸠儿！
現(xiàn)在我們獲取了每個(gè)朋友住址的經(jīng)緯度钥屈，那么請(qǐng)開始你的表演悟民。
很簡(jiǎn)單，只需要調(diào)用我們?cè)绾玫妮喿覤iKMeans
#球面距離計(jì)算函數(shù)：球面余弦距離(向量夾角*地球半徑)
#求球面上兩向量vecA,vecB的距離
def distSLC(vecA,vecB):
a=sin(vecA[0,1]*pi/180)*sin(vecB[0,1]*pi/180) #由于抓取的經(jīng)緯度為角度篷就，需要通過(guò) *pi/180來(lái)轉(zhuǎn)換為弧度
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(r'D:\DM\python\data\MLiA_SourceCode\machinelearninginaction\Ch10\places.txt').readlines():
lineArr=line.strip().split('\t')
dataList.append([float(lineArr[4]),float(lineArr[3])]) #讀取地址的經(jīng)緯度
dataMat=mat(dataList)
center,clusterAss=biKMeans(dataMat,k,distMeas=distSLC) #將地址按經(jīng)緯度聚類
#作圖
fig=plt.figure(figsize=(10,8))
rect=[0.1,0.1,0.8,0.8] #用于設(shè)置坐標(biāo)軸刻度，[]中前兩個(gè)值表示左邊起始位置竭业，后兩個(gè)值對(duì)應(yīng)坐標(biāo)長(zhǎng)度
scatterMarkers=['^','o','h','8','p','d','v','s','>','<'] #用于設(shè)置散點(diǎn)圖點(diǎn)的形狀
axprops=dict(xticks=[],yticks=[])
ax0=fig.add_axes(rect,label='ax0',**axprops)
#讀圖智润，并將圖片顯示在設(shè)定好的坐標(biāo)軸中
imgP=plt.imread(r'D:\DM\python\data\MLiA_SourceCode\machinelearninginaction\Ch10\Portland.png')
ax0.imshow(imgP)
#將地址按聚類結(jié)果作散點(diǎn)圖
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)
#標(biāo)出質(zhì)心
ax1.scatter(center[:,0].flatten().A[0],center[:,1].flatten().A[0],marker='+',s=300)
plt.show()
好了，我們寫的這個(gè)“上車點(diǎn)規(guī)劃器”就可以使用了未辆，只需要輸入聚類簇?cái)?shù)量即可 窟绷，下面看看幾個(gè)測(cè)試結(jié)果：
但是到底聚為幾類比較合理呢？這個(gè)合理是指即不讓你的朋友跑太遠(yuǎn)去坐車咐柜，也不用你跑太多地方去接人兼蜈。
05 如何挑選最佳的聚類簇?cái)?shù)
不論是KMeans算法還是BiKMeans算法，都仍有一個(gè)缺點(diǎn)沒有解決：簇?cái)?shù)量k由用戶指定拙友，用戶指定的k不一定是最佳的簇?cái)?shù)量为狸。這也是上一節(jié)我們沒有解決的問(wèn)題。
怎么辦呢遗契？
使用輪廓系數(shù)辐棒，最佳k值由輪廓系數(shù)S(i)決定，取輪廓系數(shù)最大的分類數(shù)。
什么是輪廓系數(shù)涉瘾，輪廓系數(shù)S(i)用于衡量聚類效果知态，取值在-1~1之間，越接近1表示聚類效果越好立叛，公式如下负敏，
其中a(i)=i點(diǎn)到同簇各點(diǎn)距離均值；b(i)=min(i點(diǎn)到某個(gè)非同簇各點(diǎn)均值)秘蛇，即i點(diǎn)到其他簇質(zhì)心距離的最小值其做，整體數(shù)據(jù)集的輪廓系數(shù)是各點(diǎn)輪廓系數(shù)的平均值
基于此，我們可以寫一個(gè)計(jì)算聚類結(jié)果輪廓系數(shù)的函數(shù)赁还，然后看看輪廓系數(shù)與聚類數(shù)量k的關(guān)系妖泄，從而找到最佳的聚類數(shù)量k。
輪廓系數(shù)計(jì)算函數(shù)
def outlineOfCluster(filename,maxClusterNum,distMeas=distEclud):
dataList=[]
for line in open(filename).readlines():
lineArr=line.strip().split('\t')
dataList.append([float(lineArr[4]),float(lineArr[3])]) #讀取地址的經(jīng)緯度
dataMat=mat(dataList)
#遍歷2~nn個(gè)質(zhì)心艘策，求不同數(shù)量的簇的輪廓系數(shù)蹈胡，設(shè)置輪廓系數(shù)最大的簇?cái)?shù)量為k值
sk={}
for k in range(2,maxClusterNum+1):
center,clusterAss=biKMeans(dataMat,k,distMeas)
outline=[] #代表聚類為k個(gè)簇時(shí)各點(diǎn)的輪廓系數(shù)列表
for i in range(k): #遍歷各簇i
ptsInCluster=clusterAss[nonzero(clusterAss[:,0].A==i)[0],:]
for j in range(len(ptsInCluster)): #遍歷i簇的各點(diǎn)j
ptsInClusterNotJ=vstack([ptsInCluster[:j,:],ptsInCluster[j+1:,:]])
ajn=[]
for n in range(len(ptsInClusterNotJ)): #遍歷i簇非j的點(diǎn)
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簇質(zhì)心
bjm.append(distSLC(centerWithoutI[m],ptsInCluster[j]))
bj=min(bjm)
sj=(bj-aj)/max(bj,aj) #i簇中j點(diǎn)的輪廓系數(shù)
outline.append(sj) #將i簇中各點(diǎn)的輪廓系數(shù)保存在outline中，outline用于存儲(chǔ)聚類為k個(gè)簇時(shí)各點(diǎn)的輪廓系數(shù)朋蔫，在遍歷k時(shí)需要重置
sk[k]=nanmean(outline) #聚類為k個(gè)簇時(shí)的輪廓系數(shù)是各點(diǎn)輪廓系數(shù)的均值
return sk
現(xiàn)在我們可以接近上車點(diǎn)規(guī)劃問(wèn)題的遺留問(wèn)題了罚渐，規(guī)劃幾個(gè)上車點(diǎn)最合理？
sk=outlineOfCluster(r'D:\DM\python\data\MLiA_SourceCode\machinelearninginaction\Ch10\places.txt',35,distMeas=distSLC)
plt.figure(figsize=(10,5))
plt.plot(list(sk.keys()),list(sk.values()),linewidth=3)
plt.title('聚類簇?cái)?shù)-效果',fontsize=18,fontweight='bold')
plt.xlabel('最終聚類簇?cái)?shù)量',fontsize=14)
plt.ylabel('輪廓系數(shù)',fontsize=14)
可以看到驯妄，聚類簇?cái)?shù)在達(dá)到10個(gè)簇之后荷并，輪廓系數(shù)的增量就變得非常小了，因此從性價(jià)比方面考慮青扔，選擇聚類為10個(gè)簇的性價(jià)比是最高的源织，即規(guī)劃10個(gè)上車點(diǎn)就好。
06 總結(jié)
本文實(shí)踐了KMeans聚類算法微猖，分成以下幾個(gè)部分谈息，
- KMeans算法理論和代碼實(shí)現(xiàn)
- 改進(jìn)，BiKMeans算法理論和代碼實(shí)現(xiàn)
- 實(shí)例励两，上車點(diǎn)規(guī)劃
- 抉擇黎茎，如何挑選最佳的聚類簇?cái)?shù)？
下期我們將實(shí)踐Apriori算法当悔，這是一種挖掘關(guān)聯(lián)規(guī)則得到頻繁項(xiàng)集的算法傅瞻，比如可以預(yù)測(cè)你買了A商品后會(huì)買B商品，敬請(qǐng)期待~
07 參考
《機(jī)器學(xué)習(xí)實(shí)戰(zhàn)》 Peter Harrington Chapter10
