作者：naughty
本文讨论了使用Python对单变量数据的异常点分析的几种方法。原文链接：https://my.oschina.net/taogang/blog/279402
大数据时代，数据的异常分析被广泛的用于各个场合。今天我们就来看一看其中的一种场景，对于单变量数据集的异常检测。
所谓单变量，就是指数据集中只有一个变化的值，下面我们来看看今天我们要分析的的数据，点击这里数据文件下载数据文件。
分析数据的第一步是要加载文件， 本文使用了 numpy ，pandas ，scikit learn 等常见的数据分析要用到的 Python 库。
importnumpy asnp
importpandas aspddf = pd.read_csv( "farequote.csv")
Pandas是一个常用的数据分析的Python库，提供对数据的加载，清洗，抽取，变形等操作。Pandas依赖numpy，numpy提供了基于列/多维数组（List/N-D
Array ）的数据结构的操作。许多科学计算和数据分析的库都依赖于 numpy 。
df 是 Pandas 中常用的数据类型 dataframe，dataframe 类似一个数据库的表，使用 df.head()可以得到数据的头几行，以便了解数据的概貌。
该数据结构中，第一列式 Pandas 添加的索引，第一行是每一列数据的名字，除了第一列，每一列数据可以看成是一个变量，所以该数据集共有三个变量，时间（time ）、航空公司名称 ( airline)、响应时间（responsetime）。我们可以这样理解，该数据集记录了一段时间内，各个航空公司飞机延误的时间。我们希望通过分析找出是否存在异常的情况。
注意，我们是要分析单变量，所以所有的分析都是基于某一个航空公司的数据，所以就需要对该数据集做一个查询，找出要分析的航空公司。首先要知道有哪些航空公司，使用 np.unique(df.airline)可以找到所有的航空公司代码，类似 SQL 的 Unique 命令。
array([ 'AAL', 'ACA', 'AMX', 'ASA', 'AWE', 'BAW', 'DAL', 'EGF', 'FFT',
'JAL', 'JBU', 'JZA', 'KLM', 'NKS', 'SWA', 'SWR', 'TRS', 'UAL', 'VRD'], dtype= '|S3')
查询某个航空公司的数据使用 dataframe 的 query 方法，类似 SQL 的 select。Query 返回的结果仍然是一个 dataframe 对象。
dd = df.query( 'airline=="KLM"') # 得到法航的数据
我们先了解一下数据的大致信息，使用 describe 方法
dd.responsetime.describe()
得到如下的结果：
count 1724.000000
mean 1500.613766
std 100.085320
min 1209.766800
25% 1434.084625
50% 499.135000
75% 1567.831025
max 1818.774100
Name: responsetime, dtype: float64
该结果返回了数据集 responsetime 维度上的主要统计指标，个数，均值，方差，最大最小值等等，也可以调用单独的方法例如 min(）， mean(）等来获得某一个指标。
基于标准差得异常检测
下面我们就可以开始异常点的分析了，对于单变量的异常点分析，最容易想到的就是基于标准差（Standard Deviation）的方法了。我们假定数据的正态分布的，利用概率密度函数，我们知道
95.449974面积在平均数左右两个标准差的范围内
99.730020%的面积在平均数左右三个标准差的范围内
99.993666的面积在平均数左右三个标准差的范围内
所以我们95%也就是大概两个标准差为门限，凡是落在门限外的都认为是异常点。代码如下
defa1(dataframe, threshold=.95):d = dataframe[ 'responsetime'] dataframe[ 'isAnomaly'] = d d.quantile(threshold) returndataframe printa1(dd)
运行以上程序我们得到如下结果
_time airline responsetime isAnomaly
202013- 02- 01T23: 57: 59.000- 0700KLM 1481.4945False
762013- 02- 01T23: 52: 34.000- 0700KLM 1400.9050False
1242013- 02- 01T23: 47: 10.000- 0700KLM 1501.4313False
2032013- 02- 01T23: 39: 08.000- 0700KLM 1278.9509False
2812013- 02- 01T23: 32: 27.000- 0700KLM 1386.4157False
3362013- 02- 01T23: 26: 09.000- 0700KLM 1629.9589False
3642013- 02- 01T23: 23: 52.000- 0700KLM 1482.5900False
4482013- 02- 01T23: 16: 08.000- 0700KLM 1553.4988False
5112013- 02- 01T23: 10: 39.000- 0700KLM 1555.1894False
5162013- 02- 01T23: 10: 08.000- 0700KLM 1720.7862True
5532013- 02- 01T23: 06: 29.000- 0700KLM 1306.6489False
5932013- 02- 01T23: 03: 03.000- 0700KLM 1481.7081False
6092013- 02- 01T23: 01: 29.000- 0700KLM 1521.0253False
6662013- 02- 01T22: 56: 04.000- 0700KLM 1675.2222True
... ... ... ...
结果数据集上多了一列 isAnomaly 用来标记每一行记录是否是异常点，我们看到已经有一些点被标记为异常点了。
我们看看程序的详细内容：
方法 a 1定义了一个异常检测的函数
dataframe[‘responsetime’] 等价于 dataframe.responsetime ,该操作取出 responsetime 这一列的值
d.quantile(threshold)用正态分布假定返回位于95%的点的值，大于该值得点都落在正态分布95%之外
d d.quantile(threshold)是一个数组操作，返回的新数组是 responsetime 和 threshold 的比较结果， [False , False , True ,… … False]
然后通过 dataframe 的赋值操作增加一个新的列，标记所有的异常点。
数据可视化往往是数据分析的最后一步，我们看看结果如何：
importmatplotlib.pyplot aspltda = a1(dd)fig = plt.figure()ax1 = fig.add_subplot( 2, 1, 1)ax2 = fig.add_subplot( 2, 1, 2)ax1.plot(da[ 'responsetime'])ax2.plot(da[ 'isAnomaly'])
这异常点也太多了，用99%在试试：
现在似乎好一点，然而我们知道，对于数据集的正态分布的假定往往是不成立的，假如数据分布在大小两头，那么这样的异常检测就很难奏效了。我们看看其他一些改进的方法。
基于 ZSCORE 的异常检测
zscore 的计算如下
sd 是标准差， X 是均值。一般建议门限值取为3.5
代码如下：
defa2(dataframe, threshold=3.5):d = dataframe[ 'responsetime'] zscore = (d - d.mean()) / d.std() dataframe[ 'isAnomaly'] = zscore.abs() threshold
returndataframe
另外还有一种增强的 zscore 算法，基于 MAD 。 MAD 的定义是
其中 X 是中位数。
增强的 zscore 算法如下：
defa3(dataframe, threshold=3.5):dd = dataframe[ 'responsetime'] MAD = (dd - dd.median()).abs().median() zscore = ((dd - dd.median()) * 0.6475/ MAD).abs() dataframe[ 'isAnomaly'] = zscore threshold
returndataframe
用 zscore 算法得到：
调整门限为3得到
如果换一组数据 AAL ，结果会怎么样呢？
我们发现有一段时间，所有的响应都很慢，我们想要把这些点都标记为异常，可能么？
基于 KMEAN 聚集的异常检测
通常基于 KMEAN 的聚集算法并不适用于异常点检测，以为聚集算法总是试图平衡每一个聚集中的点的数目，所以对于少数的异常点，聚集非常不好用，但是我们这个例子中，异常点都聚在一起，所以应该可以使用。
首先，为了看清聚集，我们使用时间序列的常用分析方法，增加一个维度，该维度是每一个点得前一个点得响应时间。
preresponse = 0
newcol = []newcol.append( 0)
forindex, row indd.iterrows():
ifpreresponse != 0: newcol.append(preresponse) preresponse = row.responsetimedd[ "t0"] = newcolplt.scatter(dd.t0, dd.responsetime)
我们利用 iterrows 来循环数据，把前一个点的响应时间增加到当前点，第一个点的该值为0，命名该列为 t0。然后用 scatter plot 把它画出来。
上面是法航 KLM 的数据，其中最左边的点是一个无效的点，因为前一个点的响应时间不知道所以填了0，分析时应该过滤该店。
对于 AAL ，我们可以清楚的看到两个聚集：
其中右上方的聚集，也就是点数目比较少得聚集就是我们希望检测到的异常点得集合。
我们看看如何使用 KMEAN 算法来检测吧：
defa4(dataframe, threshold = .9):# add one dimention of previous responsepreresponse = 0newcol = [] newcol.append( 0)
forindex, row indataframe.iterrows():
ifpreresponse != 0: newcol.append(preresponse) preresponse = row.responsetime dataframe[ "t0"] = newcol # remove first row as there is no previous event for timedd = dataframe.drop(dataframe.head( 1).index) clf = cluster.KMeans(n_clusters= 2) X = np.array(dd[[ 'responsetime', 't0']]) cls = clf.fit_predict(X) freq = itemfreq(cls) (A, B) = (freq[ 0, 1], freq[ 1, 1]) t = abs(A - B) / max(A, B)
ift threshold : # "Anomaly Detected!"index = freq[ 0, 0]
ifA B : index = freq[ 1, 0] dd[ 'isAnomaly'] = (cls == index)
else: # "No Anomaly Point"dd[ 'isAnomaly'] = Falsereturndd
其核心代码是以下这几行：
clf = cluster.KMeans(n_clusters= 2)X = np.array(dd[[ 'responsetime', 't0']])cls = clf.fit_predict(X)
cluster.KMeans 返回一个预测模型，我们假定有两个聚集。你可以试着加大聚集的数量，结果没什么影响。
dd[[ ‘ responsetime ‘,’ t 0’ ]] 返回一个2 * n 的数组，并赋值给 X ，用于聚集计算。
fit_pridict 方法是对 X 做聚集运算，并计算每一个点对应的聚集编号。
freq = itemfreq(cls)
itemfreq 返回聚集结果中每一个聚集的发生频率，如果其中一个比另一个显著地多，我们则认为那个少得是异常点聚集。
用该方法可以把所有聚集里的点标记为异常点。
这里我用红色标记结果让大家看的清楚一点，注意因为是 line chart ，连个竖线间的都是异常点。
总结
除了上述的算法，还有其它一些相关的算法，大家如果对背后的数据知识有兴趣的话，可以参考这篇相关介绍。
单变量的异常检测算法相对比较简单，但是要做到精准检测就更难，因为掌握的信息更少。另外boxplot也经常被用于异常检测，他和基于方差的异常检测是一致的，只不过用图形让大家一目了然的获得结果，大家有兴趣可以了解一下。
点击阅读原文，查看更多 Python 教程和资源
