模型评估与模型优化

第6章 模型评估与模型优化

6.1模型评估一般方法

对模型的性能或优劣进行评估,通过评价指标来衡量。模型评估是机器学习任务中非常重要的一环。不同的机器学习任务有着不同的评价指标,同时同一种机器学习任务也有着不同的评价指标,每个指标的着重点不一样。如分类(classification)、回归(regression)、聚类(clustering)、SVM、推荐(recommendation)等。并且很多指标可以对多种不同的机器学习模型进行评价,如精确率-召回率(precision-recall),可以用在分类、推荐等。像分类、回归、SVM都是监督式机器学习,本节的重点便是监督式机器学习的一些评价指标。
在实际情况中,我们会用不同的度量去评估我们的模型,而度量的选择,完全取决于模型的类型和模型以后要做的事。下面我们就会学习到一些用于评价模型的常用度量和图表以及它们各自的使用场景。

6.1.1分类评价指标
分类是指对给定的数据记录预测该记录所属的类别。并且类别空间已知。它包括二分类与多分类,二分类便是指只有两种类别,其输出不是0就是1;或是与否等,也可以称为“负”(negative)与正(positive)两种类别,一般在实际计算中,将其映射到“0”-“1” class中;而多分类则指类别数超过两种。下面主要根据二分类的评价指标进行讲解,不过同时它们也可扩展到多分类任务中。下面对分类中一些常用的评价指标进行介绍。
1)混淆矩阵(Confusion Matrix)
混淆矩阵显示了分类模型相对数据的真实输出(目标值)的正确预测和不正确预测数目。矩阵为NxN,其中N为目标值(类)数目。这类模型的性能通常使用矩阵中的数据评估。下表为两个类别(阳性和阴性)的2x2混淆矩阵。

常用术语:
阳性 (P, positive)
阴性 (N, Negative)
真阳性 (TP, true positive):正确的肯定。又称:命中 (hit)
真阴性 (TN, true negative):正确的否定。又称:正确拒绝 (correct rejection)
伪阳性 (FP, false positive):错误的肯定,又称:假警报 (false alarm)、第二型错误
伪阴性 (FN, false negative):错误的否定,又称:未命中(miss)、第一型错误
灵敏度(sensitivity)或真阳性率(TPR, true positive rate):sensitivity = TP/(TP + FN);
灵敏度又称为召回率(Recall)=TP/(TP+FN) = 1 - FN/P
特异性(真阴性率)specificity = TN /(FP + TN)
准确度 (ACC, accuracy):ACC = (TP + TN) / (P + N)
预测正确的数占样本数的比例。
F1评分:
精度和灵敏度的调和平均数。
F1 = 2 precision * recall / (precision+recall) = 2TP/(2TP+FP+FN)

示例:

混淆矩阵 目标
阳性 阴性
模型 阳性 70 20 阳性预测值 0.78
阴性 30 80 阴性预测值 0.73
灵敏度 特异度 准确度 = 0.75
0.7 0.8

2)ROC曲线
判定方法:ROC曲线应尽量偏离参考线。
原理:ROC全称为Receiver Operation Characteristic Curve,ROC曲线其实就是从混淆矩阵衍生出来的图形,其横坐标为1-Specificity,纵坐标为Sensitivity。

上面那条曲线就是ROC曲线,随着阈值的减小,更多的值归于正类,敏感度和1-特异度也相应增加,所以ROC曲线呈递增趋势。而那条45度线是一条参照线,也就是说ROC曲线要与这条曲线比较。
简单的说,如果我们不用模型,直接随机把客户分类,我们得到的曲线就是那条参照线,然而我们使用了模型进行预测,就应该比随机的要好,所以ROC曲线要尽量远离参照线,越远,我们的模型预测效果越好。
3)AUC(ROC曲线下面积)
判定方法:AUC应该大于0.5.
原理:ROC曲线是根据与那条参照线进行比较来判断模型的好坏,但这只是一种直觉上的定性分析,如果我们需要精确一些,就要用到AUC,也就是ROC曲线下面积。

看上图,参考线的面积是0.5,ROC曲线与它偏离越大,ROC曲线就越往左上方靠拢,它下面的面积(AUC)也就越大,这里面积是0.869。我们可以根据AUC的值与0.5相比,来评估一个分类模型的预测效果。
6.1.2回归评估指标
均方差(MSE,Mean Squared Error):
(∑(prec-act)**2)/n (prec为预测值,act为实际值,n为总样本数)
均方根差(RMSE,Root Mean Squared Error):
就是MSE开根号
平均绝对值误差(MAE,Mean Absolute Error):
(∑|prec-act|)/n
6.1.3评估指标在Sciktlearn中对应的函数
评估分类模型:

评估回归模型:

6.2模型评估与优化

对模型评估和优化时,我们常常使用训练集-验证集二划分验证(Hold-out validation)、K折交叉验证(k-Cross-validation)、超参数调优(hyperparameter tuning)等方法。三者从不同的层次对机器学习模型进行校验。Hold-out validation是评估机器学习模型泛化性能一个经典且常用的方法,通过这种方法,我们将把数据集划分为训练集和测试集,前者用于模型训练,后者用于性能的评估。
6.2.1Hold-out validation方法
使用Hold-out validation还可以采用更好的一种划分方法,把数据划分为三个部分,训练集、验证集和测试集,训练集用于模型的拟合,模型在验证集上的性能表现作为模型选择的标准,模型评估在测试集上。Hold-out validation的一般流程图如下:

6.2.2 K-fold交叉验证
一种更好,但是计算量更大的交叉验证方法是K-fold交叉验证。如同Holdout方法,K-fold交叉验证也依赖于训练数据的若干个相互独立子集。主要的区别在于K-fold交叉验证一开始就随机把数据分割成K个不相连的子集,成为folds(一般称作K折交叉验证,K的取值有5、10或者20)。每次留一份数据作为测试集,其余数据用于训练模型。
当每一份数据都轮转一遍之后,将预测的结果整合,并与目标变量的真实值比较来计算准确率。K-fold交叉验证的图形展示如下:

以下我们通过实例来进一步说明K-折交叉验证的使用方法及其优越性。这里我们使用UCI数据集中一个有关威斯康星乳腺癌数据集。
1)读取数据

import pandas as pd

df = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/wdbc.data', header=None)
df.head()

部分数据:

2)对标签数据转换为整数,目前是字符M(恶性)或B(良性)

from sklearn.preprocessing import LabelEncoder
X = df.loc[:, 2:].values
y = df.loc[:, 1].values
le = LabelEncoder()
y = le.fit_transform(y)
le.transform(['M', 'B'])

转换结果如下:
array([1, 0])
3)把数据划分为训练数据和测试数据,划分比例为7:3

from sklearn.cross_validation import train_test_split

X_train, X_test, y_train, y_test = \
train_test_split(X, y, test_size=0.20, random_state=1)

4)使用PCA对原数据集进行降维
这里我们采用管道或流水线方式,把对数据的标准化等组装在一起。

from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import Pipeline

pipe_lr = Pipeline([('scl', StandardScaler()),
('pca', PCA(n_components=2)),
('clf', LogisticRegression(random_state=1))])

5)在训练集上进行k-折交叉验证
这里我们取k=10,采用sklearn中StratifiedKFold函数。

import numpy as np
from sklearn.cross_validation import StratifiedKFold

kfold = StratifiedKFold(y=y_train,
n_folds=10,
random_state=1)

scores = []
for k, (train, test) in enumerate(kfold):
pipe_lr.fit(X_train[train], y_train[train])
score = pipe_lr.score(X_train[test], y_train[test])
scores.append(score)
print('Fold: %s, Class dist.: %s, Acc: %.3f' % (k+1, np.bincount(y_train[train]), score))

print('\nCV accuracy: %.3f +/- %.3f' % (np.mean(scores), np.std(scores)))

运行结果如下:
Fold: 1, Class dist.: [256 153], Acc: 0.891
Fold: 2, Class dist.: [256 153], Acc: 0.978
Fold: 3, Class dist.: [256 153], Acc: 0.978
Fold: 4, Class dist.: [256 153], Acc: 0.913
Fold: 5, Class dist.: [256 153], Acc: 0.935
Fold: 6, Class dist.: [257 153], Acc: 0.978
Fold: 7, Class dist.: [257 153], Acc: 0.933
Fold: 8, Class dist.: [257 153], Acc: 0.956
Fold: 9, Class dist.: [257 153], Acc: 0.978
Fold: 10, Class dist.: [257 153], Acc: 0.956

CV accuracy: 0.950 +/- 0.029
准确率达到95%,说明结果还不错。
6)改进的k-折交叉验证
k-折交叉验证比较耗资源,而且抽取数据是随机的,这就可能导致有些分发中的类别比例不很恰当,是否有更好方法?以下我们采用分层k-折交叉验证法,使用scikit-learn中cross_val_score。

from sklearn.cross_validation import cross_val_score

scores = cross_val_score(estimator=pipe_lr,
X=X_train,
y=y_train,
cv=10,
n_jobs=1)
print('CV accuracy scores: %s' % scores)
print('CV accuracy: %.3f +/- %.3f' % (np.mean(scores), np.std(scores)))

运算结果如下:
CV accuracy scores: [ 0.89130435 0.97826087 0.97826087 0.91304348 0.93478261 0.97777778 0.93333333 0.95555556 0.97777778 0.95555556]
CV accuracy: 0.950 +/- 0.029

n_jobs参数很特别,使用这个参数,我们可以把K轮的交叉验证分布到几个CPU块上,如果n_jobs=2就是2块,如果n_jobs=-1,则系统将利用所有CPU进行计算。

6.3 利用验证曲线查看是否过拟合或欠拟合

验证曲线反应参数与性能指标之间的关系图,通过验证曲线图有利于定位过拟合或欠拟合等问题,是提高模型性能的有效工具。以下我们以建立在上个数据集的逻辑斯谛模型中的正则化参数C与准确率间的关系为例,使用validation_curve函数,该函数对分类算法,默认使用分层K折交叉验证来评估模型性能,使用scikit-learn来绘制验证曲线。

import matplotlib.font_manager as fm
myfont = fm.FontProperties(fname='/home/hadoop/anaconda3/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/simhei.ttf')

from sklearn.learning_curve import validation_curve

param_range = [0.001, 0.01, 0.1, 1.0, 10.0, 100.0]
train_scores, test_scores = validation_curve(
estimator=pipe_lr,
X=X_train,
y=y_train,
param_name='clf__C',
param_range=param_range,
cv=10)

train_mean = np.mean(train_scores, axis=1)
train_std = np.std(train_scores, axis=1)
test_mean = np.mean(test_scores, axis=1)
test_std = np.std(test_scores, axis=1)

plt.plot(param_range, train_mean,
color='blue', marker='o',
markersize=5, label='训练准确率')

plt.fill_between(param_range, train_mean + train_std,
train_mean - train_std, alpha=0.15,
color='blue')

plt.plot(param_range, test_mean,
color='green', linestyle='--',
marker='s', markersize=5,
label='验证准确率')

plt.fill_between(param_range,
test_mean + test_std,
test_mean - test_std,
alpha=0.15, color='green')

plt.grid()
plt.xscale('log')
plt.legend(loc='lower right',prop=myfont)
plt.xlabel('参数C',fontproperties=myfont,size=12)
plt.ylabel('准确率',fontproperties=myfont,size=12)
plt.ylim([0.8, 1.0])
plt.tight_layout()
plt.show()


从以上图形不难看出,参数C的最优点是0.1附近,如果加大正则化强度(C较小),会导致一定程度的欠拟合;如果降低正则化强度,可能导致模型的过拟合。

6.4 利用参数网格搜索调优模型

在机器学习中,有两类参数:一种是通过训练数据得到的参数,如回归模型、神经网络中系统;一种是需要单独进行优化的参数,如正则化系统、迭代次数、树的深度等,这类参数称为调优参数,又称为超参。
这里我们将介绍一种功能强大的超参数优化方法,即网格搜索(grid search)方法,通过这种方法自动选择最优参数组合,并选择最优模型。以下利用Scikit learn中GridSearchCV函数来实现网格搜索。

from sklearn.grid_search import GridSearchCV
from sklearn.svm import SVC

pipe_svc = Pipeline([('scl', StandardScaler()),
('clf', SVC(random_state=1))])

param_range = [0.0001, 0.001, 0.01, 0.1, 1.0, 10.0, 100.0, 1000.0]

param_grid = [{'clf__C': param_range,'clf__gamma': param_range, 'clf__kernel': ['linear','rbf']}]
gs = GridSearchCV(estimator=pipe_svc,
param_grid=param_grid,
scoring='accuracy',
cv=10,
n_jobs=-1)
gs = gs.fit(X_train, y_train)
print(gs.best_score_)
print(gs.best_params_)

运行结果:
0.978021978021978
{'clf__C': 0.1, 'clf__gamma': 0.0001, 'clf__kernel': 'linear'}

由此可知,最优参数组合是:clf__C=0.1 , clf__kernel='linear', clf__gamma=0.0001
这个组合对应模型性能达到近98%。
此外,还有随机搜索(randomized search)、嵌套交叉验证选择模型等方法。大家可作为练习去练习。这里就不展开来说了。

数据探索和预处理

1、整体思路:

理解问题:查看每个变量并且根据他们的意义和对问题的重要性进行哲学分析。
单因素研究:只关注因变量(‘SalePrice’),并且进行更深入的了解。
多因素研究:分析因变量和自变量之间的关系。、
基础清洗:清洗数据集并且对缺失数据,异常值和分类数据进行一些处理。
检验假设:检查数据是否和多元分析方法的假设达到一致。
开始之前,先导入需要的库及文件:

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from scipy.stats import norm
from sklearn.preprocessing import StandardScaler
from scipy import stats
import warnings
warnings.filterwarnings('ignore')
%matplotlib inline

#bring in the six packs
df_train = pd.read_csv('./data/house_train.csv')

查看特征

df_train.columns
Index(['Id', 'MSSubClass', 'MSZoning', 'LotFrontage', 'LotArea', 'Street', 'Alley', 'LotShape', 'LandContour', 'Utilities', 'LotConfig','LandSlope', 'Neighborhood', 'Condition1', 'Condition2', 'BldgType','HouseStyle', 'OverallQual', 'OverallCond', 'YearBuilt', 'YearRemodAdd','RoofStyle', 'RoofMatl', 'Exterior1st', 'Exterior2nd', 'MasVnrType','MasVnrArea', 'ExterQual', 'ExterCond', 'Foundation', 'BsmtQual','BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinSF1','BsmtFinType2', 'BsmtFinSF2', 'BsmtUnfSF', 'TotalBsmtSF', 'Heating','HeatingQC', 'CentralAir', 'Electrical', '1stFlrSF', '2ndFlrSF','LowQualFinSF', 'GrLivArea', 'BsmtFullBath', 'BsmtHalfBath', 'FullBath','HalfBath', 'BedroomAbvGr', 'KitchenAbvGr', 'KitchenQual','TotRmsAbvGrd', 'Functional', 'Fireplaces', 'FireplaceQu', 'GarageType','GarageYrBlt', 'GarageFinish', 'GarageCars', 'GarageArea', 'GarageQual','GarageCond', 'PavedDrive', 'WoodDeckSF', 'OpenPorchSF','EnclosedPorch', '3SsnPorch', 'ScreenPorch', 'PoolArea', 'PoolQC','Fence', 'MiscFeature', 'MiscVal', 'MoSold', 'YrSold', 'SaleType','SaleCondition', 'SalePrice'], dtype='object')

2、准备工作

为了了解我们的数据,我们可以分析每个变量并且尝试理解他们的意义和与该问题的相关程度。
首先建立一个 Excel 电子表格,有如下目录:

变量 – 变量名。
类型 – 该变量的类型。这一栏只有两个可能值,“数据” 或 “类别”。 “数据” 是指该变量的值是数字,“类别” 指该变量的值是类别标签。
划分 – 指示变量划分. 我们定义了三种划分:建筑,空间,位置。
期望 – 我们希望该变量对房价的影响程度。我们使用类别标签 “高”,“中” 和 “低” 作为可能值。
结论 – 我们得出的该变量的重要性的结论。在大概浏览数据之后,我们认为这一栏和 “期望” 的值基本一致。
评论 – 我们看到的所有一般性评论。
我们首先阅读了每一个变量的描述文件,同时思考这三个问题:
我们买房子的时候会考虑这个因素吗?
如果考虑的话,这个因素的重要程度如何?
这个因素带来的信息在其他因素中出现过吗?
我们根据以上内容填好了电子表格,并且仔细观察了 “高期望” 的变量。然后绘制了这些变量和房价之间的散点图,填在了 “结论” 那一栏,也正巧就是对我们的期望值的校正。

我们总结出了四个对该问题起到至关重要的作用的变量:
OverallQual
YearBuilt.
TotalBsmtSF.
GrLivArea.
最重要的事情——分析 “房价”

df_train['SalePrice'].describe()
count 1460.000000
mean 180921.195890
std 79442.502883
min 34900.000000
25% 129975.000000
50% 163000.000000
75% 214000.000000
max 755000.000000
Name: SalePrice, dtype: float64

绘制直方图

sns.distplot(df_train['SalePrice']);


用seaborn绘图,得出结论:
-偏离正态分布
-数据正偏
-有峰值

数据偏度和峰度度量:

print("Skewness: %f" % df_train['SalePrice'].skew())
print("Kurtosis: %f" % df_train['SalePrice'].kurt())
Skewness: 1.882876
Kurtosis: 6.536282

偏度(Skewness)是描述某变量取值分布对称性的统计量。如果偏度=0,和正态分布的偏度相同;Skewness>0,长尾巴拖在右边;Skewness<0,长尾巴拖在左边,Skewness越大,分布形态偏移程度越大。 峰度(Kurtosis)是描述某变量所有取值分布形态陡缓程度的统计量。它是和正态分布相比较的:Kurtosis=0,与正态分布的陡缓程度相同;Kurtosis>0,比正态分布的高峰更加陡峭;反之亦然。

3、与数值类型的关系

(1)grlivarea/saleprice的散点图绘制
首先观察一下saleprice和数值型变量grlivarea之间的关系。我们还是可以观察到整体的趋势的,随着grlivarea的增大,saleprice有变高的趋势。存在线性关系!

var = 'GrLivArea'
data = pd.concat([df_train['SalePrice'], df_train[var]], axis=1)
data.plot.scatter(x=var, y='SalePrice', ylim=(0,800000));


可以看出’’SalePrice’和’'GrLivArea' 关系很密切,并且基本呈线性关系。

var = 'TotalBsmtSF'
data = pd.concat([df_train['SalePrice'], df_train[var]], axis=1)
data.plot.scatter(x=var, y='SalePrice', ylim=(0,800000));

TotalBsmtSF' 和 'SalePrice'关系也很密切,从图中可以看出基本呈指数分布,但从最左侧的点可以看出特定情况下’TotalBsmtSF' 对'SalePrice' 没有产生影响。

4、与类别类型的关系

(1)overallqual/saleprice的箱形图
这里有一个问题,为什么我们不用散点图了呢?那么先画一个散点图看下效果。如下图:

#box plot overallqual/saleprice
var = 'OverallQual'
data = pd.concat([df_train['SalePrice'], df_train[var]], axis=1)
f, ax = plt.subplots(figsize=(8, 6))
fig = sns.boxplot(x=var, y="SalePrice", data=data)
fig.axis(ymin=0, ymax=800000);


可见也是有相互关系的,如果用箱形图表示的话,会更加清晰:
(2)YearBuilt/saleprice的箱形图
两个变量之间的关系没有很强的趋势性,但是可以看出建筑时间较短的房屋价格更高。

var = 'YearBuilt'
data = pd.concat([df_train['SalePrice'], df_train[var]], axis=1)
f, ax = plt.subplots(figsize=(16, 8))
fig = sns.boxplot(x=var, y="SalePrice", data=data)
fig.axis(ymin=0, ymax=800000);
plt.xticks(rotation=90);

两个变量之间的关系没有很强的趋势性,但是可以看出建筑时间较短的房屋价格更高。

总结:

'GrLivArea' 和 'TotalBsmtSF' 与 'SalePrice'似乎线性相关,并且都是正相关。 对于 'TotalBsmtSF',线性关系的斜率十分的高。
'OverallQual' 和 'YearBuilt' 与 'SalePrice'也有关系。'OverallQual'的相关性更强, 箱型图显示了随着整体质量的增长,房价的增长趋势。

我们只分析了四个变量,但是还有许多其他变量我们也应该分析,这里的技巧在于选择正确的特征(特征选择)而不是定义他们之间的复杂关系(特征工程)。它们之间的复杂关系可以通过模型学习。

5、客观分析

上面的分析过于主观,我们应该更客观的去分析。
主要有三方面的分析:
-相关矩阵热图
-SalePrice相关矩阵热图
-最相关变量之间的散点图

(1)相关矩阵热图

#correlation matrix
corrmat = df_train.corr()
f, ax = plt.subplots(figsize=(12, 9))
sns.heatmap(corrmat, vmax=.8, square=True);


热图是观察特征和特征、特征和label之间关系的一种快速的方式。首先,最吸引注意的是两个红色方块。第一处指向TotalBsmtSF和1stFlrSF这两个变量,第二处指向Garage相关的变量。这两块显示了TotalBsmtSF和1stFlrSF的相关性、GarageCar和GarageArea的相关性很大。事实上,这种关联性很强,暗示它们存在多重共线性。我们能够推断当变量间具有多重共线性,那么它们给出的信息也基本上是一样的。Heatmaps恰恰可以检测多重共线性的情况并决定选择哪些特征,所以是一个非常重要的工具。

另外观察特征和SalePrice之间的关系,可以看到GrLivArea、TotalBsmtSF、OverallQual,也看到了一些其他的变量应该纳入我们考虑。

(2)SalePrice相关矩阵热图
筛选与SalePrice相关性强的特征,重点观察。

#saleprice correlation matrix
k = 10 #number of variables for heatmap
cols = corrmat.nlargest(k, 'SalePrice')['SalePrice'].index
cm = np.corrcoef(df_train[cols].values.T)
sns.set(font_scale=1.25)
hm = sns.heatmap(cm, cbar=True, annot=True, square=True, fmt='.2f', annot_kws={'size': 10}, yticklabels=cols.values, xticklabels=cols.values)
plt.show()

从图中可以看出:
'OverallQual', 'GrLivArea' 以及 'TotalBsmtSF' 与 'SalePrice'有很强的相关性。
'GarageCars' 和 'GarageArea' 也是相关性比较强的变量. 车库中存储的车的数量是由车库的面积决定的,它们就像双胞胎,所以不需要专门区分'GarageCars' 和 'GarageArea' ,所以我们只需要其中的一个变量。这里我们选择了'GarageCars'因为它与'SalePrice' 的相关性更高一些。
'TotalBsmtSF' 和 '1stFloor' 与上述情况相同,我们选择 'TotalBsmtSF' 。
'FullBath'几乎不需要考虑。
'TotRmsAbvGrd' 和 'GrLivArea'也是变量中的双胞胎。
'YearBuilt' 和 'SalePrice'相关性似乎不强。
(3).'SalePrice' 和相关变量之间的散点图

#scatterplot
sns.set()
cols = ['SalePrice', 'OverallQual', 'GrLivArea', 'GarageCars', 'TotalBsmtSF', 'FullBath', 'YearBuilt']
sns.pairplot(df_train[cols], size = 2.5)
plt.show();


尽管我们已经知道了一些主要特征,这一丰富的散点图给了我们一个关于变量关系的合理想法。其中,'TotalBsmtSF' 和 'GrLiveArea'之间的散点图是很有意思的。我们可以看出这幅图中,一些点组成了线,就像边界一样。大部分点都分布在那条线下面,这也是可以解释的。地下室面积和地上居住面积可以相等,但是一般情况下不会希望有一个比地上居住面积还大的地下室。
'SalePrice' 和'YearBuilt' 之间的散点图也值得我们思考。在“点云”的底部,我们可以观察到一个几乎呈指数函数的分布。我们也可以看到“点云”的上端也基本呈同样的分布趋势。并且可以注意到,近几年的点有超过这个上端的趋势。

6、缺失数据

关于缺失数据需要思考的重要问题:
这一缺失数据的普遍性如何?
缺失数据是随机的还是有律可循?
这些问题的答案是很重要的,因为缺失数据意味着样本大小的缩减,这会阻止我们的分析进程。除此之外,以实质性的角度来说,我们需要保证对缺失数据的处理不会出现偏离或隐藏任何难以忽视的真相。

#missing data
total = df_train.isnull().sum().sort_values(ascending=False)
percent = (df_train.isnull().sum()/df_train.isnull().count()).sort_values(ascending=False)
missing_data = pd.concat([total, percent], axis=1, keys=['Total', 'Percent'])
missing_data.head(20)


(1)当超过15%的数据都缺失的时候,我们应该删掉相关变量且假设该变量并不存在。根据这一条,一系列变量都应该删掉,例如'PoolQC', 'MiscFeature', 'Alley'等等,这些变量都不是很重要,因为他们基本都不是我们买房子时会考虑的因素。
(2)'GarageX' 变量群的缺失数据量都相同,由于关于车库的最重要的信息都可以由'GarageCars' 表达,并且这些数据只占缺失数据的5%,我们也会删除上述的'GarageX' 变量群。同样的逻辑也适用于 'BsmtX' 变量群。
(3)对于 'MasVnrArea' 和 'MasVnrType',我们可以认为这些因素并不重要。除此之外,他们和'YearBuilt' 以及 'OverallQual'都有很强的关联性,而这两个变量我们已经考虑过了。所以删除 'MasVnrArea'和 'MasVnrType'并不会丢失信息。

(4)最后,由于'Electrical'中只有一个损失的观察值,所以我们删除这个观察值,但是保留这一变量。

#dealing with missing data
df_train = df_train.drop((missing_data[missing_data['Total'] > 1]).index,1)
df_train = df_train.drop(df_train.loc[df_train['Electrical'].isnull()].index)
df_train.isnull().sum().max() #just checking that there's no missing data missing...

总结:处理空值,我们删掉所有有较多空值的特征,除了特征Electrical,因为只有一个缺失值,所以删除含该空值的样本即可。

7、异常值处理

异常值也是我们应该注意的东西。因为异常值能明显的影响我们的模型,并且是一个有价值的信息来源,帮助我们对特定行为有更多的见解。异常值是一个复杂的主题,并且值得研究。在这里,我们将用SalePrice的标准差和一系列的散点来进行快速分析。

7.1单因素分析
这里的关键在于如何建立阈值,定义一个观察值为异常值。
我们对数据进行正态化,意味着把数据值转换成均值为0,方差为1的数据。

#standardizing data
saleprice_scaled = StandardScaler().fit_transform(df_train['SalePrice'][:,np.newaxis]);
low_range = saleprice_scaled[saleprice_scaled[:,0].argsort()][:10]
high_range= saleprice_scaled[saleprice_scaled[:,0].argsort()][-10:]
print('outer range (low) of the distribution:')
print(low_range)
print('\nouter range (high) of the distribution:')
print(high_range)
outer range (low) of the distribution:
[[-1.83820775]
[-1.83303414]
[-1.80044422]
[-1.78282123]
[-1.77400974]
[-1.62295562]
[-1.6166617 ]
[-1.58519209]
[-1.58519209]
[-1.57269236]]

outer range (high) of the distribution:
[[ 3.82758058]
[ 4.0395221 ]
[ 4.49473628]
[ 4.70872962]
[ 4.728631 ]
[ 5.06034585]
[ 5.42191907]
[ 5.58987866]
[ 7.10041987]
[ 7.22629831]]

进行正态化后,可以看出:
低范围的值都比较相似并且在0附近分布。
高范围的值离0很远,并且7点几的值远在正常范围之外。
现在看,我们不把任何值作为异常值,但是我们应该注意这两个大于7的值。
7.2双变量分析

#bivariate analysis saleprice/grlivarea
var = 'GrLivArea'
data = pd.concat([df_train['SalePrice'], df_train[var]], axis=1)
data.plot.scatter(x=var, y='SalePrice', ylim=(0,800000));


这是之前GrLivArea和SalePrice关系散点图:
(1)有两个离群的'GrLivArea' 值很高的数据,我们可以推测出现这种情况的原因。或许他们代表了农业地区,也就解释了低价。 这两个点很明显不能代表典型样例,所以我们将它们定义为异常值并删除。

(2)图中顶部的两个点是七点几的观测值,他们虽然看起来像特殊情况,但是他们依然符合整体趋势,所以我们将其保留下来。

删除异常点:

df_train.sort_values(by = 'GrLivArea', ascending = False)[:2]
df_train = df_train.drop(df_train[df_train['Id'] == 1299].index)
df_train = df_train.drop(df_train[df_train['Id'] == 524].index)

(3). 'TotalBsmtSF'和'SalePrice'双变量分析

#bivariate analysis saleprice/grlivarea
var = 'TotalBsmtSF'
data = pd.concat([df_train['SalePrice'], df_train[var]], axis=1)
data.plot.scatter(x=var, y='SalePrice', ylim=(0,800000));


我们可以感觉到试图消除一些观察(例如,TotalBsmtSF> 3000),但是我们认为这是不值得的。 我们可以保留它们,所以我们不会做任何事情。

8、核心部分

核心部分
“房价”到底是谁?
这个问题的答案,需要我们验证根据数据基础进行多元分析的假设。
我们已经进行了数据清洗,并且发现了“SalePrice”的很多信息,现在我们要更进一步理解‘SalePrice’如何遵循统计假设,可以让我们应用多元技术。

应该测量4个假设量:
(1)正态性:它的重要在于很多统计检验是基于正态分布的,在房价预测的问题中,我们只检查了单变量的正态性。但是单变量正态性不能确保多变量的正态性,但是会其帮助作用。
(2)同方差性:假设在预测变量的范围因变量表现出同等水平的方差
(3)线性:通过观察散点图,看是否为线性关系,如果不是,需要数据转换,但是大多数情况下都是满足线性关系的。
(4)相关错误缺失
1). ‘SalePrice’
正态性:

应主要关注以下两点:
直方图 - 峰度和偏度。
正态概率图 - 数据分布应紧密跟随代表正态分布的对角线。

#histogram and normal probability plot
sns.distplot(df_train['SalePrice'], fit=norm);
fig = plt.figure()
res = stats.probplot(df_train['SalePrice'], plot=plt)

可以看出,房价分布不是正态的,显示了峰值,正偏度。但是,信息没有丢失,简单的数据变换就可以解决这个问题,为了防止正偏,log变换作用很好。

#applying log transformation
df_train['SalePrice'] = np.log(df_train['SalePrice'])

#transformed histogram and normal probability plot
sns.distplot(df_train['SalePrice'], fit=norm);
fig = plt.figure()
res = stats.probplot(df_train['SalePrice'], plot=plt)

从图中可以看出:
显示出了偏度
大量为0的观察值(没有地下室的房屋)
含0的数据无法进行对数变换,可用log(x+1)

2). GrLivArea
绘制直方图和正态概率曲线图:

#histogram and normal probability plot
sns.distplot(df_train['GrLivArea'], fit=norm);
fig = plt.figure()
res = stats.probplot(df_train['GrLivArea'], plot=plt)

进行对数变换:

#data transformation
df_train['GrLivArea'] = np.log(df_train['GrLivArea'])

#transformed histogram and normal probability plot
sns.distplot(df_train['GrLivArea'], fit=norm);
fig = plt.figure()
res = stats.probplot(df_train['GrLivArea'], plot=plt)

3). 'TotalBsmtSF'

绘制直方图和正态概率曲线图:

#histogram and normal probability plot
sns.distplot(df_train['TotalBsmtSF'], fit=norm);
fig = plt.figure()
res = stats.probplot(df_train['TotalBsmtSF'], plot=plt)

从图中可以看出:
 显示出了偏度
 大量为 0 的观察值(没有地下室的房屋)
 含 0 的数据无法进行对数变换

我们建立了一个变量,可以得到有没有地下室的影响值(二值变量),我们选择忽略零值,只对非零值进行对数变换。这样我们既可以变换数据,也不会损失有没有地下室的影响。

#create column for new variable (one is enough because it's a binary categorical feature)
#if area>0 it gets 1, for area==0 it gets 0
df_train['HasBsmt'] = pd.Series(len(df_train['TotalBsmtSF']), index=df_train.index)
df_train['HasBsmt'] = 0
df_train.loc[df_train['TotalBsmtSF']>0,'HasBsmt'] = 1

对该变量进行转换

#transform data
df_train.loc[df_train['HasBsmt']==1,'TotalBsmtSF'] = np.log(df_train['TotalBsmtSF'])

#histogram and normal probability plot
sns.distplot(df_train[df_train['TotalBsmtSF']>0]['TotalBsmtSF'], fit=norm);
fig = plt.figure()
res = stats.probplot(df_train[df_train['TotalBsmtSF']>0]['TotalBsmtSF'], plot=plt)

9、同方差性

最好的测量两个变量的同方差性的方法就是图像。
1). SalePrice 和 GrLivArea 同方差性
绘制散点图:

#scatter plot
plt.scatter(df_train['GrLivArea'], df_train['SalePrice']);

2).'SalePrice' with 'TotalBsmtSF'同方差性

绘制散点图:

#scatter plot
plt.scatter(df_train[df_train['TotalBsmtSF']>0]['TotalBsmtSF'], df_train[df_train['TotalBsmtSF']>0]['SalePrice']);


可以看出'SalePrice'在整个 'TotalBsmtSF'变量范围内显示出了同等级别的变化。
这就是正态化的作用,我们仅仅需要保证特征的正态化,就能解决同方差的问题。

10、虚拟变量

对类别变量进行虚拟编码:

df_train['SaleCondition'].unique()
array(['Normal', 'Abnorml', 'Partial', 'AdjLand', 'Alloca', 'Family'], dtype=object)

#convert categorical variable into dummy
df_train = pd.get_dummies(df_train)

df_train

1457 rows × 222 columns
[/cceN_python]

结论

整个方案中,我们使用了很多《多元数据分析》中提出的方法。我们对变量进行了哲学分析,不仅对'SalePrice'进行了单独分析,还结合了相关程度最高的变量进行分析。我们处理了缺失数据和异常值,我们验证了一些基础统计假设,并且将类别变量转换为虚拟变量。

但问题还没有结束,我们还需要预测房价的变化趋势,房价预测是否适合线性回归正则化的方法?是否适合组合方法?或者一些其他的方法?这篇文章解决了数据预处理的问题,预处理后,我们选择正确的模型算法进行拟合。

参考以下网站:
https://www.kaggle.com/pmarcelino/comprehensive-data-exploration-with-python
http://www.jianshu.com/p/defcbbfdb324

Scikit-learn简介

第4章 Scikit-learn简介

自2007年发布以来,Scikit-learn已经成为Python重要的机器学习库了。Scikit-Learn简称Sklearn,支持包括分类、回归、降维和聚类四大机器学习算法。还包含了特征提取、数据处理和模型评估三大模块。
Sklearn是Scipy的扩展,建立在NumPy和matplotlib库的基础上。利用这几大模块的优势,可以大大提高机器学习的效率。
Sklearn拥有着完善的文档,上手容易,具有着丰富的API,在学术界颇受欢迎。Sklearn已经封装了大量的机器学习算法,同时Sklearn内置了大量数据集,节省了获取和整理数据集的时间
1、 Sklearn安装
推荐使用Anaconda科学计算环境,里面已经内置了NumPy、SciPy、sklearn等模块,直接可用。或者使用conda进行包管理。

conda install scikit-learn

当然也可使用pip安装,在安装前,需要先安装如下依赖库:
Python(>=2.7 or >=3.5)
NumPy
SciPy
安装以后,我们可以参考安装的版本。

import sklearn
sklearn.__version__
'0.18.1'

2、 Sklearn基础
2.1 估计器(Estimator)
估计器,很多时候可以直接理解成分类器,主要包含两个函数:
fit():训练算法,设置内部参数。接收训练集和类别两个参数。
predict():预测测试集类别,参数为测试集。
大多数scikit-learn估计器接收和输出的数据格式均为numpy数组或类似格式。
2.2 转换器(Transformer)
转换器用于数据预处理和数据转换,主要是三个方法:
fit():训练算法,设置内部参数。
transform():数据转换。
fit_transform():合并fit和transform两个方法。
2.3 流水线(Pipeline)
流水线的功能:
跟踪记录各步骤的操作(以方便地重现实验结果),对各步骤进行一个封装,确保代码的复杂程度不至于超出掌控范围。
基本使用方法:
流水线的输入为一连串的数据挖掘步骤,其中最后一步必须是估计器,前几步是转换器。输入的数据集经过转换器的处理后,输出的结果作为下一步的输入。最后,用位于流水线最后一步的估计器对数据进行分类。
每一步都用元组( ‘名称’,步骤)来表示。现在来创建流水线。
使用示例:

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import MinMaxScaler
from sklearn.neighbors import KNeighborsClassifier

scaling_pipeline = Pipeline([('scale', MinMaxScaler()),('predict', KNeighborsClassifier())])

2.4 预处理
要在from sklearn.preprocessing包下。
规范化:
MinMaxScaler :最大最小值规范化
Normalizer :使每条数据各特征值的和为1
StandardScaler :为使各特征的均值为0,方差为1
编码:
LabelEncoder :把字符串类型的数据转化为整型
OneHotEncoder :特征用一个二进制数字来表示
Binarizer :为将数值型特征的二值化
MultiLabelBinarizer:多标签二值化
pandas.get_dummies: 对离散型特征进行one-hot编码,以下是对get_dummies使用示例:
离散特征的编码分为两种情况:
1、离散特征的取值之间没有大小的意义,比如color:[red,blue],那么就使用one-hot编码
2、离散特征的取值有大小的意义,比如size:[X,XL,XXL],那么就使用数值的映射{X:1,XL:2,XXL:3}
使用pandas可以很方便的对离散型特征进行one-hot编码
(1)生成一个矩阵

import pandas as pd
df = pd.DataFrame([
['green', 'M', 10.1, 'class1'],
['red', 'L', 13.5, 'class2'],
['blue', 'XL', 15.3, 'class1']])

df.columns = ['color', 'size', 'prize', 'class label']
df

(2)把size列用数值化,对应关系为:

size_mapping = { 'XL': 3,'L': 2, 'M': 1}
size_mapping = {
'XL': 3,
'L': 2,
'M': 1}
df['size'] = df['size'].map(size_mapping)
df


(3)对列class label数值化

class_mapping = {label:idx for idx,label in enumerate(set(df['class label']))}
df['class label'] = df['class label'].map(class_mapping)
df


(4)使用get_dummies进行one-hot编码

pd.get_dummies(df)


2.5 特征
sklearn常见的转换功能:


归一化使不同规格的数据转换到同一规格。无量刚好与标准化(Normalizer)的区别:
简单来说,无量纲化是依照特征矩阵的列处理数据,标准化是依照特征矩阵的行处理数据。标准化的前提是样本各特征值服从正态分布,标准化后将其转换成标准正态分布。标准化的公式类似于标准化,不同的是样本均值和样本标准差改为特征值均值和特征值标准差。
2.5.1 特征抽取
包:from sklearn.feature_extraction
特征抽取是机器学习任务最为重要的一个环节,一般而言,它对最终结果的影响要高过机器学习算法本身。只有先把现实用特征表示出来,才能借助机器学习的力量找到问题的答案。特征选择的另一个优点在于:降低真实世界的复杂度,模型比现实更容易操纵。
一般最常使用的特征抽取技术都是高度针对具体领域的,对于特定的领域,如图像处理,在过去一段时间已经开发了各种特征抽取的技术,但这些技术在其他领域的应用却非常有限。
DictVectorizer: 将dict类型的list数据,转换成numpy array
FeatureHasher : 特征哈希,相当于一种降维技巧
image:图像相关的特征抽取
text: 文本相关的特征抽取
text.CountVectorizer:将文本转换为每个词出现的个数的向量
text.TfidfVectorizer:将文本转换为tfidf值的向量
text.HashingVectorizer:文本的特征哈希
示例

CountVectorize只数出现个数:

TfidfVectorizer:个数+归一化(不包括idf:

HashingVectorizer:

2.5.2 特征选择
包:sklearn.feature_selection
特征选择的原因如下:
(1)降低复杂度
(2)降低噪音
(3)增加模型可读性
VarianceThreshold: 删除特征值的方差达不到最低标准的特征
SelectKBest: 返回k个最佳特征
SelectPercentile: 返回表现最佳的前r%个特征
单个特征和某一类别之间相关性的计算方法有很多。最常用的有卡方检验(χ2)。其他方法还有互信息和信息熵。
2.6 降维
包:sklearn.decomposition
主成分分析算法(Principal Component Analysis, PCA)的目的是找到能用较少信息描述数据集的特征组合。它意在发现彼此之间没有相关性、能够描述数据集的特征,确切说这些特征的方差跟整体方差没有多大差距,这样的特征也被称为主成分。这也就意味着,借助这种方法,就能通过更少的特征捕获到数据集的大部分信息。
2.7 组合
包:sklearn.ensemble
组合技术即通过聚集多个分类器的预测来提高分类准确率。
常用的组合分类器方法:
(1)通过处理训练数据集。即通过某种抽样分布,对原始数据进行再抽样,得到多个训练集。常用的方法有装袋(bagging)和提升(boosting)。
(2)通过处理输入特征。即通过选择输入特征的子集形成每个训练集。适用于有大量冗余特征的数据集。随机森林(Random forest)就是一种处理输入特征的组合方法。
(3)通过处理类标号。适用于多分类的情况,将类标号随机划分成两个不相交的子集,再把问题变为二分类问题,重复构建多次模型,进行分类投票。

BaggingClassifier: Bagging分类器组合
BaggingRegressor: Bagging回归器组合
AdaBoostClassifier: AdaBoost分类器组合
AdaBoostRegressor: AdaBoost回归器组合
GradientBoostingClassifier:GradientBoosting分类器组合
GradientBoostingRegressor: GradientBoosting回归器组合
ExtraTreeClassifier:ExtraTree分类器组合
ExtraTreeRegressor: ExtraTree回归器组合
RandomTreeClassifier:随机森林分类器组合
RandomTreeRegressor: 随机森林回归器组合
使用举例

from sklearn.ensemble import AdaBoostClassifier
from sklearn.tree import DecisionTreeClassifier
AdaBoostClassifier(DecisionTreeClassifier(max_depth=1),
algorithm="SAMME",
n_estimators=200)

解释
装袋(bagging):根据均匀概率分布从数据集中重复抽样(有放回),每个自助样本集和原数据集一样大,每个自助样本集含有原数据集大约60%的数据。训练k个分类器,测试样本被指派到得票最高的类。
提升(boosting):通过给样本设置不同的权值,每轮迭代调整权值。不同提升算法之间的差别,一般是(1)如何更新样本的权值,(2)如何组合每个分类器的预测。其中Adaboost中,样本权值是增加那些被错误分类的样本的权值,分类器C_i的重要性依赖于它的错误率。
Boosting主要关注降低偏差,因此Boosting能基于泛化性能相当弱的学习器构建出很强的集成;Bagging主要关注降低方差,因此它在不剪枝的决策树、神经网络等学习器上效用更为明显。偏差指的是算法的期望预测与真实预测之间的偏差程度,反应了模型本身的拟合能力;方差度量了同等大小的训练集的变动导致学习性能的变化,刻画了数据扰动所导致的影响。

2.8 模型评估(度量)
包:sklearn.metrics
sklearn.metrics包含评分方法、性能度量、成对度量和距离计算。
分类结果度量
参数大多是y_true和y_pred。
accuracy_score:分类准确度
condusion_matrix :分类混淆矩阵
classification_report:分类报告
precision_recall_fscore_support:计算精确度、召回率、f、支持率
2.9 交叉验证
包:sklearn.cross_validation

KFold:K-Fold交叉验证迭代器。接收元素个数、fold数、是否清洗
LeaveOneOut:LeaveOneOut交叉验证迭代器
LeavePOut:LeavePOut交叉验证迭代器
LeaveOneLableOut:LeaveOneLableOut交叉验证迭代器
LeavePLabelOut:LeavePLabelOut交叉验证迭代器
LeaveOneOut(n) 相当于 KFold(n, n_folds=n) 相当于LeavePOut(n, p=1)。
LeaveP和LeaveOne差别在于leave的个数,也就是测试集的尺寸。LeavePLabel和LeaveOneLabel差别在于leave的Label的种类的个数。
LeavePLabel这种设计是针对可能存在第三方的Label,比如我们的数据是一些季度的数据。那么很自然的一个想法就是把1,2,3个季度的数据当做训练集,第4个季度的数据当做测试集。这个时候只要输入每个样本对应的季度Label,就可以实现这样的功能。
以下是实验代码:

import numpy as np
import sklearn
import sklearn.cross_validation as cross_validation
X = np.array([[1, 2], [3, 4], [5, 6], [7, 8],[9, 10]])
y = np.array([1, 2, 1, 2, 3])
def show_cross_val(method):
if method == "lolo":
labels = np.array(["summer", "winter", "summer", "winter", "spring"])
cv = cross_validation.LeaveOneLabelOut(labels)
elif method == 'lplo':
labels = np.array(["summer", "winter", "summer", "winter", "spring"])
cv = cross_validation.LeavePLabelOut(labels,p=2)
elif method == 'loo':
cv = cross_validation.LeaveOneOut(n=len(y))
elif method == 'lpo':
cv = cross_validation.LeavePOut(n=len(y),p=3)
for train_index, test_index in cv:
print("TRAIN:", train_index, "TEST:", test_index)
X_train, X_test = X[train_index], X[test_index]
y_train, y_test = y[train_index], y[test_index]
print ("X_train: ",X_train )
print ("y_train: ", y_train )
print ("X_test: ",X_test )
print ("y_test: ",y_test )
if __name__ == '__main__':
show_cross_val("lpo")

常用方法
• train_test_split:分离训练集和测试集(不是K-Fold)
• cross_val_score:交叉验证评分,可以指认cv为上面的类的实例
• cross_val_predict:交叉验证的预测
2.10 网格搜索
包:sklearn.grid_search
网格搜索最佳参数

GridSearchCV:搜索指定参数网格中的最佳参数
ParameterGrid:参数网格
ParameterSampler:用给定分布生成参数的生成器
RandomizedSearchCV:超参的随机搜索
通过best_estimator_.get_params()方法,获取最佳参数。
2.11 多分类、多标签分类
包:sklearn.multiclass
OneVsRestClassifier:1-rest多分类(多标签)策略
OneVsOneClassifier:1-1多分类策略
OutputCodeClassifier:1个类用一个二进制码表示
示例代码

from sklearn import metrics
from sklearn import cross_validation
from sklearn.svm import SVC
from sklearn.multiclass import OneVsRestClassifier
from sklearn.preprocessing import MultiLabelBinarizer
import numpy as np
from numpy import random
X=np.arange(15).reshape(5,3)
y=np.arange(5)
Y_1 = np.arange(5)
random.shuffle(Y_1)
Y_2 = np.arange(5)
random.shuffle(Y_2)
Y = np.c_[Y_1,Y_2]
def multiclassSVM():
X_train, X_test, y_train, y_test = cross_validation.train_test_split(X, y, test_size=0.2,random_state=0)
model = OneVsRestClassifier(SVC())
model.fit(X_train, y_train)
predicted = model.predict(X_test)
print (predicted)
def multilabelSVM():
Y_enc = MultiLabelBinarizer().fit_transform(Y)
X_train, X_test, Y_train, Y_test = cross_validation.train_test_split(X, Y_enc, test_size=0.2, random_state=0)
model = OneVsRestClassifier(SVC())
model.fit(X_train, Y_train)
predicted = model.predict(X_test)
print (predicted)
if __name__ == '__main__':
multiclassSVM()

上面的代码测试了svm在OneVsRestClassifier的包装下,分别处理多分类和多标签的情况。特别注意,在多标签的情况下,输入必须是二值化的。所以需要MultiLabelBinarizer()先处理。
3、 创建自己的转换器
在特征抽取的时候,经常会发现自己的一些数据预处理的方法,sklearn里可能没有实现,但若直接在数据上改,又容易将代码弄得混乱,难以重现实验。这个时候最好自己创建一个转换器,在后面将这个转换器放到pipeline里,统一管理。
如果我们想接收一个numpy数组,根据其均值将其离散化,任何高于均值的特征值替换为1,小于或等于均值的替换为0。
代码实现:

from sklearn.base import TransformerMixin
from sklearn.utils import as_float_array

class MeanDiscrete(TransformerMixin):

#计算出数据集的均值,用内部变量保存该值。
def fit(self, X, y=None):
X = as_float_array(X)
self.mean = np.mean(X, axis=0)
#返回self,确保在转换器中能够进行链式调用(例如调用transformer.fit(X).transform(X))
return self

def transform(self, X):
X = as_float_array(X)
assert X.shape[1] == self.mean.shape[0]
return X > self.mean

参考博客:
http://www.jianshu.com/p/516f009c0875