采用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')