kmeans算法用Python怎么实现

1、从Kmeans说起

Kmeans是一个非常基础的聚类算法，使用了迭代的思想，关于其原理这里不说了。下面说一下如何在matlab中使用kmeans算法。

创建7个二维的数据点：

复制代码 代码如下:

x=[randn(3,2)*.4randn(4,2)*.5+ones(4,1)*[4 4]]

使用kmeans函数：

复制代码 代码如下:

class = kmeans(x, 2)

x是数据点，x的每一行代表一个数据；2指定要有2个中心点，也就是聚类结果要有2个簇。 class将是一个具有70个元素的列向量，这些元素依次对应70个数据点，元素值代表着其对应的数据点所处的分类号。某次运行后，class的值是：

复制代码 代码如下:

2

2

2

1

1

1

1

这说明x的前三个数据点属于簇2，而后四个数据点属于簇1。 kmeans函数也悉袭枝可以像下面这样使用：

复制代码 代码如下:

>>[class, C, sumd, D] = kmeans(x, 2)

class =

2

2

2

1

1

1

1

C =

4.06294.0845

-0.13410.1201

sumd =

1.2017

0.2939

D =

34.37270.0184

29.56440.1858

36.35110.0898

0.1247 37.4801

0.7537 24.0659

0.1979 36.7666

0.1256 36.2149

class依旧代表着每个数据点的分类C包含最终的中心点，一行代表一个中心点；sumd代表着每个中心点与所属簇内各个数据点的距离之和；D的

每一行也对应一个数据点，行中的数值依次是该数据点与各个中心点之间的距离，Kmeans默认使用的距离是欧几里得距离（参考资料[3]）的平方值。

kmeans函数使用的距离，也可以是曼哈顿距离（L1-距离），以及其他类型的距离，可以通过添加参数指定。

kmeans有几个缺点（这在很多资料上都有说明）：

1、最终簇的类别数目（即中心点或者说种子点的数目）k并不一定能事先知道，所以如何选一个合适的k的值是一个问题。

2、最睁敏开始的种子点的选择的好坏会影响到聚类结果。

3、对噪声和离群点敏感。

4、等等。

2、kmeans++算法的基本思路

kmeans++算法的主要工作体现在种子点的选择上，基本原则是使得各个种子点之间的距离尽可能的大，但是又得排除噪声的影响。 以下为基本思路：

1、从输入的数据点集合（要求有k个聚类）中随机选择一个点作为第一个聚类中心

2、对于数据集中的每一个点x，计算它与最近聚类中心(指已选择的聚类中心)的距离D(x)

3、选择一个新的数据点作为新的聚类中心，选择的原则是：D(x)较大的点，被选取作为聚类中心的概率较大

4、重复2和3直到k个聚类中心被选出来

5、利用这k个初始的聚类中心来运行标准的k-means算法

假定数据点集合X有n个数据点，依次用X(1)、X(2)、……、X(n)表示，那么，在第2步中依次计算每个数据点与最近的种子点（聚类中心）的

距离，依次得到D(1)、D(2)、……、D(n)构成的集合D。在D中，为了避免噪声，不能直接选取值最大的元素，应该选择值较大的元素，然后将其对应

的数据点作为种子点。

如何选择值较大的元素呢，下面是一种思路（暂未找到最初的来源，在资料[2]等地方均有提及，笔者换了一种让自己更好理解的说法）：

把集合D中的每个元素D(x)想象为一根线L(x)，线的长度就是元素的值。将这些线依次按照L(1)、L(2)、……、L(n)的顺序连接起来，组成长

线L。L(1)、L(2)、……、L(n)称为L的子线。根据概率的相关知识，如果我们在L上随机选择一个点，那么这个点所在的子线很有可能是比较长的子

线，而这个子线对应的数据点就可以作为种子点。下文中kmeans++的两种实现均是这个原理。

3、python版本的kmeans++

在http://rosettacode.org/wiki/K-means%2B%2B_clustering 中能找到禅孝多种编程语言版本的Kmeans++实现。下面的内容是基于python的实现（中文注释是笔者添加的）：

复制代码 代码如下:

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)

'''

points是数据点，nclusters是给定的簇类数目

cluster_centers包含初始化的nclusters个中心点，开始都是对象->(0,0,0)

'''

def kpp(points, cluster_centers):

cluster_centers[0] = copy(choice(points)) #随机选取第一个中心点

d = [0.0 for _ in xrange(len(points))] #列表，长度为len(points)，保存每个点离最近的中心点的距离

for i in xrange(1, len(cluster_centers)): # i=1...len(c_c)-1

sum = 0

for j, p in enumerate(points):

d[j] = nearest_cluster_center(p, cluster_centers[:i])[1] #第j个数据点p与各个中心点距离的最小值

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]

'''

points是数据点，nclusters是给定的簇类数目

'''

def lloyd(points, nclusters):

cluster_centers = [Point() for _ in xrange(nclusters)] #根据指定的中心点个数，初始化中心点，均为(0,0,0)

# call k++ init

kpp(points, cluster_centers) #选择初始种子点

# 下面是kmeans

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()

上述代码实现的算法是针对二维数据的，所以Point对象有三个属性，分别是在x轴上的值、在y轴上的值、以及所属的簇的标识。函数lloyd是

kmeans++算法的整体实现，其先是通过kpp函数选取合适的种子点，然后对数据集实行kmeans算法进行聚类。kpp函数的实现完全符合上述

kmeans++的基本思路的2、3、4步。

在聚类分析中，K-均值聚类算法（k-means algorithm）是无监督分类中的一种基本方法，其也称为C-均值算法，其基本思想是：通过迭代的方法，逐次更新各聚类中心的值，直至得到最好的聚类结果。假设要把绝渗友样本集分为c个类别，算法如下：（1）适当选择c个类的初始中心；（2）在第k次迭代中，对任意一个样本，求其到c个中心的距离，将该样本归到距离最短的中心所在的类，（3）利用均值等方法更新该类的中心值；（4）并槐对于所有的c个聚类中心，如果利用（2）（3）的迭代法更新后，值保持不变，则迭代结束，否则继续迭代。下面介绍作者编写的一个分两类的程序，可以把其作为函数调用。%% function [samp1,samp2]=kmeans(samp)作为调用喊丛函数时去掉注释符samp=[11.1506 6.7222 2.3139 5.9018 11.0827 5.7459 13.2174 13.8243 4.8005 0.9370 12.3576]%样本集[l0 l]=size(samp)%%利用均值把样本分为两类，再将每类的均值作为聚类中心th0=mean(samp)n1=0n2=0c1=0.0c1=double(c1)c2=c1for i=1:lif samp(i)<th0c1=c1+samp(i)n1=n1+1elsec2=c2+samp(i)n2=n2+1endendc1=c1/n1c2=c2/n2%初始聚类中心t=0cl1=c1cl2=c2c11=c1c22=c2%聚类中心while t==0samp1=zeros(1,l)samp2=samp1n1=1n2=1for i=1:lif abs(samp(i)-c11)<abs(samp(i)-c22)samp1(n1)=samp(i)cl1=cl1+samp(i)n1=n1+1c11=cl1/n1elsesamp2(n2)=samp(i)cl2=cl2+samp(i)n2=n2+1c22=cl2/n2endendif c11==c1 &&c22==c2t=1endcl1=c11cl2=c22c1=c11c2=c22end %samp1,samp2为聚类的结果。初始中心值这里采用均值的办法，也可以根据问题的性质，用经验的方法来确定，或者将样本集随机分成c类，计算每类的均值。k-均值算法需要事先知道分类的数量，这是其不足之处。

---