我们经常接触到的聚类分析,一般都是数值聚类,一种常见的做法是同时提取 N 种特征,将它们放在一起组成一个 N 维向量,从而得到一个从原始数据集合到 N 维向量空间的映射——你总是需要显式地或者隐式地完成这样一个过程,然后基于某种规则进行分类,在该规则下,同组分类具有最大的相似性。
假设我们提取到原始数据的集合为(x1, x2, …, xn),并且每个xi为d维的向量,K-means聚类的目的就是,在给定分类组数k(k ≤ n)值的条件下,将原始数据分成k类
S = {S1, S2, …, Sk},在数值模型上,即对以下表达式求最小值:
\underset{\mathbf{S}} {\operatorname{arg\,min}} \sum_{i=1}^{k} \sum_{\mathbf x_j \in S_i} \left\| \mathbf x_j - \boldsymbol\mu_i \right\|^2
这里μi 表示分类Si 的平均值。
那么在计算机编程中,其又是如何实现的呢?其算法步骤一般如下:
1、从D中随机取k个元素,作为k个簇的各自的中心。
2、分别计算剩下的元素到k个簇中心的相异度,将这些元素分别划归到相异度最低的簇。相异度一般用欧氏距离或者马氏距离来运算,距离越小,说明二者之间越相似。
3、根据聚类结果,重新计算k个簇各自的中心,计算方法是取簇中所有元素各自维度的算术平均数。
4、将D中全部元素按照新的中心重新聚类。
5、重复第4步,直到聚类结果不再变化。
6、将结果输出。
具体的过程参考文章:http://www.cnblogs.com/moondark/archive/2012/03/08/2385770.html
在数据挖掘中,K-Means算法是一种cluster analysis的算法,其主要是来计算数据聚集的算法,主要通过不断地取离种子点最近均值的算法。
问题
K-Means算法主要解决的问题如下图所示。我们可以看到,在图的左边有一些点,我们用肉眼可以看出来有四个点群,但是我们怎么通过计算机程序找出这几个点群来呢?于是就出现了我们的K-Means算法(Wikipedia链接)
K-Means要解决的问题
算法概要
这个算法其实很简单,如下图所示:
从上图中,我们可以看到,A,B,C,D,E是五个在图中点。而灰色的点是我们的种子点,也就是我们用来找点群的点。有两个种子点,所以K=2。
然后,K-Means的算法如下:
- 随机在图中取K(这里K=2)个种子点。
- 然后对图中的所有点求到这K个种子点的距离,假如点Pi离种子点Si最近,那么Pi属于Si点群。(上图中,我们可以看到A,B属于上面的种子点,C,D,E属于下面中部的种子点)
- 接下来,我们要移动种子点到属于他的“点群”的中心。(见图上的第三步)
- 然后重复第2)和第3)步,直到,种子点没有移动(我们可以看到图中的第四步上面的种子点聚合了A,B,C,下面的种子点聚合了D,E)。
这个算法很简单,但是有些细节我要提一下,求距离的公式我不说了,大家有初中毕业水平的人都应该知道怎么算的。我重点想说一下“求点群中心的算法”。
求距离即相异度的算法
1)Minkowski Distance公式——λ可以随意取值,可以是负数,也可以是正数,或是无穷大。
2)Euclidean Distance公式——也就是第一个公式λ=2的情况
3)CityBlock Distance公式——也就是第一个公式λ=1的情况
这三个公式的求中心点有一些不一样的地方,我们看下图(对于第一个λ在0-1之间)。
(1)Minkowski Distance (2)Euclidean Distance (3) CityBlock Distance
上面这几个图的大意是他们是怎么个逼近中心的,第一个图以星形的方式,第二个图以同心圆的方式,第三个图以菱形的方式。
K-Means的演示
如果你以”K Means Demo“为关键字到Google里查你可以查到很多演示。这里推荐一个演示:http://home.dei.polimi.it/matteucc/Clustering/tutorial_html/AppletKM.html
操作是,鼠标左键是初始化点,右键初始化“种子点”,然后勾选“Show History”可以看到一步一步的迭代。
注:这个演示的链接也有一个不错的K Means Tutorial。
K-Means++算法
K-Means主要有两个最重大的缺陷——都和初始值有关:
- K是事先给定的,这个K值的选定是非常难以估计的。很多时候,事先并不知道给定的数据集应该分成多少个类别才最合适。(ISODATA算法通过类的自动合并和分裂,得到较为合理的类型数目K)
- K-Means算法需要用初始随机种子点来搞,这个随机种子点太重要,不同的随机种子点会有得到完全不同的结果。(K-Means++算法可以用来解决这个问题,其可以有效地选择初始点)
我在这里重点说一下K-Means++算法步骤:
- 先从我们的数据库随机挑个随机点当“种子点”。
- 对于每个点,我们都计算其和最近的一个“种子点”的距离D(x)并保存在一个数组里,然后把这些距离加起来得到Sum(D(x))。
- 然后,再取一个随机值,用权重的方式来取计算下一个“种子点”。这个算法的实现是,先取一个能落在Sum(D(x))中的随机值Random,然后用Random -= D(x),直到其<=0,此时的点就是下一个“种子点”。
- 重复第(2)和第(3)步直到所有的K个种子点都被选出来。
- 进行K-Means算法。
相关的代码你可以在这里找到“implement the K-means++ algorithm”(墙)另,Apache的通用数据学库也实现了这一算法
K-Means算法应用
看到这里,你会说,K-Means算法看来很简单,而且好像就是在玩坐标点,没什么真实用处。而且,这个算法缺陷很多,还不如人工呢。是的,前面的例子只是玩二维坐标点,的确没什么意思。但是你想一下下面的几个问题:
1)如果不是二维的,是多维的,如5维的,那么,就只能用计算机来计算了。
2)二维坐标点的X,Y 坐标,其实是一种向量,是一种数学抽象。现实世界中很多属性是可以抽象成向量的,比如,我们的年龄,我们的喜好,我们的商品,等等,能抽象成向量的目的就是可以让计算机知道某两个属性间的距离。如:我们认为,18岁的人离24岁的人的距离要比离12岁的距离要近,鞋子这个商品离衣服这个商品的距离要比电脑要近,等等。
只要能把现实世界的物体的属性抽象成向量,就可以用K-Means算法来归类了。
在《k均值聚类(K-means)》 这篇文章中举了一个很不错的应用例子,作者用亚洲15支足球队的2005年到1010年的战绩做了一个向量表,然后用K-Means把球队归类,得出了下面的结果,呵呵。
- 亚洲一流:日本,韩国,伊朗,沙特
- 亚洲二流:乌兹别克斯坦,巴林,朝鲜
- 亚洲三流:中国,伊拉克,卡塔尔,阿联酋,泰国,越南,阿曼,印尼
其实,这样的业务例子还有很多,比如,分析一个公司的客户分类,这样可以对不同的客户使用不同的商业策略,或是电子商务中分析商品相似度,归类商品,从而可以使用一些不同的销售策略,等等。
最后给一个挺好的算法的幻灯片:http://www.cs.cmu.edu/~guestrin/Class/10701-S07/Slides/clustering.pdf
下面是我的MATLAB代码,效果不是特别好,可能需要改进,因为找中心的方式有好多种,这种方式可能会出现最终找不到最好的中心点的情况:
主函数:
clc;
clear all;
close all;
data1=ones(30,1);%产生数据
data2=ones(30,1)*50;
data3=ones(30,1)*2500;
data=[data1;data2;data3];
% [u re]=myKMeans(data,3);
[myerror,label,diedai]=mymykmeans(data,4);
%%%%%%%%%%%%%%%%
function [myerror,label,diedai]=mymykmeans(data,K)
[m,n]=size(data);%求data的大小,m代表样本的个数,n代表数据特征数
zhongxin=zeros(K,n);%存储中心点
juli=zeros(m,K);%存储各个样本到中心点的距离,行是所有的样本,列是每个样本与中心点的距离
juli1=zeros(m,K);
label=zeros(m,1);%存储各个样本的标签
diedaimax=500;%设置最大迭代次数
for i=1:K
zhongxin(i,:)=data(floor(rand*m),:);%随机产生K个聚类中心
end
% for i=1:n
% ma(i)=max(data(:,i)); %每一维最大的数每一列的最大值
% mi(i)=min(data(:,i)); %每一维最小的数
% for j=1:K
% zhongxin(j,i)=ma(i)+(mi(i)-ma(i))*rand(); %随机初始化,不过还是在每一维[min max]中初始化好些
% end
% end
%%下面开始进行迭代
for diedai=1:diedaimax
for i=1:m
for j=1:K
juli1(i,j)=sqrt(sum((data(i,:)-zhongxin(j,:)).^2));%一个样本与所有聚类中心的距离
end%一行代表一个样本,K列代表与K个聚类中心的距离
end
for i=1:m
[julisort,zuobiao]=sort(juli1(i,:));%将距离按照由小到大排序
label(i,1)=zuobiao(1,1);
% for k=1:K
% if zuobiao(1,1)==k%如果最小的距离的标签是k
% label(i,1)=k;%则第k个样本的标签就是k
% end
% end
end
sumaver=zeros(K,n);%初始化中心点的特征向量sumaver
geshu=zeros(K,1);%为求中心点的特征向量的平均值而设定的计数器
%计算聚类中心
for i=1:m
sumaver(label(i,1),:)=sumaver(label(i,1),:)+data(i,:);
geshu(label(i,1),1)=geshu(label(i,1),1)+1;
end
for k=1:K
sumaver(k,:)=sumaver(k,:)./geshu(k,1);
zhongxin(k,:)=sumaver(k,:);%更新聚类中心
end
myerror=0;
for i=1:m
for k=1:K%如果距离不再变化,即中心点不再变化,或者达到了最大的迭代次数,则停止迭代
myerror=myerror+sum((juli1(i,k)-juli(i,k)).^2);
end
end
juli=juli1;
juli1=zeros(m,K);
if myerror==0%如果所有的中心点不再移动
break;
end
end
另外,MATLAB中有自带的kmeans函数,可以直接使用:
K-means聚类算法采用的是将N*P的矩阵X划分为K个类,使得类内对象之间的距离最大,而类之间的距离最小。
使用方法:
Idx=Kmeans(X,K)
[Idx,C]=Kmeans(X,K)
[Idc,C,sumD]=Kmeans(X,K)
[Idx,C,sumD,D]=Kmeans(X,K)
各输入输出参数介绍:
X---N*P的数据矩阵
K---表示将X划分为几类,为整数
Idx---N*1的向量,存储的是每个点的聚类标号
C---K*P的矩阵,存储的是K个聚类质心位置
sumD---1*K的和向量,存储的是类内所有点与该类质心点距离之和
D---N*K的矩阵,存储的是每个点与所有质心的距离
[┈]=Kmeans(┈,’Param1’,’Val1’,’Param2’,’Val2’,┈)
其中参数Param1、Param2等,主要可以设置为如下:
1、’Distance’---距离测度
‘sqEuclidean’---欧氏距离
‘cityblock’---绝对误差和,又称L1
‘cosine’---针对向量
‘correlation’---针对有时序关系的值
‘Hamming’---只针对二进制数据
2、’Start’---初始质心位置选择方法
‘sample’---从X中随机选取K个质心点
‘uniform’---根据X的分布范围均匀的随机生成K个质心
‘cluster’---初始聚类阶段随机选取10%的X的子样本(此方法初始使用’sample’方法)
Matrix提供一K*P的矩阵,作为初始质心位置集合
3、’Replicates’---聚类重复次数,为整数
使用案例:
data=
5.0 3.5 1.3 0.3 -1
5.5 2.6 4.4 1.2 0
6.7 3.1 5.6 2.4 1
5.0 3.3 1.4 0.2 -1
5.9 3.0 5.1 1.8 1
5.8 2.6 4.0 1.2 0
[Idx,C,sumD,D]=Kmeans(data,3,’dist’,’sqEuclidean’,’rep’,4)
运行结果:
Idx =
1
2
3
1
3
2
C =
5.0000 3.4000 1.3500 0.2500 -1.0000
5.6500 2.6000 4.2000 1.2000 0
6.3000 3.0500 5.3500 2.1000 1.0000
sumD =
0.0300
0.1250
0.6300
D =
0.0150 11.4525 25.5350
12.0950 0.0625 3.5550
29.6650 5.7525 0.3150
0.0150 10.7525 24.9650
21.4350 2.3925 0.3150
10.2050 0.0625 4.0850
另外这个链接里面讲kmeans和kmeans++都非常好,还有MATLAB和Python的代码:http://www.jb51.NET/article/49395.htm,看完了这个里面的内容,kmeans和kmeans++ 的MATLAB和python实现就都懂了,强力推荐
下面是python中实现kmeans++的程序,在之前给的连接里也有:
from math import pi, sin, cos
from collections import namedtuple
from random import random, choice
from copy import copy
try:
import psyco
psyco.full()
except ImportError:
pass
FLOAT_MAX = 1e100
class Point:
__slots__ = ["x", "y", "group"]
def __init__(self, x=0.0, y=0.0, group=0):
self.x, self.y, self.group = x, y, group
def generate_points(npoints, radius):
points = [Point() for _ in xrange(npoints)]
# note: this is not a uniform 2-d distribution
for p in points:
r = random() * radius
ang = random() * 2 * pi
p.x = r * cos(ang)
p.y = r * sin(ang)
return points
def nearest_cluster_center(point, cluster_centers):
"""Distance and index of the closest cluster center"""
def sqr_distance_2D(a, b):
return (a.x - b.x) ** 2 + (a.y - b.y) ** 2
min_index = point.group
min_dist = FLOAT_MAX
for i, cc in enumerate(cluster_centers):
d = sqr_distance_2D(cc, point)
if min_dist > d:
min_dist = d
min_index = i
return (min_index, min_dist)
def kpp(points, cluster_centers):
cluster_centers[0] = copy(choice(points))
d = [0.0 for _ in xrange(len(points))]
for i in xrange(1, len(cluster_centers)):
sum = 0
for j, p in enumerate(points):
d[j] = nearest_cluster_center(p, cluster_centers[:i])[1]
sum += d[j]
sum *= random()
for j, di in enumerate(d):
sum -= di
if sum > 0:
continue
cluster_centers[i] = copy(points[j])
break
for p in points:
p.group = nearest_cluster_center(p, cluster_centers)[0]
def lloyd(points, nclusters):
cluster_centers = [Point() for _ in xrange(nclusters)]
# call k++ init
kpp(points, cluster_centers)
lenpts10 = len(points) >> 10
changed = 0
while True:
# group element for centroids are used as counters
for cc in cluster_centers:
cc.x = 0
cc.y = 0
cc.group = 0
for p in points:
cluster_centers[p.group].group += 1
cluster_centers[p.group].x += p.x
cluster_centers[p.group].y += p.y
for cc in cluster_centers:
cc.x /= cc.group
cc.y /= cc.group
# find closest centroid of each PointPtr
changed = 0
for p in points:
min_i = nearest_cluster_center(p, cluster_centers)[0]
if min_i != p.group:
changed += 1
p.group = min_i
# stop when 99.9% of points are good
if changed <= lenpts10:
break
for i, cc in enumerate(cluster_centers):
cc.group = i
return cluster_centers
def print_eps(points, cluster_centers, W=400, H=400):
Color = namedtuple("Color", "r g b");
colors = []
for i in xrange(len(cluster_centers)):
colors.append(Color((3 * (i + 1) % 11) / 11.0,
(7 * i % 11) / 11.0,
(9 * i % 11) / 11.0))
max_x = max_y = -FLOAT_MAX
min_x = min_y = FLOAT_MAX
for p in points:
if max_x < p.x: max_x = p.x
if min_x > p.x: min_x = p.x
if max_y < p.y: max_y = p.y
if min_y > p.y: min_y = p.y
scale = min(W / (max_x - min_x),
H / (max_y - min_y))
cx = (max_x + min_x) / 2
cy = (max_y + min_y) / 2
print "%%!PS-Adobe-3.0\n%%%%BoundingBox: -5 -5 %d %d" % (W + 10, H + 10)
print ("/l {rlineto} def /m {rmoveto} def\n" +
"/c { .25 sub exch .25 sub exch .5 0 360 arc fill } def\n" +
"/s { moveto -2 0 m 2 2 l 2 -2 l -2 -2 l closepath " +
" gsave 1 setgray fill grestore gsave 3 setlinewidth" +
" 1 setgray stroke grestore 0 setgray stroke }def")
for i, cc in enumerate(cluster_centers):
print ("%g %g %g setrgbcolor" %
(colors[i].r, colors[i].g, colors[i].b))
for p in points:
if p.group != i:
continue
print ("%.3f %.3f c" % ((p.x - cx) * scale + W / 2,
(p.y - cy) * scale + H / 2))
print ("\n0 setgray %g %g s" % ((cc.x - cx) * scale + W / 2,
(cc.y - cy) * scale + H / 2))
print "\n%%%%EOF"
def main():
npoints = 30000
k = 7 # # clusters
points = generate_points(npoints, 10)
cluster_centers = lloyd(points, k)
print_eps(points, cluster_centers)
main()
另外,python中,Scikit-learn 中有一个 K-Means implementation that uses k-means++ by default.
可以直接调用python中的包去用kmeans,这个网站里面的就是http://my.oschina.Net/u/175377/blog/84420,可以找到各种分类的工具包,非常好,里面的kmeans算法就是用Kmeans++去初始化开始的中心点的。学习的时候可以自己去尝试编一下,但是了解了算法的具体实现以后还是调别人的更加方便简单高效。
python代码:
from sklearn import cluster
import numpy as np
import random
#def genMatrix(rows,cols): #产生一个随机的二维数组
# matrix =[[random.uniform(0,1) for col in range(cols)] for row in range(rows)]
# return matrix
#matrix=genMatrix(2,3) #自定义的数组产生,不过array更好用
#m=np.array([random.uniform(0,1) for i in range(15)])
#m=m.reshape(3,5,1)
#m.shape
#m.shape[0]
#m.shape[1]
#m.shape[2]
data1=np.array([random.uniform(0,1) for i in range(0,30)])
data2=np.array([random.uniform(10,20) for i in range(0,30)])
data3=np.array([random.uniform(20,30) for i in range(0,30)])
data=np.array([data1,data2,data3])
k_means = cluster.KMeans(3)#函数在使用的时候是每一行代表一个数据
k_means.fit(data)
print k_means.labels_