[机器学习]K-means原理与源码实现

K-means算法的主要思想就是以空间中的K个点为中心进行聚类，对最靠近它的对象进行归类。通过迭代的方法不断的更新各聚类中心的值，直到最好的聚类结果。

---

K的取值：
确定聚类数K没有最佳的方法，通常需要根据具体的问题由人工进行选择。非监督聚类没有比较直接的聚类评估方法，但是可以从簇内的稠密程度和簇间的离散程度来评估聚类的效果。最常见的方法有轮廓系数Silhouette Coefficient和Calinski-Harabaz Index。其中Calinski-Harabaz Index计算直接简单，得到的结果越大则聚类效果越好。计算公式如下：
$$s(K)=(\frac{tr(B_k)}{tr(W_k)}\cdot \frac{m-k}{k-1})$$

其中：m为训练集样本数，k为类别数。Bk为类别之间的协方差矩阵，Wk为内部数据之间的协方差矩阵。tr为矩阵的迹。
也就是说内部数据的协方差越小越好，类别之间的协方差越大越好，这样对应的Calinski-Harabaz Index分数也就越高

---

主要步骤：

1. 在N个数据中，随机挑选K个数据（也就是最后聚类为K类）做为聚类的初始中心。
2. 分别计算每个数据点到这K个中心点的欧式距离，离哪个中心点最近就分配到哪个簇中。
3. 重新计算这K个簇数据的坐标均值，将新的均值作为聚类的中心。
4. 重复2和3步骤，直到簇中心的坐标不再变换或者达到规定的迭代次数，形成最终的K个聚类。
   
   

---

源码实现：
有两个py文件。kmeans.py 和 kmeanstest.py。源码在kmeans.py中，kmeanstest.py为实例（后期会给出sklean运行代码）。

# %load kmeans.py

#导入模块
import numpy as np
import matplotlib.pyplot as plt
from math import sqrt

#计算欧式距离
def eucDistance(vec1,vec2):
    return sqrt(sum(pow(vec2-vec1,2)))

#初始聚类中心选择
def initCentroids(dataSet,k):
    numSamples,dim = dataSet.shape
    centroids = np.zeros((k,dim))
    for i in range(k):
        index = int(np.random.uniform(0,numSamples))
        centroids[i,:] = dataSet[index,:]
    return centroids

#K-means聚类算法，迭代
def kmeanss(dataSet,k):
    numSamples = dataSet.shape[0]
    clusterAssement = np.mat(np.zeros((numSamples,2)))
    clusterChanged = True
    #  初始化聚类中心
    centroids = initCentroids(dataSet,k)
    while clusterChanged:
        clusterChanged = False
        for i in range(numSamples):
            minDist = 100000.0
            minIndex = 0
            # 找到哪个与哪个中心最近
            for j in range(k):
                distance = eucDistance(centroids[j,:],dataSet[i,:])
                if distance<minDist:
                    minDist = distance
                    minIndex = j
              # 更新簇
            clusterAssement[i,:] = minIndex,minDist**2
            if clusterAssement[i,0]!=minIndex:
                clusterChanged = True
         # 坐标均值更新簇中心
        for j in range(k):
            pointsInCluster = dataSet[np.nonzero(clusterAssement[:0].A==j)[0]]
            centroids[j,:] = np.mean(pointsInCluster,axis=0)
    print('Congratulations,cluster complete!')
    return centroids,clusterAssement

#聚类结果显示
def showCluster(dataSet,k,centroids,clusterAssement):
    numSamples,dim = dataSet.shape
    mark = ['or','ob','og','ok','^r','+r','<r','pr']
    if k>len(mark):
        print('Sorry!')
        return 1
    for i in np.xrange(numSamples):
        markIndex = int(clusterAssement[i,0])
        plt.plot(centroids[i,0],centroids[i,1],mark[i],markersize=12)
    plt.show()

有如下80个二维数据的样本，存在testSet的文档中，经过数据预处理和分析得知数据集共有4个类，所以确定K=4。

# %load kmeanstest.py

#导入模块
import kmeans
import  numpy as np
import matplotlib.pyplot as plt
from math import sqrt

#从文件加载数据集
dataSet=[]
fileIn = open('./testSet.txt')
for line in fileIn.readlines():
    lineArr = line.strip().split('\t')
    dataSet.append([float(lineArr[0]),float(lineArr[1])])

#调用k-means进行数据聚类
dataSet = np.mat(dataSet)
k = 4
centroids,clusterAssement = kmeans.kmeanss(dataSet,k)

#显示结果
kmeans.showCluster(dataSet,centroids,clusterAssement)

使用k-means算法聚类结果如下：

---

用sklearn实现Kmeans：

随机生成二维多类数据，样本大致分为4类，调用sklearn中的聚类函数进行分析不同k值对聚类结果的影响以及评估方法。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets.samples_generator import make_blobs
from sklearn.cluster import KMeans

X, y = make_blobs(n_samples=1000, n_features=2, centers=[[-1,-1], [0,0], [1,1], [2,2]],\
                 cluster_std = [0.4, 0.2, 0.2, 0.2], random_state=9)
plt.scatter(X[:, 0], X[:, 1], marker='+')
plt.show()

K=2：

y_pred = KMeans(n_clusters=2, random_state=9).fit_predict(X)
plt.scatter(X[:, 0], X[:, 1], c=y_pred)
plt.show()

调用Calinski-Harabaz Index来评估聚类分数:

from sklearn import metrics
metrics.calinski_harabaz_score(X, y_pred)

K=3:

y_pred = KMeans(n_clusters=3, random_state=9).fit_predict(X)
plt.scatter(X[:, 0], X[:, 1], c=y_pred)
plt.show()
from sklearn import metrics
metrics.calinski_harabaz_score(X, y_pred)

K=4:

y_pred = KMeans(n_clusters=4, random_state=9).fit_predict(X)
plt.scatter(X[:, 0], X[:, 1], c=y_pred)
plt.show()
from sklearn import metrics
metrics.calinski_harabaz_score(X, y_pred)

可以发现当K=4的时候得分最高，在往后分值就会有所下降，所以现在的K值比较合适。

K=5:

y_pred = KMeans(n_clusters=5, random_state=9).fit_predict(X)
plt.scatter(X[:, 0], X[:, 1], c=y_pred)
plt.show()
from sklearn import metrics
metrics.calinski_harabaz_score(X, y_pred)