用Python做单变量数据集的异常点分析
2014年6月12日星期四
大数据时代，数据的异常分析被广泛的用于各个场合。 今天我们就来看一看其中的一种场景，对于单变量数据集的异常检测。
所谓单变量，就是指数据集中只有一个变化的值，下面我们来看看今天我们要分析的的数据，点击这里数据文件下载数据文件。
分析数据的第一步是要加载文件， 本文使用了numpy，pandas，scikit learn等常见的数据分析要用到的Python库。
import numpy as np
import pandas as pd
df = 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% 1499.135000 75% 1567.831025 max 1818.774100 Name: responsetime, dtype: float64
该结果返回了数据集responsetime维度上的主要统计指标，个数，均值，方差，最大最小值等等，也可以调用单独的方法例如min（），mean（）等来获得某一个指标。
基于标准差得异常检测
下面我们就可以开始异常点的分析了，对于单变量的异常点分析，最容易想到的就是基于标准差（Standard Deviation）的方法了。我们假定数据的正态分布的，利用概率密度函数，我们知道
- 95.449974面积在平均数左右两个标准差的范围内
- 99.730020%的面积在平均数左右三个标准差的范围内
- 99.993666的面积在平均数左右三个标准差的范围内
所以我们95%也就是大概两个标准差为门限，凡是落在门限外的都认为是异常点。代码如下
def a1(dataframe, threshold=.95):
d = dataframe['responsetime']
dataframe['isAnomaly'] = d > d.quantile(threshold)
return dataframe
print a1(dd)
运行以上程序我们得到如下结果
_time airline responsetime isAnomaly 20 2013-02-01T23:57:59.000-0700 KLM 1481.4945 False 76 2013-02-01T23:52:34.000-0700 KLM 1400.9050 False 124 2013-02-01T23:47:10.000-0700 KLM 1501.4313 False 203 2013-02-01T23:39:08.000-0700 KLM 1278.9509 False 281 2013-02-01T23:32:27.000-0700 KLM 1386.4157 False 336 2013-02-01T23:26:09.000-0700 KLM 1629.9589 False 364 2013-02-01T23:23:52.000-0700 KLM 1482.5900 False 448 2013-02-01T23:16:08.000-0700 KLM 1553.4988 False 511 2013-02-01T23:10:39.000-0700 KLM 1555.1894 False 516 2013-02-01T23:10:08.000-0700 KLM 1720.7862 True 553 2013-02-01T23:06:29.000-0700 KLM 1306.6489 False 593 2013-02-01T23:03:03.000-0700 KLM 1481.7081 False 609 2013-02-01T23:01:29.000-0700 KLM 1521.0253 False 666 2013-02-01T22:56:04.000-0700 KLM 1675.2222 True ... ... ... ...
结果数据集上多了一列isAnomaly用来标记每一行记录是否是异常点，我们看到已经有一些点被标记为异常点了。
我们看看程序的详细内容：
- 方法a1定义了一个异常检测的函数
- dataframe['responsetime']等价于dataframe.responsetime,该操作取出responsetime这一列的值
- d.quantile(threshold)用正态分布假定返回位于95%的点的值，大于该值得点都落在正态分布95%之外
- d > d.quantile(threshold)是一个数组操作，返回的新数组是responsetime和threshold的比较结果，[False,False,True,... ... False]
- 然后通过dataframe的赋值操作增加一个新的列，标记所有的异常点。
数据可视化往往是数据分析的最后一步，我们看看结果如何：
import matplotlib.pyplot as plt da = 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
代码如下：
def a2(dataframe, threshold=3.5):
d = dataframe['responsetime']
zscore = (d - d.mean())/d.std()
dataframe['isAnomaly'] = zscore.abs() > threshold
return dataframe
另外还有一种增强的zscore算法，基于MAD。MAD的定义是
其中X是中位数。
增强的zscore算法如下：
def a3(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
return dataframe
用zscore算法得到：
调整门限为3得到
如果换一组数据AAL，结果会怎么样呢?
我们发现有一段时间，所有的响应都很慢，我们想要把这些点都标记为异常，可能么？
基于KMEAN聚集的异常检测
通常基于KMEAN的聚集算法并不适用于异常点检测，以为聚集算法总是试图平衡每一个聚集中的点的数目，所以对于少数的异常点，聚集非常不好用，但是我们这个例子中，异常点都聚在一起，所以应该可以使用。
首先，为了看清聚集，我们使用时间序列的常用分析方法，增加一个维度，该维度是每一个点得前一个点得响应时间。
preresponse = 0
newcol = []
newcol.append(0)
for index, row in dd.iterrows():
if preresponse != 0:
newcol.append(preresponse)
preresponse = row.responsetime
dd["t0"] = newcol
plt.scatter(dd.t0,dd.responsetime)
我们利用iterrows来循环数据，把前一个点的响应时间增加到当前点，第一个点的该值为0，命名该列为t0。然后用scatter plot把它画出来。
上面是法航KLM的数据，其中最左边的点是一个无效的点，因为前一个点的响应时间不知道所以填了0，分析时应该过滤该店。
对于AAL，我们可以清楚的看到两个聚集：
其中右上方的聚集，也就是点数目比较少得聚集就是我们希望检测到的异常点得集合。
我们看看如何使用KMEAN算法来检测吧：
def a4(dataframe, threshold = .9):
## add one dimention of previous response
preresponse = 0
newcol = []
newcol.append(0)
for index, row in dataframe.iterrows():
if preresponse != 0:
newcol.append(preresponse)
preresponse = row.responsetime
dataframe["t0"] = newcol
## remove first row as there is no previous event for time
dd = 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)
if t > threshold :
## "Anomaly Detected!"
index = freq[0,0]
if A > B :
index = freq[1,0]
dd['isAnomaly'] = (cls == index)
else :
## "No Anomaly Point"
dd['isAnomaly'] = False
return dd
其核心代码是以下这几行：
clf = cluster.KMeans(n_clusters=2) X=np.array(dd[['responsetime','t0']]) cls = clf.fit_predict(X)
cluster.KMeans返回一个预测模型，我们假定有两个聚集。你可以试着加大聚集的数量，结果没什么影响。
dd[['responsetime','t0']]返回一个2*n的数组，并赋值给X，用于聚集计算。
fit_pridict方法是对X做聚集运算,并计算每一个点对应的聚集编号。
freq = itemfreq(cls)
itemfreq返回聚集结果中每一个聚集的发生频率，如果其中一个比另一个显著地多，我们则认为那个少得是异常点聚集。
用该方法可以把所有聚集里的点标记为异常点。
这里我用红色标记结果让大家看的清楚一点，注意因为是line chart，连个竖线间的都是异常点。
总结
除了上述的算法，还有其它一些相关的算法，大家如果对背后的数据知识有兴趣的话，可以参考这篇相关介绍。
单变量的异常检测算法相对比较简单，但是要做到精准检测就更难，因为掌握的信息更少。另外boxplot也经常被用于异常检测，他和基于方差的异常检测是一致的，只不过用图形让大家一目了然的获得结果，大家有兴趣可以了解一下。
标签： Data, Data Analysis, Python
使用Python抓取欧洲足球联赛数据
背景
Web Scraping
在大数据时代，一切都要用数据来说话，大数据处理的过程一般需要经过以下的几个步骤
- 数据的采集和获取
- 数据的清洗，抽取，变形和装载
- 数据的分析，探索和预测
- 数据的展现
其中首先要做的就是获取数据，并提炼出有效地数据，为下一步的分析做好准备。
数据的来源多种多样，以为我本身是足球爱好者，而世界杯就要来了，所以我就想提取欧洲联赛的数据来做一个分析。许多的网站都提供了详细的足球数据，例如：
这些网站都提供了详细的足球数据，然而为了进一步的分析，我们希望数据以格式化的形式存储，那么如何把这些网站提供的网页数据转换成格式化的数据呢？这就要用到Web scraping的技术了。简单地说，Web Scraping就是从网站抽取信息， 通常利用程序来模拟人浏览网页的过程，发送http请求，从http响应中获得结果。
Web Scraping 注意事项
在抓取数据之前，要注意以下几点：
- 阅读网站有关数据的条款和约束条件，搞清楚数据的拥有权和使用限制
- 友好而礼貌，使用计算机发送请求的速度飞人类阅读可比，不要发送非常密集的大量请求以免造成服务器压力过大
- 因为网站经常会调整网页的结构，所以你之前写的Scraping代码，并不总是能够工作，可能需要经常调整
- 因为从网站抓取的数据可能存在不一致的情况，所以很有可能需要手工调整
Python Web Scraping 相关的库
Python提供了很便利的Web Scraping基础，有很多支持的库。这里列出一小部分
- BeautifulSoup http://www.crummy.com/software/BeautifulSoup/
- Scrapy http://scrapy.org/
- webscraping https://code.google.com/p/webscraping/
当然也不一定要用Python或者不一定要自己写代码，推荐关注import.io
Web Scraping 代码
下面，我们就一步步地用Python，从腾讯体育来抓取欧洲联赛13/14赛季的数据。
首先要安装Beautifulsoup
pip install beautifulsoup4
我们先从球员的数据开始抓取。
球员数据的Web请求是http://soccerdata.sports.qq.com/playerSearch.aspx?lega=epl&pn=2 ，返回的内容如下图所示：
该web服务有两个参数，lega表示是哪一个联赛，pn表示的是分页的页数。
首先我们先做一些初始化的准备工作
from urllib2 import urlopen import urlparse import bs4 BASE_URL = "http://soccerdata.sports.qq.com" PLAYER_LIST_QUERY = "/playerSearch.aspx?lega=%s&pn=%d" league = ['epl','seri','bund','liga','fran','scot','holl','belg'] page_number_limit = 100 player_fields = ['league_cn','img','name_cn','name','team','age','position_cn','nation','birth','query','id','teamid','league']
urlopen,urlparse,bs4是我们将要使用的Python库。
BASE_URL,PLAYER_LIST_QUERY,league,page_number_limit和player_fields是我们会用到的一些常量。
下面是抓取球员数据的具体代码：
def get_players(baseurl):
html = urlopen(baseurl).read()
soup = bs4.BeautifulSoup(html, "lxml")
players = [ dd for dd in soup.select('.searchResult tr') if dd.contents[1].name != 'th']
result = []
for player in players:
record = []
link = ''
query = []
for item in player.contents:
if type(item) is bs4.element.Tag:
if not item.string and item.img:
record.append(item.img['src'])
else :
record.append(item.string and item.string.strip() or 'na')
try:
o = urlparse.urlparse(item.a['href']).query
if len(link) == 0:
link = o
query = dict([(k,v[0]) for k,v in urlparse.parse_qs(o).items()])
except:
pass
if len(record) != 10:
for i in range(0, 10 - len(record)):
record.append('na')
record.append(unicode(link,'utf-8'))
record.append(unicode(query["id"],'utf-8'))
record.append(unicode(query["teamid"],'utf-8'))
record.append(unicode(query["lega"],'utf-8'))
result.append(record)
return result
result = []
for url in [ BASE_URL + PLAYER_LIST_QUERY % (l,n) for l in league for n in range(page_number_limit) ]:
result = result + get_players(url)
我们来看看抓取球员数据的详细过程：
首先我们定义了一个get_players方法，该方法会返回某一请求页面上所有球员的数据。为了得到所有的数据，我们通过一个for循环，因为要循环各个联赛，每个联赛又有多个分页，一般情况下是需要一个双重循环的：
for i in league:
for j in range(0, 100):
url = BASE_URL + PLAYER_LIST_QUERY % (l,n)
## send request to url and do scraping
Python的list comprehension可以很方便的通过构造一个列表的方式来减少循环的层次。
另外Python还有一个很方便的语法来合并连个列表： list = list1 + list2
好我们再看看如何使用BeautifulSoup来抓取网页中我们需要的内容。
首先调用urlopen读取对应url的内容，通常是一个html，用该html构造一个beautifulsoup对象。
beautifulsoup对象支持很多查找功能，也支持类似css的selector。通常如果有一个DOM对象是<xx class='cc'>,我们使用以下方式来查找：
obj = soup.find("xx","cc")
另外一种常见的方式就是通过CSS的selector方式，在上述代码中，我们选择class=searchResult元素里面，所有的tr元素，过滤掉th也就是表头元素。
for dd in soup.select('.searchResult tr') if dd.contents[1].name != 'th'
对于每一行记录tr，生成一条球员记录，并存放在一个列表中。所以我们就循环tr的内容tr.contents,获得对应的field内容。
对于每一个tr的content，我们先检查其类型是不是一个Tag,对于Tag类型有几种情况，一种是包含img的情况，我们需要取出球员的头像图片的网址。
另一种是包含了一个链接，指向其他数据内容
所以在代码中要分别处理这些不同的情况。
对于一个Tag对象，Tag.x可以获得他的子对象，Tag['x']可以获得Tag的attribute的值。
所以用item.img['src']可以获得item的子元素img的src属性。
对已包含链接的情况，我们通过urlparse来获取查询url中的参数。这里我们利用了dict comprehension的把查询参数放入一个dict中，然后添加到列表中。
dict([(k,v[0]) for k,v in urlparse.parse_qs(o).items()])
对于其它情况，我们使用Python 的and or表达式以确保当Tag的内容为空时，我们写入‘na’，该表达式类似C/C++或Java中的三元操作符 X ? A : B
然后有一段代码判断当前记录的长度是否大于10，不大于10则用空值填充，目的是避免一些不一致的地方。
if len(record) != 10:
for i in range(0, 10 - len(record)):
record.append('na')
最后，我们把query中的一些相关的参数如球员的id，球队的id，所在的联赛代码等加入到列表。
record.append(unicode(link,'utf-8')) record.append(unicode(query["id"],'utf-8')) record.append(unicode(query["teamid"],'utf-8')) record.append(unicode(query["lega"],'utf-8'))
最后我们把本页面所有球员的列表放入一个列表返回。
好了，现在我们拥有了一个包含所有球员的信息的列表，我们需要把它存下来，以进一步的处理，分析。通常，csv格式是一个常见的选择。
import csv
def write_csv(filename, content, header = None):
file = open(filename, "wb")
file.write('\xEF\xBB\xBF')
writer = csv.writer(file, delimiter=',')
if header:
writer.writerow(header)
for row in content:
encoderow = [dd.encode('utf8') for dd in row]
writer.writerow(encoderow)
write_csv('players.csv',result,player_fields)
这里需要注意的就是关于encode的问题。因为我们使用的时utf-8的编码方式，在csv的文件头，需要写入\xEF\xBB\xBF，详见这篇文章
好了现在大功告成，抓取的csv如下图：
因为之前我们还抓取了球员本赛季的比赛详情，所以我们可以进一步的抓取所有球员每一场比赛的记录
抓取的代码如下
def get_player_match(url):
html = urlopen(url).read()
soup = bs4.BeautifulSoup(html, "lxml")
matches = [ dd for dd in soup.select('.shtdm tr') if dd.contents[1].name != 'th']
records = []
for item in [ dd for dd in matches if len(dd.contents) > 11]: ## filter out the personal part
record = []
for match in [ dd for dd in item.contents if type(dd) is bs4.element.Tag]:
if match.string:
record.append(match.string)
else:
for d in [ dd for dd in match.contents if type(dd) is bs4.element.Tag]:
query = dict([(k,v[0]) for k,v in urlparse.parse_qs(d['href']).items()])
record.append('teamid' in query and query['teamid'] or query['id'])
record.append(d.string and d.string or 'na')
records.append(record)
return records[1:] ##remove the first record as the header
def get_players_match(playerlist, baseurl = BASE_URL + '/player.aspx?'):
result = []
for item in playerlist:
url = baseurl + item[10]
print url
result = result + get_player_match(url)
return result
match_fields = ['date_cn','homeid','homename_cn','matchid','score','awayid','awayname_cn','league_cn','firstteam','playtime','goal','assist','shoot','run','corner','offside','foul','violation','yellowcard','redcard','save']
write_csv('m.csv',get_players_match(result),match_fields)
抓取的过程和之前类似。
下一步做什么
现在我们拥有了详细的欧洲联赛的数据，那么下一步要怎么做呢，我推荐大家把数据导入BI工具来做进一步的分析。有两个比较好的选择：
Tableau在数据可视化领域可谓无出其右，Tableau Public完全免费，用数据可视化来驱动数据的探索和分析，拥有非常好的用户体验
当然你也可以用Excel。 另外大家如果有什么好的免费的数据分析的平台，欢迎交流。
订阅：
博文 (Atom)