特征工程

特征筛选的方法主要包括：Filter（过滤法）、Wrapper（封装法）、Embedded（嵌入法）

filter: 过滤法

特征选择方法一：去掉取值变化小的特征（Removing features with low variance）

方法虽然简单但是不太好用，可以把它作为特征选择的预处理，先去掉那些取值变化小的特征
如果机器资源充足，并且希望尽量保留所有信息，可以把阈值设置得比较高，或者只过滤离散型特征只有一个取值的特征。
离散型变量：95%的实例的该特征取值都是1，那就可以认为这个特征作用不大。 如果100%都是1，那这个特征就没意义了。
连续型变量：抛弃那些方差小于某个阈值的特征。

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)

特征选择实现方法二：单变量特征选择

单变量特征选择方法独立的衡量每个特征与响应变量之间的关系，单变量特征选择能够对每一个特征进行测试，
衡量该特征和响应变量之间的关系，根据得分扔掉不好的特征。该方法简单，易于运行，易于理解，
通常对于理解数据有较好的效果（但对特征优化、提高泛化能力来说不一定有效）
单变量特征选择可以用于理解数据、数据的结构、特点，也可以用于排除不相关特征，但是它不能发现冗余特征。

　　1.Pearson相关系数（Pearson Correlation） -- 主要用于连续型特征的筛选，只对线性关系敏感

import numpy as np
from scipy.stats import pearsonr

np.random.seed(2019)
size=1000
x = np.random.normal(0, 1, size)
# 计算两变量间的相关系数
print("Lower noise {}".format(pearsonr(x, x + np.random.normal(0, 1, size))))
print("Higher noise {}".format(pearsonr(x, x + np.random.normal(0, 10, size))))

　　2.互信息和最大信息系数（Mutual information and maximal information coefficient）

熵H(Y)与条件熵H(Y|X)之间的差称为互信息 -- 只能用于离散型特征的选择，对离散化的方式很敏感
由于互信息法并不方便直接用于特征选择，因此引入了最大信息系数。最大信息数据首先寻找一种最优的离散方式，然后把互信息取值转换成一种度量方式，取值区间为[0,1]。

x = np.random.normal(0,10,300)
z = x *x
pearsonr(x,z)
# 输出-0.1
from minepy import MINE
m = MINE()
m.compute_score(x, z)
print(m.mic())

　　3.距离相关系数（Distance correlation） -- 为了克服Pearson相关系数的弱点而生的。

from scipy.spatial.distance import pdist, squareform
import numpy as np

from numbapro import jit, float32

def distcorr(X, Y):
    """ Compute the distance correlation function

    >>> a = [1,2,3,4,5]
    >>> b = np.array([1,2,9,4,4])
    >>> distcorr(a, b)
    0.762676242417
    """
    X = np.atleast_1d(X)
    Y = np.atleast_1d(Y)
    if np.prod(X.shape) == len(X):
        X = X[:, None]
    if np.prod(Y.shape) == len(Y):
        Y = Y[:, None]
    X = np.atleast_2d(X)
    Y = np.atleast_2d(Y)
    n = X.shape[0]
    if Y.shape[0] != X.shape[0]:
        raise ValueError('Number of samples must match')
    a = squareform(pdist(X))
    b = squareform(pdist(Y))
    A = a - a.mean(axis=0)[None, :] - a.mean(axis=1)[:, None] + a.mean()
    B = b - b.mean(axis=0)[None, :] - b.mean(axis=1)[:, None] + b.mean()

    dcov2_xy = (A * B).sum()/float(n * n)
    dcov2_xx = (A * A).sum()/float(n * n)
    dcov2_yy = (B * B).sum()/float(n * n)
    dcor = np.sqrt(dcov2_xy)/np.sqrt(np.sqrt(dcov2_xx) * np.sqrt(dcov2_yy))
    return dcor

　　距离相关系数的另一种算法，链接

import numpy as np

def dist(x, y):
    #1d only
    return np.abs(x[:, None] - y)
    
def d_n(x):
    d = dist(x, x)
    dn = d - d.mean(0) - d.mean(1)[:,None] + d.mean()
    return dn

def dcov_all(x, y):
    dnx = d_n(x)
    dny = d_n(y)
    
    denom = np.product(dnx.shape)
    dc = (dnx * dny).sum() / denom
    dvx = (dnx**2).sum() / denom
    dvy = (dny**2).sum() / denom
    dr = dc / (np.sqrt(dvx) * np.sqrt(dvy))
    return dc, dr, dvx, dvy


import matplotlib.pyplot as plt
fig = plt.figure()
for case in range(1,5):
    np.random.seed(9854673)
    x = np.linspace(-1,1, 501)
    if case == 1:
        y = - x**2 + 0.2 * np.random.rand(len(x))
    elif case == 2:
        y = np.cos(x*2*np.pi) + 0.1 * np.random.rand(len(x))
    elif case == 3:
        x = np.sin(x*2*np.pi) + 0.0 * np.random.rand(len(x))  #circle
    elif case == 4:
        x = np.sin(x*1.5*np.pi) + 0.1 * np.random.rand(len(x))  #bretzel
    dc, dr, dvx, dvy = dcov_all(x, y)
    print dc, dr, dvx, dvy
    
    ax = fig.add_subplot(2,2, case)
    #ax.set_xlim(-1, 1)
    ax.plot(x, y, '.')
    yl = ax.get_ylim()
    ax.text(-0.95, yl[0] + 0.9 * np.diff(yl), 'dr=%4.2f' % dr)

plt.show()

　　4.基于学习模型的特征排序（Model based ranking）

这种方法的思路是直接使用你要用的机器学习算法，针对每个单独的特征和响应变量建立预测模型。如果特征与响应变量之间的关系是非线性的，则有许多替代方案，例如基于树的方法（决策树，随机森林）、或者扩展的线性模型等。基于树的方法是最简单的方法之一，因为他们可以很好地模拟非线性关系，不需要太多的调整。
但是要避免的主要是过度拟合，因此树的深度应该相对较小，并且应该应用交叉验证。

from sklearn.datasets import load_boston
from sklearn.model_selection  import train_test_split, cross_val_score, ShuffleSplit
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestRegressor
#from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error

#Load boston housing dataset as an example
boston = load_boston()
#print(boston.DESCR)
#x_train,x_test,y_train,y_test = train_test_split(boston.data,boston.target,random_state=33,test_size=0.25)
X = boston["data"]
Y = boston["target"]
names = boston["feature_names"]

rf = RandomForestRegressor(n_estimators=20, max_depth=4)

scores = []
# 使用每个特征单独训练模型，并获取每个模型的评分来作为特征选择的依据。
for i in range(X.shape[1]):
     score = cross_val_score(rf, X[:, i:i+1], Y, scoring="r2",
                              cv=ShuffleSplit(len(X), 3, .3))
     scores.append((round(np.mean(score), 3), names[i]))
print(sorted(scores, reverse=True))


from sklearn.svm import SVR
l_svr = SVR(kernel='linear')  #poly, rbf
l_svr.fit(x_train,y_train)
l_svr.score(x_test,y_test)

from sklearn.neighbors import KNeighborsRegressor
knn = KNeighborsRegressor(weights="uniform")
knn.fit(x_train,y_train)
knn.score(x_test,y_test)

from sklearn.tree import DecisionTreeRegressor
dt = DecisionTreeRegressor()
dt.fit(x_train,y_train)
dt.score(x_test,y_test)

from sklearn.ensemble import RandomForestRegressor
rfr = RandomForestRegressor()
rfr.fit(x_train,y_train)
rfr.score(x_test,y_test)

from sklearn.ensemble import ExtraTreesRegressor
etr = ExtraTreesRegressor()
etr.fit(x_train,y_train)
etr.score(x_test,y_test)

from sklearn.ensemble import GradientBoostingRegressor
gbr = GradientBoostingRegressor()
gbr.fit(x_train,y_train)
gbr.score(x_test,y_test)

　　5.卡方检验 -- 只适用于分类问题中离散型特征筛选

卡方值描述两个事件的独立性或者描述实际观察值与期望值的偏离程度。卡方值越大，表名实际观察值与期望值偏离越大，也说明两个事件的相互独立性越弱。

#导入sklearn库中的SelectKBest和chi2
from sklearn.feature_selection import SelectKBest ,chi2
#选择相关性最高的前5个特征
X_chi2 = SelectKBest(chi2, k=5).fit_transform(X, y)
X_chi2.shape

　　

特征选择实现方法三：线性模型与正则化

当所有特征在相同尺度上时，最重要的特征应该在模型中具有最高系数，而与输出变量不相关的特征应该具有接近零的系数值。即使使用简单的线性回归模型，当数据不是很嘈杂（或者有大量数据与特征数量相比）并且特征（相对）独立时，这种方法也能很好地工作。
正则化就是把额外的约束或者惩罚项加到已有模型（损失函数）上，以防止过拟合并提高泛化能力。
当用在线性模型上时，L1正则化和L2正则化也称为Lasso和Ridge。
Lasso能够挑出一些优质特征，同时让其他特征的系数趋于0。当如需要减少特征数的时候它很有用，但是对于数据理解来说不是很好用。
Ridge将回归系数均匀的分摊到各个关联变量上，L2正则化对于特征理解来说更加有用
多元线性回归方程演变成求θ。每个特征都有对应的权重系数coef，特征的权重系数的正负值代表特征与目标值是正相关还是负相关，特征的权重系数的绝对值代表重要性。

#获取boston数据
from sklearn.linear_model import LinearRegression
boston=load_boston()
x=boston.data
y=boston.target
#过滤掉异常值
x=x[y<50]
y=y[y<50]
reg=LinearRegression()
reg.fit(x,y)
#求排序后的coef
coefSort=reg.coef_.argsort()
#featureNameSort: 按对标记值的影响，从小到大的各特征值名称
#featureCoefSore：按对标记值的影响，从小到大的coef_
featureNameSort=boston.feature_names[coefSort]
featureCoefSore=reg.coef_[coefSort]
print("featureNameSort:", featureNameSort)
print("featureCoefSore:", featureCoefSore)

#A helper method for pretty-printing linear models
def pretty_print_linear(coefs, names = None, sort = False):
    if names == None:
        names = ["X%s" % x for x in range(len(coefs))]
    lst = zip(coefs, names)
    if sort:
        lst = sorted(lst,  key = lambda x:-np.abs(x[0]))
    return " + ".join("%s * %s" % (round(coef, 3), name)
                                   for coef, name in lst)

# lasso回归
from sklearn.linear_model import Lasso
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_boston

boston = load_boston()
scaler = StandardScaler()
X = scaler.fit_transform(boston["data"])
Y = boston["target"]
names = boston["feature_names"]

lasso = Lasso(alpha=.3)
lasso.fit(X, Y)

print("Lasso model: {}".format(
      pretty_print_linear(lasso.coef_, names, sort = True)))

# 岭回归
from sklearn.linear_model import Ridge
from sklearn.metrics import r2_score
size = 100

#We run the method 10 times with different random seeds
for i in range(10):
    print("Random seed {}".format(i))
    np.random.seed(seed=i)
    X_seed = np.random.normal(0, 1, size)
    X1 = X_seed + np.random.normal(0, .1, size)
    X2 = X_seed + np.random.normal(0, .1, size)
    X3 = X_seed + np.random.normal(0, .1, size)
    Y = X1 + X2 + X3 + np.random.normal(0, 1, size)
    X = np.array([X1, X2, X3]).T


    lr = LinearRegression()
    lr.fit(X,Y)
    print("Linear model: {}".format(pretty_print_linear(lr.coef_)))

    ridge = Ridge(alpha=10)
    ridge.fit(X,Y)
    print("Ridge model: {}".format(pretty_print_linear(ridge.coef_)))

　　

特征选择实现方法四：随机森林选择

　　1.平均不纯度减少（mean decrease impurity）
当训练决策树的时候，可以计算出每个特征减少了多少树的不纯度。对于一个决策树森林来说，可以算出每个特征平均减少了多少不纯度，并把它平均减少的不纯度作为特征选择的标准。

from sklearn.datasets import load_boston
from sklearn.ensemble import RandomForestRegressor
import numpy as np

#Load boston housing dataset as an example
boston = load_boston()
X = boston["data"]
Y = boston["target"]
names = boston["feature_names"]
# 训练随机森林模型，并通过feature_importances_属性获取每个特征的重要性分数。rf = RandomForestRegressor()
rf.fit(X, Y)
print("Features sorted by their score:")
print(sorted(zip(map(lambda x: round(x, 4), rf.feature_importances_), names),
             reverse=True))

　　2.平均精确度减少（mean decrease accuracy）

通过直接度量每个特征对模型精确率的影响来进行特征选择。
主要思路是打乱每个特征的特征值顺序，并且度量顺序变动对模型的精确率的影响。
对于不重要的变量来说，打乱顺序对模型的精确率影响不会太大。
对于重要的变量来说，打乱顺序就会降低模型的精确率。

from sklearn.model_selection import ShuffleSplit
from sklearn.metrics import r2_score
from collections import defaultdict
X = boston["data"]
Y = boston["target"]
rf = RandomForestRegressor()
scores = defaultdict(list)
#crossvalidate the scores on a number of different random splits of the data
for train_idx, test_idx in ShuffleSplit(len(X), 100, .3):
    X_train, X_test = X[train_idx], X[test_idx]
    Y_train, Y_test = Y[train_idx], Y[test_idx]
    # 使用修改前的原始特征训练模型，其acc作为后续混洗特征值后的对比标准。
    r = rf.fit(X_train, Y_train)
    acc = r2_score(Y_test, rf.predict(X_test))
    # 遍历每一列特征
    for i in range(X.shape[1]):
        X_t = X_test.copy()
        # 对这一列特征进行混洗，交互了一列特征内部的值的顺序
        np.random.shuffle(X_t[:, i])
        shuff_acc = r2_score(Y_test, rf.predict(X_t))
        # 混洗某个特征值后，计算平均精确度减少程度。scores[names[i]].append((acc-shuff_acc)/acc)
print("Features sorted by their score:")
print(sorted([(round(np.mean(score), 4), feat) for feat, score in scores.items()], reverse=True))

　　

特征选择实现方法五：顶层特征选择

　　1.稳定性选择（Stability selection）
它的主要思想是在不同的数据子集和特征子集上运行特征选择算法，不断的重复，最终汇总特征选择结果。比如可以统计某个特征被认为是重要特征的频率（被选为重要特征的次数除以它所在的子集被测试的次数）。
理想情况下，重要特征的得分会接近100%。稍微弱一点的特征得分会是非0的数，而最无用的特征得分将会接近于0。

from sklearn.linear_model import RandomizedLasso
from sklearn.datasets import load_boston
boston = load_boston()
#using the Boston housing data.
#Data gets scaled automatically by sklearn's implementation
X = boston["data"]
Y = boston["target"]
names = boston["feature_names"]
rlasso = RandomizedLasso(alpha=0.025)
rlasso.fit(X, Y)
print("Features sorted by their score:")
print(sorted(zip(map(lambda x: round(x, 4), rlasso.scores_), names),
             reverse=True))

　　2.递归特征消除（Recursive feature elimination，RFE）

递归特征消除的主要思想是反复的构建模型（如SVM或者回归模型）然后选出最好的（或者最差的）的特征（可以根据系数来选），把选出来的特征放到一遍，然后在剩余的特征上重复这个过程，直到所有特征都遍历了。
这个过程中特征被消除的次序就是特征的排序。因此，这是一种寻找最优特征子集的贪心算法。
RFE的稳定性很大程度上取决于在迭代的时候底层用哪种模型。
假如RFE采用的普通的回归，没有经过正则化的回归是不稳定的，那么RFE就是不稳定的。
假如RFE采用的是Ridge，而用Ridge正则化的回归是稳定的，那么RFE就是稳定的。

from sklearn.feature_selection import RFE
from sklearn.linear_model import LinearRegression
boston = load_boston()
X = boston["data"]
Y = boston["target"]
names = boston["feature_names"]
#use linear regression as the model
lr = LinearRegression()
#rank all features, i.e continue the elimination until the last one
rfe = RFE(lr, n_features_to_select=1)
rfe.fit(X,Y)
print("Features sorted by their rank:")
print(sorted(zip(map(lambda x: round(x, 4), rfe.ranking_), names)))

　　正则化的线性模型可用于特征理解和特征选择。相比起L1正则化，L2正则化的表现更加稳定，L2正则化对于数据的理解来说很合适。

　　由于响应变量和特征之间往往是非线性关系，可以采用basis expansion的方式将特征转换到一个更加合适的空间当中，在此基础上再考虑运用简单的线性模型。

　　随机森林是一种非常流行的特征选择方法，它易于使用。但它有两个主要问题：

• 重要的特征有可能得分很低（关联特征问题）
• 这种方法对特征变量类别多的特征越有利（偏向问题）

　　当选择最优特征以提升模型性能的时候，可以采用交叉验证的方法来验证某种方法是否比其他方法要好。
　　当用特征选择的方法来理解数据的时候要留心，特征选择模型的稳定性非常重要，稳定性差的模型很容易就会导致错误的结论。

 参考资料：

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

特征工程系列：特征筛选的原理与实现（上）

特征工程系列：特征筛选的原理与实现（下）

特征选择/筛选方法总结