数据算法之kmeans聚类

一、聚类算法

聚类属于无监督学习，是数据挖掘十大经典算法之一 。

---

二、k-means聚类算法简介

1、k-means聚类算法的逻辑

• a. 给定一组数据集，先确定好需要分类个数k，最大归类迭代次数maxit，分类点最小移动值cutoff；

• b. 对数据集进行归一化处理；

  
  归一化公式
• c. 初始化一个k行 X n列的矩阵centroid，n为分类数据集的特征个数，随机在数据集中抽取k组数据特征填入其中
• d. 遍历整个数据集，计算其与每个矩阵centroid的欧式距离，并将其赋值给距离最短的一个
• e. 一轮遍历之后重新计算每个类的中心，中心点的值为同一类所有点的平均值。并与原矩阵centroid计算欧式距离，若中心点偏移不超过cutoff值，则聚类结束，最终的centroid为每个类的中心点。否则用重新计算的中心循环过程d、e，直到maxit归0。

python实现代码如下：

# -*- coding: utf-8 -*-
"""
Created on Tue May 29 14:27:54 2018

@author: junzhu
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

#import data
def import_data():
    with pd.ExcelFile('users.xlsx') as x:
        df = pd.read_excel(x,sheetname = 'Sheet1')
    return df

#data standardization
def data_standardization(dataset):
    min_col = np.amin(dataset[:,1:],axis = 0)
#    max_col = np.amax(dataset[:,1:],axis = 0)
    ptp_col = np.ptp(dataset[:,1:],axis = 0)
    dataset[:,1:] = (dataset[:,1:]-min_col) / ptp_col
    return dataset

#add one column to the end
def add_column(dataset):
    add_col = np.zeros((len(dataset),1))
    added_dataset = np.concatenate((dataset,add_col),axis = 1)
    return added_dataset

#calculate,classify,update
def kmeans(k,maxit,cutoff,dataset_fin,centroid):
    old_centroid = np.array([])
    while not stop_signal(old_centroid,centroid,cutoff,maxit):
        #calculate distance and classify
        maxit -= 1
        row_num = dataset_fin.shape[0]
        for i in range(row_num):
            mindis = np.inf
            sample_data = dataset_fin[i,1:-1]
            for j in range(k):
                centroid_data = centroid[j,:-1]
                distance = o_distance(sample_data,centroid_data)
                if distance < mindis:
                    mindis = distance
                    dataset_fin[i,-1] = j + 1
        
        #update centroid
        old_centroid = centroid.copy()
        for n in range(1,k+1):
            new_tempt = dataset_fin[np.nonzero(dataset_fin[:,-1] == n)[0],1:-1]
            centroid[np.nonzero(centroid[:,-1] == n)[0],:-1] = np.mean(new_tempt,axis = 0)
    
    return dataset_fin,centroid


def o_distance(x1,x2):
    distance = np.sqrt(np.sum(np.power((x1-x2),2)))
    return distance
    
def stop_signal(old_centroid,centroid,cutoff,maxit = 20):  
    if len(old_centroid) == 0:
        return False
    elif (o_distance(old_centroid[:,:-1],centroid[:,:-1]) < cutoff):
        return True
    elif maxit == 0:
        return True
    else:
        return False
  

#################################################################
def run_main():
    
    #set-up k,cutoff,maxiteration
    k = 3
    cutoff = 0.02
    maxit = 10
    
    #import data
    dataset = import_data()
    row_num,column_num = dataset.shape
    
    #data standardization
    dataset = dataset.as_matrix()
    dataset_std = data_standardization(dataset)
    
    #add one column to the end
    dataset_fin = add_column(dataset_std)
    
    #initialize k-matrix
    centroid = np.zeros((k,column_num))
    centroid[:,-1] = np.arange(1,k+1)
    for i in range(k):
        centroid[i,:-1] = dataset_fin[np.random.choice(np.arange(row_num)),1:-1]
    
    #finally result
    dataset_kmeans , centroid_kmeans = kmeans(k,maxit,cutoff,dataset_fin,centroid)
    
    print('dataset_kmeans:\n',dataset_kmeans)
    

if __name__ == '__main__':
    run_main()

以上是k-means聚类的思路和代码。

2、k-means聚类的优点是：

1）原理比较简单，实现也是很容易，收敛速度快。
2）聚类效果较优。
3）算法的可解释度比较强。
4）主要需要调参的参数仅仅是簇数k。

3、k-means聚类的缺点有：

1）K值的选取不好把握。
2）如果各隐含类别的数据不平衡，比如各隐含类别的数据量严重失衡，或者各隐含类别的方差不同，则聚类效果不佳。
3）采用迭代方法，得到的结果只是局部最优。
4）对噪音和异常点比较的敏感。(改进1：离群点检测的LOF算法，通过去除离群点后再聚类，可以减少离群点和孤立点对于聚类效果的影响；改进2：改成求点的中位数，这种聚类方式即K-Mediods聚类（K中值）)。

针对k-means聚类的缺点，有一个改进的方法是采用二分k-means聚类。

---

三、二分k-means聚类

二分k-means是一种层次聚类方法，简单说就是一个不断分裂,择优录用的过程。
它的一个原则是：因为聚类的误差平方和能够衡量聚类性能，该值越小表示数据点月接近于它们的质心，聚类效果就越好。所以我们就需要对误差平方和最大的簇进行再一次的划分，因为误差平方和越大，表示该簇聚类越不好，越有可能是多个簇被当成一个簇了，所以我们首先需要对这个簇进行划分。

1、二分k-means聚类的逻辑

• a. 给定一组数据集，先确定好需要分类个数k；

• b. 对数据集进行归一化处理；

  
  归一化公式
• c. 初始化一个m行 X 2列的矩阵clusterAssment，m为待分类数据集的样本个数。clusterAssment[0]存放分类编号，clusterAssment[1]存放距离差的平方和
• d. 初始化列表centList，用于存放中心点值
• e. 首先利用k-means法将原数据分成2部分：A,B，并记录
• f. 然后将A,B再利用k-means法各分成2部分：A1,A2,B1,B2，计算A1,A2,B的距离差平方x和计算B1,B2,A的距离差平方y，选择x，y中小的那个值，例如x最小，则本轮后的分类为A1,A2,B三类，并记录
• g. 重复f步骤， 利用k-means法将A1,A2,B各分成2部分，并择优录用，直到最终分类个数达到k则停止。

python实现代码如下：

#encoding=utf-8

from numpy import*
import matplotlib.pyplot as plt

def getDistance(vecA,vecB):
    '''
    本函数用于计算欧氏距离
    :param vecA: 向量A
    :param vecB: 向量B
    :return:欧氏距离
    '''
    return sqrt(sum(power(vecA-vecB,2)))

def randCent(dataSet,k):
    '''
    本函数用于生成k个随机质心
    :param dataSet: 数据集，具有矩阵形式
    :param k:指定的质心个数
    :return:随机质心，具有矩阵形式
    '''
    n = shape(dataSet)[1] # 获取特征数目
    centRoids = mat(zeros((k,n)))
    for j in range(n):
        minJ = min(dataSet[:,j]) # 获取每个特征的最小值
        rangeJ = float(max(dataSet[:,j]-minJ)) # 获取每个特征的范围
        centRoids[:,j] = minJ + rangeJ*random.rand(k,1) # numpy下的rand表示随机生成k*1的随机数矩阵，范围0-1
    return centRoids

#k-means算法，同上面的k-means算法在代码上略有不同
def kMeans(dataSet,k,disMens = getDistance,createCent = randCent):
    '''
    本函数用于k均值聚类
    :param dataSet: 数据集，要求有矩阵形式
    :param k: 指定聚类的个数
    :param disMens: 求解距离的方式，除欧式距离还可以定义其他距离计算方式
    :param createCent: 生成随机质心方式
    :return:随机质心，簇索引和误差距离矩阵
    '''
    m = shape(dataSet)[0]
    clusterAssment = mat(zeros((m,2))) # 要为每个样本建立一个簇索引和相对的误差，所以需要m行的矩阵，m就是样本数
    centRoids = createCent(dataSet,k) # 生成随机质心
    clusterChanged = True
    while clusterChanged:
        clusterChanged = False
        for i in range(m): # 遍历所有样本
            minDist = inf;minIndex = -1 # 初始化最小值
            for j in range(k): # 遍历所有质心
                disJI = disMens(centRoids[j,:],dataSet[i,:])
                if disJI < minDist:
                    minDist = disJI;minIndex = j # 找出距离当前样本最近的那个质心
            if clusterAssment[i,0] != minIndex: # 更新当前样本点所属于的质心
                clusterChanged = True # 如果当前样本点不属于当前与之距离最小的质心，则说明簇分配结果仍需要改变
                clusterAssment[i,:] = minIndex,minDist**2
        for cent in range(k):
            ptsInClust = dataSet[nonzero(clusterAssment[:,0].A == cent)[0]]
            # nonzero 返回的是矩阵中所有非零元素的坐标，坐标的行数与列数个存在一个数组或矩阵当中
            # 矩阵支持检查元素的操作，所有可以写成matrix == int这种形式，返回的一个布尔型矩阵，代表矩阵相应位置有无此元素
            # 这里指寻找当前质心下所聚类的样本
            centRoids[cent,:] = mean(ptsInClust,axis = 0) # 更新当前的质心为所有样本的平均值，axis = 0代表对列求平均值
    return centRoids,clusterAssment

def binKMeans(dataSet, k, distMeas = getDistance):
    '''
    本函数用于二分k均值算法
    :param dataSet: 数据集，要求有矩阵形式
    :param k: 指定聚类个数
    :param distMeas: 求解距离的方式
    :return:质心，簇索引和误差距离矩阵
    '''
    m = shape(dataSet)[0]
    clusterAssment = mat(zeros((m,2)))
    centRoids0 = mean(dataSet,axis = 0).tolist()[0] # 初始化一个簇，只有一个质心，分量就是就是所有特征的均值
    # 注意，tolist函数用于将矩阵转化为一个列表，此列表为嵌套列表
    centList = [centRoids0]
    for j in range(m): # 遍历所有样本，计算所有样本与当前质心的距离作为误差
        clusterAssment[j,1] = distMeas(mat(centRoids0),dataSet[j,:])**2

    while (len(centList) < k): # 循环条件为当前质心数目还不够指定数目
        lowestSSE = inf  #初始化距离差平方和为一个无限大的数

        for i in range(len(centList)): # 遍历所有质心
            ptsCurrCluster = dataSet[nonzero(clusterAssment[:,0].A == i)[0],:] # 搜索到当前质心所聚类的样本
            centroidsMat,splitClusterAss = kMeans(ptsCurrCluster,2,distMeas) # 将当前分割成两个簇
            sseSplit = sum(splitClusterAss[:,1]) # 计算分裂簇后的SSE
            sseNotSplit = sum(clusterAssment[nonzero(clusterAssment[:,0].A != i)[0],1])

            # 计算分裂之前的SSE
            if (sseSplit + sseNotSplit) < lowestSSE: # 如果分裂之后的SSE小，则更新
                bestCent2Split = i
                bestNewCents = centroidsMat
                bestClustAss = splitClusterAss.copy()
                lowestSSE = sseSplit + sseNotSplit


        #重新编制簇的编号，凡是分裂后编号为1的簇，编号为质心列表长度，编号为0的簇，编号为最佳分裂质心的编号，以此更新
        bestClustAss[nonzero(bestClustAss[:,0].A == 1)[0],0] = len(centList)
        bestClustAss[nonzero(bestClustAss[:,0].A == 0)[0],0] = bestCent2Split

        centList[bestCent2Split] = bestNewCents[0,:].tolist()[0] # 添加分裂的质心到质心列表中
        centList.append(bestNewCents[1,:].tolist()[0])

        clusterAssment[nonzero(clusterAssment[:,0].A == bestCent2Split)[0],:] = bestClustAss

    return mat(centList),clusterAssment

dataSet = random.randint(10,high = 20,size = (30,3))
k = 3

centList,clusterAssment = binKMeans(dataSet, k, distMeas=getDistance)

---

四、密度DBscan聚类

和传统的K-Means算法相比，DBSCAN最大的不同就是不需要输入类别数k，当然它最大的优势是可以发现任意形状的聚类簇，而不是像K-Means，一般仅仅使用于凸的样本集聚类。同时它在聚类的同时还可以找出异常点，这点和BIRCH算法类似。
一般来说，如果数据集是稠密的，并且数据集不是凸的，那么用DBSCAN会比K-Means聚类效果好很多。
如果数据集不是稠密的，则不推荐用DBSCAN来聚类。

1、密度DBscan聚类的逻辑

• a. 定义邻域半径e和样本数量阈值Minpts（在领域半径内的最少点数量）；
• b. 找到满足条件a的点集合，称为核心对象x；
• c. 选取其中1个核心对象x：
  1. 查找其核心对象x领域半径内的点y
  2. 对每一个y，检查其是否也是核心对象，如果yes，对新的核心对象查找其领域半径内的点z
  3. 循环1、2，直到找到所有的核心对象及其领域半径内的点，所有找到的核心对象和点都标注为同一类。（类似于不断发展下线）
• d. 在原始样本集中去掉过程c的结果，
• e. 对剩余的核心对象重复过程c

2、密度DBscan聚类的优点：

1） 可以对任意形状的稠密数据集进行聚类，相对的，K-Means之类的聚类算法一般只适用于凸数据集。
2） 可以在聚类的同时发现异常点，对数据集中的异常点不敏感。
3） 聚类结果没有偏倚，相对的，K-Means之类的聚类算法初始值对聚类结果有很大影响。

3、密度DBscan聚类的缺点：

1）如果样本集的密度不均匀、聚类间距差相差很大时，聚类质量较差，这时用DBSCAN聚类一般不适合。
2） 如果样本集较大时，聚类收敛时间较长，此时可以对搜索最近邻时建立的KD树或者球树进行规模限制来改进。
3） 调参相对于传统的K-Means之类的聚类算法稍复杂，主要需要对距离阈值ϵ，邻域样本数阈值MinPts联合调参，不同的参数组合对最后的聚类效果有较大影响。

python实现代码如下：

import math
import numpy as np
import pylab as pl

 #数据集：每三个是一组分别是西瓜的编号，密度，含糖量
data = """
1,0.697,0.46,2,0.774,0.376,3,0.634,0.264,4,0.608,0.318,5,0.556,0.215,
6,0.403,0.237,7,0.481,0.149,8,0.437,0.211,9,0.666,0.091,10,0.243,0.267,
11,0.245,0.057,12,0.343,0.099,13,0.639,0.161,14,0.657,0.198,15,0.36,0.37,
16,0.593,0.042,17,0.719,0.103,18,0.359,0.188,19,0.339,0.241,20,0.282,0.257,
21,0.748,0.232,22,0.714,0.346,23,0.483,0.312,24,0.478,0.437,25,0.525,0.369,
26,0.751,0.489,27,0.532,0.472,28,0.473,0.376,29,0.725,0.445,30,0.446,0.459"""

#数据处理 dataset是30个样本（密度，含糖量）的列表
a = data.split(',')
dataset = [(float(a[i]), float(a[i+1])) for i in range(1, len(a)-1, 3)]

#计算欧几里得距离,a,b分别为两个元组
def dist(a, b):
    return math.sqrt(math.pow(a[0]-b[0], 2)+math.pow(a[1]-b[1], 2))

#算法模型
def DBSCAN(D, e, Minpts):
    #初始化核心对象集合T,聚类个数k,聚类集合C, 未访问集合P,
    T = set(); k = 0; C = []; P = set(D)
    
    #遍历样本集，寻找核心对象：
    #核心对象：某点x，在x周围e范围内有不少于Minpts个点，则x是核心对象
    for d in D:
        if len([ i for i in D if dist(d, i) <= e]) >= Minpts:
            T.add(d)

    #开始聚类
    while len(T):
        P_old = P

        #随机选择1个核心对象赋给o
        o = list(T)[np.random.randint(0, len(T))]

        #在原数据集中删除选中的核心对象o
        P = P - set(o)
        Q = []
        Q.append(o)
        
        while len(Q):          
            q = Q[0]
            #在原数据集中收集核心对象周围满足e条件的点
            Nq = [i for i in D if dist(q, i) <= e]

            #如果满足条件的点数量大于等于Minpts
            if len(Nq) >= Minpts:
               
                # S是P和Nq的交集                 
                S = P & set(Nq)

                # 把S添加到Q中，目的是循环S中所有点，看其中是否也有核心对象
                # 如果也有核心对象，则所有找到的核心对象及其周围点都是一类的
                Q += (list(S))
                                          
                P = P - S
                
            Q.remove(q)
        
        k += 1
        Ck = list(P_old - P)
        T = T - set(Ck)
        C.append(Ck)    
    
    return C

# 画图
def draw(C):
   colValue = ['r', 'y', 'g', 'b', 'c', 'k', 'm']
   for i in range(len(C)):
       coo_X = []    #x坐标列表
       coo_Y = []    #y坐标列表
       for j in range(len(C[i])):
           coo_X.append(C[i][j][0])
           coo_Y.append(C[i][j][1])
       pl.scatter(coo_X, coo_Y, marker='x', color=colValue[i%len(colValue)], label=i)

   pl.legend(loc='upper right')
   pl.show()

C = DBSCAN(dataset, 0.11, 5)
draw(C)

以上内容参考、借鉴了以下资料:
1、简单聚类方法K-means方法的实现
2、《机器学习实战》二分-kMeans算法（二分K均值聚类）
3、DBSCAN密度聚类算法