from numpy import *
def load_data(file_name):
data=[]
fr=open(file_name)
for line in fr.readlines():
cur_line=line.strip().split('\t')
flt_line=map(float,cur_line)
data.append(flt_line)
return data
def distance(a,b):
'''计算俩向量的欧式距离'''
return sqrt(sum(power(a-b,2)))
def rand_cent(data,k):
n=shape(data)[1]
#n数据的列数
centroids=mat(zeros((k,n)))
#初始质点矩阵,k为质点数
for j in range(n):
min_j=min(data[:,j])
#选出第j列中最小的数
range_j=float(max(data[:,j])-min_j)
#用j列中最大的减去最小的
centroids[:,j]=min_j+range_j*random.rand(k,1)
#获得随机质点,质点在数据范围之内
return centroids
def k_means(data,k,distance=distance,create_cent=rand_cent):
'''k均值聚类算法'''
m=shape(data)[0]
cluster_assment=mat(zeros((m,2)))
#cluster_assment第一列为质点,第二列为距离
centroids=create_cent(data,k)
cluster_changed=True
while cluster_changed:
cluster_changed=False
for i in range(m):
min_distance=inf;min_index=-1
#初始最小距离为正无穷,最小索引为-1
for j in range(k):
#j 相当于质点脚标
#k为4 ,则就有4个质点,j就是每个质点的脚标
dist_ji=distance(centroids[j,:],data[i,:])
#计算每个随机质点与各数据的距离
if dist_ji<min_distance:
#若随机质点与数据的距离 ji 小于最小距离,则最小距离就继承ji
#最小索引为j
min_distance=dist_ji;min_index=j
if cluster_assment[i,0]!=min_index:cluster_changed=True
#cluster_assment 第i行第1个值是质点的脚标,就是这行数据属于哪个质点往哪聚
#若这行数据的质点不等于上面得出的质点,则 cluster_changed为真\
#这个while还要在循环一下,循环到判断不成立为止
cluster_assment[i,:]=min_index,min_distance**2
#让第i行数据的质点等于上面得到的质点,第二个参数是距离的平方
print centroids
#质点是在变换的
for cent in range(k):
pts_in_clust=data[nonzero(cluster_assment[:,0].A==cent)[0]]
#pts_in_clust是数据中质点等于cent且非零的那一行元素所构成的矩阵
centroids[cent,:]=mean(pts_in_clust,axis=0)
#更新质点为(x,y);x=pts_in_clust第一列的平均数;y=pts_in_clust第二列的平均数
#这一步也就是k均值聚类的核心,所谓的均值
#就是说质点是接近于这一类数据的平均数,想想圆心就明白了
#各数据到质点的距离要短,每一集群的点都是到这一质点距离最短的点
return centroids,cluster_assment
def bisection_kmeans(data,k):
'''二分k-均值聚类算法'''
m=shape(data)[0]
cluster_assment=mat(zeros((m,2)))
#初始聚类估值矩阵
centroid0=mean(data,axis=0).tolist()[0]
#初始质点为所有数据的平均数.仅一个点
#数据有两列,正好两个平均值构成一个二维点
#tolist()转换成列表
cent_list=[centroid0]
#质点列表
for j in range(m):
#各数据点到预估质点的欧式距离的平方
#平方是为了拉开差距,有利于二分。使小于1的越小,大于1的越大
cluster_assment[j,1]=distance(mat(centroid0),data[j,:])**2
while (len(cent_list)<k):
#不应该自定义k
lowest_sse=inf
#最低sse
#sse误差平方和
for i in range(len(cent_list)):
pts_in_curr_cluster=data[nonzero(cluster_assment[:,0].A==i)[0],:]
#取出非零下预估的质点为i的数据
centroid_mat,split_clust_ass=k_means(pts_in_curr_cluster,2,distance)
#再根据k均值聚类法把上面取出的数据再分下,k设为2就是为了二分
#因为是二分所以centroid_mat有两个质点
sse_after_split=sum(split_clust_ass[:,1])
#划分后sse
#split_clust_ass[:,1]这一列是各数据点到各自质点的距离
#怎么才是最佳聚类,就是让这个所有的距离和最小
#以每个数据点为质点距离和为0,但就不是聚类的
#所以最佳的聚类是要质点的数量与距离和之间的平衡点
sse_before_split=sum(cluster_assment[nonzero(cluster_assment[:,0].A!=i)[0],1])
#划分前sse
#取出估值中质点不等于i的距离值,然后求和
print 'sse after and before split',sse_after_split,sse_before_split
if (sse_after_split+sse_before_split)<lowest_sse:
#以此来判断是否需要再分下去
#若划分前后的数据sse之和小于最低sse
#则更新最佳划分质点,最佳质点矩阵,最佳预估结果,最低sse
best_cent_to_split=i
best_new_cents=centroid_mat.tolist()
best_clust_ass=split_clust_ass.copy()
lowest_sse=sse_after_split+sse_before_split
best_clust_ass[nonzero(best_clust_ass[:,0].A==1)[0],0]=len(cent_list)
#最佳预估结果的中质点等于1的更改为质点len(cent_list)
#因为质点列表程增长状态(二分法)
best_clust_ass[nonzero(best_clust_ass[:,0].A==0)[0],0]=best_cent_to_split
#最佳预估结果中质点等于0 的更改为上面得出的最佳质点,因为0是默认值
print 'the best_cent_to_split is :',best_cent_to_split
print 'the len of best_clust_ass is :',len(best_clust_ass)
cent_list[best_cent_to_split]=best_new_cents[0]
#更新质点列表
cent_list.append(best_new_cents[1])
cluster_assment[nonzero(cluster_assment[:,0].A==best_cent_to_split)[0],:]=best_clust_ass
return cent_list,cluster_assment
def test():
data=load_data('testSet2.txt')
data=mat(data)
cent,clus=k_means(data,3)
cent_list,cluster_assment=bisection_kmeans(data,3)
s= mat(cent_list)
#print s
#print 'centroids',cent
# print 'cluster',clus
import matplotlib.pyplot as plt
fig=plt.figure()
ax=fig.add_subplot(111)
x=[float(i) for i in data[:,0]]
y=[float(i) for i in data[:,1]]
ax.scatter(x,y)
x1=s[:,0]
y1=s[:,-1]
#print x1,y1
x2=cent[:,0]
y2=cent[:,1]
ax.plot(x1,y1,'+',c='red')
ax.plot(x2,y2,'*',c='green')
plt.show()
test()
图中红色十字为二分法为质点,绿星味均值法质点。两者能得到相同的结果。
样本数据:
3.275154 2.957587
-3.344465 2.603513
0.355083 -3.376585
1.852435 3.547351
-2.078973 2.552013
-0.993756 -0.884433
2.682252 4.007573
-3.087776 2.878713
-1.565978 -1.256985
2.441611 0.444826
-0.659487 3.111284
-0.459601 -2.618005
2.177680 2.387793
-2.920969 2.917485
-0.028814 -4.168078
3.625746 2.119041
-3.912363 1.325108
-0.551694 -2.814223
2.855808 3.483301
-3.594448 2.856651
0.421993 -2.372646
1.650821 3.407572
-2.082902 3.384412
-0.718809 -2.492514
4.513623 3.841029
-4.822011 4.607049
-0.656297 -1.449872
1.919901 4.439368
-3.287749 3.918836
-1.576936 -2.977622
3.598143 1.975970
-3.977329 4.900932
-1.791080 -2.184517
3.914654 3.559303
-1.910108 4.166946
-1.226597 -3.317889
1.148946 3.345138
-2.113864 3.548172
0.845762 -3.589788
2.629062 3.535831
-1.640717 2.990517
-1.881012 -2.485405
4.606999 3.510312
-4.366462 4.023316
0.765015 -3.001270
3.121904 2.173988
-4.025139 4.652310
-0.559558 -3.840539
4.376754 4.863579
-1.874308 4.032237
-0.089337 -3.026809
3.997787 2.518662
-3.082978 2.884822
0.845235 -3.454465
1.327224 3.358778
-2.889949 3.596178
-0.966018 -2.839827
2.960769 3.079555
-3.275518 1.577068
0.639276 -3.412840