K中心点算法之PAM

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Oct 22 20:31:32 2017

@author: LPS
"""

import numpy as np
import pandas as pd
import copy

df = np.loadtxt('waveform.txt',delimiter=',')    # 载入waveform数据集，22列，最后一列为标签0，1，2
s= np.array(df)
print(s.shape)
print(s[0:10])

data0 = s[s[:,s.shape[1]-1]==0][:100]   # 取标签为0的前100个样本
data1 = s[s[:,s.shape[1]-1]==1][:100]   # 取标签为1的前100个样本
data2 = s[s[:,s.shape[1]-1]==2][:100]   # 取标签为2的前100个样本

data = np.array([data0,data1,data2])   
data = data.reshape(-1,22)    

def dis(data_a, data_b):
    return np.sqrt(np.sum(np.square(data_a - data_b), axis=1))    # 返回欧氏距离

def kmeans_wave(n=10, k=3, data=data):
    data_new = copy.deepcopy(data)  # 前21列存放数据，不可变。最后1列即第22列存放标签，标签列随着每次迭代而更新。
    data_now = copy.deepcopy(data)  # data_now用于存放中间过程的数据

    center_point = np.random.choice(300,3,replace=False)
    center = data_new[center_point,:20]   # 随机形成的3个中心，维度为（3，21）


    distance = [[] for i in range(k)]
    distance_now = [[] for i in range(k)]  # distance_now用于存放中间过程的距离
    lost = np.ones([300,k])*float('inf')   # 初始lost为维度为（300，3）的无穷大

    for j in range(k):   # 首先完成第一次划分，即第一次根据距离划分所有点到三个类别中
        distance[j] = np.sqrt(np.sum(np.square(data_new[:,:20] - np.array(center[j])), axis=1))
    data_new[:, 21] = np.argmin(np.array(distance), axis=0)  # data_new 的最后一列，即标签列随之改变，变为距离某中心点最近的标签，例如与第0个中心点最近，则为0
 
    for i in range(n):    # 假设迭代n次

        for m in range(k):   # 每一次都要分别替换k=3个中心点，所以循环k次。这层循环结束即算出利用所有点分别替代3个中心点后产生的900个lost值

            for l in range(300):  # 替换某个中心点时都要利用全部点进行替换，所以循环300次。这层循环结束即算出利用所有点分别替换1个中心点后产生的300个lost值

                center_now = copy.deepcopy(center)   # center_now用于存放中间过程的中心点
                center_now[m] = data_now[l,:20]   # 用第l个点替换第m个中心点
                for j in range(k):  # 计算暂时替换1个中心点后的距离值
                    distance_now[j] = np.sqrt(np.sum(np.square(data_now[:,:20] - np.array(center_now[j])), axis=1))
                data_now[:, 21] = np.argmin(np.array(distance), axis=0)  # data_now的标签列更新，注意data_now时中间过程，所以这里不能选择更新data_new的标签列

                lost[l, m] = (dis(data_now[:, :20], center_now[data_now[:, 21].astype(int)]) \
                       - dis(data_now[:, :20], center[data_new[:, 21].astype(int)])).sum()   # 这里很好理解lost的维度为什么为300*3了。lost[l,m]的值代表用第l个点替换第m个中心点的损失值

        if np.min(lost) <0:   # lost意味替换代价，选择代价最小的来完成替换
            index = np.where(np.min(lost) == lost)  # 即找到min(lost)对应的替换组合
            index_l = index[0][0]   # index_l指将要替代某个中心点的候选点
            index_m = index[1][0]   # index_m指将要被替代的某个中心点，即用index_l来替代index_m

        center[index_m] = data_now[index_l,:20]       #更新聚类中心

        for j in range(k):
            distance[j] = np.sqrt(np.sum(np.square(data_now[:, :20] - np.array(center[j])), axis=1))
        data_new[:, 21] = np.argmin(np.array(distance), axis=0)  # 更新参考矩阵,至此data_new的标签列得以更新，即完成了一次迭代

    return data_new  # 最后返回data_new，其最后一列即为最终聚好的标签


if __name__ == '__main__':
    data_new = kmeans_wave(10,3,data)
    print(data_new.shape)
    print(np.mean(data[:,21] == data_new[:,21]))   # 验证划分准确度