现在有一个存有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。
效果图如下:

请问有数据集的文件嘛
@1654007348475: 链接:https://pan.baidu.com/s/1h_LQQrTSxYcVKUc2RpqCKg
提取码:2022
imgP = plt.imread(‘Portland.png’)
问下’Portland.png’ 是怎么生成的?
@1664528592856: 图是原来就有的,只不过把坐标标出来了而已,这个案例来自科班的机器学习课程
return acos(a+b)*6371.0
a+b可能为负报错:
return acos(a + b) * 6371.0
ValueError: math domain error
# 计算Kmeans聚类算法的SSE值
def calcSSE(clusterAssment):
matrix = sum(clusterAssment[:, -1])
# squeeze the matrix into a vector
return matrix[0, 0]
# 参考:https://blog.csdn.net/Nafisa06/article/details/105648820
def distSLC2(vacA, vacB):
lat1 = vacA[0, 1]
lon1 = vacA[0, 0]
lat2 = vacB[0, 1]
lon2 = vacB[0, 0]
lat1 = lat1 * np.pi / 180
lon1 = lon1 * np.pi / 180
lat2 = lat2 * np.pi / 180
lon2 = lon2 * np.pi / 180
dis = np.arccos(round((np.sin(lat1) * np.sin(lat2)) + (np.cos(lat1) * np.cos(lat2) * np.cos(lon2 – lon1)), 10)) * 6371.004
return dis