汽车是现代社会不可或缺的交通工具。借助汽车带来的便利,人们可以打破距离的限制,更好地享受生活. Vast14向三个行业的出版商、市场和搜索引擎提供数据,这三个行业为汽车、房地产和休闲、住宿和旅游. 本案例提供一份来自 Vast 的二手汽车数据集. 数据集共有 417 个样本,包含 6 个特征。这些数据由用户手工录入,容易受到人为失误的影响,比如用户在错误的字段中填写信息,无意中发生错误或者输入胖手指值等.

案例将首先对二手汽车价格数据集进行预处理，并拟合线性回归模型，考察衡量汽车的属性和汽车价格之间的关系。然后再使用可视化，统计，LOF等方法检测数据集中是否含有异常值。

离群点（Outliers），简单而言就是离其余数据点非常远的数据点。它们会极大的影响后续的分析结果，甚至产生有误导的分析结果。

Vast向3个行业的出版商、市场和搜索引擎提供数据，这三个行业包括汽车、房地产和休闲、住宿和旅游。Vast的系统通过白标签集成，并在一些非常受欢迎的消费应用程序（Southwest GetAway Finder，AOL Travel，Yahoo! Travel，Car and Driver等等）中提供搜索结果、产品建议和特别优惠。

Vast的汽车数据是由成千上万的二手汽车卖家提供，并且向市场公布。由于这些数据是由用户手工录入，因此容易受到人为失误的影响，比如用户在错误的字段中提交值，或者无意中发生错误或胖手指值。对于8岁的车辆，里程表读数为100,000英里。直觉告诉我们，100,000美元是大多数小型车的不寻常的价格。一个上市的42000美元是合理的，比如2013年的凯迪拉克ATS豪华版，而对另一个（例如1997年的别克莱斯布雷）来说，这可能是意料之外的。

图1： 2013年的凯迪拉克ATS豪华版与1997年的别克莱斯布雷汽车信息

检测离群值有利于纠正错误，向用户提供卓越和合适的产品。

1. 数据集描述¶

数据集包括训练数据集和测试数据集：

accord_sedan_testing.csv
accord_sedan_training.csv

测试数据集accord_sedan_training.csv包含了417辆本田雅阁轿车的信息列表。各个特征如下：

特征名称	含义	类型	取值示例
$price$	价格	int	14995
$mileage$	已行驶英里数	int	67697
$year$	上市年份	int	2006
$trim$	档次	str	ex;：高配款且带皮革内饰 ex：高配款 lx：低配款
$engine$	引擎缸数	int	4 Cyl：4缸 6 Cyl：6缸
$transmission$	换挡方式	str	Manual：手动挡 Automatic：自动挡

3. 导入数据集并切分¶

import pandas as pd # 读取数据表并进行基于DataFrame结构的操作
import numpy as np
import seaborn as sns
from matplotlib import pyplot as plt
from sklearn import preprocessing
from sklearn.feature_extraction import DictVectorizer
from sklearn.linear_model import LinearRegression
%matplotlib inline

import warnings
warnings.filterwarnings('ignore') # 不显示warning信息
pd.options.display.width = 900 # Dataframe 显示宽度设置

画图配色方案¶

绿色 #5dbe80
蓝色 #2d9ed8
紫色 #a290c4
橙色 #efab40
红色 #EE5150

train = pd.read_csv('./input/accord_sedan_training.csv')
print 'The shape is',train.shape
train.head(7)

The shape is (417, 6)

# 将 train 划分为 x 和 y
train_x = train.drop('price',1)
train_y = train['price']

4. 特征提取并构建线性回归模型¶

使用sklearn.feature_extraction中的DictVectorizer将名义型变量实现One-Hot编码，得到一个10维向量空间，包括：

engine=4 Cyl
engine=6 Cyl
mileage
price
transmission=Automatic
transmission=Manual
trim=ex
trim=exl
trim=lx
year

注意：在异常值检测时需要利用到$~price~$特征。

# One-Hot编码方法一：使用`DictVectorizer`实现
dv = DictVectorizer()
dv.fit(train.T.to_dict().values()) # one-hot 编码
print 'Dimension is',len(dv.feature_names_)
dv.feature_names_

Dimension is 10

['engine=4 Cyl',
 'engine=6 Cyl',
 'mileage',
 'price',
 'transmission=Automatic',
 'transmission=Manual',
 'trim=ex',
 'trim=exl',
 'trim=lx',
 'year']

也可以使用pandas的get_dummies函数实现one-hot编码：

# One-Hot编码方法二：使用`pandas`的`get_dummies`函数实现
nomial_var = ['engine','trim','transmission']
multi_dummies = [] # 存储三个 DataFrame
train_x_dummies = train_x[['mileage','year']]
for col in nomial_var:
    dummies = pd.get_dummies(train_x[col], prefix = col)
    train_x_dummies = pd.concat([train_x_dummies,dummies],axis=1) # 将编码结果与非编码特征水平拼接起来
train_x_dummies.head()

# 构建线性回归模型
train_x_array = dv.transform(train_x.T.to_dict().values())
# train_x_array = train_x_dummies.values # 也可以使用get_dummies得到的结果

LR = LinearRegression().fit(train_x_array,train_y)
' + '.join([format(LR.intercept_, '0.2f')] 
           + map(lambda (f,c): "(%0.2f %s)" % (c, f), 
               zip(dv.feature_names_, LR.coef_)))

'12084.24 + (-337.20 engine=4 Cyl) + (337.20 engine=6 Cyl) + (-0.05 mileage) + (0.00 price) + (420.68 transmission=Automatic) + (-420.67 transmission=Manual) + (208.93 trim=ex) + (674.60 trim=exl) + (-883.53 trim=lx) + (2.23 year)'

即拟合得到的模型为：

$price≈12084.24 − 337.20(engine=4Cyl) + 337.20(engine=6Cyl) - 0.05(mileage) + 420.68(transmission=Automatic) - 420.67(transmission=Manual) + 208.93(trim=ex) + 674.60(trim=exl) - 883.53(trim=lx) + 2.23(year)$

5. 离群值检测¶

由测试集拟合得到的模型，我们可以预测测试集中的价格，计算每个样本的绝对误差，并得出

pred_y = LR.predict(dv.transform(train.T.to_dict().values()))
train_error = abs(pred_y - train_y) # 计算绝对误差
np.percentile(train_error,[75,90,95,99]) # 计算绝对误差数据的百分位数

array([ 1391.7140425 ,  2200.1916772 ,  2626.94049427,  3857.45969199])

我们也可以画出 train_error 的盒图来观察离群值。 boxplot

sns.boxplot(x = train_error,palette = "Set2")

<matplotlib.axes._subplots.AxesSubplot at 0x7f8e19335290>

可以看到测试集中存在部分离群值。

在本案例中，我们设定置信水平为0.95，即认为超过95%百分位数的train_error为离群值。下面我们在二维空间中画出正常值（蓝色）与离群值（红色）：

outlierIndex = train_error >= np.percentile(train_error, 95)
inlierIndex = train_error < np.percentile(train_error, 95)

# 得到train_error最大的index值，即极端离群值
most_severe = train_error[outlierIndex].idxmax() 

fig = plt.figure(figsize=(7,7))
indexes = [inlierIndex, outlierIndex, most_severe]
color = ['#2d9ed8','#EE5150','#a290c4']
label = ['normal points', 'outliers', 'extreme outliers']
for i,c,l in zip(indexes,color,label):
    plt.scatter(train['mileage'][i], 
                train_y[i], 
                c=c,
                marker='^',
                label=l)
plt.legend(loc = 'upper right',
           frameon=True,
           framealpha=1,
           fontsize = 12)
plt.xlabel('$mileage$')
plt.ylabel('$price$')
plt.grid('on')
sns.set_style('dark')

我们来看看离群值的数量有多少？

outlierIndex.value_counts()

False    396
True      21
Name: price, dtype: int64

上图结果也符合我们的经验理解，二手车的行驶公里数越高，它卖出去的价格就应该越低，所以对于处在右上和左下区域的点可能是一些离群值（对于同一款车而言）。比如左下区域的点，一些行驶里程数低，价格也比较低的车辆，有可能该车辆是事故车辆或者有损坏，而右上区域的离群值有可能是真实的离群值，相对来讲不容易有合理的解释，可能是输入失误或者胖手指输入造成。

本案例中的数据只有400多条，如果数据再多一些，则检测的结果会更加可靠。

6. 标准化对离群值检测的影响¶

通常情况下，为了避免不同尺度的影响。我们在进行线性回归模型拟合之前，需要对各个特征进行标准化。常见的标准化有z-score标准化、0-1标准化等，这里我们选择z-score标准化来观察标准化对离群值检测的影响。

# 利用 preprocessing.scale函数将特征标准化
columns = train_x_dummies.columns
train_x_zscore = pd.DataFrame(preprocessing.scale(train_x_dummies),columns = columns)
#train_y_zscore = pd.DataFrame(preprocessing.scale(pd.DataFrame(train_y,columns=['price'])),columns = ['price'])

# 线性模型拟合
LR_zscore = LinearRegression().fit(train_x_zscore.values,train_y)
' + '.join([format(LR_zscore.intercept_, '0.2f')] 
           + map(lambda (f,c): "(%0.2f %s)" % (c, f), 
               zip(dv.feature_names_, LR_zscore.coef_)))

'12084.24 + (-1353.40 engine=4 Cyl) + (-0.00 engine=6 Cyl) + (-166.91 mileage) + (166.91 price) + (225.03 transmission=Automatic) + (138.41 transmission=Manual) + (-274.18 trim=ex) + (116.65 trim=exl) + (-116.65 trim=lx)'

pred_y_zscore = LR_zscore.predict(train_x_zscore)
train_error_zscore = abs(pred_y_zscore - train_y) # 计算绝对误差
np.percentile(train_error_zscore,[75,90,95,99]) # 计算绝对误差数据的百分位数

array([ 1391.71708207,  2200.19426726,  2626.93763764,  3857.46054116])

outlierIndex_zscore = train_error_zscore >= np.percentile(train_error_zscore, 95)
inlierIndex_zscore = train_error_zscore < np.percentile(train_error_zscore, 95)
diff = (outlierIndex_zscore != outlierIndex)  # diff 用于存储标准化前后的离群值检测结果不同的index
diff.value_counts()

False    417
Name: price, dtype: int64

# 画出标准化前后的检测差异点
fig = plt.figure(figsize=(7,7))

# rep_inlierIndex为标准化前后都为正常值的index
rep_inlierIndex = (inlierIndex == inlierIndex_zscore)

indexes = [rep_inlierIndex, outlierIndex, outlierIndex_zscore]
color = ['#2d9ed8','#EE5150','#a290c4']
markers = ['^','<','>']
label = ['inliers', 'outliers before z-score', 'outliers after z-score']
for i,c,m,l in zip(indexes,color,markers,label):
    plt.scatter(train['mileage'][i], 
                train_y[i], 
                c=c,
                marker=m,
                label=l)
plt.xlabel('$mileage$')
plt.ylabel('$price$')
plt.grid('on')
plt.legend(loc = 'upper right',
           frameon=True,
           framealpha=1,
           fontsize = 12)
sns.set_style('dark')

从结果可以看到，绝大多数样本的检测结果一致。有两个样本存在差别，其中一个样本在标准化之前会被检测为离群值，另外一个样本在标准化之后会被检测为离群值。虽然在本例中，标准化前后的检测效果差异不是很大，我们仍然建议在线性建模之前对特征进行标准化。

7. 测试集的验证¶

我们先以$mileage$为横坐标，$price$为纵坐标画出训练集和测试集的所有样本点。

test = pd.read_csv('./input/accord_sedan_testing.csv')

datasets = [train,test]
color = ['#2d9ed8','#EE5150']
label = ['training set', 'testing set']
fig = plt.figure(figsize=(7,7))
for i,c,l in zip(range(len(datasets)),color,label):
    plt.scatter(datasets[i]['mileage'], 
                datasets[i]['price'], 
                c=c,
                marker='^',
                label=l)
plt.xlabel('$mileage$')
plt.ylabel('$price$')
plt.grid('on')
plt.legend(loc = 'upper right',
           frameon=True,
           framealpha=1,
           fontsize = 12)
sns.set_style('dark')

我们来看看利用在训练集上训练得到的模型在测试集上的泛化效果：

pred_y_test = LR.predict(dv.transform(test.T.to_dict().values()))
test_error = abs(pred_y_test - test['price'])

# 使用分布图观察测试集误差
fig = plt.figure(figsize=(7,7))
sns.distplot(test_error,kde=False)
plt.xlabel('$test\_error$')
plt.ylabel('$count$')
plt.grid('on')

# 找出极端离群值
most_severe_test = test_error.idxmax()
test.iloc[most_severe_test]

price                2612
mileage             59308
year                 2006
trim                   ex
engine              6 Cyl
transmission    Automatic
Name: 49, dtype: object

从分布图中可以看到，我们的模型对测试集上其中一个样本的预测表现非常差。该样本是一个极端的离群样本。该车是一个6缸高配版的车，并且其已行驶英里数只有~60,000英里左右，但是其卖出的价格才\$2612。

根据经验，我们猜测这个离群样本出现的两种可能：

在网站里填写时出错；
该车辆有车体的损伤或者有汽车所有权问题（偷来的或者劫来的）

8. 在测试集上使用LOF进行离群值检测¶

test = pd.read_csv('./input/accord_sedan_testing.csv')
fig = plt.figure(figsize=(7,7))
plt.scatter(test['mileage'], 
            test['price'], 
            c='#EE5150',
            marker='^',
            label='testing set')
plt.xlabel('$mileage$')
plt.ylabel('$price$')
plt.grid('on')
plt.legend(loc = 'upper right',
           frameon=True,
           framealpha=1,
           fontsize = 12)
sns.set_style('dark')

可达密度计算公式¶

$$lrd_k(x) = \left( \frac{1}{k} \sum_{y \in N_k(x)} rd_k(x, y) \right)^{-1}$$

其中 $rd_k(x, y)$ 样本点$x$到样本点$y$的第$k$可达距离

LOF因子计算公式¶

$$lof_{k}(x) = \frac{1}{k} \sum_{y \in N_k(x)} \frac{lrd_k(y)}{lrd_k(x)}$$

LOF算法的一般流程可以描述为：

初始化$k$，用于后续计算第$k$距离；
计算每个样本点与其他点的距离，并对其排序；
计算每个样本点的第$k$距离，第$k$领域；
计算每个样本点的可达密度以及LOF值；
对所有样本点的LOF值进行排序，与1作比较，越大于1，越可能是离群值。

data = test[['mileage','price']]

# 导入sklearn用于计算最近邻相关数据的方法 NearestNeighbors
from sklearn.neighbors import NearestNeighbors
import numpy as np

# 定义函数计算第k可达距离
def k_Distance(data, k):
    neigh = NearestNeighbors(k)
    model = neigh.fit(data)
    
    nums = data.shape[0]

    k_distance = []
    neighbor_info = []
    
    for index in xrange(nums):
        
        # K个距离
        dist = neigh.kneighbors([data[index]], n_neighbors = k + 1)
        # 最大的dist
        k_distance.append(dist[0][-1][-1])
        
        # neighbor为k个邻居的索引，dists为相应的距离
        dists, neighbor = neigh.radius_neighbors([data[index]], radius = k_distance[index])
        
        # 排除自身
        mask = neighbor[0] != index     
        neighbor_info.append([neighbor[0][mask], dists[0][mask]])
    return k_distance, neighbor_info

# 定义函数计算局部可达密度
def reach_density(data, k_distance, neighbor_info):
    
    nums = data.shape[0]
    density_list = []
    
    for index in xrange(nums):
        
        neighbors, dists = neighbor_info[index]
        nums_neigh = len(neighbors)
    
        sum_dist = 0
        
        for item in xrange(nums_neigh):
            
            k_dist_o = k_distance[neighbors[item]] 
    
            direct_dist = dists[item]
        
            reach_dist_p_o = max(k_dist_o, direct_dist)
        
            sum_dist += reach_dist_p_o
           
        density_list.append(nums_neigh/sum_dist)
    
    return density_list

# 定义函数计算LOF因子
def cal_lof(index, data, neighbor_info, lrd_list):
    
    point_p = data[index]
    neighbors, _ = neighbor_info[index]
    
    nums_neigh = len(neighbors)
    sum_density = 0
    
    for item in xrange(nums_neigh):
        
        sum_density += lrd_list[item]

    return sum_density/(nums_neigh*lrd_list[index])

dists, neighbor_info = k_Distance(data.values, 2)

lrd_list = reach_density(data, dists, neighbor_info)

nums = data.shape[0]

lof_list = []

for index in xrange(nums):
    
    lof = cal_lof(index, data.values, neighbor_info, lrd_list)
    lof_list.append(lof)

boolean_array = [item > 5 for item in lof_list]

indicy = []
for key, value in enumerate(boolean_array):
    if value:
        indicy.append(key)

print indicy

[49, 52]

画出离群值：

fig = plt.figure(figsize=(7,7))
for i in data.index:
    if i not in indicy:
        plt.scatter(data.iloc[i]['mileage'], 
                data.iloc[i]['price'], 
                c='#2d9ed8',
                s=50,
                marker='^',
                label='inliers')
    else:
        plt.scatter(data.iloc[i]['mileage'], 
                data.iloc[i]['price'], 
                c='#EE5150',
                s=50,
                marker='^',
                label='outliers')
        plt.xlabel('$mileage$')
        plt.ylabel('$price$')
        plt.grid('on')

另外一种比较乱的写法¶

test_2d = test[['mileage','price']]
from sklearn.neighbors import NearestNeighbors
import numpy as np
neigh = NearestNeighbors(5) # 默认为欧式距离
model = neigh.fit(test_2d)

data = test_2d
# dist为每个样本点与第k距离邻域内的点的距离（包括自身）,neighbor为第k距离邻域点的编号（包括自身）
dist, neighbor=neigh.kneighbors(test_2d,n_neighbors=6)

k_distance_p = np.max(dist,axis=1)

nums = data.shape[0]
lrdk_p = []
lof = []
for p_index in xrange(nums):
    rdk_po = []
    neighbor_p = neighbor[p_index][neighbor[p_index]!=p_index]
    for o_index in neighbor_p:
        rdk_po.append(max(k_distance_p[o_index],int(dist[p_index][neighbor[p_index]==o_index])))
    lrdk_p.append(float(len(neighbor_p))/sum(rdk_po))

for p_index in xrange(nums):
    lrdk_o=[]
    neighbor_p = neighbor[p_index][neighbor[p_index]!=p_index]
    for o_index in neighbor_p:
        lrdk_o.append(lrdk_p[o_index])
    lof.append(float(sum(lrdk_o))/(len(neighbor_p)*(lrdk_p[p_index])))    

fig = plt.figure(figsize=(7,7))

for index,size in zip(xrange(nums),lof):
    if index in indicy:
        plt.scatter(data['mileage'][index], 
                    data['price'][index], 
                    s=np.exp(lof[index])*50,
                    c='#efab40',
                    alpha=0.6,
                    marker='o')
        plt.text(data['mileage'][index]-np.exp(lof[index])*50, 
                 data['price'][index]-np.exp(lof[index])*50,
                str(round(lof[index],2)))
    else:
        plt.scatter(data['mileage'][index], 
                    data['price'][index], 
                    s=np.exp(lof[index])*50,
                    c='#5dbe80',
                    alpha=0.6,
                    marker='o')
        plt.text(data['mileage'][index]-np.exp(lof[index])*50, 
                 data['price'][index]-np.exp(lof[index])*50,
                str(round(lof[index],2)),
                fontsize=7)

plt.xlabel('mileage')
plt.ylabel('price')
plt.grid('off')

	mileage	year	engine_4 Cyl	trim_ex	trim_lx	transmission_Automatic	transmission_Manual
0	67697	2006	1.0	1.0	0.0	0.0	1.0
1	73738	2006	1.0	1.0	0.0	0.0	1.0
2	80313	2006	1.0	0.0	1.0	1.0	0.0
3	86096	2006	1.0	0.0	1.0	1.0	0.0
4	79607	2006	1.0	0.0	1.0	1.0	0.0