k-means聚类算法的离群点检测

采用k-means算法进行离群点检测，并包含通过SSE算法找到最佳聚类数、PCA降维方法。

import pandas as pd
import numpy as np
import matplotlib.dates as md
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans

data1= pd.read_csv("设备1.csv",encoding = 'ANSI',engine='python' )

#找到合适的降维特征数
X = data1.values
#标准化
X_std = StandardScaler().fit_transform(X)
mean_vec = np.mean(X_std, axis=0)
#特征值和特征向量的提取
cov_mat = np.cov(X_std.T)## 
eig_vals, eig_vecs = np.linalg.eig(cov_mat)
#特征值与特征向量一一对应
eig_pairs = [ (np.abs(eig_vals[i]),eig_vecs[:,i]) for i in range(len(eig_vals))]
#特征值从高到低排序
eig_pairs.sort(key = lambda x: x[0], reverse= True)
#特征值归一化
tot = sum(eig_vals)
var_exp = [(i/tot)*100 for i in sorted(eig_vals, reverse=True)] # Individual explained variance
cum_var_exp = np.cumsum(var_exp) # Cumulative explained variance
K=0
#这里我设置保存90的主成分，得到的k值为5，或者通过图像观察选择k值，输出保留成分比例。
for a in cum_var_exp:
    K = K + 1
    if a>=90:
        break 
#特征值的可视化
plt.figure(figsize=(10, 5))
plt.bar(range(len(var_exp)), var_exp, alpha=0.3, align='center', label='individual explained variance', color = 'g')
plt.step(range(len(cum_var_exp)), cum_var_exp, where='mid',label='cumulative explained variance')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal components')
plt.legend(loc='best')
plt.show()

#降维到5个特征数 K=5
X_std = StandardScaler().fit_transform(X)
data = pd.DataFrame(X_std)

pca = PCA(n_components=K)
data = pca.fit_transform(data)
#对降维后的5个特征的数据进行k-means异常检测
#通过SSE方法获取最佳的聚类数k
#求二阶导数最大的拐点为聚类数，保存模型bestmodel
SSE = []
SSE_d1 = [] #sse的一阶导数
SSE_d2 = [] #sse的二阶导数
for k in range(1, 20):
    kmeans_model = KMeans(n_clusters=k )
    kmeans_model.fit(data)
    SSE.append(kmeans_model.inertia_)  # 保存每一个k值的SSE值
X = range(1,20)
plt.xlabel('k')
plt.ylabel('SSE')
plt.plot(X, SSE, 'o-')
print(SSE)
plt.show()

由图可知，取手肘部差不多k=4或K=5即可。

i=0
b=[]

# 为了可视化，只取PCA降维的2个特征
#标准化
scaler = StandardScaler()
np_scaled = scaler.fit_transform(data)
data = pd.DataFrame(np_scaled)
#聚类数为观察得到的4
kmeans1 = KMeans(n_clusters=4 ).fit(data)
df['cluster'] = kmeans1.predict(data)
df.index = data.index
#降维后的两个特征
df['principal_feature1'] = data[0]
df['principal_feature2'] = data[1]
df['cluster'].value_counts()
#计算距离的函数
def getDistanceByPoint(data, model):
   distance = pd.Series()
   for i in range(0,len(data)):
       Xa = np.array(data.loc[i])
       Xb = model.cluster_centers_[model.labels_[i]-1]
       distance.set_value(i, np.linalg.norm(Xa-Xb))
   return distance
 #设置异常数据比例
outliers_fraction = 0.004
#  计算每个点与其最近的质心点之间的距离，最大的距离被认为是异常的。
distance = getDistanceByPoint(data,kmeans1)
number_of_outliers = int(outliers_fraction*len(distance))
threshold = distance.nlargest(number_of_outliers).min()
# 异常检测结果anomaly1包含了上述方法（0：正常，1：异常）
df['anomaly1'] = (distance >= threshold).astype(int)

#使用聚类视图可视化异常点
fig, ax = plt.subplots(figsize=(10,6))
colors = {0:'blue', 1:'red'}
ax.scatter(df['principal_feature1'], df['principal_feature2'], c=df["anomaly1"].apply(lambda x: colors[x]))
plt.xlabel('principal feature1')
plt.ylabel('principal feature2')
plt.show();

#保存异常数据
data=df.loc[df['anomaly1']==1]
data.to_csv("ckmianamoly.csv",na_rep="NULL" , encoding='ANSI')