发表自话题：k曲线怎么看

本文开始，将启动新的系列文章，探索一个有趣的话题：异常检测。这是数据科学领域中的核心问题之一，也是一个炙手可热的话题。

本文使用之前的系列文章中用过的物联网数据集，继续通过Python语言，借助数理统计方法、聚类算法或者专门的异常检测算法，进行异常点检测，并分析效果，试图找到场景合适、效率高、准确度高的异常检测方法。这里先列出找到的14个异常检测算法：

箱体图法

正态分布法

DBSCAN算法

局部离群因子LOF检测法

孤立森林Isolation Forest算法

One Class SVM

Robust covariance协方差法

Feature bagging方法

Histogram-Based OutlierDetection方法

Angle-Based OutlierDetection方法

主成分分析方法

K近邻KNN方法

最小协方差法

Z-Score方法

系列文章将逐步探索上述算法在目标物联网数据集上的使用效果。

说到离群点和异常点，翻遍数据科学、可视化、数理统计等领域的著作或者学术文章，也没个精确的定义。不过，感谢语言本身的魅力，它完美的解释了这个问题。不管是中文的离群点和异常点，还是英文的Outlier和Anomaly都非常形象地描述了这个问题。

其实，从哲学的角度分析，离群点也好，异常点也好，都是相对的。当评判的标准，或者群和正常数据的数量或者数值发生变化，离群点和异常点都可能有所不同。接下来，就进入正题，逐一尝试这些异常检测方法。

箱体图法

基于箱体图，详细内容参见链接：Python这波神操作你打几分？物联网数据 — 描述性分析(2)，可以定义某个数值型变量中的异常点和极端异常点，它们的判断表达式(其中Q3和Q1为上下四分位，1.5为可调参数)：

按照如上标准，异常点分别为427个，479个，0个，1009个，0个。绘图如下：

修改标准1.5为2，异常点分别为300个，133个，0个，776个，0个。绘图如下：

修改标准1.5为3，异常点分别为5个，0个，0个，38个，0个。绘图如下：

修改标准1.5为4，异常点分别为2个，0个，0个，3个，0个。绘图如下：

箱体图的优点是计算快，也能明显识别异常点。但不能显示异常点的值和ID，而且只能计算一维数据。

本文不再给出异常值的列表，请使用如下代码自行体验：

import pandas as pdimport matplotlib.pyplot as plt# 导入数据dataset1 =pd.read_csv('../input/dataset2.csv',index_col='date',parse_dates=['date'])def print_outliers(dataset1,ColumnName,whis):
    sunspots=dataset1[ColumnName]# 计算下四分位数和上四分位Q1 = sunspots.(q = 0.25)
    Q3 = sunspots.(q = 0.75)# 基于1.5倍的四分位差计算上下须对应的值low_whisker = Q1 - whis*(Q3 - Q1)
    up_whisker = Q3 + whis*(Q3 - Q1)
    cnts = dataset1[(dataset1[ColumnName]>= up_whisker) | (dataset1[ColumnName] 'markerfacecolor':'red', 'markeredgecolor':'red', 'markersize':4},# 指定均值点的标记符号(菱形)、填充色和大小meanprops = {'marker':'D','markerfacecolor':'black', 'markersize':4},medianprops = {'linestyle':'--','color':'orange'}, # 指定中位数的标记符号(虚线)和颜色labels = [''] # 去除箱线图的x轴刻度值)def Boxplot_anomalies(dataset1,ColumnName,whis):
    print_outliers(dataset1, ColumnName,whis)
    draw_Boxplot(dataset1, ColumnName,whis)
whis=4plt.figure(figsize=(30, 20))
plt.subplot(231)
ColumnName='Light'Boxplot_anomalies(dataset1,ColumnName,whis)
plt.subplot(232)
ColumnName='Temperature'Boxplot_anomalies(dataset1,ColumnName,whis)
plt.subplot(233)
ColumnName='Humidity'Boxplot_anomalies(dataset1,ColumnName,whis)
plt.subplot(234)
ColumnName='CO2'Boxplot_anomalies(dataset1,ColumnName,whis)
plt.subplot(235)
ColumnName='HumidityRatio'Boxplot_anomalies(dataset1,ColumnName,whis)
plt.savefig('../input/Boxplot_anomalies4.png')

正态分布法

根据正态分布的定义可知，数据点落在偏离均值正负1倍标准差(即sigma值)内的概率为68.2%；数据点落在偏离均值正负2倍标准差内的概率为95.4%；数据点落在偏离均值正负3倍标准差内的概率为99.6%。

所以，换个角度思考上文提到的概率值，如果数据点落在偏离均值正负2倍标准差之外的概率就不足5%，它属于小概率事件，即认为这样的数据点为异常点。同理，如果数据点落在偏离均值正负3倍标准差之外的概率将会更小，可以认为这些数据点为极端异常点。为使读者直观地理解文中提到的概率值，可以查看标准正态分布的概率密度图，如下图所示：

进一步，基于上图的结论，可以按照下表中的判断条件，识别出数值型变量的异常点和极端异常点，如下表所示：

利用正态分布的知识点，结合pyplot子模块中的plot函数绘制折线图和散点图，并借助于两条水平参考线识别异常值或极端异常值。

上图为温度曲线及其上下2西格玛控制线和3西格玛控制线，异常点159个，极端异常点0个。

上图为CO2曲线及其上下2西格玛控制线和3西格玛控制线，异常点897个，极端异常点64个。

正态分布法的优点是计算快，但不能显示异常点的值和ID，而且只能计算一维数据。

从原理分析，不论过程变量还是离散的测量点，按照正态分布法都能找到异常检测点。即便这些点仅仅是正常的极值，也很可能被识别为异常点(参照上述图中的温度曲线异常检测结果)。

从适用场景情况来看，

对于过程变量来说，仪表故障、传输错误、瞬间极值这些情况很容易被检测出来，但要对结果进行甄别，防止正常的瞬间极值被识别为异常点。

对于集散的测量点，比如零件尺寸数据、产品质量数据，使用此方法可以非常容易识别极值和异常，进而进行质量改善。事实上，基于数理统计的质量管理方法，六西格玛和SPC都含有类似的工具。

其它的物联网测量参数的数据不再赘述。请使用如下代码自行体验：

import pandas as pdimport matplotlib.pyplot as plt# 读入外部数据dataset2 =pd.read_csv('../input/dataset2.csv',index_col='date',parse_dates=['date'])# 绘制单条折线图，并在折线图的基础上添加点图plt.figure(figsize=(30, 20))
plt.plot(dataset2.index, # x轴数据dataset2["CO2"], # y轴数据linestyle = '-', # 设置折线类型linewidth = 2, # 设置线条宽度color = 'steelblue', # 设置折线颜色marker = 'o', # 往折线图中添加圆点markersize = 4, # 设置点的大小markeredgecolor='black', # 设置点的边框色markerfacecolor='black') # 设置点的填充色# 添加上下界的水平参考线(便于判断异常点，如下判断极端异常点，只需将2改为3)plt.axhline(y = dataset2["CO2"].mean() - 2* dataset2["CO2"].std(), linestyle = '--', color = 'gray')
plt.axhline(y = dataset2["CO2"].mean() + 2* dataset2["CO2"].std(), linestyle = '--', color = 'gray')
plt.axhline(y = dataset2["CO2"].mean() - 3* dataset2["CO2"].std(), linestyle = '--', color = 'red')
plt.axhline(y = dataset2["CO2"].mean() + 3* dataset2["CO2"].std(), linestyle = '--', color = 'red')# 显示图形plt.savefig('../input/hist_anomalies.png')# 导入模块，用于日期刻度的修改(因为默认格式下的日期刻度标签并不是很友好)import matplotlib as mpl# 获取图的坐标信息ax = plt.gca()# 设置日期的显示格式date_format =mpl.dates.DateFormatter("%m-%d")
ax.xaxis.set_major_formatter(date_format)# 设置x轴每个刻度的间隔天数xlocator = mpl.ticker.MultipleLocator(7)
ax.xaxis.set_major_locator(xlocator)# 为了避免x轴刻度标签的紧凑，将刻度标签旋转45度plt.xticks(rotation=45)# 计算判断异常点和极端异常点的临界值outlier_ll =dataset2["CO2"].mean() - 2* dataset2["CO2"].std()
outlier_ul = dataset2["CO2"].mean() + 2* dataset2["CO2"].std()
extreme_outlier_ll = dataset2["CO2"].mean() - 3* dataset2["CO2"].std()
extreme_outlier_ul = dataset2["CO2"].mean() + 3* dataset2["CO2"].std()# 寻找异常点cnts1 =dataset2[(dataset2["CO2"] > outlier_ul) | (dataset2["CO2"] < extreme_outlier_ll)]print(cnts1["CO2"])# 寻找极端异常点cnts2 =dataset2[(dataset2["CO2"] > extreme_outlier_ul) | (dataset2["CO2"] < extreme_outlier_ll)]print(cnts2["CO2"])

DBSCAN算法

DBSCAN(Density-BasedSpatial Clustering of Applications with Noise，具有噪声的基于密度的聚类方法)是一种很典型的密度聚类算法。

关于DBSCAN的一些简单介绍以及用法，参见文章什么？ MindSphere可以这样玩？(12)—— 异常检测

使用物联网数据进行异常检测的结果如下图：

左上图为未经异常检测的原始数据的散点图。

中上图为设置半径200，最小分组为6的结果。

右上图为设置半径150，最小分组为6的结果。

左下图为设置半径100，最小分组为6的结果。

中下图为设置半径50，最小分组为6的结果。

右下图为设置半径10，最小分组为3的结果。

由六个图结合可知，半径50，最小分组为6可以很好地识别异常点。同时，从图上方的数据可以得出DBSCAN的计算效率非常高，2600个数据点小于0.1秒。

此图源代码如下；

import matplotlib.pyplotas pltimport  pandas as pdimport timefrom sklearn.cluster import KMeansfrom sklearn.cluster import DBSCAN
plt.figure(figsize=(30, 20))
dataset1 = pd.read_csv('../input/dataset1.csv',index_col='date',parse_dates=['date'])
X=dataset1.loc[:,['CO2','Light']]
plt.subplot(231)
plt.scatter(X['CO2'], X['Light'])# DBSCANt0 = time.time()
dbscan = DBSCAN(eps=200, min_samples=6).fit(X)
t = time.time() - t0
plt.subplot(232)
plt.scatter(X['CO2'], X['Light'], c=dbscan.labels_)
plt.title('time : %f' % t)# DBSCANt0 = time.time()
dbscan = DBSCAN(eps=150, min_samples=6).fit(X)
t = time.time() - t0
plt.subplot(233)
plt.scatter(X['CO2'], X['Light'], c=dbscan.labels_)
plt.title('time : %f' % t)# DBSCANt0 = time.time()
dbscan = DBSCAN(eps=100, min_samples=6).fit(X)
t = time.time() - t0
plt.subplot(234)
plt.scatter(X['CO2'], X['Light'], c=dbscan.labels_)
plt.title('time : %f' % t)# DBSCANt0 = time.time()
dbscan = DBSCAN(eps=50, min_samples=6).fit(X)
t = time.time() - t0
plt.subplot(235)
plt.scatter(X['CO2'], X['Light'], c=dbscan.labels_)
plt.title('time : %f' % t)# DBSCANt0 = time.time()
dbscan = DBSCAN(eps=10, min_samples=3).fit(X)
t = time.time() - t0
plt.subplot(236)
plt.scatter(X['CO2'], X['Light'], c=dbscan.labels_)
plt.title('time : %f' % t)
plt.savefig('../input/kmeans_DBSCAN.png')

DBSCAN是MartinEster, Hans-PeterKriegel等人于1996年提出的一种基于密度的空间的数据聚类方法，该算法是最常用的一种聚类方法。该算法将具有足够密度区域作为距离中心，不断生长该区域,算法基于一个事实：一个聚类可以由其中的任何核心对象唯一确定。该算法利用基于密度的聚类的概念，即要求聚类空间中的一定区域内所包含对象(点或其他空间对象)的数目不小于某一给定阈值。该方法能在具有噪声的空间数据库中发现任意形状的簇，可将密度足够大的相邻区域连接，能有效处理异常数据，主要用于对空间数据的聚类，优缺点总结如下:

该算法的优点：

(1)聚类速度快且能够有效处理噪声点和发现任意形状的空间聚类；

(2)与K-MEANS比较起来，不需要输入要划分的聚类个数；

(3)聚类簇的形状没有偏倚；

(4)可以在需要时输入过滤噪声的参数。

该算法的缺点：

(1)当数据量增大时，要求较大的内存支持I/O消耗也很大；

(2)当空间聚类的密度不均匀、聚类间距差相差很大时，聚类质量较差，因为这种情况下参数MinPts和Eps选取困难。

(3)算法聚类效果依赖与距离公式选取，实际应用中常用欧式距离，对于高维数据，存在“维数灾难”。

该算法基本概念：

(1)Eps邻域：给定对象半径Eps内的邻域称为该对象的Eps邻域;

(2)核心点(core point)：如果对象的Eps邻域至少包含最小数目MinPts的对象，则称该对象为核心对象;

(3)边界点(edge point)：边界点不是核心点，但落在某个核心点的邻域内;

(4)噪音点(outlier point)：既不是核心点，也不是边界点的任何点;

(5)直接密度可达(directly density-reachable)：给定一个对象集合内，与核心点的距离在Eps邻域内，则是直接密度可达的;

(6)密度可达(density-reachable)：如果在一个对象集合内，有一组点是依次密度直接可达，则该组内的首尾两点关于Eps和MinPts密度可达的;

(7)密度相连(density-connected)：给定一个对象集合内，如果有两个点是关于Eps和MinPts密度可达，说明这两个点是关于Eps和MinPts密度相连的。

“直接密度可达”和“密度可达”概念描述：根据前文基本概念的描述，由于有标记的各点­M、P、O和R的Eps近邻均包含3个以上的点，因此它们都是核对象；M­是从P“直接密度可达”；而Q则是从­M“直接密度可达”；基于上述结果，Q是从P“密度可达”；但P从Q无法“密度可达”(非对称)。类似地，S和R从O是“密度可达”的；O、R和S均是“密度相连”(对称)的。

参数设置：

DBSCAN共包括3个输入数据：数据集D，给定点在邻域内成为核心对象的最小邻域点数：MinPts,邻域半径：Eps，其中Eps和MinPts需要根据具体应用人为设定。

(1) Eps的值可以使用绘制k-距离曲线(k-distance graph)方法得当，在k-距离曲线图明显拐点位置为对应较好的参数。若参数设置过小，大部分数据不能聚类；若参数设置过大，多个簇和大部分对象会归并到同一个簇中。

K-距离：K距离的定义在DBSCAN算法原文中给出了详细解说，给定K邻域参数k,对于数据中的每个点，计算对应的第k个最近邻域距离，并将数据集所有点对应的最近邻域距离按照降序方式排序，称这幅图为排序的k距离图，选择该图中第一个谷值点位置对应的k距离值设定为Eps。一般将k值设为4。

(2) MinPts的选取有一个指导性的原则(a rule of thumb)，MinPts≥dim+1,其中dim表示待聚类数据的维度。MinPts设置为1是不合理的，因为设置为1，则每个独立点都是一个簇，MinPts≤2时，与层次距离最近邻域结果相同，因此，MinPts必须选择大于等于3的值。若该值选取过小，则稀疏簇中结果由于密度小于MinPts，从而被认为是边界点儿不被用于在类的进一步扩展；若该值过大，则密度较大的两个邻近簇可能被合并为同一簇。因此，该值是否设置适当会对聚类结果造成较大影响。

该算法时间复杂度：

(1)DBSCAN的基本时间复杂度是O(N*找出Eps领域中的点所需要的时间), N是点的个数。最坏情况下时间复杂度是O(N2)

(2)在低维空间数据1中,有一些数据结构如KD树，使得可以有效的检索特定点给定距离内的所有点，时间复杂度可以降低到O(NlogN)

该算法空间复杂度：

低维和高维数据中，其空间都是O(N)，对于每个点它只需要维持少量数据，即簇标号和每个点的标识(核心点或边界点或噪音点)

原文链接：https://blog.csdn.net/zhouxianen1987/article/details/68945844

一句话总结：使用该算法的关键在于参数的调整，即如何设置合适的参数实现最佳的检测效果。

局部离群因子LOF检测法

LocalOutlier Factor(LOF)是基于密度的经典算法(Breuning et.al. 2000), 文章发表于 SIGMOD 2000, 到目前已经有 3000+ 的引用。在 LOF 之前的异常检测算法大多是基于统计方法的，或者是借用了一些聚类算法用于异常点的识别(比如，DBSCAN，OPTICS)。但是，基于统计的异常检测算法通常需要假设数据服从特定的概率分布，这个假设往往是不成立的。而聚类的方法通常只能给出 0/1 的判断(即：是不是异常点)，不能量化每个数据点的异常程度。相比较而言，基于密度的LOF算法要更简单、直观。它不需要对数据的分布做太多要求，还能量化每个数据点的异常程度(outlierness)。

LOF算法详细介绍参见https://zhuanlan.zhihu.com/p/37753692

分别调整距离K为8，10，20得出的异常检测结果如下三个图，其中红色圆圈中的点为异常点，圆圈越大异常程度越大。

如上三个图形可知，是否能彻底检测出所有的异常点，需要认真调整参数。

算法实现的代码如下：

import numpy as npimport pandas as pdimport matplotlib.pyplot as pltdef localoutlierfactor(data, predict,k):from sklearn.neighbors import LocalOutlierFactor
    clf = LocalOutlierFactor(n_neighbors=k + 1, algorithm='auto', contamination=0.1, n_jobs=-1)
    clf.fit(data)# 记录 k 邻域距离predict['k distances'] = clf.kneighbors(predict)[0].max(axis=1)# 记录 LOF 离群因子，做相反数处理predict['local outlier factor'] =-clf._decision_function(predict.iloc[:, :-1])return predictdef plot_lof(result, method):import matplotlib.pyplot as plt
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号plt.figure(figsize=(8, 4)).add_subplot(111)
    plt.scatter(result[result['local outlier factor'] > method].index,
                result[result['local outlier factor'] > method]['local outlier factor'], c='red', s=50,marker='.', alpha=None,label='离群点')
    plt.scatter(result[result['local outlier factor'] method].sort_values(by='local outlier factor')
    inliers = predict[predict['local outlier factor']

标签组：[聚类] [层次聚类方法] [dbscan] 

上一篇：干完这7道题，包你学会从对数幅频特性曲线图中求开环增益K！

下一篇：means聚类到底应该聚为几类？