Python K-均值聚类实验

        利用某地2021年、2022年的3月、6月、9月、12月份的逐日平均气温T、相对湿度RH、位势高度Z850的再分析资料构建训练样本X=[T,RH,Z850]T。利用K-Means算法实现对基于这三个要素构成的特征空间的聚类分析。

1、利用Scikit-Learn实现K-Means聚类

        取K=2，3，4，可视化聚类效果

代码如下：

import numpy as np
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

data = np.load('data_experiment3.npy')
n_clusters = [4]

def fit_plot_kmean_model(n_clusters,X):
    kmean = KMeans(n_clusters=n_clusters)
    kmean.fit_predict(X)
    markers = ['o','^','*','s']
    colors= ['#440453',  '#30688D',  '#32B67B',  '#FDE73A']
    labels = kmean.labels_
    centers = kmean.cluster_centers_
    plt.title("k={}".format(n_clusters))
    for c in range(n_clusters):
        cluster = X[labels==c]
        ax.scatter(cluster[:,0],cluster[:,1],cluster[:,2],s=20,c=colors[c])
    ax.scatter(centers[:,0],centers[:,1],centers[:,2],marker='o',facecolors='none', edgecolors='red',s=300)
    for i,c in enumerate(centers):
        ax.scatter(c[0],c[1],c[2],marker='x',s=50,c=colors[i])
    ax.set_xlabel('T')
    ax.set_ylabel('RH')
    ax.set_zlabel('Z850')


for i,c in enumerate(n_clusters):
    fig = plt.figure(figsize=(8, 8),dpi=144)
    ax = fig.add_subplot(111, projection='3d')
    fit_plot_kmean_model(c,data)
plt.show()

二、用Numpy或PyTorch实现K-means 2.0

        先从训练集中随机采样𝐾样本作为初始簇中心点𝝁1 , 𝝁2, … , 𝝁k，然后寻找距离样本𝒙 𝑖 最近的簇中心的索引，最后选择代价函数最小的𝑐 1 , 𝑐 2 , … , 𝑐 𝑚 , 𝝁1 , 𝝁2, … , 𝝁k。代价函数主要使用矩阵的范数来计算。

代码如下：

import numpy as np
import matplotlib.pyplot as plt

class KMeans:
    def __init__(self,k):
        self.k = k

    def train(self,X,MAXITER = 100, TOL = 1e-3):
        centroids = np.random.rand(self.k,X.shape[1])
        centroidsold = centroids.copy()
        for iter_ in range(MAXITER):
            dist = np.linalg.norm(X - centroids[0,:],axis=1).reshape(-1,1)
            for class_ in range(1,self.k):
                dist = np.append(dist,np.linalg.norm(X - centroids[class_,:],axis=1).reshape(-1,1),axis=1)
            classes = np.argmin(dist,axis=1)
            # update position
            for class_ in set(classes):
                centroids[class_,:] = np.mean(X[classes == class_,:],axis=0)
            if np.linalg.norm(centroids - centroidsold) < TOL:
                break
                print('Centroid converged')
        self.centroids = centroids
    
    def predict(self,X):
        dist = np.linalg.norm(X - self.centroids[0,:],axis=1).reshape(-1,1)
        for class_ in range(1,self.k):
            dist = np.append(dist,np.linalg.norm(X - self.centroids[class_,:],axis=1).reshape(-1,1),axis=1)
        classes = np.argmin(dist,axis=1)
        return classes


X= -0.5 + np.random.rand(100,3)
X1 = 0.5 + np.random.rand(33,3)
X2 = 2 + np.random.rand(33,3)
X[33:66, :] = X1
X[67:, :] = X2


from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize = (8,5))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(X[:,0],X[:,1],X[:,2])

kmeans = KMeans(3)
kmeans.train(X)
kmeans.centroids
classes = kmeans.predict(X)
classes
fig = plt.figure(figsize = (8,5))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(X[classes == 0,0],X[classes == 0,1],X[classes == 0,2])
ax.scatter(X[classes == 1,0],X[classes == 1,1],X[classes == 1,2])
ax.scatter(X[classes == 2,0],X[classes == 2,1],X[classes == 2,2])
plt.show()