今天手寫實現了k-means，目的是加深對這個演算法原理的理解，有不足的地方請多指教。

ris鳶尾花資料集包含3個不同品種的鳶尾花（Setosa，Versicolour，and Virginica）資料，花瓣和萼片長度，儲存在一個150*4的 numpy.ndarry中

150行4列，150行指150多花，4列分別是Sepal Length，Sepal Width, Petal Length and Petal Width

使用pandas官方demo

本程式碼使用k-means實現對鶯尾花種類class的區分，最後進行了視覺化

匯入庫

from numpy import *
import matplotlib.pyplot as plt
import pandas as pd
 

# Load dataset
url = "https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data"
names = ['sepal-length', 'sepal-width', 'petal-length', 'petal-width', 'class']
dataset = pd.read_csv(url, names=names)

dataset['class'][dataset['class']=='Iris-setosa']=0
dataset['class'][dataset['class']=='Iris-versicolor']=1
dataset['class'][dataset['class']=='Iris-virginica']=2
#對類別進行編碼，3個類別分別賦值0，1，2
 

#算距離
def distEclud(vecA, vecB):                  #兩個向量間歐式距離
    return sqrt(sum(power(vecA - vecB, 2))) #la.norm(vecA-vecB)

#初始化聚類中心：通過在區間範圍隨機產生的值作為新的中心點
def randCent(dataSet, k):
    #獲取特徵維度
    n = shape(dataSet)[1]   
    #建立聚類中心0矩陣 k x n
    centroids = mat(zeros((k,n)))
    #遍歷n維特徵         
    for j in range(n):     
        #第j維特徵屬性值min   ,1x1矩陣                 
        minJ = min(dataSet[:,j])        
        #區間值max-min，float數值    
        rangeJ = float(max(dataSet[:,j]) - minJ)   
        #第j維，每次隨機生成k箇中心
        centroids[:,j] = mat(minJ + rangeJ * random.rand(k,1))
    return centroids
 

def randChosenCent(dataSet,k):
    # 樣本數
    m=shape(dataSet)[0]
    # 初始化列表
    centroidsIndex=[]
    #生成類似於樣本索引的列表
    dataIndex=list(range(m))
    for i in range(k):
        #生成隨機數
        randIndex=random.randint(0,len(dataIndex))
        #將隨機產生的樣本的索引放入centroidsIndex
        centroidsIndex.append(dataIndex[randIndex])
        #刪除已經被抽中的樣本
        del dataIndex[randIndex]
    #根據索引獲取樣本
    centroids = dataSet.iloc[centroidsIndex]
    return mat(centroids)

def kMeans(dataSet, k):
    # 樣本總數
    m = shape(dataSet)[0]
    #分配樣本到最近的簇：存[簇序號,距離的平方]
    # m行  2 列
    clusterAssment = mat(zeros((m,2)))

    #step1:
    #通過隨機產生的樣本點初始化聚類中心
    centroids = randChosenCent(dataSet, k)
    print('最初的中心=',centroids)

    #標誌位，如果迭代前後樣本分類發生變化值為Tree，否則為False
    clusterChanged = True
    #檢視迭代次數
    iterTime=0
    #所有樣本分配結果不再改變，迭代終止
    while clusterChanged:   
        clusterChanged = False        
        #step2:分配到最近的聚類中心對應的簇中
        for i in range(m):
            #初始定義距離為無窮大
            minDist = inf;
            #初始化索引值
            minIndex = -1
            # 計算每個樣本與k箇中心點距離
            for j in range(k):
                #計算第i個樣本到第j箇中心點的距離
                distJI = distEclud(centroids[j,:],dataSet.values[i,:])
                #判斷距離是否為最小
                if distJI < minDist:
                    #更新獲取到最小距離
                    minDist = distJI
                    #獲取對應的簇序號
                    minIndex = j
            #樣本上次分配結果跟本次不一樣，標誌位clusterChanged置True
            if clusterAssment[i,0] != minIndex:
                clusterChanged = True
            clusterAssment[i,:] = minIndex,minDist**2 #分配樣本到最近的簇
        iterTime+=1
        sse=sum(clusterAssment[:,1])
        print('the SSE of %d'%iterTime + 'th iteration is %f'%sse)
        #step3:更新聚類中心
        for cent in range(k):#樣本分配結束後，重新計算聚類中心
            #獲取該簇所有的樣本點
            ptsInClust = dataSet.iloc[nonzero(clusterAssment[:,0].A==cent)[0]]
            #更新聚類中心：axis=0沿列方向求均值。
            centroids[cent,:] = mean(ptsInClust, axis=0) 
    return centroids, clusterAssment

def kMeansSSE(dataSet,k,distMeas=distEclud, createCent=randChosenCent):
    m = shape(dataSet)[0]
    #分配樣本到最近的簇：存[簇序號,距離的平方]
    clusterAssment=mat(zeros((m,2)))
    #step1:#初始化聚類中心
    centroids = createCent(dataSet, k)
    print('initial centroids=',centroids)
    sseOld=0
    sseNew=inf
    iterTime=0 #檢視迭代次數
    while(abs(sseNew-sseOld)>0.0001):
        sseOld=sseNew
        #step2:將樣本分配到最近的質心對應的簇中
        for i in range(m):
            minDist=inf;minIndex=-1
            for j in range(k):
                #計算第i個樣本與第j個質心之間的距離
                distJI=distMeas(centroids[j,:],dataSet.values[i,:])
                #獲取到第i樣本最近的質心的距離,及對應簇序號
                if distJI<minDist:
                    minDist=distJI;minIndex=j
            clusterAssment[i,:]=minIndex,minDist**2 #分配樣本到最近的簇
        iterTime+=1
        sseNew=sum(clusterAssment[:,1])
        print('the SSE of %d'%iterTime + 'th iteration is %f'%sseNew)
        #step3:更新聚類中心
        for cent in range(k):
            #樣本分配結束後，重新計算聚類中心
            ptsInClust=dataSet[nonzero(clusterAssment[:,0].A==cent)[0]]
            #按列取平均,mean()對array型別
            centroids[cent,:] = mean(ptsInClust, axis=0)
    return centroids, clusterAssment

#2維資料聚類效果顯示
def datashow(dataSet,k,centroids,clusterAssment):  #二維空間顯示聚類結果
    from matplotlib import pyplot as plt
    num,dim=shape(dataSet)  #樣本數num ,維數dim
    
    if dim!=2:
        print('sorry,the dimension of your dataset is not 2!')
        return 1
    marksamples=['or','ob','og','ok','^r','^b','<g'] #樣本圖形標記
    if k>len(marksamples):
        print('sorry,your k is too large,please add length of the marksample!')
        return 1
        #繪所有樣本
    for i in range(num):
        markindex=int(clusterAssment[i,0])#矩陣形式轉為int值, 簇序號
        #特徵維對應座標軸x,y；樣本圖形標記及大小
        plt.plot(dataSet.iat[i,0],dataSet.iat[i,1],marksamples[markindex],markersize=6)

    #繪中心點            
    markcentroids=['o','*','^']#聚類中心圖形標記
    label=['0','1','2']
    c=['yellow','pink','red']
    for i in range(k):
        plt.plot(centroids[i,0],centroids[i,1],markcentroids[i],markersize=15,label=label[i],c=c[i])
        plt.legend(loc = 'upper left')
    plt.xlabel('sepal length')  
    plt.ylabel('sepal width') 
   
    plt.title('k-means cluster result') #標題        
    plt.show()

#畫出實際影象
def trgartshow(dataSet,k,labels):
    from matplotlib import pyplot as plt

    num,dim=shape(dataSet)
    label=['0','1','2']
    marksamples=['ob','or','og','ok','^r','^b','<g']
    # 通過迴圈的方式，完成分組散點圖的繪製
    for i in range(num): 
        plt.plot(datamat.iat[i,0],datamat.iat[i,1],marksamples[int(labels.iat[i,0])],markersize=6 )
    for i in range(0,num,50): 
        plt.plot(datamat.iat[i,0],datamat.iat[i,1],marksamples[int(labels.iat[i,0])],markersize=6,label=label[int(labels.iat[i,0])] )
    plt.legend(loc = 'upper left')
    # 新增軸標籤和標題

    plt.xlabel('sepal length')  
    plt.ylabel('sepal width') 
   
    plt.title('iris true result') #標題 

    # 顯示圖形
    plt.show()   
    #label=labels.iat[i,0]

#聚類前，繪製原始的樣本點
def originalDatashow(dataSet):
        #樣本的個數和特徵維數
    num,dim=shape(dataSet)
    marksamples=['ob'] #樣本圖形標記
    for i in range(num):
        plt.plot(datamat.iat[i,0],datamat.iat[i,1],marksamples[0],markersize=5)
    plt.title('original dataset')
    plt.xlabel('sepal length')  
    plt.ylabel('sepal width') #標題
    plt.show()

if __name__=='__main__':
#=====kmeans聚類
    # # #獲取樣本資料
    datamat=dataset.loc[:, ['sepal-length','sepal-width']]
    #真實的標籤
    labels=dataset.loc[:, ['class']]
    # #原始資料顯示
    originalDatashow(datamat)

    # #*****kmeans聚類
    k=3 #使用者定義聚類數
    mycentroids,clusterAssment=kMeans(datamat,k)
    #mycentroids,clusterAssment=kMeansSSE(datamat,k)

    #繪圖顯示
    datashow(datamat,k,mycentroids,clusterAssment)
    trgartshow(datamat,3,labels)

最初的中心= [[ 5.8  2.7]
 [ 6.7  3.3]
 [ 6.3  3.3]]
the SSE of 1th iteration is 89.870000
the SSE of 2th iteration is 56.076254
the SSE of 3th iteration is 51.738059
the SSE of 4th iteration is 50.664805
the SSE of 5th iteration is 49.569934
the SSE of 6th iteration is 48.696179
the SSE of 7th iteration is 46.824378
the SSE of 8th iteration is 45.085428
the SSE of 9th iteration is 44.384855
the SSE of 10th iteration is 43.591498
the SSE of 11th iteration is 41.904928
the SSE of 12th iteration is 39.066514
the SSE of 13th iteration is 38.316500
the SSE of 14th iteration is 37.912536
the SSE of 15th iteration is 37.423306
the SSE of 16th iteration is 37.136261
the SSE of 17th iteration is 37.123702

總結：

為了在平面上進行視覺化，只選取了4個特徵中的2個特徵，可能用3個特徵會更好，之前在tf的官網上看過minists資料集在tensorboard上的三維畫圖，挺漂亮的。還不會用matlibplot畫3維散點，echarts應該可以。有時間可以嘗試一下。