第十一章 K-Means(K均值)演算法模型實現（中）

來自專欄 數據模型的那些事兒

python2 代碼實現：

from numpy import *import numpydef loadDataSet(fileName): #general function to parse tab -delimited floats dataMat = [] #assume last column is target value fr = open(fileName) for line in fr.readlines(): curLine = line.strip().split( ) fltLine = map(float,curLine) #map all elements to float() dataMat.append(fltLine) return dataMatdef distEclud(vecA, vecB): return sqrt(sum(power(vecA - vecB, 2))) #la.norm(vecA-vecB)def randCent(dataSet, k): n = shape(dataSet)[1] centroids = mat(zeros((k,n)))#create centroid mat for j in range(n):#create random cluster centers, within bounds of each dimension minJ = min(dataSet[:,j]) rangeJ = float(max(dataSet[:,j]) - minJ) centroids[:,j] = mat(minJ + rangeJ * random.rand(k,1)) return centroidsdef kMeans(dataSet, k, distMeas=distEclud, createCent=randCent): m = shape(dataSet)[0] clusterAssment = mat(zeros((m,2)))#create mat to assign data points #to a centroid, also holds SE of each point centroids = createCent(dataSet, k) clusterChanged = True while clusterChanged: clusterChanged = False for i in range(m):#for each data point assign it to the closest centroid minDist = inf; minIndex = -1 for j in range(k): distJI = distMeas(centroids[j,:],dataSet[i,:]) if distJI < minDist: minDist = distJI; minIndex = j if clusterAssment[i,0] != minIndex: clusterChanged = True clusterAssment[i,:] = minIndex,minDist**2 print centroids for cent in range(k):#recalculate centroids ptsInClust = dataSet[nonzero(clusterAssment[:,0].A==cent)[0]]#get all the point in this cluster centroids[cent,:] = mean(ptsInClust, axis=0) #assign centroid to mean return centroids, clusterAssment#以上代碼基本跟上篇python3的代碼一致，除了2跟3語法不同的語句外

二分K均值演算法

為了克服K均值演算法收斂於局部最小值的問題，提出了二分K均值演算法。

演算法思想

該演算法首先將所有點作為一個簇，然後將該簇一分為2，之後選擇其中一個簇繼續進行劃分，劃分規則是按照最大化SSE（目標函數）的值。

主要步驟：

• 將所有點看成一個簇
• 計算每一個簇的總誤差
• 在給定的簇上進行K均值聚類，計算將簇一分為二的總誤差
• 選擇使得誤差最小的那個簇進行再次劃分
• 重複步驟2，直到簇的個數滿足要求

具體實現

def biKmeans(dataSet, k, distMeas=distEclud): m = shape(dataSet)[0] clusterAssment = mat(zeros((m,2)))#創建一個矩陣存儲每個點的簇分配結果及平方誤差 centroid0 = mean(dataSet, axis=0).tolist()[0]#計算整個數據集的質心 centList =[centroid0] #create a list with one centroid#使用一個列表來保留所有的質心 for j in range(m):#calc initial Error 遍曆數據集中所有點 clusterAssment[j,1] = distMeas(mat(centroid0), dataSet[j,:])**2#計算每個點到質心的誤差值 while (len(centList) < k):#該循環會不停對簇進行劃分，直到得到想要的簇數目為止，為此需要比較劃分前後的sse lowestSSE = inf#開始將最小SSE設為無窮大 for i in range(len(centList)):#遍歷簇列表centList中的每個簇來決定最佳的簇進行劃分 ptsInCurrCluster = dataSet[nonzero(clusterAssment[:,0].A==i)[0],:]#get the data points currently in cluster i 對每個簇，對該簇中的所有點看成一個小的數據集ptsInCurrCluster centroidMat, splitClustAss = kMeans(ptsInCurrCluster, 2, distMeas)#將ptsInCurrCluster輸入到函數kmeans()（k=2）中進行處理生成2個質心簇，並給出每個簇的誤差值#誤差與剩餘數據集的誤差之和將作為本次劃分的誤差 sseSplit = sum(splitClustAss[:,1])#compare the SSE to the currrent minimum sseNotSplit = sum(clusterAssment[nonzero(clusterAssment[:,0].A!=i)[0],1]) print "sseSplit, and notSplit: ",sseSplit,sseNotSplit if (sseSplit + sseNotSplit) < lowestSSE:#如果該劃分的sse值最小，則本次劃分保存...一旦決定了要劃分的簇，就要執行實際劃分操作，即將要劃分的簇中所有點的簇分配結果進行修改即可。當使用KMEANS()函數並簇數為2時，得到兩個編號0與1的結果簇，需要將這些簇編號修改改為劃分簇與新加簇的編號，該過程通過2個數組過濾器完成... bestCentToSplit = i bestNewCents = centroidMat bestClustAss = splitClustAss.copy() lowestSSE = sseSplit + sseNotSplit# bestClustAss[nonzero(bestClustAss[:,0].A == 1)[0],0] = len(centList) #change 1 to 3,4, or whatever bestClustAss[nonzero(bestClustAss[:,0].A == 0)[0],0] = bestCentToSplit print the bestCentToSplit is: ,bestCentToSplit print the len of bestClustAss is: , len(bestClustAss) centList[bestCentToSplit] = bestNewCents[0,:].tolist()[0]#replace a centroid with two best centroids #新的簇分配結果被更新 centList.append(bestNewCents[1,:].tolist()[0])#新的質心添加到centlist中 clusterAssment[nonzero(clusterAssment[:,0].A == bestCentToSplit)[0],:]= bestClustAss#reassign new clusters, and SSE return mat(centList), clusterAssment #while循環結束後，同kmeans（）函數一樣，函數返回質心列表與簇分配結果

下面幾個函數是將一個API字典里的地址轉換為經度與緯度，然後對給出的地址坐標進行聚類，最後畫出所有點以及簇中心，並看看聚類結果到底如何。

import urllibimport jsondef geoGrab(stAddress, city): #函數geoGrab()從雅虎返回一個字典 apiStem = http://where.yahooapis.com/geocode? #create a dict and constants for the goecoder params = {} params[flags] = J#JSON return type返回類型為json格式 params[appid] = aaa0VN6k params[location] = %s %s % (stAddress, city) url_params = urllib.urlencode(params)#urllib的urlencode()函數將創建的字典轉換為可通過URL進行傳遞的字元串格式 yahooApi = apiStem + url_params print yahooApi #列印輸出的URLc=urllib.urlopen(yahooApi)#打開url return json.loads(c.read())#讀取返回值由於返回值是json格式，所以可以使用json的python模塊來將其解碼為一個字典，一旦返回了解碼後的字典，也就意味著你成功的對一個地址進行了地理編碼。from time import sleepdef massPlaceFind(fileName):#massplacefind()函數將所有這些封裝起來並且將相關信息保存到文件中 fw = open(places.txt, w)#打開一個文本文件 for line in open(fileName).readlines(): line = line.strip() lineArr = line.split( )#tab分割的文件 retDict = geoGrab(lineArr[1], lineArr[2])#獲取第2列第3列結果，並輸入到geoGrab函數中 if retDict[ResultSet][Error] == 0:#檢查geoGrab()的輸出字典判斷有沒有錯誤，如果沒有錯誤，就可以從字典中讀取經緯度， lat = float(retDict[ResultSet][Results][0][latitude]) lng = float(retDict[ResultSet][Results][0][longitude]) print "%s %f %f" % (lineArr[0], lat, lng) fw.write(%s %f %f
% (line, lat, lng))#這些值被添加到原來的原來對應的行上，同時寫入到新的文件中 else: print "error fetching" #如果有錯誤，就不需要抽取經緯度 sleep(1)#利用sleep()函數將massPlaceFind()函數推遲1秒，為了保護不要在短時間內過於頻繁的調用API,因為如果頻繁調用，請求可能會被封掉，所以推遲一下比較好。 fw.close()對地理坐標進行聚類def distSLC(vecA, vecB):#Spherical Law of Cosines a = sin(vecA[0,1]*pi/180) * sin(vecB[0,1]*pi/180) b = cos(vecA[0,1]*pi/180) * cos(vecB[0,1]*pi/180) * cos(pi * (vecB[0,0]-vecA[0,0]) /180)#使用球面餘弦定理來計算2個經緯度之間的距離，因為兩級與赤道同時走相同的距離對經緯度的變化不同 return arccos(a + b)*6371.0 #pi is imported with numpy#返回地球表面兩點之間的距離

簇聚類及繪圖函數import matplotlibimport matplotlib.pyplot as pltdef clusterClubs(numClust=3): datList = [] for line in open(C:/Users/HZF/Desktop/python數據挖掘十大演算法實現/11.csv).readlines(): lineArr = line.strip().split(,) #print lineArr #np.array fltLine = map(float,lineArr) #print fltLine datList.append(fltLine) #print dataList #datList=numpy.array(datList) #print datList datMat = mat(datList) myCentroids, clustAssing= kMeans(datMat, numClust)#使用kmeans演算法聚類，這裡也可以調用bikmeans聚類 #print myCentroids, clustAssing fig = plt.figure() rect=[0.1,0.1,0.8,0.8] scatterMarkers=[s, o, ^, 8, p, d, v, h, >, <]#標記類型 axprops = dict(xticks=[], yticks=[]) ax0=fig.add_axes(rect, label=ax0, **axprops) plt.savefig(glp.png, dpi = 75) imgp=plt.imread(portland.png)#imread()函數基於一副圖像來創建矩陣 ax0.imshow(imgP)#繪製該矩陣 ax1=fig.add_axes(rect, label=ax1, frameon=False) for i in range(numClust):#遍歷每一個簇 ptsInCurrCluster = datMat[nonzero(clustAssing[:,0].A==i)[0],:] markerStyle = scatterMarkers[i % len(scatterMarkers)]#使用索引i%len(scatterMarkers)選擇標記形狀，當有更多簇時，可以循環使用這些標記 ax1.scatter(ptsInCurrCluster[:,0].flatten().A[0], ptsInCurrCluster[:,1].flatten().A[0], marker=markerStyle, s=90)#一一畫出來 ax1.scatter(myCentroids[:,0].flatten().A[0], myCentroids[:,1].flatten().A[0], marker=+, s=300)#使用+字標記表示簇中心 plt.show()#顯示在圖中if __name__==__main__: #主函數調用 clusterClubs(numClust=3)

這篇是對kmeans演算法的python2實現，下篇主要寫應用及參考文獻。會儘快更新！

（待續）

推薦閱讀：