（干货）结合Scikit-learn介绍几种常用的特征选择方法

系统版本：window 7 （64bit）
python版本：python 3.5
我的GitHub：https://github.com/weepon
写在前面：前段时间正好用到特征选择的知识，有幸读到这篇文章，本文也主要参考这篇文章写成，但与原文章有不同之处：第一、纠正了原始文章中的代码错误，使其能用python3.5正常运行；第二、增加了一些新的特征选择内容，使其更加完善。本文中所有代码均放在Github

特征选择(排序)对于数据科学家、机器学习从业者来说非常重要。好的特征选择能够提升模型的性能，更能帮助我们理解数据的特点、底层结构，这对进一步改善模型、算法都有着重要作用。

特征选择主要有两个功能：

• 减少特征数量、降维，使模型泛化能力更强，减少过拟合
• 增强对特征和特征值之间的理解

拿到数据集，一个特征选择方法，往往很难同时完成这两个目的。通常情况下，我们经常不管三七二十一，选择一种自己最熟悉或者最方便的特征选择方法（往往目的是降维，而忽略了对特征和数据理解的目的）。

在许多机器学习相关的书里，很难找到关于特征选择的内容，因为特征选择要解决的问题往往被视为机器学习的一种副作用，一般不会单独拿出来讨论。

本文将结合Scikit-learn提供的例子介绍几种常用的特征选择方法，它们各自的优缺点和问题。源代码

1 去掉取值变化小的特征 (Removing features with low variance)

这应该是最简单的特征选择方法了：假设某特征的特征值只有0和1，并且在所有输入样本中，95%的实例的该特征取值都是1，那就可以认为这个特征作用不大。如果100%都是1，那这个特征就没意义了。当特征值都是离散型变量的时候这种方法才能用，如果是连续型变量，就需要将连续变量离散化之后才能用，而且实际当中，一般不太会有95%以上都取某个值的特征存在，所以这种方法虽然简单但是不太好用。可以把它作为特征选择的预处理，先去掉那些取值变化小的特征，然后再从接下来提到的的特征选择方法中选择合适的进行进一步的特征选择。

scikit-learn中的VarianceThreshold 是特征选择中的一项基本方法。它会移除所有方差不满足阈值的特征。默认设置下，它将移除所有方差为0的特征，即那些在所有样本中数值完全相同的特征。假设我们有一个带有布尔特征的数据集，我们要移除那些超过80%的数据都为1或0的特征。布尔特征是伯努利随机变量，该类变量的方差为：

$VAR[X]=p(1-p)$
我们可以使用阈值 0.8 * (1 - 0.8):

from sklearn.feature_selection import VarianceThreshold
X = [[0, 0, 1], [0, 1, 0], [1, 0, 0], [0, 1, 1], [0, 1, 0], [0, 1, 1]]
sel = VarianceThreshold(threshold=(.8 * (1 - .8)))
sel.fit_transform(X)

array([[0, 1],
   [1, 0],
   [0, 0],
   [1, 1],
   [1, 0],
   [1, 1]])

果然, VarianceThreshold 移除了第一列特征，因为第一列中特征值为0的概率达到了 p = 5/6 > 0.8 。

2 单变量特征选择 (Univariate feature selection)

单变量特征选择能够对每一个特征进行测试，衡量该特征和响应变量之间的关系，根据得分扔掉不好的特征。对于回归和分类问题可以采用卡方检验等方式对特征进行测试。

这种方法比较简单，易于运行，易于理解，通常对于理解数据有较好的效果（但对特征优化、提高泛化能力来说不一定有效）；这种方法有许多改进的版本、变种。

2.1 Pearson相关系数 Pearson Correlation

皮尔森相关系数是一种最简单的，能帮助理解特征和响应变量之间关系的方法，该方法衡量的是变量之间的线性相关性，结果的取值区间为[-1，1]，-1表示完全的负相关(这个变量下降，那个就会上升)，+1表示完全的正相关，0表示没有线性相关。

Pearson Correlation速度快、易于计算，经常在拿到数据(经过清洗和特征提取之后的)之后第一时间就执行。Scipy的pearsonr方法能够同时计算相关系数和p-value

import numpy as np  
from scipy.stats import pearsonr  
np.random.seed(0)  
size = 300  
x = np.random.normal(0, 1, size)  
print( "Lower noise", pearsonr(x, x + np.random.normal(0, 1, size)))  
print( "Higher noise", pearsonr(x, x + np.random.normal(0, 10, size))) 

Lower noise (0.71824836862138408, 7.3240173129983507e-49)
Higher noise (0.057964292079338155, 0.31700993885324752)

这个例子中，我们比较了变量在加入噪音之前和之后的差异。当噪音比较小的时候，相关性很强，p-value很低。
Scikit-learn提供的f_regrssion方法能够批量计算特征的p-value，非常方便，参考sklearn的pipeline

Pearson相关系数的一个明显缺陷是，作为特征排序机制，他只对线性关系敏感。如果关系是非线性的，即便两个变量具有一一对应的关系，Pearson相关性也可能会接近0。

x1 = np.random.uniform(-1, 1, 100000)  
print( pearsonr(x1, x1**2)[0])

-0.00218133842327

更多类似的例子参考sample plots。另外，如果仅仅根据相关系数这个值来判断的话，有时候会具有很强的误导性，如Anscombe’s quartet，最好把数据可视化出来，以免得出错误的结论。

2.2 互信息和最大信息系数 (Mutual information and maximal information coefficient (MIC) )

互信息公式
想把互信息直接用于特征选择其实不是太方便：1、它不属于度量方式，也没有办法归一化，在不同数据及上的结果无法做比较；2、对于连续变量的计算不是很方便（X和Y都是集合，x，y都是离散的取值），通常变量需要先离散化，而互信息的结果对离散化的方式很敏感。

最大信息系数克服了这两个问题。它首先寻找一种最优的离散化方式，然后把互信息取值转换成一种度量方式，取值区间在[0，1]。minepy提供了MIC功能。

反过头来看y=x^2这个例子，MIC算出来的互信息值为1(最大的取值)。(建议用pip install minepy, 之前试过其他方式安装都失败了。。。)

from minepy import MINE  
m = MINE()  
x = np.random.uniform(-1, 1, 10000)  
m.compute_score(x, x**2)  
print( m.mic())  

1.0000

MIC的统计能力遭到了一些质疑，当零假设不成立时，MIC的统计就会受到影响。在有的数据集上不存在这个问题，但有的数据集上就存在这个问题。

2.3 距离相关系数 (Distance correlation)

距离相关系数是为了克服Pearson相关系数的弱点而生的。在x和x^2这个例子中，即便Pearson相关系数是0，我们也不能断定这两个变量是独立的（有可能是非线性相关）；但如果距离相关系数是0，那么我们就可以说这两个变量是独立的。

R的energy包里提供了距离相关系数的实现，另外这是Python gist的实现。

#R-code
> x = runif (1000, -1, 1)
> dcor(x, x**2)
[1] 0.4943864 

尽管有MIC和距离相关系数在了，但当变量之间的关系接近线性相关的时候，Pearson相关系数仍然是不可替代的。第一、Pearson相关系数计算速度快，这在处理大规模数据的时候很重要。第二、Pearson相关系数的取值区间是[-1，1]，而MIC和距离相关系数都是[0，1]。这个特点使得Pearson相关系数能够表征更丰富的关系，符号表示关系的正负，绝对值能够表示强度。当然，Pearson相关性有效的前提是两个变量的变化关系是单调的。

2.4 基于学习模型的特征排序 (Model based ranking)

这种方法的思路是直接使用你要用的机器学习算法，针对每个单独的特征和响应变量建立预测模型。其实Pearson相关系数等价于线性回归里的标准化回归系数。假如某个特征和响应变量之间的关系是非线性的，可以用基于树的方法（决策树、随机森林）、或者扩展的线性模型等。基于树的方法比较易于使用，因为他们对非线性关系的建模比较好，并且不需要太多的调试。但要注意过拟合问题，因此树的深度最好不要太大，再就是运用交叉验证。

在波士顿房价数据集上使用sklearn的随机森林回归给出一个单变量选择的例子：

from sklearn.cross_validation import cross_val_score, ShuffleSplit  
from sklearn.datasets import load_boston  
from sklearn.ensemble import RandomForestRegressor  

#Load boston housing dataset as an example  
boston = load_boston()  
X = boston["data"]  
Y = boston["target"]  
names = boston["feature_names"]  

rf = RandomForestRegressor(n_estimators=20, max_depth=4)  
scores = []  
for i in