1. 特征选择------sklearn代码

1.1 特征选择------方差法

忽略warning错误

import warnings

warnings.filterwarnings("ignore")

# 方差法

from sklearn.feature_selection import VarianceThreshold

X = [[0, 2, 0, 3], [0, 1, 4, 3], [0, 1, 1, 3]]

selector = VarianceThreshold()

X_new = selector.fit_transform(X)

#

print('X_new:\n',selector.fit_transform(X))

print('get_params:\n',selector.get_params())

print('get_support:\n',selector.get_support())

print('inverse_transform:\n',selector.inverse_transform(X_new))

运行结果

X_new:

[[2 0]

[1 4]

[1 1]]

get_params:

{'threshold': 0.0}

get_support:

[False True True False]

inverse_transform:

[[0 2 0 0]

[0 1 4 0]

[0 1 1 0]]

1.2 特征选择------单变量特征选择 (卡方,F分布,互信息)

1.2.1 什么是卡方分布

Compute chi-squared stats between each non-negative feature and class.This score can be used to select the n_features features with the highest values for the test chi-squared statistic from X, which must contain only non-negative features.

检测每个自变量与因变量的相关性，使用此函数“除去”最可能独立于类的特征。要优于方差检验

1.2.2 时间复杂度

O(n_classes * n_features)

1.2.3 sklearn中卡方的两个函数

* SelectKBest(score_func，k=)。保留评分最高得分的 K 个特征;

* SelectPercentile(score_func，percentile=)。保留最高得分的百分比特征;

1.2.4 score_func的选择

1.2.3 的两个函数 默认：f_classif 函数。这里想使用卡方分布法，选择函数chi2.

对于回归:

f_regression：相关系数，计算每个变量与目标变量的相关系数，然后计算出F值和P值；

mutual_info_regression：互信息;

对于分类:

chi2：卡方检验；

f_classif：方差分析，计算方差分析(ANOVA)的F值 (组间均方 / 组内均方)；

mutual_info_classif：互信息；

注：chi2 , mutual_info_classif , mutual_info_regression 可以保持数据的稀疏性。

注： 关于f_classif和f_regression的了解。点此链接

1.2.5 两个属性

scores_ : array-like, shape=(n_features,)，Scores of features.

pvalues_ : array-like, shape=(n_features,)，p-values of feature scores, None if score_func returned only scores.

from sklearn.datasets import load_iris

from sklearn.feature_selection import SelectKBest,SelectPercentile

from sklearn.feature_selection import chi2,f_classif,mutual_info_classif

iris = load_iris()

X, y = iris.data, iris.target

print('old_data_shape:\n',X.shape)

# SelectKBest 保留评分最高的 K 个特征

X_SelectKBest = SelectKBest(chi2, k=2)

X_new_SelectKBest = X_SelectKBest.fit_transform(X, y)

print('new_SelectKBest_data_shape:\n',X_new_SelectKBest.shape)

# SelectPercentile 保留最高得分百分比的特征

# 默认：f_classif 函数

# 对于回归: f_regression , mutual_info_regression

# 对于分类: chi2 , f_classif , mutual_info_classif

X_SelectPercentile = SelectPercentile(chi2,percentile=50)

X_new_SelectPercentile = X_SelectPercentile.fit_transform(X, y)

print('new_SelectPercentile_data_shape:\n',X_new_SelectPercentile.shape)

# Attributes

# scores_ : array-like, shape=(n_features,)

# Scores of features.

# pvalues_ : array-like, shape=(n_features,)

# p-values of feature scores, None if `score_func` returned only scores.

print('KBest scores： \n',X_SelectKBest.scores_)

print('Percentile scores：\n',X_SelectPercentile.scores_)

print('KBest pvalues_： \n',X_SelectKBest.pvalues_)

print('Percentile pvalues_：\n',X_SelectPercentile.pvalues_)

# 查看 选择的那两列

# 对比发现 SelectKBest 和 SelectPercentile 挑选的都是最后两列

print('old_data_10:\n',X[0:10,:])

print('X_new_SelectKBest:\n',X_new_SelectKBest[0:10,:])

print('X_new_SelectPercentile:\n',X_new_SelectPercentile[0:10,:])

运行结果

old_data_shape:

(150, 4)

new_SelectKBest_data_shape:

(150, 2)

new_SelectPercentile_data_shape:

(150, 2)

KBest scores：

[ 10.81782088 3.59449902 116.16984746 67.24482759]

Percentile scores：

[ 10.81782088 3.59449902 116.16984746 67.24482759]

KBest pvalues_：

[4.47651499e-03 1.65754167e-01 5.94344354e-26 2.50017968e-15]

Percentile pvalues_：

[4.47651499e-03 1.65754167e-01 5.94344354e-26 2.50017968e-15]

old_data_10:

[[5.1 3.5 1.4 0.2]

[4.9 3. 1.4 0.2]

[4.7 3.2 1.3 0.2]

[4.6 3.1 1.5 0.2]

[5. 3.6 1.4 0.2]

[5.4 3.9 1.7 0.4]

[4.6 3.4 1.4 0.3]

[5. 3.4 1.5 0.2]

[4.4 2.9 1.4 0.2]

[4.9 3.1 1.5 0.1]]

X_new_SelectKBest:

[[1.4 0.2]

[1.4 0.2]

[1.3 0.2]

[1.5 0.2]

[1.4 0.2]

[1.7 0.4]

[1.4 0.3]

[1.5 0.2]

[1.4 0.2]

[1.5 0.1]]

X_new_SelectPercentile:

[[1.4 0.2]

[1.4 0.2]

[1.3 0.2]

[1.5 0.2]

[1.4 0.2]

[1.7 0.4]

[1.4 0.3]

[1.5 0.2]

[1.4 0.2]

[1.5 0.1]]

# 对比 不同的 score_func 有什么区别，是不是不同的函数选择的特征都是一样的？

# 说明：这里小样本，结果具有偶然性，实际开发项目时可以对比一下。

# 个人觉得，根据理论方法，选择不同的score_func会选择出不同的特征。

# 分类函数比较: chi2 , f_classif , mutual_info_classif

X_SelectPercentile_chi2 = SelectPercentile(chi2,percentile=50)

X_SelectPercentile_f_classif = SelectPercentile(f_classif,percentile=50)

X_SelectPercentile_mutual_info_classif = SelectPercentile(mutual_info_classif,percentile=50)

X_new_SelectKBest_chi2 = X_SelectPercentile_chi2.fit_transform(X, y)

X_new_SelectKBest_f_classif = X_SelectPercentile_f_classif.fit_transform(X, y)

X_new_SelectKBest_mutual_info_classif = X_SelectPercentile_mutual_info_classif.fit_transform(X, y)

print('chi2 scores： \n',X_SelectPercentile_chi2.scores_)

print('f_classif scores： \n',X_SelectPercentile_f_classif.scores_)

print('mutual_info_classif scores：\n',X_SelectPercentile_mutual_info_classif.scores_)

print('chi2 pvalues_： \n',X_SelectPercentile_chi2.pvalues_)

print('f_classif pvalues_： \n',X_SelectPercentile_f_classif.pvalues_)

print('mutual_info_classif pvalues_：\n',X_SelectPercentile_mutual_info_classif.pvalues_)

运行结果

chi2 scores：

[ 10.81782088 3.59449902 116.16984746 67.24482759]

f_classif scores：

[ 119.26450218 47.3644614 1179.0343277 959.32440573]

mutual_info_classif scores：

[0.47980703 0.26374852 0.98914392 0.97626377]

chi2 pvalues_：

[4.47651499e-03 1.65754167e-01 5.94344354e-26 2.50017968e-15]

f_classif pvalues_：

[1.66966919e-31 1.32791652e-16 3.05197580e-91 4.37695696e-85]

mutual_info_classif pvalues_：

None

# 互信息

from minepy import MINE

import pandas as pd

def my_mine(x_data,y):

# 定义一个列表，将互信息的值放在列表中

mic_score = []

m = MINE()

# 取出 x 的每一列

for column in x_data.columns:

m.compute_score(x_data[column],y)

mic_score.append(m.mic())

# 转为 dataframe 并将其 列名 一一对应

mic_score = pd.DataFrame(mic_score )

# 更改 索引

mic_score.index = x_data.columns

# 更改列名

mic_score = mic_score.rename(columns = {0:'MICSCORE'})

return mic_score

from sklearn.datasets import load_iris

iris = load_iris()

x ,y= pd.DataFrame(iris.data).rename(columns={0:'a',1:'b',2:'c',3:'d'}), iris.target

my_mine(x,y)

# 说明：

# sklearn.feature_selection 中封装的互信息，和此处有一点的差异，但是差异不大

# 推荐使用 sklearn 官网的sklearn.feature_selection 库

# [0.50056966 0.24913013 0.99109155 0.99213811]

运行结果

MICSCORE

a 0.642196

b 0.401504

c 0.918296

d 0.918296

1.3 特征选择------Pearson相关系数 (Pearson Correlation)

皮尔森相关系数是一种最简单的，能帮助理解特征和响应变量之间关系的方法，该方法衡量的是变量之间的线性相关性，结果的取值区间为[-1，1]，-1表示完全的负相关，+1表示完全的正相关，0表示没有线性相关。 Pearson相关系数的一个明显缺陷是，作为特征排序机制，只对线性关系敏感。如果关系是非线性的，即便两个变量具有一一对应的关系，Pearson相关性也可能会接近0。

# Scipy的 pearsonr 方法

import pandas as pd

from pandas import Series

from scipy.stats import pearsonr

def my_pearsonr(x_data,y_data):

# 遍历 x_data的每一列

sorce_p_value = []

for column in x_data.columns:

sorce = pearsonr(x_data[column], y_data)

sorce_p_value.append(list(sorce))

# 转化为 dataframe

sorce_p_value = pd.DataFrame(sorce_p_value)

# 更改 索引

sorce_p_value.index = x_data.columns

# 更改列名

sorce_p_value = sorce_p_value.rename(columns = {0:'sorce',1:'p-value'})

return sorce_p_value

from sklearn.datasets import load_iris

x ,y= pd.DataFrame(iris.data).rename(columns={0:'a',1:'b',2:'c',3:'d'}), iris.target

sorce_p_value = my_pearsonr(x,y)

sorce_p_value

运行结果

sorce p-value

a 0.782561 2.890478e-32

b -0.419446 9.159985e-08

c 0.949043 4.155478e-76

d 0.956464 4.775002e-81

# 相关系数法对非线性不明感

import numpy as np

x = np.random.uniform(-1, 1, 100000)

print(pearsonr(x, x**2))

运行结果

(-0.006741987892572058, 0.03300673581074622)

1.4 特征选择------距离相关系数 (Distance correlation)

# 预备知识

# 维度改变

## atleast_xd 支持将输入数据直接视为 x维。这里的 x 可以表示：1，2，3。

from scipy.spatial.distance import pdist, squareform

import numpy as np

np.atleast_1d([1])

np.atleast_2d([1])

np.atleast_3d([1])

print('np.atleast_1d:',np.atleast_1d([1]))

print('np.atleast_2d:',np.atleast_2d([1]))

print('np.atleast_3d:',np.atleast_3d([1]))

# 欧几里德距离

# Pairwise distance between pairs of object(Pdist函数用于各种距离的生成)

X = [[1,2,3,4,5],

[2,4,6,8,10],

[3,6,9,12,15],

[1,2,3,4,5]]

x_pdist=pdist(X)

# pdist(x) 计算m*n的数据矩阵中对象之间的欧几里得距离

# 得到的是一个长度为m(m-1)/2的距离向量，距离是按顺序排列的(2,1)，(3,1)…….(m,1),(3,2)……..(m,2)………(m,m-1)

print('欧氏距离:\n',x_pdist)

#矩阵中的(i,j),对应于原始数据集中的i列和j列之间的欧氏距离。

print('矩阵:\n',squareform(x_pdist))

# 连乘操作

print(np.prod([1,2])) # 1*2

print(np.prod([[1,2],[3,4]])) # (1*2) * (3*4)

print(np.prod([[1,2],[3,4]],axis=0)) # (1*3) * (2*4)

print(np.prod([[1,2],[3,4]],axis=1)) # (1*2) (3*4)

运行结果

np.atleast_1d: [1]

np.atleast_2d: [[1]]

np.atleast_3d: [[[1]]]

欧氏距离:

[ 7.41619849 14.83239697 0. 7.41619849 7.41619849 14.83239697]

矩阵:

[[ 0. 7.41619849 14.83239697 0. ]

[ 7.41619849 0. 7.41619849 7.41619849]

[14.83239697 7.41619849 0. 14.83239697]

[ 0. 7.41619849 14.83239697 0. ]]

2

24

[3 8]

[ 2 12]

from scipy.spatial.distance import pdist, squareform

import numpy as np

def my_distcorr(X, Y):

""" Compute the distance correlation function

>>> a = [1,2,3,4,5]

>>> b = np.array([2,4,6,8,10])

>>> 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

a = [1,2,3,4,5]

b = np.array([2,4,6,8,10])

print(my_distcorr(a, b))

# 查看 Pearson相关系数 接近零，其 距离相关系数

my_distcorr(a, [i*i for i in a])

运行结果

1.0

0.9869160440537483

1.5 特征选择------基于模型的特征排序 (Model based ranking)

from sklearn.cross_validation import cross_val_score, ShuffleSplit

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"]

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", # 注意X[:, i]和X[:, i:i+1]的区别

cv=ShuffleSplit(len(X), 5, .3),

# cv=5,

n_jobs=-1)

scores.append((format(np.mean(score), '.3f'), names[i]))

print(sorted(scores, reverse=True))

运行结果

[('0.665', 'LSTAT'), ('0.545', 'RM'), ('0.394', 'NOX'), ('0.327', 'INDUS'), ('0.320', 'TAX'), ('0.318', 'PTRATIO'), ('0.201', 'CRIM'), ('0.179', 'ZN'), ('0.161', 'RAD'), ('0.124', 'DIS'), ('0.114', 'B'), ('0.078', 'AGE'), ('0.031', 'CHAS')]

1.6 特征选择------递归特征消除 (Recursive Feature Elimination)

属性：

support_：选择的特征，布尔类型。

ranking_：特征的排名位置。 选择的特征被指定为等级1。越差的特征，等级数越高。

# 属性：

# support_：选择的特征，布尔类型。

# ranking_：特征的排名位置。 选择的特征被指定为等级1。越差的特征，等级数越高。

# 不采用交叉验证

from sklearn.datasets import make_friedman1

from sklearn.feature_selection import RFE

from sklearn.svm import SVR

X, y = make_friedman1(n_samples=50, n_features=10, random_state=0)

# 可以选择不同的基模型

estimator = SVR(kernel="linear")

#参数estimator为基模型，参数n_features_to_select为选择的特征个数

selector = RFE(estimator, n_features_to_select=5)

selector = selector.fit(X, y)

print('selector.support_:\n',selector.support_)

print('selector.ranking_:\n',selector.ranking_)

print(selector)

运行结果

selector.support_:

[ True True True True True False False False False False]

selector.ranking_:

[1 1 1 1 1 6 4 3 2 5]

RFE(estimator=SVR(C=1.0, cache_size=200, coef0=0.0, degree=3, epsilon=0.1, gamma='auto',

kernel='linear', max_iter=-1, shrinking=True, tol=0.001, verbose=False),

n_features_to_select=5, step=1, verbose=0)

# # 采用交叉验证，

from sklearn.datasets import make_friedman1

from sklearn.feature_selection import RFECV

from sklearn.svm import SVR

X, y = make_friedman1(n_samples=50, n_features=10, random_state=0)

estimator = SVR(kernel="linear")

# 参数：

# step： 如果大于或等于1，那么`step`对应于(整数)每次迭代时要删除的要素数。

# 如果在(0.0,1.0)内，那么`step`对应于百分比(向下舍入)要在每次迭代时删除的要素。

selector = RFECV(estimator, step=1, cv=6)

selector = selector.fit(X, y)

print('selector.support_:\n',selector.support_)

print('selector.ranking_:\n',selector.ranking_)

运行结果

selector.support_:

[ True True True True True False False False False False]

selector.ranking_:

[1 1 1 1 1 6 4 3 2 5]

import matplotlib.pyplot as plt

from sklearn.svm import SVC

from sklearn.cross_validation import StratifiedKFold

from sklearn.feature_selection import RFECV

from sklearn.datasets import make_classification

# Build a classification task using 3 informative features

X, y = make_classification(n_samples=1000, n_features=25, n_informative=3,

n_redundant=2, n_repeated=0, n_classes=8,

n_clusters_per_class=1, random_state=0)

# Create the RFE object and compute a cross-validated score.

svc = SVC(kernel="linear")

# The "accuracy" scoring is proportional to the number of correct

# classifications

rfecv = RFECV(estimator=svc, step=1, cv=StratifiedKFold(y, 2),

scoring='accuracy')

rfecv.fit(X, y)

print("Optimal number of features : %d" % rfecv.n_features_)

# Plot number of features VS. cross-validation scores

plt.figure()

plt.xlabel("Number of features selected")

plt.ylabel("Cross validation score (nb of correct classifications)")

plt.plot(range(1, len(rfecv.grid_scores_) + 1), rfecv.grid_scores_)

plt.show()

运行结果

Optimal number of features : 3

1.7 特征选择------基于L1的特征选择 (L1-based feature selection)

使用L1范数作为惩罚项的线性模型(Linear models)会得到稀疏解：大部分特征对应的系数为0。当你希望减少特征的维度以用于其它分类器时，可以通过 feature_selection.SelectFromModel 来选择不为0的系数。

特别指出，常用于此目的的稀疏预测模型有 linear_model.Lasso(回归),linear_model.LogisticRegression 和 svm.LinearSVC(分类).

SelectFromModel 作为meta-transformer，能够用于拟合后任何拥有coef_或feature_importances_ 属性的预测模型。 如果特征对应的coef_ 或 feature_importances_ 值低于设定的阈值threshold，那么这些特征将被移除。除了手动设置阈值，也可通过字符串参数调用内置的启发式算法(heuristics)来设置阈值，包括：平均值(“mean”), 中位数(“median”)以及他们与浮点数的乘积，如”0.1*mean”。

# svm.LinearSVC

from sklearn.svm import LinearSVC

from sklearn.datasets import load_iris

from sklearn.feature_selection import SelectFromModel

iris = load_iris()

X, y = iris.data, iris.target

lsvc = LinearSVC(C=0.01, penalty="l1", dual=False).fit(X, y)

model = SelectFromModel(lsvc, prefit=True)

X_new = model.transform(X)

print('X_shape:',X.shape)

print('X_new.shape:',X_new.shape)

运行结果

X_shape: (150, 4)

X_new.shape: (150, 3)

## 带L1和L2惩罚项的逻辑回归作为基模型的特征选择

# 对于SVM和逻辑回归，参数C控制稀疏性：C越小，被选中的特征越少。对于Lasso，参数alpha越大，被选中的特征越少。

from sklearn.feature_selection import SelectFromModel

from sklearn.linear_model import LogisticRegression as LR

X_new = SelectFromModel(LR(C=0.1)).fit_transform(iris.data, iris.target)

print('X_new.shape:',X_new.shape)

运行结果

X_new.shape: (150, 2)

# Use SelectFromModel meta-transformer along with Lasso to select the best couple of features from the Boston dataset.

import matplotlib.pyplot as plt

import numpy as np

from sklearn.datasets import load_boston

from sklearn.feature_selection import SelectFromModel

from sklearn.linear_model import LassoCV

# Load the boston dataset.

boston = load_boston()

X, y = boston['data'], boston['target']

# We use the base estimator LassoCV since the L1 norm promotes sparsity of features.

clf = LassoCV()

# Set a minimum threshold of 0.25

sfm = SelectFromModel(clf, threshold=0.25)

sfm.fit(X, y)

n_features = sfm.transform(X).shape[1]

print('X_shape:',X.shape)

print('new_X_shape:',sfm.transform(X).shape)

# Reset the threshold till the number of features equals two.

# Note that the attribute can be set directly instead of repeatedly

# fitting the metatransformer.

while n_features > 2:

sfm.threshold += 0.1

X_transform = sfm.transform(X)

n_features = X_transform.shape[1]

# Plot the selected two features from X.

plt.title(

"Features selected from Boston using SelectFromModel with "

"threshold %0.3f." % sfm.threshold)

feature1 = X_transform[:, 0]

feature2 = X_transform[:, 1]

plt.plot(feature1, feature2, 'r.')

plt.xlabel("Feature number 1")

plt.ylabel("Feature number 2")

plt.ylim([np.min(feature2), np.max(feature2)])

plt.show()

运行结果

X_shape: (506, 13)

new_X_shape: (506, 5)

1.8 特征选择------ 随机稀疏模型 (Randomized sparse models)

基于L1的稀疏模型的局限在于，当面对一组互相关的特征时，它们只会选择其中一项特征。为了减轻该问题的影响,可以使用随机化技术，通过多次重新估计稀疏模型来扰乱设计矩阵，或通过多次下采样数据来统计一个给定的回归量被选中的次数。

RandomizedLasso 实现了使用这项策略的Lasso，RandomizedLogisticRegression 使用逻辑回归，适用于分类任务。要得到整个迭代过程的稳定分数，你可以使用 lasso_stability_path。

注意到对于非零特征的检测，要使随机稀疏模型比标准F统计量更有效， 那么模型的参考标准需要是稀疏的，换句话说，非零特征应当只占一小部分。

import matplotlib.pyplot as plt

import numpy as np

from scipy import linalg

from sklearn.linear_model import (RandomizedLasso, lasso_stability_path,

LassoLarsCV)

from sklearn.feature_selection import f_regression

from sklearn.preprocessing import StandardScaler

from sklearn.metrics import auc, precision_recall_curve

from sklearn.ensemble import ExtraTreesRegressor

from sklearn.utils.extmath import pinvh

def mutual_incoherence(X_relevant, X_irelevant):

"""Mutual incoherence, as defined by formula (26a) of [Wainwright2006].

"""

projector = np.dot(np.dot(X_irelevant.T, X_relevant),

pinvh(np.dot(X_relevant.T, X_relevant)))

return np.max(np.abs(projector).sum(axis=1))

for conditioning in (1, 1e-4):

###########################################################################

# Simulate regression data with a correlated design

n_features = 501

n_relevant_features = 3

noise_level = .2

coef_min = .2

# The Donoho-Tanner phase transition is around n_samples=25: below we

# will completely fail to recover in the well-conditioned case

n_samples = 25

block_size = n_relevant_features

rng = np.random.RandomState(42)

# The coefficients of our model

coef = np.zeros(n_features)

coef[:n_relevant_features] = coef_min + rng.rand(n_relevant_features)

# The correlation of our design: variables correlated by blocs of 3

corr = np.zeros((n_features, n_features))

for i in range(0, n_features, block_size):

corr[i:i + block_size, i:i + block_size] = 1 - conditioning

corr.flat[::n_features + 1] = 1

corr = linalg.cholesky(corr)

# Our design

X = rng.normal(size=(n_samples, n_features))

X = np.dot(X, corr)

# Keep [Wainwright2006] (26c) constant

X[:n_relevant_features] /= np.abs(

linalg.svdvals(X[:n_relevant_features])).max()

X = StandardScaler().fit_transform(X.copy())

# The output variable

y = np.dot(X, coef)

y /= np.std(y)

# We scale the added noise as a function of the average correlation

# between the design and the output variable

y += noise_level * rng.normal(size=n_samples)

mi = mutual_incoherence(X[:, :n_relevant_features],

X[:, n_relevant_features:])

###########################################################################

# Plot stability selection path, using a high eps for early stopping

# of the path, to save computation time

alpha_grid, scores_path = lasso_stability_path(X, y, random_state=42,

eps=0.05)

plt.figure()

# We plot the path as a function of alpha/alpha_max to the power 1/3: the

# power 1/3 scales the path less brutally than the log, and enables to

# see the progression along the path

hg = plt.plot(alpha_grid[1:] ** .333, scores_path[coef != 0].T[1:], 'r')

hb = plt.plot(alpha_grid[1:] ** .333, scores_path[coef == 0].T[1:], 'k')

ymin, ymax = plt.ylim()

plt.xlabel(r'$(\alpha / \alpha_{max})^{1/3}$')

plt.ylabel('Stability score: proportion of times selected')

plt.title('Stability Scores Path - Mutual incoherence: %.1f' % mi)

plt.axis('tight')

plt.legend((hg[0], hb[0]), ('relevant features', 'irrelevant features'),

loc='best')

###########################################################################

# Plot the estimated stability scores for a given alpha

# Use 6-fold cross-validation rather than the default 3-fold: it leads to

# a better choice of alpha:

# Stop the user warnings outputs- they are not necessary for the example

# as it is specifically set up to be challenging.

with warnings.catch_warnings():

warnings.simplefilter('ignore', UserWarning)

lars_cv = LassoLarsCV(cv=6).fit(X, y)

# Run the RandomizedLasso: we use a paths going down to .1*alpha_max

# to avoid exploring the regime in which very noisy variables enter

# the model

alphas = np.linspace(lars_cv.alphas_[0], .1 * lars_cv.alphas_[0], 6)

clf = RandomizedLasso(alpha=alphas, random_state=42).fit(X, y)

trees = ExtraTreesRegressor(100).fit(X, y)

# Compare with F-score

F, _ = f_regression(X, y)

plt.figure()

for name, score in [('F-test', F),

('Stability selection', clf.scores_),

('Lasso coefs', np.abs(lars_cv.coef_)),

('Trees', trees.feature_importances_),

]:

precision, recall, thresholds = precision_recall_curve(coef != 0,

score)

plt.semilogy(np.maximum(score / np.max(score), 1e-4),

label="%s. AUC: %.3f" % (name, auc(recall, precision)))

plt.plot(np.where(coef != 0)[0], [2e-4] * n_relevant_features, 'mo',

label="Ground truth")

plt.xlabel("Features")

plt.ylabel("Score")

# Plot only the 100 first coefficients

plt.xlim(0, 100)

plt.legend(loc='best')

plt.title('Feature selection scores - Mutual incoherence: %.1f'

% mi)

plt.show()

1.9 特征选择------ 基于树的特征选择 (Tree-based feature selection)

基于树的预测模型,能够用来计算特征的重要程度，因此能用来去除不相关的特征:

import numpy as np

import matplotlib.pyplot as plt

from sklearn.datasets import make_classification

from sklearn.ensemble import ExtraTreesClassifier

# Build a classification task using 3 informative features

X, y = make_classification(n_samples=1000,

n_features=10,

n_informative=3,

n_redundant=0,

n_repeated=0,

n_classes=2,

random_state=0,

shuffle=False)

# Build a forest and compute the feature importances

forest = ExtraTreesClassifier(n_estimators=250,

random_state=0)

forest.fit(X, y)

importances = forest.feature_importances_

std = np.std([tree.feature_importances_ for tree in forest.estimators_],

axis=0)

indices = np.argsort(importances)[::-1]

# Print the feature ranking

print("Feature ranking:")

for f in range(X.shape[1]):

print("%d. feature %d (%f)" % (f + 1, indices[f], importances[indices[f]]))

# Plot the feature importances of the forest

plt.figure()

plt.title("Feature importances")

plt.bar(range(X.shape[1]), importances[indices],

color="r", yerr=std[indices], align="center")

plt.xticks(range(X.shape[1]), indices)

plt.xlim([-1, X.shape[1]])

plt.show()

运行结果

Feature ranking:

1. feature 1 (0.295902)

2. feature 2 (0.208351)

3. feature 0 (0.177632)

4. feature 3 (0.047121)

5. feature 6 (0.046303)

6. feature 8 (0.046013)

7. feature 7 (0.045575)

8. feature 4 (0.044614)

9. feature 9 (0.044577)

10. feature 5 (0.043912)

2019061208.png

# 用于人脸识别数据

from time import time

import matplotlib.pyplot as plt

from sklearn.datasets import fetch_olivetti_faces

from sklearn.ensemble import ExtraTreesClassifier

# Number of cores to use to perform parallel fitting of the forest model

n_jobs = 1

# Load the faces dataset

data = fetch_olivetti_faces()

X = data.images.reshape((len(data.images), -1))

y = data.target

mask = y < 5 # Limit to 5 classes

X = X[mask]

y = y[mask]

# Build a forest and compute the pixel importances

print("Fitting ExtraTreesClassifier on faces data with %d cores..." % n_jobs)

t0 = time()

forest = ExtraTreesClassifier(n_estimators=1000,

max_features=128,

n_jobs=n_jobs,

random_state=0)

forest.fit(X, y)

print("done in %0.3fs" % (time() - t0))

importances = forest.feature_importances_

importances = importances.reshape(data.images[0].shape)

# Plot pixel importances

plt.matshow(importances, cmap=plt.cm.hot)

plt.title("Pixel importances with forests of trees")

plt.show()