现在有一个存有70个地址和城市名的文本，而没有这些地点的距离信息，而我们想要对这些地点进行聚类，找到每个簇的质心地点，从而可以安排合理的行程，即不同簇中的地点之间选择交通工具抵达，而位于同一个簇内的地点之间可以采取步行的方法抵达。使用Kmeans算法可以为我们找到一种更加经济而且高效的出行方式。
数据说明：东经取经度的正值(Longitude)，西经取经度负值(-Longitude)，北纬取90-纬度值(90-Latitude),南纬取90+纬度值(90+Latitude)，则经过上述处理过后的两点被计为(MLonA, MLatA)和(MLonB, MLatB)。地点对应的经纬度经过以上处理过后，保存在places.txt文件和places.csv中的最后两列。
使用球面余弦定理来计算两个经纬度之间的实际距离
地球平均半径6371千米。如果我们假设地球半径为R。设A点的经纬度为(LonA, LatA)，B点的经纬度为(LonB, LatB)，两点距离的如下公式：
- a = sin(LatA*Pi/180)*sin(LatB*Pi/180)
- b = cos(LatA*Pi/180)*cos(LatB*Pi/180)*cos((MLonA-MLonB)*Pi/180)
- Distance = R* arccos(a + b)*Pi/180
# 计算根据经纬度计算两点之间的球面距离
def distSLC(vacA, vacB):
a = sin(vacA[0, 1]*np.pi/180) * sin(vacB[0, 1]*np.pi/180)
b = cos(vacA[0, 1]*np.pi/180) * cos(vacB[0, 1]*np.pi/180)*cos(np.pi*(vacB[0, 0]-vacA[0, 0])/180)
return acos(a+b)*6371.0
K值的确定
本例使用拐点法，就是在不同的k值下计算簇内离差平方和，然后通过可视化的方法找到”拐点”所对应的k值，J为Kmeans算法的目标函数，随着簇数量的增加，簇中的样本量会越来越少，进而导致目标函数J的值也会越来越小，通过可视化方法，重点关注的是斜率的变化，当斜率由大突然变小时，并且之后的斜率变化缓慢，则认为突然变化的点就是寻找的目标点，因为继续随着簇数K的增加，聚类效果不再有大的变化。
完整代码及具体实现
import numpy as np
import pandas as pd
from math import cos, sin, acos
import matplotlib.pyplot as plt
# 加载数据
def loaddATEsET(filename):
data = pd.read_csv(filename, header=None)
dataArr = data.iloc[:, -2:].values
return np.mat(dataArr)
# 计算根据经纬度计算两点之间的球面距离
def distSLC(vacA, vacB):
a = sin(vacA[0, 1]*np.pi/180) * sin(vacB[0, 1]*np.pi/180)
b = cos(vacA[0, 1]*np.pi/180) * cos(vacB[0, 1]*np.pi/180)*cos(np.pi*(vacB[0, 0]-vacA[0, 0])/180)
return acos(a+b)*6371.0
# 随机初始化K个质心(质心满足数据边界之内)
def randCent(dataSet, k):
n = np.shape(dataSet)[1]
centroids = np.mat(np.zeros((k, n)))
for j in range(n):
minJ = np.min(dataSet[:, j])
rangeJ = float(max(dataSet[:, j])-minJ)
centroids[:, j] = minJ+rangeJ*np.random.rand(k, 1)
return centroids
# Kmeans聚类算法
def kMeans(dataSet, k, distMeas=distSLC, createCent=randCent):
m = np.shape(dataSet)[0]
clusterAssment = np.mat(np.zeros((m, 2)))
centroids = createCent(dataSet, k)
clusterChanged = True
while clusterChanged:
clusterChanged = False
for i in range(m):
minDist = np.inf
minIndex = -1
for j in range(k):
distJI = distMeas(centroids[j, :], dataSet[i, :])
if distJI < minDist:
minDist = distJI
minIndex = j
if clusterAssment[i, 0] != minIndex:
clusterChanged = True
clusterAssment[i, :] = minIndex, minDist**2
# print(centroids)
for cent in range(k):
ptsInClust = dataSet[np.nonzero(clusterAssment[:, 0].A == cent)[0]]
if len(ptsInClust) != 0:
centroids[cent, :] = np.mean(ptsInClust, axis=0)
return centroids, clusterAssment
# 计算Kmeans聚类算法的SSE值
def calcSSE(clusterAssment):
return sum(clusterAssment[:, -1])
# 对地图坐标进行聚类，并在地图上显示聚类结果
def plotKMeans(dataSet, clusterAssment, centroids, k):
fig = plt.figure()
rect = [0.1, 0.1, 0.8, 0.8]
scatterMarkers = ['s', 'o', '^', '8', 'p', 'd', 'v', 'h', '>', '<']
axprops = dict(xticks=[], yticks=[])
ax0 = fig.add_axes(rect, label='ax0', **axprops)
imgP = plt.imread('Portland.png')
ax0.imshow(imgP)
ax1 = fig.add_axes(rect, label='ax1', frameon=False)
for i in range(k):
ptsInCurrCluster = dataSet[np.nonzero(clusterAssment[:, 0].A == i)[0], :]
markerStyle = scatterMarkers[i % len(scatterMarkers)]
ax1.scatter(ptsInCurrCluster[:, 0].flatten().A[0],
ptsInCurrCluster[:, 1].flatten().A[0],
marker=markerStyle, s=90)
ax1.scatter(centroids[:, 0].flatten().A[0], centroids[:, 1].flatten().A[0],
marker='+', s=300)
plt.show()
# plt.savefig('fig.png', bbox_inches='tight')
dataSet = loaddATEsET("places.csv")
for i in range(7):
centroids, clusterAssment = kMeans(dataSet, i+2)
sseValue = calcSSE(clusterAssment)
print("当k值为：%d时，对应的SSE值为：%f" % (i+2, sseValue))
centroids, clusterAssment = kMeans(dataSet, 3)
plotKMeans(dataSet, clusterAssment, centroids, 3)
运行结果与分析
当k值分别为k=2,3,4,5,6,7,8时，计算出对应的SSE值如下：
| K取值 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
| SSE值 | 3131.80 | 1819.44 | 1649.26 | 1077.58 | 949.05 | 762.87 | 665.26 |
折线图如下：
结论：
由上图可见，本例的K值确定为3。
效果图如下：
