西瓜书9.10：实现自动确定聚类数目的k_means算法

问题：

试设计一个能自动确定聚类数的改进k均值算法，编程实现并在西瓜数据集4.0上运行。

数据集：

西瓜数据集4.0

数据集描述：

该数据集共有30个样本，每个样本有密度和含糖度两个特征。

思路：

如何确定k的取值：

这里希望每一类别内部样本点距离较小而不同聚类之间的误差较大，所以用：

作为损失函数，其值越小越好，希望取得局部最优点，当值在k时增大，则k=k-1时最优

如何选择较好的k个均值：

随机选取的问题：

若是均值随机选，则很可能寻找到两个距离很近的均值点，在迭代过程中，这两个均值点会逐渐靠近最后重合，使得最后聚类数目减少。

算法思想:

首先随机选取L个均值点，接着根据这L个均值点对样本聚类，减去那些聚类内样本较小的点。接着在剩下的L1个店中随机选取一个均值点，并且寻找距离该点最远的均值点，依次迭代直至选取到K个均值点。

算法流程：

k_means算法：

1. 初始化：首先选取k个均值点
2. 将样本点划分到距离自己最近的均值点所在类别
3. 根据聚类结果更新均值点
4. 重复步骤2.3直至均值点不再改变
5. 输出聚类结果

结果：

通过对损失函数计算，得到当k=3时，聚类结果最优，聚类结果以及图形化展示如下：

源码 :

损失函数计算：

clu_unique=np.unique(cluster)    Di=[]    D=len(data)    E=0.0    for i in range(len(mean)):        Di.append(len(np.where(cluster==i)[0]))    for i in range(D):        E+=np.linalg.norm(data[i]-mean[cluster[i][0]],ord=2)    for i in range(len(mean)):        for j in range(len(mean)):            E+=np.linalg.norm(mean[i]-mean[j],ord=2)    E-=np.log(len(mean)/D)    return E

优化版本K个均值值选取：

def find_k_means(data,K):    L=int(K*np.log(K))    if L<K:        L=K    np.random.seed(int(time.time()))    r_index=random_unique(0,data.shape[0],L)    mean=data[r_index]#随机选取L个中心    cluster = classify(data,mean)  # 纪录每个样本所属类别    remove_index=remove_center(cluster,data.shape[0],L)#删除以该中心开始聚类数目最少的中心点    new_mean=[]#新的中心点    for i in range(mean.shape[0]):        if (i not in remove_index) and (mean[i].tolist() not in new_mean):            new_mean.append(mean[i].tolist())    if len(new_mean)>K:#new_mean 里面元素不同        k_mean=[]        r_i=np.random.randint(0,len(new_mean),1)[0]        old_v=new_mean[r_i]        k_mean.append(old_v)        while(len(k_mean)<K):            dis=[]            for i in range(len(new_mean)):                dis.append([new_mean[i],np.linalg.norm(np.array(new_mean[i])-np.array(old_v),ord=2)])            max=np.max(np.array(dis)[:,1])            max_index=np.where(np.array(dis)[:,1]==max)[0]            while(dis[max_index[0]][0] in k_mean):                dis.pop(max_index[0])                max = np.max(np.array(dis)[:, 1])                max_index = np.where(np.array(dis)[:, 1] == max)[0]            old_v=dis[max_index[0]][0]#上一个距离最远的样本            k_mean.append(old_v)    else:        k_mean=new_mean    return np.array(k_mean)

K均值算法：

oldE=100color = ['green', 'red', 'purple', 'pink', 'yellow','green','brown','tan','seashell','salmon']mark=['^','o','*','.','#']old_cluster=[]for k in np.arange(2,len(data),1):    k_mean=find_k_means(data,k)    old_mean = np.zeros(shape=k_mean.shape)    cluster = np.zeros(shape=k_mean.shape)    while (not ((old_mean == k_mean).all())):#直到mean值不再改变 达到最优        cluster = classify(data, k_mean)        old_mean = k_mean        k_mean = update_means(data, cluster)    E = loss(data, cluster, k_mean)    if oldE<E:        print("在k="+str(k-1)+"次达到最优")        for i in np.unique(old_cluster):            index = np.where(old_cluster == i)            x = data[index, 0]            y = data[index, 1]            plt.scatter(x, y, color=color[i%10], marker=mark[i%5])        plt.show()        break    oldE=E    old_cluster=cluster    print("聚类结果：" + str(cluster.reshape((cluster.shape[1], cluster.shape[0]))))    print("聚类损失：" + str(E))