上课，老师让做作业，kmeans分类，将iris进行聚类分类，于是就做了这个作业。很简单，我就将其设置为分三类，重点是我在选择初始center的时候，尝试了使用随机选取和使用大招选取两种方式，随机选择初始点的效果不如放大招。这个大招是这样的：先随机选取一个作为center1，再选取一个距离这个center1最远的点作为center2，再选择一个距离center1和center2最远的点作为center3。

在训练过程中更新质心的时候，我不是选择online模式，来一个点更新一次质心，而是选择了mini_batch的模式，先送进来mini_batch个样本，分别将其划归到对应的中心处，然后再进行更新。说这么多屁话，不如把代码放上来，代码如下：

import numpy as np
import csv
import random

features = np.loadtxt('iris.csv',delimiter=',',usecols=(1,2,3,4))  #read features
z_min, z_max = features.min(axis=0), features.max(axis=0)          #features normalized
features = (features - z_min)/(z_max - z_min)

csv_file = open('iris.csv')   #transform string to num label
csv_reader_lines = csv.reader(csv_file)
classes_list = []
for i in csv_reader_lines:
    classes_list.append(i[-1])
labels = []
for i in classes_list:
    if i=='setosa':
        labels.append(0)
    elif i=='versicolor':
        labels.append(1)
    else:
        labels.append(2)

labels = np.array(labels)
labels = labels.reshape((150,1))   # transformm list to numpy type

data_index = np.arange(features.shape[0])
np.random.shuffle(data_index)

train_input = features[ data_index[0:120] ]
train_label = labels[ data_index[0:120] ]

test_input = features[ data_index[120:150] ]
test_label = labels[ data_index[120:150] ]

train_length = 120
K = 3
center_1_pos = random.randint(0,train_length)
center1 = train_input[ center_1_pos ]
# center1 = train_input[0]
# center2 = train_input[1]
# center3 = train_input[2]
# print(center1)
# print(center2)
# print(center3)

biggest_distance = 0.0
center_2_pos = 0

for i in range(train_length):#选择center2
    dist = np.sum(pow( (center1 - train_input[i]),2 ))
    if dist > biggest_distance:
        biggest_distance = dist
        center_2_pos = i

center2 = train_input[center_2_pos]


biggest_distance = 0.0
center_3_pos = 0

for i in range(train_length):#选择center3
    dist = np.sum(pow( (center1 - train_input[i]), 2 )) + np.sum(pow( (center2 - train_input[i]) , 2))
    if dist > biggest_distance:
        biggest_distance = dist
        center_3_pos = i

center3 = train_input[center_3_pos]
mini_batch = 20

for epoch in range(10):#在整个数据集上训练10次
    for i in range(6):
        belong1 = []
        belong2 = []
        belong3 = []
        for j in range(mini_batch):#mini_batch
            temp_index = mini_batch * i + j
            belong = 1
            dist_1 = np.sum(pow( ( center1 - train_input[mini_batch*i+j] ),2 ))
            temp_dist = dist_1

            dist_2 = np.sum(pow((center2 - train_input[mini_batch * i + j]), 2))
            dist_3 = np.sum(pow((center3 - train_input[mini_batch * i + j]), 2))

            if(dist_2 < temp_dist):
                temp_dist = dist_2
                belong = 2
            if(dist_3 < temp_dist):
                belong = 3

            if belong==1:
                belong1.append( temp_index )
            elif belong == 2:
                belong2.append(temp_index)
            else:
                belong3.append(temp_index)


        for k in belong1:
            center1 = center1 + train_input[k]
        center1 = center1 / (1 + len(belong1))

        for k in belong2:
            center2 = center2 + train_input[k]
        center2 = center2 / (1 + len(belong2))

        for k in belong3:
            center3 = center3 + train_input[k]
        center3 = center3 / (1 + len(belong3))

    b_1=[]
    b_2=[]
    b_3=[]
    for l in range(test_input.shape[0]):#在测试机上进行测试
        belong = 1
        dist_1 = np.sum(pow((center1 - test_input[l]), 2))
        temp_dist = dist_1

        dist_2 = np.sum(pow((center2 - test_input[ l ]), 2))
        dist_3 = np.sum(pow((center3 - test_input[ l ]), 2))

        if (dist_2 < temp_dist):
            temp_dist = dist_2
            belong = 2
        if (dist_3 < temp_dist):
            belong = 3



        if belong == 1:
            b_1.append(test_label[l][0])
        elif belong == 2:
            b_2.append(test_label[l][0])
        else:
            b_3.append(test_label[l][0])
    print()
    print('epoch : {} / 10' .format(epoch+1))
    print('center1: ',b_1)
    print('center2',b_2)
    print('center3: ',b_3)

下面是我运行程序的结果：

epoch : 1 / 10
center1:  [2, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1]
center2 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
center3:  [2, 2, 2, 2, 2]

epoch : 2 / 10
center1:  [2, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1]
center2 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
center3:  [2, 2, 2, 2, 2]

epoch : 3 / 10
center1:  [2, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1]
center2 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
center3:  [2, 2, 2, 2, 2]

epoch : 4 / 10
center1:  [2, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1]
center2 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
center3:  [2, 2, 2, 2, 2]

epoch : 5 / 10
center1:  [2, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1]
center2 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
center3:  [2, 2, 2, 2, 2]

epoch : 6 / 10
center1:  [2, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1]
center2 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
center3:  [2, 2, 2, 2, 2]

epoch : 7 / 10
center1:  [2, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1]
center2 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
center3:  [2, 2, 2, 2, 2]


     小数据集上结果还是挺好的。

iris数据集的一个下载链接如下：链接: https://pan.baidu.com/s/1PpIwncqbtQbEuGKSxKyBMg  密码: 7ins

注意：我在读取iris的数据时候，将第一行的属性名称删除了，要不然处理起来麻烦。

哎，做这个小作业也挺有意思的，学习了Nummpy数据的归一化还，csv文件的numpy读取，以及如何将csv文件中的字符类型转化为数字标签。

首先介绍一下（yuan）鸢尾花数据集，该数据集测量了所有150个样本的4个特征单位都是cm，分别是：

1.sepal length（花萼长度）

2.sepal width（花萼宽度）

3.petal length（花瓣长度）

4.petal width（花瓣宽度）

对应一个150行4列的矩阵，本文思路是通过PCA压缩数据将其从150行四列变成150行3列，4维变3维（每个特征对应一个维度）

在网上看了很多，k-means和PCA的结果都是二维的，我想用三维图来展示分类的的效果图就自己将代码改了一下。发个博客记录一下代码，思想大概就是将鸢尾花4维的数据转为三维后再使用k-means分类，因为提前知道了这些数据是来自三类的花，所以提前给出K-means要分类的个数为3，算法会提前随机给3个类似于重心的东西，然后离重心近的会被分到该簇。之后每个簇的重心会随着其簇内样本的增加不断改变。重心大概就是结果图中X所在的位置。
#coding=utf-8
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.datasets import load_iris
from mpl_toolkits.mplot3d import Axes3D
from sklearn.cluster import KMeans
import numpy as np





##计算欧式距离
def distEuclid(x, y):
    return np.sqrt(np.sum((x - y) ** 2))


##随机产生n个dim维度的数据 (这里为了展示结果 dim取2或者3)
def genDataset(n, dim):
    data = []
    while len(data) < n:
        p = np.around(np.random.rand(dim) * size, decimals=2)
        data.append(p)
    return data


## 初始化簇中心点 一开始随机从样本中选择k个 当做各类簇的中心
def initCentroid(data, k):
    num, dim = data.shape
    centpoint = np.zeros((k, dim))
    l = [x for x in range(num)]
    np.random.shuffle(l)
    for i in range(k):
        index = int(l[i])
        centpoint[i] = data[index]
    return centpoint


##进行KMeans分类
def KMeans(data, k):
    ##样本个数
    num = np.shape(data)[0]

    ##记录各样本 簇信息 0:属于哪个簇 1:距离该簇中心点距离
    cluster = np.zeros((num, 2))
    cluster[:, 0] = -1

    ##记录是否有样本改变簇分类
    change = True
    ##初始化各簇中心点
    cp = initCentroid(data, k)

    while change:
        change = False

        ##遍历每一个样本
        for i in range(num):
            minDist = 9999.9
            minIndex = -1

            ##计算该样本距离每一个簇中心点的距离 找到距离最近的中心点
            for j in range(k):
                dis = distEuclid(cp[j], data[i])
                if dis < minDist:
                    minDist = dis
                    minIndex = j

            ##如果找到的簇中心点非当前簇 则改变该样本的簇分类
            if cluster[i, 0] != minIndex:
                change = True
                cluster[i, :] = minIndex, minDist

        ## 根据样本重新分类  计算新的簇中心点
        for j in range(k):
            pointincluster = data[[x for x in range(num) if cluster[x, 0] == j]]
            cp[j] = np.mean(pointincluster, axis=0)

    print("finish!")
    return cp, cluster


##展示结果  各类簇使用不同的颜色  中心点使用X表示
def Show(data, k, cp, cluster):
    num, dim = data.shape
    color = ['b', 'r', 'g', 'c', 'y', 'm', 'k']
    ##二维图
    if dim == 2:
        for i in range(num):
            mark = int(cluster[i, 0])
            plt.plot(data[i, 0], data[i, 1], color[mark] + 'o')

        for i in range(k):
            plt.plot(cp[i, 0], cp[i, 1], color[i] + 'x')
    ##三维图
    elif dim == 3:
        ax = plt.subplot(122, projection='3d')
        for i in range(num):
            mark = int(cluster[i, 0])
            ax.scatter(data[i, 0], data[i, 1], data[i, 2], c=color[mark])
            ax.set_title('k-means result')
        for i in range(k):
            ax.scatter(cp[i, 0], cp[i, 1], cp[i, 2], c=color[i], marker='x')

    plt.show()
if __name__ == "__main__":
        data = load_iris()#以字典形式加载鸢尾花数据集
        y = data.target #使用y表示数据集中的标签
        x = data.data #使用x表示数据集中的属性数据
        #使用PCA 算法，设置降维后主成分数目为 2
        #print(x,'\n', y)
        #print(x)
        #print(type(x))
        size = 20  ##取值范围
        pca = PCA(n_components=3)
        #对原始数据进行降维，保存在 reduced_X 中
        reduced_X = pca.fit_transform(x)
        print('降维后的数据为:\n',reduced_X)#降维后的数据
        print('各主成分方差解释度为：',pca.explained_variance_ratio_)#方差解释度
        print('主成分对应的载荷矩阵为',pca.components_)
        red_x, red_y , red_z = [], [],[]
        blue_x, blue_y , blue_z= [], [],[]
        green_x, green_y, green_z=[],[],[]

        for i in range(len(reduced_X)):
            #标签为0时，3维标签数据保存到列表red_x,red_y,redz中
            if y[i] == 0:
                red_x.append(reduced_X[i][0])
                red_y.append(reduced_X[i][1])
                red_z.append(reduced_X[i][2])#
            elif y[i] == 1:
                blue_x.append(reduced_X[i][0])
                blue_y.append(reduced_X[i][1])
                blue_z.append(reduced_X[i][2])#
            else:
                green_x.append(reduced_X[i][0])
                green_y.append(reduced_X[i][1])
                green_z.append(reduced_X[i][2])#
        X = reduced_X[:, :3]  # #表示我们取特征空间中的3个维度
       #print(X.shape)
       #print(X)
       #print(type(X))
        num = 50 ##点个数
        k=3 ##分类个数
        data = X
        cp,cluster = KMeans(data,k)
        ax=plt.figure().add_subplot(121,projection='3d')
        ax.scatter(red_x, red_y,red_z,c='r', marker='o')
        ax.scatter(blue_x, blue_y,blue_z,c='b', marker='o')
        ax.scatter(green_x, green_y,green_z,c='g', marker='o')#散点图中用s，其余图用markersize可调节散点的大小
        ax.set_title('PCA result')
        Show(data,k,cp,cluster)

本文使用sklearn的KMeans方法对无标签数据集，天气数据集minute_weather进行处理，同时尝试自己编写Kmeans算法进行数据处理，对比两者差异，检查自定义算法可用性。
数据集minute_weather：

链接：https://pan.baidu.com/s/1Ko6YK2xJNiDRsq2befcYsQ

提取码：wwww
说明
本次实验所使用的minute_weather数据集并不“干净”，所以先进行了数据清洗，其中“垃圾”数据包括空值数据，整列或整行为0的数据，最终由于数据过大无法有效处理，只使用其中1000行数据进行室验。
本次通过聚类后得出的数据可以看出，自己编写的kmeans算法与sklearn自带的kmeans算法所得聚类情况相差不大，由可视化图展示的结果来看，数据的聚类情况并不是很明显，分类之间无明显间隔，当然，这也与我们的可视化方式有关，若我采用其他两个参数作为[x,y]进行展示，或许会更加明显。
代码实现
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn.cluster import KMeans

path = 'D:\myData\code\python\data\\minute_weather.csv'
pdData = pd.read_csv(path)
pdData.head()




  

    

      

      
rowID
      
hpwren_timestamp
      
air_pressure
      
air_temp
      
avg_wind_direction
      
avg_wind_speed
      
max_wind_direction
      
max_wind_speed
      
min_wind_direction
      
min_wind_speed
      
rain_accumulation
      
rain_duration
      
relative_humidity
    
  
  

    

      
0
      
0
      
2011-09-10 00:00:49
      
912.3
      
64.76
      
97.0
      
1.2
      
106.0
      
1.6
      
85.0
      
1.0
      
NaN
      
NaN
      
60.5
    
    

      
1
      
1
      
2011-09-10 00:01:49
      
912.3
      
63.86
      
161.0
      
0.8
      
215.0
      
1.5
      
43.0
      
0.2
      
0.0
      
0.0
      
39.9
    
    

      
2
      
2
      
2011-09-10 00:02:49
      
912.3
      
64.22
      
77.0
      
0.7
      
143.0
      
1.2
      
324.0
      
0.3
      
0.0
      
0.0
      
43.0
    
    

      
3
      
3
      
2011-09-10 00:03:49
      
912.3
      
64.40
      
89.0
      
1.2
      
112.0
      
1.6
      
12.0
      
0.7
      
0.0
      
0.0
      
49.5
    
    

      
4
      
4
      
2011-09-10 00:04:49
      
912.3
      
64.40
      
185.0
      
0.4
      
260.0
      
1.0
      
100.0
      
0.1
      
0.0
      
0.0
      
58.8
    
  


数据预处理
pdData.iloc[:,2:].head()




  

    

      

      
air_pressure
      
air_temp
      
avg_wind_direction
      
avg_wind_speed
      
max_wind_direction
      
max_wind_speed
      
min_wind_direction
      
min_wind_speed
      
rain_accumulation
      
rain_duration
      
relative_humidity
    
  
  

    

      
0
      
912.3
      
64.76
      
97.0
      
1.2
      
106.0
      
1.6
      
85.0
      
1.0
      
NaN
      
NaN
      
60.5
    
    

      
1
      
912.3
      
63.86
      
161.0
      
0.8
      
215.0
      
1.5
      
43.0
      
0.2
      
0.0
      
0.0
      
39.9
    
    

      
2
      
912.3
      
64.22
      
77.0
      
0.7
      
143.0
      
1.2
      
324.0
      
0.3
      
0.0
      
0.0
      
43.0
    
    

      
3
      
912.3
      
64.40
      
89.0
      
1.2
      
112.0
      
1.6
      
12.0
      
0.7
      
0.0
      
0.0
      
49.5
    
    

      
4
      
912.3
      
64.40
      
185.0
      
0.4
      
260.0
      
1.0
      
100.0
      
0.1
      
0.0
      
0.0
      
58.8
    
  


weatherdata =pdData.iloc[:,2:]
weatherdata.head()




  

    

      

      
air_pressure
      
air_temp
      
avg_wind_direction
      
avg_wind_speed
      
max_wind_direction
      
max_wind_speed
      
min_wind_direction
      
min_wind_speed
      
rain_accumulation
      
rain_duration
      
relative_humidity
    
  
  

    

      
0
      
912.3
      
64.76
      
97.0
      
1.2
      
106.0
      
1.6
      
85.0
      
1.0
      
NaN
      
NaN
      
60.5
    
    

      
1
      
912.3
      
63.86
      
161.0
      
0.8
      
215.0
      
1.5
      
43.0
      
0.2
      
0.0
      
0.0
      
39.9
    
    

      
2
      
912.3
      
64.22
      
77.0
      
0.7
      
143.0
      
1.2
      
324.0
      
0.3
      
0.0
      
0.0
      
43.0
    
    

      
3
      
912.3
      
64.40
      
89.0
      
1.2
      
112.0
      
1.6
      
12.0
      
0.7
      
0.0
      
0.0
      
49.5
    
    

      
4
      
912.3
      
64.40
      
185.0
      
0.4
      
260.0
      
1.0
      
100.0
      
0.1
      
0.0
      
0.0
      
58.8
    
  


#清洗NUN等异常数据
#去除为空的记录
print("数据的形状为：", weatherdata.shape)
exp1 = weatherdata.notnull()
#exp2 = weatherdata["rain_duration"].notnull()
#exp = exp1 & exp2
weatherdata_notnull = weatherdata.loc[exp1,:]
print("删除缺失记录后数据的形状为：", weatherdata_notnull.shape)
weatherdata_notnull.head()

ValueError: Cannot index with multidimensional key

更换一个方法进行清洗：
#清洗NUN等异常数据
#去除为空的记录
print("数据的形状为：", weatherdata.shape)
#删除表中任何含有NaN的行
weatherdata_notnull = weatherdata.dropna(axis=0, how='any')
print("删除缺失记录后数据的形状为：", weatherdata_notnull.shape)
weatherdata_notnull.head()

数据的形状为： (1587257, 11)
删除缺失记录后数据的形状为： (1586823, 11)




  

    

      

      
air_pressure
      
air_temp
      
avg_wind_direction
      
avg_wind_speed
      
max_wind_direction
      
max_wind_speed
      
min_wind_direction
      
min_wind_speed
      
rain_accumulation
      
rain_duration
      
relative_humidity
    
  
  

    

      
1
      
912.3
      
63.86
      
161.0
      
0.8
      
215.0
      
1.5
      
43.0
      
0.2
      
0.0
      
0.0
      
39.9
    
    

      
2
      
912.3
      
64.22
      
77.0
      
0.7
      
143.0
      
1.2
      
324.0
      
0.3
      
0.0
      
0.0
      
43.0
    
    

      
3
      
912.3
      
64.40
      
89.0
      
1.2
      
112.0
      
1.6
      
12.0
      
0.7
      
0.0
      
0.0
      
49.5
    
    

      
4
      
912.3
      
64.40
      
185.0
      
0.4
      
260.0
      
1.0
      
100.0
      
0.1
      
0.0
      
0.0
      
58.8
    
    

      
5
      
912.3
      
63.50
      
76.0
      
2.5
      
92.0
      
3.0
      
61.0
      
2.0
      
0.0
      
0.0
      
62.6
    
  


#数据太大，处理时间过长，我们提取其中一部分处理
#获取前1000条数据
weatherdata_ture = weatherdata_notnull.iloc[:1000,:]
print("使用的数据形状：",weatherdata_ture.shape)

使用的数据形状： (1000, 11)

sklearn的k-means聚类
X = weatherdata_ture.values
model = KMeans(n_clusters=12)#构造聚类器，12类
model.fit(X)#聚类

KMeans(algorithm='auto', copy_x=True, init='k-means++', max_iter=300,
       n_clusters=12, n_init=10, n_jobs=None, precompute_distances='auto',
       random_state=None, tol=0.0001, verbose=0)

label_pred = model.labels_

可视化
#绘制k-means结果
x0 = X[label_pred == 0]
x1 = X[label_pred == 1]
x2 = X[label_pred == 2]
x3 = X[label_pred == 3]
x4 = X[label_pred == 4]
x5 = X[label_pred == 5]
x6 = X[label_pred == 6]
x7 = X[label_pred == 7]
x8 = X[label_pred == 8]
x9 = X[label_pred == 9]
x10 = X[label_pred ==10]
x11 = X[label_pred == 11]

plt.scatter(x0[:, 0], x0[:, 1], c = "red", marker='o', label='label0')
plt.scatter(x1[:, 0], x1[:, 1], c = "green", marker='*', label='label1')
plt.scatter(x2[:, 0], x2[:, 1], c = "blue", marker='+', label='label2')
plt.scatter(x3[:, 0], x3[:, 1], c = "red", marker='o', label='label3')
plt.scatter(x4[:, 0], x4[:, 1], c = "green", marker='*', label='label4')
plt.scatter(x5[:, 0], x5[:, 1], c = "blue", marker='+', label='label5')
plt.scatter(x6[:, 0], x6[:, 1], c = "red", marker='o', label='label6')
plt.scatter(x7[:, 0], x7[:, 1], c = "green", marker='*', label='label7')
plt.scatter(x8[:, 0], x8[:, 1], c = "blue", marker='+', label='label8')
plt.scatter(x9[:, 0], x9[:, 1], c = "red", marker='o', label='label9')
plt.scatter(x10[:, 0], x10[:, 1], c = "green", marker='*', label='label10')
plt.scatter(x11[:, 0], x11[:, 1], c = "blue", marker='+', label='label11')

plt.xlabel('petal length')
plt.ylabel('petal width')
plt.legend(loc=2)
plt.show()


自定义kmeans算法
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from random import sample


#导入数据
path = 'D:\myData\code\python\data\\minute_weather.csv'
pdData = pd.read_csv(path)
weatherdata =pdData.iloc[:,2:]
#删除表中任何含有NaN的行
weatherdata_notnull = weatherdata.dropna(axis=0, how='any')
#获取前1000条数据
weatherdata_ture = weatherdata_notnull.iloc[:1000,:]

X = weatherdata_ture.values


n = 12 #表示n个聚类
u = sample(X.tolist(),n) #随机选择n个X中的向量作为聚类中心
max_iter = 0 #记录迭代次数
u_before = u
 
while max_iter<12:
    #簇分配过程
    c = []
    #print(u_before,u)
    for i in range(len(X)):
        min = 1000
        index = 0
        for j in range(n):
            vec1 = X[i]
            vec2 = u[j]
            dist = np.sqrt(np.sum(np.square(vec1 - vec2)))
            if dist<min:
                min = dist
                index =j
        c.append(index)
        # print(i,"------",c[i])
    #移动聚类中心
    for j in range(n):
        sum = np.zeros(11)  # 定义n为零向量 此处n为11
        count = 0  # 统计不同类别的个数
        for i in range(len(X)):
            if c[i]==j:
                sum = sum+X[i]
                count = count+1
        u[j] = sum/count
 
    print(max_iter,"------------",u)
    #设置迭代次数
    max_iter = max_iter + 1
print(np.array(c))
label_pred = np.array(c)

#绘制k-means结果
x0 = X[label_pred == 0]
x1 = X[label_pred == 1]
x2 = X[label_pred == 2]
x3 = X[label_pred == 3]
x4 = X[label_pred == 4]
x5 = X[label_pred == 5]
x6 = X[label_pred == 6]
x7 = X[label_pred == 7]
x8 = X[label_pred == 8]
x9 = X[label_pred == 9]
x10 = X[label_pred ==10]
x11 = X[label_pred == 11]

plt.scatter(x0[:, 0], x0[:, 1], c = "red", marker='o', label='label0')
plt.scatter(x1[:, 0], x1[:, 1], c = "green", marker='*', label='label1')
plt.scatter(x2[:, 0], x2[:, 1], c = "blue", marker='+', label='label2')
plt.scatter(x3[:, 0], x3[:, 1], c = "red", marker='o', label='label3')
plt.scatter(x4[:, 0], x4[:, 1], c = "green", marker='*', label='label4')
plt.scatter(x5[:, 0], x5[:, 1], c = "blue", marker='+', label='label5')
plt.scatter(x6[:, 0], x6[:, 1], c = "red", marker='o', label='label6')
plt.scatter(x7[:, 0], x7[:, 1], c = "green", marker='*', label='label7')
plt.scatter(x8[:, 0], x8[:, 1], c = "blue", marker='+', label='label8')
plt.scatter(x9[:, 0], x9[:, 1], c = "red", marker='o', label='label9')
plt.scatter(x10[:, 0], x10[:, 1], c = "green", marker='*', label='label10')
plt.scatter(x11[:, 0], x11[:, 1], c = "blue", marker='+', label='label11')

plt.xlabel('petal length')
plt.ylabel('petal width')
plt.legend(loc=2)
plt.show()


可以看出与我们sklearn分类布局基本一致
以下为kmeans算法的错误尝试：
1.取k个随机质心：
def randCent(dataSet, k):
  n = np.shape(dataSet)[1]
  centroids = np.mat(np.zeros((k,n)))  # 每个质心有n个坐标值，总共要k个质心
  for j in range(n):
    minJ = min(dataSet[:,j])
    maxJ = max(dataSet[:,j])
    rangeJ = float(maxJ - minJ)
    centroids[:,j] = minJ + rangeJ * np.random.rand(k, 1)
  return centroids

2.计算距离（欧几里得）：
def distEclud(vecA, vecB):
  #return math.sqrt(sum(np.power(vecA - vecB, 2))) # 求两个向量之间的距离
  return np.sqrt(sum((vecA - vecB) ** 2)) # 求两个向量之间的距离

  File "<tokenize>", line 5
    return np.sqrt(sum((A - B) ** 2)) # 求两个向量之间的距离
    ^
IndentationError: unindent does not match any outer indentation level

3.归类并计算蔟中心：
def kMeans(dataSet, k, distMeans =distEclud, createCent = randCent):
  m = np.shape(dataSet)[0]
  clusterAssment = np.mat(np.zeros((m,2)))  # 用于存放该样本属于哪类及质心距离
  # clusterAssment第一列存放该数据所属的中心点，第二列是该数据到中心点的距离
  centroids = createCent(dataSet, k)
  clusterChanged = True  # 用来判断聚类是否已经收敛
  while clusterChanged:
    clusterChanged = False;
    for i in range(m): # 把每一个数据点划分到离它最近的中心点
      minDist = np.inf #表示无穷大,若两点间的距离小于MinDist;
      minIndex = -1;
      for j in range(k):
        distJI = distMeans(centroids[j,:], dataSet[i,:])
        if distJI < minDist:
          minDist = distJI; minIndex = j # 如果第i个数据点到第j个中心点更近，则将i归属为j
      if clusterAssment[i,0] != minIndex: clusterChanged = True; # 如果分配发生变化，则需要继续迭代
      clusterAssment[i,:] = minIndex,minDist**2  # 并将第i个数据点的分配情况存入字典
    print(centroids)
    for cent in range(k):  # 重新计算中心点
      ptsInClust = dataSet[nonzero(clusterAssment[:,0].A == cent)[0]]  # 去第一列等于cent的所有列
      centroids[cent,:] = mean(ptsInClust, axis = 0) # 算出这些数据的中心点
  return centroids, clusterAssment

带入数据：
import math
print("X.shape: ",np.shape(X)[0])
myCentroids,clustAssing = kMeans(X,12)
print("myCentroids: \n",myCentroids)
print("clustAssing: \n",clustAssing)

LinAlgError: Last 2 dimensions of the array must be square

总结：
Kmeans作为一种发现数据内在结构的聚类技术，在我们以往的实验数据处理过程中经常被使用到，它也被称为无监督学习。通过实验对看means的使用情况可以看出，Kmeans算法的优缺点同样很明显。
优点：

1）原理简单，容易实现

2）可解释度较强

3）使用广泛
缺点：

1）K值很难确定，我们几次实验都是预先知道K的取值。

2）容易局部最优

3）对噪音和异常点敏感

4）需样本存在均值（限定数据种类）

5）聚类效果依赖于聚类中心的初始化

6）计算量大，由于存在大量迭代计算的原因导致kmean算法对大型数据处理不友好，如本次实验15W条数据，等待5分钟依然没有结果，无奈放弃使用小数据进行实验。

用户分类是运营、营销过程中的很基础的一环，合理地对用户进行分类，分别实施合适的运营手段，以提高效率。

本文数据集来自于《python数据分析与挖掘实战（第二版）》，有兴趣的同学买来看看。
先介绍一下RFM和Kmeans:

RFM模型：

R，recently，表示用户最近付费/使用/登录等行为的时间

F，frequency，表示用户付费/使用等行为的频次/频率

M，money，表示用户在一个时间段内的消费金额，如果是免费产品，像抖音，应该可以用时间替换，时间就是金钱啊！
对于游戏来说，这个模型衡量你有没有氪金，玩的多不多，最近有没有玩；对于B站来说，这个模型衡量你冲大会员没有，投币没有，看得多不多；
显然，RFM模型是一个3特征模型，虽然它3个特征已经比较好地刻画出用户价值的模型，也适用于大多数场景，也好进行可视化（三维），但是基于每一种具体的业务场景，肯定有基础该场景的也重要的特征，把这些特征加进来，能更好地分类用户。
把RFM值都计算出来后，问题就来了，怎么划分？如果以每个字段（RFM）的某个特定值，如平均值、某分位数值进行划分，组合2³将用户划分成8类，三维空间也被划分为方方的8个区域；



这是经典的划分方法，后续可以根据不同的划分阈值测试调整,使收益最大化。
那这种方法的缺陷在哪呢？

1.需要人为去设定阈值，有非常丰富的经验还好，但是非常丰富的经验是需要代价的。

2.划分得太“方正”，这种按特定值的划分方式，没有考虑数据的集中程度，对特征也不敏感，比如说一个业务场景以M划分应该要分成3部分，R划分两部分，F划分3部分，这就很难去实现了。
为了弥补上述缺陷，大神们在几十年前发明了Kmeans算法。
Kmeans算法

k，表示k个中心

means，表示均值

不同于按特定值划分，Kmeans用样本到中心点的距离来分类（近墨者黑），以RFM为例,假设初始设置的k=3，即有三个聚类中心，得到的结果是这样的（灵魂画手示意图）

也因为是距离分类，所以我们要把数据消除量纲的影响，比如RFM中M的跨度很大，分成3类的时候很可能就在M一个维度上分成3类，造成对F和M不敏感。



那这些聚类中心是如何生成的呢？一般步骤如下：

1.计算每个样本到聚类中心的距离，跟哪个近就划分到哪个中心；

2.用新划分的样本群体的平均值，计算出新的聚类中心；

3.根据给定的限制条件看看是否终止（如中心的变化量，循环次数等）

4.循环步骤123
这时问题又来了，而且是好多问题：

1.k的数量，怎么确定？

2.初始k的位置，怎么确定？

3.限制条件，怎么确定？
若知后事如何，请见下文分解
航空公司客户分类模型（概览）
目的与背景
企业的资源是有限的，多角度，全方面地分析客户价值，合理分配资源，可增加收益。分析客户价值的方法很多，也根据不同的业务场景进行调整，最常见的是矩阵分析法和RFM模型。
航空公司采集统计了两年的客户数据，通过挖掘分析把客户分类，从而提供提供更精准的服务和运营手段，获得更高的收益。
分析方法和流程
（1）探索性数据分析：

着重于探索数据的分布、完整性、相关性、以及各数据的业务意义。
（2）数据清洗和处理

数据清洗：缺失值，异常值，时间范围等；

数据处理：数据标准化，特征构造和筛选等；
（3）分析建模

Kmeans客户聚类–>模型分析优化–>聚类客户群价值分析
（4）应用和局限分析
探索性数据分析
“第一次听说pandas_profiling，其实我就没拒绝”
#读取数据
air_data = pd.read_csv(r'......\data\air_data.csv')

#使用pandas_profiling生成数据探索报告
import pandas_profiling
air_EDA =pandas_profiling.ProfileReport(air_data)

#profilereport包含了pandas的describe()，head()，tail()，corr()，isnull()，values_counts()等函数，然后用直方图、热力图、表格等形式展现出来。


业务理解

由于是以RFM为基础的思路，所以在挑选数据时，先把RFM的数据挑选出来（下图红框）；其次，还要考虑其他的一些重要特征，如入会日期，他是反映忠诚度的一个指标，还有构造新的特征，平均票价，它反映客户的消费倾向和水平；


数据清洗和处理
在上一步，我们先挑选5个维度的数据：

L：FFP_DATE-入会时长； R：LAST_FIGHT_DATE-最后一次飞行时间  F：L1Y_FLIGHT_COUNT-第二年飞行次数；  M：SUM_YR_2-第二年总消费， 以及构造的新特征 A：SUM_YR_2/L1Y_FLIGHT_COUNT - 平均票价
所以，只针对这五个维度进行数据清洗和处理就可以了
#提取数据
df_LRFMA = air_data[['LOAD_TIME','FFP_DATE','LAST_TO_END','L1Y_Flight_Count','SUM_YR_2']]

#数据处理，构造平均票价字段
df_LRFMA = df_LRFMA.assign(
                           L = (
                                       (
                                       pd.to_datetime(df_LRFMA['LOAD_TIME'])
                                        - pd.to_datetime(df_LRFMA['FFP_DATE'])
                                       )/30
                                   ).map(lambda x : int(x.days)),
                           R = pd.to_numeric(df_LRFMA['LAST_TO_END']),
                           F = df_LRFMA['L1Y_Flight_Count'],
                           M = df_LRFMA['SUM_YR_2'],
                           A = df_LRFMA['SUM_YR_2']/df_LRFMA['L1Y_Flight_Count'])
#提取LRFMA数据                          
df_LRFMA = df_LRFMA[['L','R','F','M','A']]

#填充缺失值，这里的缺失值都在A上，是由0/0造成的，所以用0填充，表示平均票价是0
df_LRFMA = df_LRFMA.fillna(0)

#数据标准化，消除量纲的影响
S_LRFMA = StandardScaler().fit_transform(df_LRFMA)

分析建模
使用sklearn.cluster中的KMeans，上文提出的问题：如何确认初始值K？

sklearn 的kmeans中提供一个属性：inertia_，表示所有数据点到它所属的聚类中心的距离的总和；

随着k的增大，inertia逐渐减少，当k等于数据量n时，inertia=0

在k增大的过程中，当inertia的减少量不明显时，可以认为是k的增加把“密集”的分开了，而这是不必要的，所以k的选择要根据这个“拐点”和实际业务逻辑来确定
from sklearn.cluster import KMeans
#遍历k，找合适的k值
inertia = []
for i in range(1,15):
    model = KMeans(n_clusters=i,n_jobs=4)
    model.fit(S_LRFMA)
    inertia.append(model.inertia_)



根据数据表现，k=5或k=6都相对合理，另一方面，RFM可把用户价值分成8类，所以我选择了更接近8的6。
k= 6
model = KMeans(n_clusters=k,n_jobs=4,random_state =42)
model.fit(S_LRFMA)

kmeans_cc =model.cluster_centers_

client = pd.DataFrame(kmeans_cc.round(3),columns = ['L','R','F','M','A'])
client = pd.merge(client,group,left_index=True,right_index=True,how='left')



为各类用户群“验明正身”后，还要统计其数量，以及作进一步的描述；
#计算各群体个数
group = pd.Series(model.labels_,name ='客户数').value_counts()
#合并数据，并加上描述
client = pd.DataFrame(kmeans_cc.round(3),columns = ['L','R','F','M','A'])
client = pd.merge(client,group,left_index=True,right_index=True,how='left')

client = client.assign(用户属性 = ['VIP','流失客户','新/价值低用户','一般用户','VVIP','高消费用户'],
                       描述 = [
                            '入会时间长，飞行频率偏高，消费方面表现良好，需重点维护的VIP客户群体',
                            '各项指标都偏低，最后一次飞行时间在很久之前的用户群体',
                            '入会时间较短，最近有飞行记录，其他指标表现均最低',
                           '入会时间长，其他指标表现一般,为一般客户',
                            '各项指标表现均优秀，是公司忠实的，高价值的用户',
                            '消费频次较低，但消费金额大，需重点保持发展'])



对于kmeans的初始k位置，限制条件如何设置和确定，详细请看kmeans中的参数：tol, max_iter, k-means++
模型应用
上述模型以LRFMA这五个角度把用户分成6类，但在每一类中，还可以进行分类来实施更细致的服务。对于这个模型有以下几点的应用：

1.监控业务：

定时计算每个群体的用户数占比，如流失用户比例增大，这就需要去发掘问题的所在点；又如发现聚类点发生比较大的改变，也要去深究原因。

2.差异化营销：

一个用户的消费倾销，都刻在他所产生的数据上，根据不同的用户分群，提供不同的服务/优惠/活动，增加用户的忠诚度，提高消费比例。