数学建模及机器学习算法（一）：聚类-kmeans（Python及MATLAB实现，包括k值选取与聚类效果评估）

                 1
                #
                 -*- coding: utf-8 -*-
                 2
                import
                 math

                 3
                import
                 random

                 4
                import
                 numpy as np

                 5
                import
                 matplotlib.pyplot as plt

                 6
                 7
                class
                 Canopy:

                 8
                def
                __init__
                (self, dataset):

                 9         self.dataset = dataset
10         self.t1 = 0
11         self.t2 = 0
1213# 设置初始阈值14def setThreshold(self, t1, t2):
15if t1 > t2:
16             self.t1 = t1
17             self.t2 = t2
18else:
19print(‘t1 needs to be larger than t2!‘)
2021# 使用欧式距离进行距离的计算22def euclideanDistance(self, vec1, vec2):
23return math.sqrt(((vec1 - vec2)**2).sum())
2425# 根据当前dataset的长度随机选择一个下标26def getRandIndex(self):
27return random.randint(0, len(self.dataset) - 1)
2829def clustering(self):
30if self.t1 == 0:
31print(‘Please set the threshold.‘)
32else:
33             canopies = []  # 用于存放最终归类结果34while len(self.dataset) != 0:
35                 rand_index = self.getRandIndex()
36                 current_center = self.dataset[rand_index]  # 随机获取一个中心点，定为P点37                 current_center_list = []  # 初始化P点的canopy类容器38                 delete_list = []  # 初始化P点的删除容器39                 self.dataset = np.delete(
40                     self.dataset, rand_index, 0)  # 删除随机选择的中心点P41for datum_j in range(len(self.dataset)):
42                     datum = self.dataset[datum_j]
43                     distance = self.euclideanDistance(
44                         current_center, datum)  # 计算选取的中心点P到每个点之间的距离45if distance < self.t1:
46# 若距离小于t1，则将点归入P点的canopy类47                        current_center_list.append(datum)
48if distance < self.t2:
49                         delete_list.append(datum_j)  # 若小于t2则归入删除容器50# 根据删除容器的下标，将元素从数据集中删除51                 self.dataset = np.delete(self.dataset, delete_list, 0)
52                canopies.append((current_center, current_center_list))
53return canopies
545556def showCanopy(canopies, dataset, t1, t2):
57     fig = plt.figure()
58     sc = fig.add_subplot(111)
59     colors = [‘brown‘, ‘green‘, ‘blue‘, ‘y‘, ‘r‘, ‘tan‘, ‘dodgerblue‘, ‘deeppink‘, ‘orangered‘, ‘peru‘, ‘blue‘, ‘y‘, ‘r‘,
60‘gold‘, ‘dimgray‘, ‘darkorange‘, ‘peru‘, ‘blue‘, ‘y‘, ‘r‘, ‘cyan‘, ‘tan‘, ‘orchid‘, ‘peru‘, ‘blue‘, ‘y‘, ‘r‘, ‘sienna‘]
61     markers = [‘*‘, ‘h‘, ‘H‘, ‘+‘, ‘o‘, ‘1‘, ‘2‘, ‘3‘, ‘,‘, ‘v‘, ‘H‘, ‘+‘, ‘1‘, ‘2‘, ‘^‘,
62‘<‘, ‘>‘, ‘.‘, ‘4‘, ‘H‘, ‘+‘, ‘1‘, ‘2‘, ‘s‘, ‘p‘, ‘x‘, ‘D‘, ‘d‘, ‘|‘, ‘_‘]
63for i in range(len(canopies)):
64         canopy = canopies[i]
65         center = canopy[0]
66         components = canopy[1]
67         sc.plot(center[0], center[1], marker=markers[i],
68                 color=colors[i], markersize=10)
69         t1_circle = plt.Circle(
70             xy=(center[0], center[1]), radius=t1, color=‘dodgerblue‘, fill=False)
71         t2_circle = plt.Circle(
72             xy=(center[0], center[1]), radius=t2, color=‘skyblue‘, alpha=0.2)
73        sc.add_artist(t1_circle)
74        sc.add_artist(t2_circle)
75for component in components:
76             sc.plot(component[0], component[1],
77                     marker=markers[i], color=colors[i], markersize=1.5)
78     maxvalue = np.amax(dataset)
79     minvalue = np.amin(dataset)
80     plt.xlim(minvalue - t1, maxvalue + t1)
81     plt.ylim(minvalue - t1, maxvalue + t1)
82    plt.show()
838485if__name__ == "__main__":
86     dataset = np.random.rand(500, 2)  # 随机生成500个二维[0,1)平面点87     t1 = 0.6
88     t2 = 0.4
89     gc = Canopy(dataset)
90    gc.setThreshold(t1, t2)
91     canopies = gc.clustering()
92print(‘Get %s initial centers.‘ % len(canopies))
93 showCanopy(canopies, dataset, t1, t2)

                 1
                import
                 scipy

                 2
                from  scipy.spatial.distance import euclidean
 3from sklearn.cluster import KMeans as k_means
 4 5 dst = euclidean
 6 7 k_means_args_dict = {
 8‘n_clusters‘: 0,
 9# drastically saves convergence time10‘init‘: ‘k-means++‘,
11‘max_iter‘: 100,
12‘n_init‘: 1,
13‘verbose‘: False,
14# ‘n_jobs‘:815}
161718def gap(data, refs=None, nrefs=20, ks=range(1, 11)):
19"""20    I: NumPy array, reference matrix, number of reference boxes, number of clusters to test
21    O: Gaps NumPy array, Ks input list
2223    Give the list of k-values for which you want to compute the statistic in ks. By Gap Statistic
24    from Tibshirani, Walther.
25"""26     shape = data.shape
2728ifnot refs:
29         tops = data.max(axis=0)
30         bottoms = data.min(axis=0)
31         dists = scipy.matrix(scipy.diag(tops - bottoms))
32         rands = scipy.random.random_sample(size=(shape[0], shape[1], nrefs))
33for i in range(nrefs):
34             rands[:, :, i] = rands[:, :, i] * dists + bottoms
35else:
36         rands = refs
3738     gaps = scipy.zeros((len(ks),))
3940for (i, k) in enumerate(ks):
41         k_means_args_dict[‘n_clusters‘] = k
42         kmeans = k_means(**k_means_args_dict)
43        kmeans.fit(data)
44         (cluster_centers, point_labels) = kmeans.cluster_centers_, kmeans.labels_
4546         disp = sum(
47             [dst(data[current_row_index, :], cluster_centers[point_labels[current_row_index], :]) for current_row_index
48in range(shape[0])])
4950         refdisps = scipy.zeros((rands.shape[2],))
5152for j in range(rands.shape[2]):
53             kmeans = k_means(**k_means_args_dict)
54            kmeans.fit(rands[:, :, j])
55             (cluster_centers, point_labels) = kmeans.cluster_centers_, kmeans.labels_
56             refdisps[j] = sum(
57                 [dst(rands[current_row_index, :, j], cluster_centers[point_labels[current_row_index], :]) for58                  current_row_index in range(shape[0])])
5960# let k be the index of the array ‘gaps‘61         gaps[i] = scipy.mean(scipy.log(refdisps)) - scipy.log(disp)
6263return ks, gaps