Kmeans算法详解及MATLAB实现_使用马氏距离计算kmeans matlab

首先要来了解的一个概念就是聚类，简单地说就是把相似的东西分到一组，同 Classification (分类)不同，对于一个 classifier ，通常需要你告诉它“这个东西被分为某某类”这样一些例子，理想情况下，一个 classifier 会从它得到的训练集中进行“学习”，从而具备对未知数据进行分类的能力，这种提供训练数据的过程通常叫做 supervised learning (监督学习)，而在聚类的时候，我们并不关心某一类是什么，我们需要实现的目标只是把相似的东西聚到一起，因此，一个聚类算法通常只需要知道如何计算相似 度就可以开始工作了，因此 clustering 通常并不需要使用训练数据进行学习，这在 Machine Learning 中被称作 unsupervised learning (无监督学习)。


　　我们经常接触到的聚类分析，一般都是数值聚类，一种常见的做法是同时提取 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

　　

下面是我的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

http://www.cnblogs.com/hxsyl/p/4054583.html  这个链接里面有更加详细的使用MATLAB自带的kmeans函数的案例，可以参考。kmeans算法找中心的方式有很多种，需要综合考虑算法的效率和效果去进行选择。

另外这个链接里面讲kmeans和kmeans++都非常好，还有MATLAB和python的代码：http://www.jb51.net/article/49395.htm，看完了这个里面的内容，kmeans和kmeans++ 的MATLAB和python实现就都懂了，强力推荐

下面是python中实现kmeans++的程序，在之前给的连接里也有：

from math import pi, sin, cosfrom collections import namedtuplefrom random import random, choicefrom 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_