格式如下：
5 2 3
2 3 4 5 10 12 5 1 12 10
其中，5表示数据的数量，2表示数据的维度，3表示聚类的数量。
后面每两个实数对应一个数据点。
不过这个程序始数据维数固定为2，后来想改成4以下的任意维度，但没有改完，所以其实只能为2维。
我已经对程序做了注释。
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <conio.h>
#include <math.h>
#define SUCCESS 1
#define FAILURE 0
#define TRUE 1
#define FALSE 0
#define MAXVECTDIM 4 // 数据最大维数 (看来这个程序写了一半，维数实际上只能为2)
#define MAXPATTERN 1588 // 数据数量最大值
#define MAXCLUSTER 10 // 聚类最大值
// 聚类结构
struct aCluster {
double Center[MAXVECTDIM]// 中心/引力数据对象
int Member[MAXPATTERN] // 该聚类中数据在Pattern中的索引
int NumMembers // 数据的数量
}
struct aVector {
double Center[MAXVECTDIM]
int Size
}
static double Pattern[MAXPATTERN][MAXVECTDIM + 1]// 所有数据存放在这个数组中
// 所以的东西都搁System类里面了
class System {
private :
aCluster Cluster[MAXCLUSTER]// 聚类数组
int NumPatterns // 输入数据的数量
int SizeVector // 数据的维数
int NumClusters // 数据的聚类数
void DistributeSamples()// 根据中心聚类，重新分配数据到不同的聚类
int CalcNewClustCenters() // 重新计算新的聚类中心
double EucNorm(int, int)// 误差准则
int FindClosestCluster(int) // 查找最接近的聚类
public :
System() {}
int LoadPatterns(char * fname) // 从文件中读取数据
void InitClusters()// 初始化聚类
void RunKMeans() // 执行K-Means算法
void ShowClusters()// 显示聚类
void SaveClusters(char * fname)// 保存聚类到文件中
void ShowCenters() // 显示聚类中心数据
}
void System::ShowCenters() {
int i
printf("Cluster centers:\n")
for (i = 0i <NumClustersi++) {
Cluster[i].Member[0] = i
printf("ClusterCenter[%d]=(%f,%f)\n", i, Cluster[i].Center[0],Cluster[i].Center[1])
int b=0
}
printf("\n")
}
int System::LoadPatterns(char * fname) {
FILE * InFilePtr
int i, j
double x
if ( (InFilePtr = fopen(fname, "r")) == NULL)
return FAILURE
// 读入数据的数量，维度和聚类数的数量
fscanf(InFilePtr, "%d", &NumPatterns)
fscanf(InFilePtr, "%d", &SizeVector)
fscanf(InFilePtr, "%d", &NumClusters)
// 读入数据
for (i = 0i <NumPatternsi++) {
for (j = 0j <SizeVectorj++) {
fscanf(InFilePtr, "%lg", &x)
Pattern[i][j] = x
}
}
// 打印读入的数据
printf("Input patterns:\n")
for (i = 0i <NumPatternsi++) {
printf("Pattern[%d]=(%2.3f,%2.3f,%d,%d)\n", i, Pattern[i][0], Pattern[i][1],Pattern[i][2], Pattern[i][3])
}
printf("\n--------------------\n")
return SUCCESS
}
void System::InitClusters() {
int i, j
SizeVector=2// 看来开始数据维数固定为2，后来想改成4以下的任意维度，但没有改完
printf("Initial cluster centers:\n")
// 把前NumClusters个数据作为NumClusters个聚类的数据中心
for (i = 0i <NumClustersi++) {
Cluster[i].Member[0] = i// 记录这个数据的索引到第i个聚类中
for (j = 0j <SizeVectorj++) {
Cluster[i].Center[j] = Pattern[i][j]// 把这个数据作为数据中心
}
}
// 打印聚类的数据中心
for (i = 0i <NumClustersi++) {
printf("ClusterCenter[%d]=(%f,%f)\n", i, Cluster[i].Center[0], Cluster[i].Center[1])
}
printf("\n")
}
void System::RunKMeans() {
int converged// 是否收敛
int pass // 计算的趟数
pass = 1 // 第一趟
converged = FALSE// 没有收敛
while (converged == FALSE) { // 没有收敛的情况下反复跑
printf("PASS=%d\n", pass++)// 打印趟数
DistributeSamples()// 重新分配数据
converged = CalcNewClustCenters()// 计算新的聚类中心，如果计算结果和上次相同，认为已经收敛
ShowCenters()// 显示聚类中心
}
}
// 第p个数据到第c个聚类中心的误差准则(距离的平方)
double System::EucNorm(int p, int c) {
double dist, x// x可能是调试用，没有实际作用
int i
dist = 0
// 将每个维度的差的平方相加，得到距离的平方
for (i = 0i <SizeVectori++) {
x = (Cluster[c].Center[i] - Pattern[p][i]) *
(Cluster[c].Center[i] - Pattern[p][i])
dist += (Cluster[c].Center[i] - Pattern[p][i]) *
(Cluster[c].Center[i] - Pattern[p][i])
}
return dist
}
// 查找第pat个数据的最近聚类
int System ::FindClosestCluster(int pat) {
int i, ClustID/*最近聚类的ID*/
double MinDist/*最小误差*/, d
MinDist = 9.9e+99// 初始为一个极大的值
ClustID = -1
// 依次计算3个聚类到第pat个数据的误差，找出最小值
for (i = 0i <NumClustersi++) {
d = EucNorm(pat, i)
printf("Distance from pattern %d to cluster %d is %f\n", pat, i, sqrt(d))
if (d <MinDist) { // 如果小于前最小值，将改值作为当前最小值
MinDist = d
ClustID = i
}
}
if (ClustID <0) { // 没有找到？？！ —— 这个应该不可能发生，如果出现表示出现了不可思议的错误
printf("Aaargh")
exit(0)
}
return ClustID
}
void System ::DistributeSamples() {
int i, pat, Clustid, MemberIndex
// 将上次的记录的该聚类中的数据数量清0，重新开始分配数据
for (i = 0i <NumClustersi++) {
Cluster[i].NumMembers = 0
}
// 找出每个数据的最近聚类数据中心，并将该数据分配到该聚类
for (pat = 0pat <NumPatternspat++) {
Clustid = FindClosestCluster(pat)
printf("patern %d assigned to cluster %d\n\n", pat, Clustid)
MemberIndex = Cluster[Clustid].NumMembers // MemberIndex是当前记录的数据数量，也是新加入数据在数组中的位置
Cluster[Clustid].Member[MemberIndex] = pat// 将数据索引加入Member数组
Cluster[Clustid].NumMembers++ // 聚类中的数据数量加1
}
}
int System::CalcNewClustCenters() {
int ConvFlag, VectID, i, j, k
double tmp[MAXVECTDIM]// 临时记录新的聚类中心，因为需要比较，不能直接记录
char xs[255]
char szBuf[255]
ConvFlag = TRUE
printf("The new cluster centers are now calculated as:\n")
// 逐个计算NumClusters个聚类
for (i = 0i <NumClustersi++) {
strcpy(xs,"")
printf("Cluster Center %d (1/%d):\n",i,Cluster[i].NumMembers)
// tmp所有维度数值清0
for (j = 0j <SizeVectorj++) {
tmp[j] = 0.0
}
// 计算聚类中所有数据的所有维度数值的和，为下一步计算均值做准备
for (j = 0j <Cluster[i].NumMembersj++) {
VectID = Cluster[i].Member[j]
for (k = 0k <SizeVectork++) {
tmp[k] += Pattern[VectID][k]
sprintf(szBuf,"%d ",Pattern[VectID][k])
strcat(xs,szBuf)
}
printf("%s\n",xs)// 输出刚才计算的数据
strcpy(xs,"")
}
// 求出均值，并判断是否和上次相同
for (k = 0k <SizeVectork++) {
if(Cluster[i].NumMembers!=0)
tmp[k] = tmp[k] / Cluster[i].NumMembers
if (tmp[k] != Cluster[i].Center[k]) // 如果不同，则认为没有收敛
ConvFlag = FALSE
Cluster[i].Center[k] = tmp[k]// 用新值代替旧值
}
printf("%s,\n", xs)
}
return ConvFlag// 返回收敛情况
}
// 显示聚类中心数据
void System::ShowClusters() {
int cl
for (cl = 0cl <NumClusterscl++) {
printf("\nCLUSTER %d ==>[%f,%f]\n", cl, Cluster[cl].Center[0],
Cluster[cl].Center[1])
}
}
// 没有实现的函数
void System::SaveClusters(char * fname) {
}
#include <windows.h>
void main(int argc, char * argv[]) {
System kmeans
// 如果没有提供参数，显示用法
if (argc <2) {
printf("USAGE: KMEANS PATTERN_FILE\n")
exit(0)
}
// 参数1 为数据文件名，从中读入数据
if (kmeans.LoadPatterns(argv[1]) == FAILURE) {
printf("UNABLE TO READ PATTERN_FILE:%s\n", argv[1])
exit(0)
}
kmeans.InitClusters()
kmeans.RunKMeans()
kmeans.ShowClusters()
kmeans.ShowCenters()
return
}
补充一个Matlab实现方法：
function [cid,nr,centers] = cskmeans(x,k,nc)
% CSKMEANS K-Means clustering - general method.
%
% This implements the more general k-means algorithm, where
% HMEANS is used to find the initial partition and then each
% observation is examined for further improvements in minimizing
% the within-group sum of squares.
%
% [CID,NR,CENTERS] = CSKMEANS(X,K,NC) Performs K-means
% clustering using the data given in X.
%
% INPUTS: X is the n x d matrix of data,
% where each row indicates an observation. K indicates
% the number of desired clusters. NC is a k x d matrix for the
% initial cluster centers. If NC is not specified, then the
% centers will be randomly chosen from the observations.
%
% OUTPUTS: CID provides a set of n indexes indicating cluster
% membership for each point. NR is the number of observations
% in each cluster. CENTERS is a matrix, where each row
% corresponds to a cluster center.
%
% See also CSHMEANS
% W. L. and A. R. Martinez, 9/15/01
% Computational Statistics Toolbox
warning off
[n,d] = size(x)
if nargin <3
% Then pick some observations to be the cluster centers.
ind = ceil(n*rand(1,k))
% We will add some noise to make it interesting.
nc = x(ind,:) + randn(k,d)
end
% set up storage
% integer 1,...,k indicating cluster membership
cid = zeros(1,n)
% Make this different to get the loop started.
oldcid = ones(1,n)
% The number in each cluster.
nr = zeros(1,k)
% Set up maximum number of iterations.
maxiter = 100
iter = 1
while ~isequal(cid,oldcid) &iter <maxiter
% Implement the hmeans algorithm
% For each point, find the distance to all cluster centers
for i = 1:n
dist = sum((repmat(x(i,:),k,1)-nc).^2,2)
[m,ind] = min(dist)% assign it to this cluster center
cid(i) = ind
end
% Find the new cluster centers
for i = 1:k
% find all points in this cluster
ind = find(cid==i)
% find the centroid
nc(i,:) = mean(x(ind,:))
% Find the number in each cluster
nr(i) = length(ind)
end
iter = iter + 1
end
% Now check each observation to see if the error can be minimized some more.
% Loop through all points.
maxiter = 2
iter = 1
move = 1
while iter <maxiter &move ~= 0
move = 0
% Loop through all points.
for i = 1:n
% find the distance to all cluster centers
dist = sum((repmat(x(i,:),k,1)-nc).^2,2)
r = cid(i)% This is the cluster id for x
%%nr,nr+1
dadj = nr./(nr+1).*dist'% All adjusted distances
[m,ind] = min(dadj)% minimum should be the cluster it belongs to
if ind ~= r % if not, then move x
cid(i) = ind
ic = find(cid == ind)
nc(ind,:) = mean(x(ic,:))
move = 1
end
end
iter = iter+1
end
centers = nc
if move == 0
disp('No points were moved after the initial clustering procedure.')
else
disp('Some points were moved after the initial clustering procedure.')
end
warning on
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步。
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）