特征工程之特征选择
目录
简介
1 Filter（过滤式选择）
1.1 移除低方差特征（variance threshold）
1.2 信息增益（information gain）
1.3 单变量特征选择 (Univariate feature selection)
1.3.1 卡方检验 (chi-square test)
1.3.2 Pearson 相关系数 (Pearson Correlation)
1.3.3 费雪分数（fisher score）
1.4 Relief(Relevant Features)
2 Wrapper（包裹式选择）
2.1 递归特征消除（recursive feature elimination）
2.2 遗传算法（genetic algorithms）
2.3 拉斯维加斯方法（Las Vegas Wrapper）
3 Embedded（嵌入式选择）
3.1 L1（LASSO ）
3.2 决策树、随机森林、极限树
简介
随着科技的发展，数据量越来越大，在建立模型时，考虑的数据维度越来越广，所以建立模型前的降维越来越重要，降维的方式一般有两种，其一是用原始的维度合成新的重要维度，例如SVD和PCA，其二是在原始的维度中保留重要维度，剔除次要维度。
第一种降维方式的优点是可以简单高效的合成重要维度，缺点是合成的维度失去其现实中的可解释性。
第二种降维方式的优点是保持其原始的可解释性，缺点是计算比第一种相对复杂。
所以在图片识别，声音识别等不需要解释中间变量的模型领域经常采用第一种建模方式，在金融领域往往需要追求变量的可解释性，所以经常采用第二种降维方式。
本文主要采用第二种降维方式，也就是在原有的特征中进行特征选择。这种降维的方法其关键分为两大步：
第一步：如何遍历所有特征。
第二步：如何判断特征的重要性。
如何遍历所有特征。
在第一步中可使用前向搜索，后向搜索和双向搜索方法遍历所有特征，这三种遍历方法全部是贪心算法，最求每一步最优，不是全局最优。
1、前向搜索：首先对单特征进行遍历，找到此次遍历中最重要的特征，然后保留这个特征，遍历这个特征和其他任一特征的组合的主要性，找到第二重要的特征，保留这两个特征，遍历这个两个特征和其他任一特征的组合的重要性，找到第三个特征，保留着三个特征，以次下去，即可对全部特征进行重要性排序。
2、后向搜索：与前向搜索相反，开始在全部特征中遍历剔除一个特征，找到影响重要性最小的特征，将其剔除，然后，在剩下的n-1个特征中遍历剔除一个特征，找到影响重要性最小的特征，以此下去，即可对全部特征进行重要性排序。
3、双向搜索：同时进行向前和向后搜索，但是一定要注意，向后搜索一定不要剔除向前搜索选中的特征。
如何判断特征的重要性。
如何判断特征的重要性，有很多方法，如：信息熵、相关系数、信息价值（Information Value），具体方法下面会具体介绍。
主要方法分类
特征筛选降维的方法主要分为三大类：过滤法（Filter）、包裹法（Wrapper）、嵌入法（Embedded）。这三者的区别和具体算法下面进行具体介绍。
Filter（过滤法）
Filter的思想是，特征筛选和模型建立完全隔离开，筛选特征时，不考虑具体什么模型，只看特征对目标变量影响。
这种方法包括：移除低方差特征（variance threshold）、信息增益（information gain）、卡方检验 (chi-square test)、
Pearson 相关系数 (Pearson Correlation)、费雪分数（fisher score）
移除低方差特征（variance threshold）
其思想是，剔除方差小的特征，方差小的特征其值变化较小，认为区别力度不大，但是该方法只适用于离散型随机变量，若是连续性随机变脸需要进行woe封箱。
给出如下例子，可以看到第一个维度的取值为（0，0，1，0，0，0）其波动率非常小，所以认为这个特征应该清楚。
import numpy as np
# Variance 偏差 Threshold阈值
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=0.2)
sel.fit_transform(X)
array([[0, 1],
[1, 0],
[0, 0],
[1, 1],
[1, 0],
[1, 1]])
信息增益（information gain）
这里先给出信息熵的公式如下：
其中：D是一个集合，
对于集合D，考虑样本的属性factor1,在属性factor1上，可将集合D划分为
特征划分带来的信息增益值越大，其这个特征对结果影响越大，所以我们可以通过信息增益来判断特征的重要性，结合上面特征遍历的方法，我们就可以得到筛选特征的方法。
下面给出信息增益结合前向搜索的python例子：
#定义计算熵的函数
def ent(data):
prob1 = pd.value_counts(data) / len(data)
return sum(np.log2(prob1) * prob1 * (-1))
#定义计算信息增益的函数
def InformationGain(data,list1,str2):
e1 = data.groupby(list1).apply(lambda x:ent(x[str2])).reset_index()[0]
p1 = data.groupby(list1).size().reset_index(name="Time")['Time'] / len(data[str2])
e2 = sum(e1 * p1)
return ent(data[str2]) - e2
data1 = pd.DataFrame({'a':[1,1,1,5,3,2],
'b':[0,1,2,1,2,1],
'c':[0,1,0,2,3,4],
'd':[1,0,3,2,5,3],
'e':[1,1,3,1,2,3]})
data1[['a','b','c','d','e']]
list1=["a","c"]
label='e'
InformationGain(data1,list1,label)
list_1=['a','b','c','d']
list_2=[]
tmp_list=list_1
for factor in list_1:
factor_imp =pd.DataFrame([[i,InformationGain(data1,list_2+[i],label)] for i in tmp_list])
s1_argmax = factor_imp[factor_imp.iloc[:,1] == factor_imp.iloc[:,1].max()][0].iloc[0]
tmp_list=list(set(tmp_list).difference(set(s1_argmax)))
list_2=list_2+[s1_argmax]
print(list_2)
['d', 'b', 'a', 'c']
上面给出了，前向搜索结合信息增益筛选特征的例子。在很多情况下为了简化，只进行单变量特征选择。
单变量特征选择 (Univariate feature selection)
单变量特征选择不进行前向、后向或者双向搜索，只进行单一变量的影响排名。常用的方法有卡方检验，相关系数和费雪分数。
卡方检验 (chi-square test)
我们在讲解卡方检验以前，先来推导卡方分布。
卡方分布
设
其中：
以上分布被称为卡方分布。
证明：对于统计量
使用极坐标变换：
所以：
所以：
因为
所以：
当x趋于正无穷大时，有：
解得：
所以得到：
上面我们推导出来卡方分布，在推导的过程中，我们看到n个服从相同独立的正太分布随机变量的平方和服从卡方分布。
所以当我们构造的检验的统计量的表达形式是n个服从相同独立的正太分布随机变量的平方和时，这种检验方式称为卡方检验。
卡方统计量
卡方统计量的由来：
皮尔逊在衡量数据的分布与所选择的预期或假设分布之间的差异是，需要构造一个统计量，对这个统计量进行假设检验，皮尔逊首先考虑所有数据分布的频数
但是他发现不同数据的频数可能大于假设的分布的频率，也可能小于假设的分布的频率。所以不同数据的频率差之间会相互抵消。从而失去本来想衡量数据分布和假设分布差异的意义。
针对这一问题，皮尔逊引入了平方的概念，即把上面的和变成了：
所以上面的统计量就是最终的卡方统计量。
但是值得注意的是，这个统计量并不严格服从卡方分布。
理论次数
上面的条件也就是进行卡方检验的条件。当不符合上面的条件时，我们可以进行适当的处理，来使用卡方检验，这的处理方式不是本文重点，不展开讲述。
卡方检验
上面构造出了卡方统计量，对这个统计量进行检验的话，就是卡方检验，
卡方检验就是统计样本的实际观测值与理论推断值之间的偏离程度，实际观测值与理论推断值之间的偏离程度就决定卡方值的大小，卡方值越大，越不符合；卡方值越小，偏差越小，越趋于符合，若两个值完全相等时，卡方值就为0，表明理论值完全符合。
注意：卡方检验针对分类变量，没如果是连续变量，则要分组离散化。
卡方检验的原假设
卡方检验的用处主要有以下几种：
1、检验某个连续变量的分布是否与某种理论分布相一致。如是否符合正态分布、是否服从均匀分布、是否服从Poisson分布等。
2、检验某个分类变量各类的出现概率是否等于指定概率。如在36选7的彩票抽奖中，每个数字出现的概率是否各为1/36；掷硬币时，正反两面出现的概率是否均为0.5。
3、检验某两个分类变量是否相互独立。如吸烟(二分类变量：是、否)是否与呼吸道疾病(二分类变量：是、否)有关；产品原料种类(多分类变量)是否与产品合格(二分类变量)有关。
4、检验控制某种或某几种分类因素的作用以后，另两个分类变量是否相互独立。如在上例中，控制性别、年龄因素影响以后，吸烟是否和呼吸道疾病有关；控制产品加工工艺的影响后，产品原料类别是否与产品合格有关。
5、检验某两种方法的结果是否一致。如采用两种诊断方法对同一批人进行诊断，其诊断结果是否一致；采用两种方法对客户进行价值类别预测，预测结果是否一致。
可以看到，上面的用途中，第三种用途可以检验二分类和多分类特征和目标的关系。此时原假设
所以可以中卡方检验来进行特征筛选。下面给出具体的python例子：
from sklearn.datasets import load_iris # 鸢尾花数据集
# 用于特征选择的包
from sklearn.feature_selection import SelectKBest # 训练的套路
from sklearn.feature_selection import chi2
iris = load_iris()
x, y = iris.data, iris.target
# 参数解释：chi2 使用卡方检验，k=2 留下两个特征
selector = SelectKBest(chi2,k=2).fit(x,y)
print("scores_:",selector.scores_)
print("pvalues_:",selector.pvalues_)
print("selected index:",selector.get_support(True))
print("selected factor:",selector.get_support())
scores_: [ 10.81782088 3.59449902 116.16984746 67.24482759]
pvalues_: [4.47651499e-03 1.65754167e-01 5.94344354e-26 2.50017968e-15]
selected index: [2 3]
selected factor: [False False True True]
Pearson 相关系数 (Pearson Correlation)
这种特征选择的方式一般是模型中最长用的一种方式，需要注意的是，求相关系数时，不是用的某一特征和目标的相关系数，因为我们的模型并非线性回归，只有模型是线性回归时，求特征和目标的相关系数才有意义，在一般的模型中，我们求的是特征与特征之间的相关系数。我们会剔除一些相关性非常强的特征。
相关系数的python代码如下：
import numpy as np
from scipy.stats import pearsonr # sciPy 科学python哈哈，stats：统计学
from sklearn.datasets import load_iris # 鸢尾花数据集
iris = load_iris()
x, y = iris.data, iris.target
x_1=[a[0] for a in x]
x_2=[a[1] for a in x]
c=pearsonr(x_1, x_2)
print("相关系数:",c[0])
print("p-value:",c[1])
相关系数: -0.10936924995064935
p-value: 0.1827652152713665
费雪分数（fisher score）
我们在处理分类问题时，对于某个特征，我们希望这个特征在同类里面方差近可能的小，不同类之间方差近可能大。这样的特征肯定对分类有很大的帮助。
我们构造如下变量（fisher score）：
其中：
所以上面构造的这个fisher score越大，组内方差越小，组间的方差越大，也就是这个特征的效果越好。
相比于信息增益，fisher score 计算更简单，且跟准确。所以这种方法相对更好。
下面我们给出python例子：
import pandas as pd
import numpy as np
from sklearn.datasets import load_iris # 鸢尾花数据集
iris = load_iris()
x, y = iris.data, iris.target
data=pd.concat([pd.DataFrame(x),pd.DataFrame(y)],axis=1)
data.columns = range(5)
# # 计算fisher得分
items = list(range(4))
num_classes = len(set(y))
fisher_score = []
grouped = data.groupby([4], as_index=False)
n = [len(data[data[4] == k+1]) for k in range(num_classes)]
for i in items: # 遍历所有特征列
temp = grouped[i].agg({str(i)+'_mean': 'mean',
str(i)+'_std': 'std'}) # 已求出特征i在各类别k中的均值u_ik、方差p_ik
numerator = 0
denominator = 0
u_i = data[i].mean()
for k in range(num_classes):
n_k = n[k]
u_ik = temp.iloc[k, :][str(i)+'_mean']
p_ik = temp.iloc[k, :][str(i)+'_std']
numerator += n_k*(u_ik-u_i)**2
denominator += n_k*p_ik**2
fisher_score.append(numerator/denominator)
print("fisher score:",fisher_score)
fisher score: [1.8166050382319252, 0.8748323114550877, 21.986196311478416, 18.33205972771192]
Relief(Relevant Features)
Relief(Relevant Features)是一种著名的过滤式特征选择方法，该方法设计了一个“相关统计量”来度量特征的重要性，该统计量是一个向量，其每个分量分别对应一个初始特征，而特征子集的重要性则是由子集中每个特征对应的相关统计量分量之和来决定，于是，最终只需要指定一个阈值p,然后选择比p大的相关统计量分量对应的特征即可。
相关统计量向量
在上面的描述中，可以看出，这个相关统计量的构造及其重要，我们下面介绍其是怎么构造的。
Relief
使用Relief进行特征选择，可以处理二分类问题。我们先给出以下两个定义：
1、猜中近邻：对应一个样本
2、猜错近邻：对应一个样本
我们定义这个相关统计量向量的第j个分量为：
其中：j为样本的第j个特征，
注：其对连续的
为什么这样构造这个统计量向量。其实可以这样理解，对于属性j,当样本i与猜中近临的距离小于猜错近临的距离时，说明这个属性对分类是有正向帮助的，此时公式中
Relief-F
上面给出的是二分类问题时Relief的原理，对于多分类问题，一般使用Relief-F，其思想和方法和Relief几乎一摸一样，就是在相关统计量分量的构造上有一点一点区别，如下：
其中：
下面是python的一个例子：
# _*_ coding:utf8 _*_
import pandas as pd
import numpy as np
import numpy.linalg as la
import random
# 异常类
class ReliefError:
pass
class Relief:
def __init__(self, data_df, sample_rate, t, k):
"""
#
:param data_df: 数据框（字段为特征，行为样本）
:param sample_rate: 抽样比例
:param t: 统计量分量阈值
:param k: k近邻的个数
"""
self.__data = data_df
self.__feature = data_df.columns
self.__sample_num = int(round(len(data_df) * sample_rate))
self.__t = t
self.__k = k
# 数据处理（将离散型数据处理成连续型数据，比如字符到数值）
def get_data(self):
new_data = pd.DataFrame()
for one in self.__feature[:-1]:
col = self.__data[one]
if (str(col.tolist()[0]).split(".")[0]).isdigit() or str(col.tolist()[0]).isdigit() or (str(col.tolist()[0]).split('-')[-1]).split(".")[-1].isdigit():
new_data[one] = self.__data[one]
# print '%s 是数值型' % one
else:
# print '%s 是离散型' % one
keys = list(set(list(col)))
values = list(xrange(len(keys)))
new = dict(zip(keys, values))
new_data[one] = self.__data[one].map(new)
new_data[self.__feature[-1]] = self.__data[self.__feature[-1]]
return new_data
# 返回一个样本的k个猜中近邻和其他类的k个猜错近邻
def get_neighbors(self, row):
df = self.get_data()
row_type = row[df.columns[-1]]
right_df = df[df[df.columns[-1]] == row_type].drop(columns=[df.columns[-1]])
aim = row.drop(df.columns[-1])
f = lambda x: eulidSim(np.mat(x), np.mat(aim))
right_sim = right_df.apply(f, axis=1)
right_sim_two = right_sim.drop(right_sim.idxmin())
right = dict()
right[row_type] = (right_sim_two.sort_values().index[0:self.__k]).tolist()
# print list(right_sim_two.sort_values().index[0:self.__k])
types =list(set(df[df.columns[-1]]) - set([row_type]))
wrong = dict()
for one in types:
wrong_df = df[df[df.columns[-1]] == one].drop(columns=[df.columns[-1]])
wrong_sim = wrong_df.apply(f, axis=1)
wrong[one] = wrong_sim.sort_values().index[0:self.__k].tolist()
#print (right, wrong)
return right, wrong
# 计算特征权重
def get_weight(self, feature, index, NearHit, NearMiss):
# data = self.__data.drop(self.__feature[-1], axis=1)
data = self.__data
row = data.iloc[index]
right = 0
for one in list(NearHit.values())[0]:
nearhit = data.iloc[one]
if (str(row[feature]).split(".")[0]).isdigit() or str(row[feature]).isdigit() or (str(row[feature]).split('-')[-1]).split(".")[-1].isdigit():
max_feature = data[feature].max()
min_feature = data[feature].min()
right_one = pow(round(abs(row[feature] - nearhit[feature]) / (max_feature - min_feature), 2), 2)
else:
right_one = 0 if row[feature] == nearhit[feature] else 1
right += right_one
right_w = round(right / self.__k, 2)
wrong_w = 0
# 样本row所在的种类占样本集的比例
p_row = round(float(list(data[data.columns[-1]]).count(row[data.columns[-1]])) / len(data), 2)
for one in NearMiss.keys():
# 种类one在样本集中所占的比例
p_one = round(float(list(data[data.columns[-1]]).count(one)) / len(data), 2)
wrong_one = 0
for i in NearMiss[one]:
nearmiss = data.iloc[i]
if (str(row[feature]).split(".")[0]).isdigit() or str(row[feature]).isdigit() or (str(row[feature]).split('-')[-1]).split(".")[-1].isdigit():
max_feature = data[feature].max()
min_feature = data[feature].min()
wrong_one_one = pow(round(abs(row[feature] - nearmiss[feature]) / (max_feature - min_feature), 2), 2)
else:
wrong_one_one = 0 if row[feature] == nearmiss[feature] else 1
wrong_one += wrong_one_one
wrong = round(p_one / (1 - p_row) * wrong_one / self.__k, 2)
wrong_w += wrong
w = wrong_w - right_w
return w
# 过滤式特征选择
def reliefF(self):
sample = self.get_data()
# print sample
m, n = np.shape(self.__data) # m为行数，n为列数
score = []
sample_index = random.sample(range(0, m), self.__sample_num)
#print ('采样样本索引为 %s ' % sample_index)
num = 1
for i in sample_index: # 采样次数
one_score = dict()
row = sample.iloc[i]
NearHit, NearMiss = self.get_neighbors(row)
#print ('第 %s 次采样，样本index为 %s，其NearHit k近邻行索引为 %s ，NearMiss k近邻行索引为 %s' % (num, i, NearHit, NearMiss))
for f in self.__feature[0:-1]:
w = self.get_weight(f, i, NearHit, NearMiss)
one_score[f] = w
#print ('特征 %s 的权重为 %s.' % (f, w))
score.append(one_score)
num += 1
f_w = pd.DataFrame(score)
#print('采样各样本特征权重如下：')
#print(f_w)
print('平均特征权重如下：')
print(f_w.mean())
return f_w.mean()
# 返回最终选取的特征
def get_final(self):
f_w = pd.DataFrame(self.reliefF(), columns=['weight'])
final_feature_t = f_w[f_w['weight'] > self.__t]
print(final_feature_t)
# final_feature_k = f_w.sort_values('weight').head(self.__k)
# print final_feature_k
return final_feature_t
# 几种距离求解
def eulidSim(vecA, vecB):
return la.norm(vecA - vecB)
def cosSim(vecA, vecB):
"""
:param vecA: 行向量
:param vecB: 行向量
:return: 返回余弦相似度（范围在0-1之间）
"""
num = float(vecA * vecB.T)
denom = la.norm(vecA) * la.norm(vecB)
cosSim = 0.5 + 0.5 * (num / denom)
return cosSim
def pearsSim(vecA, vecB):
if len(vecA) < 3:
return 1.0
else:
return 0.5 + 0.5 * np.corrcoef(vecA, vecB, rowvar=0)[0][1]
if __name__ == '__main__':
iris = load_iris()
x, y = iris.data, iris.target
data=pd.concat([pd.DataFrame(x),pd.DataFrame(y)],axis=1)
data.columns = range(5)
f = Relief(data, 1, 0.2, 2)
f.reliefF()
平均特征权重如下：
0 0.047867
1 0.084133
2 0.163933
3 0.199867
dtype: float64
Wrapper（包裹法）
Wrapper（包裹法）的主要思想是选择特征时，会考虑接下来具体使用什么模型，不会像过滤法一样，将模型和特征选择完全分离开。
这种方法主要包括：递归特征消除（recursive feature elimination）、序列特征选择（sequential feature selection algorithms）、遗传算法（genetic algorithms）
递归特征消除（RFE）
这种算法的做法是，使用前向或者后向遍历特征，这里判断特征重要性，不是通过某一个值，而是对这些特征直接进行模型训练，用贪心算法找到每一步重要的特征或者剔除不重要的特征。
这种做法的缺点是效果的稳定性完全取决于模型的稳定性，像线性回归这种不稳定的模型，RFE的效果不是很稳定。
下面给出RFE的python例子：
"""
Created on 2020-09-12
@author:dingweijie
递归特征消除（RFE）与多分类里的逻辑回归
"""
from sklearn.model_selection import train_test_split ###数据集划分训练集和测试集
from sklearn.feature_selection import RFECV
from sklearn.linear_model import LogisticRegressionCV
from sklearn.model_selection import StratifiedKFold #多分类的交叉验证
from sklearn.datasets import load_iris ##数据集
import pandas as pd
import numpy as np
def notch_diff_table(actual,pred):
df=pd.DataFrame({'actual':np.ravel(actual),'pred':pred })
df=df.eval('notch=actual-pred')
df=df.eval('abs_notch=abs(actual-pred)')
df_notch=pd.DataFrame(df.notch.value_counts(normalize=True).mul(100).sort_index(axis=0)).round(2)
df_notch.notch=['%.2f%%'% notch for notch in df_notch.notch ]
df_abs_notch=pd.DataFrame(df['abs_notch'].value_counts(normalize=True).mul(100).sort_index(axis=0)).round(2)
df_abs_notch=df_abs_notch.cumsum()
df_abs_notch.abs_notch=['%.2f%%'% abs_notch for abs_notch in df_abs_notch.abs_notch ]
return df_abs_notch,df_notch
#引入数据集。以下三行替换为信评数据集。
iris = load_iris()
X_data = iris.data
y_data = iris.target
#训练集测试集模型划分。
X_train, X_test, y_train, y_test = train_test_split(X_data,y_data, test_size=.4, random_state=0)#这里是按照6:4对训练集测试集进行划分
model = LogisticRegressionCV(multi_class='multinomial')
#递归特征消除
rfe = RFECV(estimator=model, cv=StratifiedKFold(n_splits=3, random_state=1),scoring='neg_mean_squared_error')
rfe = rfe.fit(X_data, y_data)
# 和传参对应，所选择的属性的个数
print(rfe.n_features_)
# 打印的是相应位置上属性的排名
print(rfe.ranking_)
# 属性选择的一种模糊表示，选择的是true，未选择的是false
print(rfe.support_)
# 第1个属相的排名
print(rfe.ranking_[1])
# 外部估计函数的相关信息
print(rfe.estimator_)
best_parameters=rfe.estimator_.get_params()
print(best_parameters)
lr=LogisticRegressionCV(Cs=best_parameters['Cs'],
class_weight=best_parameters['class_weight'],
cv=best_parameters['cv'],
dual=best_parameters['dual'],
fit_intercept=best_parameters['fit_intercept'],
intercept_scaling=best_parameters['intercept_scaling'],
max_iter=best_parameters['max_iter'],
multi_class=best_parameters['multi_class'],
n_jobs=best_parameters['n_jobs'],
penalty=best_parameters['penalty'],
random_state=best_parameters['random_state'],
refit=best_parameters['refit'],
scoring=best_parameters['scoring'],
solver=best_parameters['solver'],
tol=best_parameters['tol'],
verbose=best_parameters['verbose'])
re = lr.fit(X_train, y_train)
y_test_pred = re.predict(X_test)
y_train_pred = re.predict(X_train)
score = 0
for y, y_pred in zip(y_test, y_test_pred):
score += 1 if y == y_pred else 0
print("accuracy of test_dataset :",score / len(y_test))
train_score = 0
for y, y_pred in zip(y_train, y_train_pred):
train_score += 1 if y == y_pred else 0
print("accuracy of train_dataset :",train_score / len(y_train))
test_df_abs_notch,test_df_notch=notch_diff_table(y_test,y_test_pred)
print('测试集极差矩阵:')
print(test_df_abs_notch)
print(test_df_notch)
train_df_abs_notch,train_df_notch=notch_diff_table(y_train,y_train_pred)
print('训练集集极差矩阵:')
print(train_df_abs_notch)
print(train_df_notch)
4
[1 1 1 1]
[ True True True True]
1
LogisticRegressionCV(Cs=10, class_weight=None, cv=None, dual=False,
fit_intercept=True, intercept_scaling=1.0, max_iter=100,
multi_class='multinomial', n_jobs=1, penalty='l2',
random_state=None, refit=True, scoring=None, solver='lbfgs',
tol=0.0001, verbose=0)
{'Cs': 10, 'class_weight': None, 'cv': None, 'dual': False, 'fit_intercept': True, 'intercept_scaling': 1.0, 'max_iter': 100, 'multi_class': 'multinomial', 'n_jobs': 1, 'penalty': 'l2', 'random_state': None, 'refit': True, 'scoring': None, 'solver': 'lbfgs', 'tol': 0.0001, 'verbose': 0}
accuracy of test_dataset : 0.9166666666666666
accuracy of train_dataset : 0.9888888888888889
测试集极差矩阵:
abs_notch
0.0 91.67%
1.0 100.00%
notch
-1 1.67%
0 91.67%
1 6.67%
训练集集极差矩阵:
abs_notch
0.0 98.89%
1.0 100.00%
notch
-1 1.11%
0 98.89%
遗传算法（genetic algorithms）
遗传算法（ Genetic Algorithm）是一种启发式的智能优化算法，其是模拟达尔文生物进化论中的自然选择和遗传学机理，实现了一种通过模拟自然进化的过程搜索最优解的方法。
因为是模拟达尔文自然进化的过程，所以其主要特点是直接对结构对象进行建模操作，不存在求导和函数连续性的限定；且采用概率化的寻优方法，不需要确定的规则就能自动获取和指导优化的搜索空间，自适应地调整搜索方向。
可以看到，遗传算法的优点很多，但是因为遗传算法是一种启发式的优化算法，所以他有点像我们初中学习几何问题时使用的辅助线，极具艺术性，在编码、适应性函数、交叉和变异时，需要些灵感性的构造，使用的好了非常简单高效的解决问题，使用不好根本解决不了问题。
遗传算法（GA）数学理论
遗传算法虽然是模拟的自然选择，但是其数学理论很简单，主要包括两部分：模式定理和积木块假设（building block hypochesis）。
模式定理：其给出了遗传算法在优化问题和搜索问题里面的数学有效性。但是此定理并没有在数学上给出其收敛性。
积木块假设：在给定的前提下，可以假设满足特定条件的遗传算法问题可以收敛到全局最优解。注意这是个假设。
模式定理
在介绍模式定理前，我们先介绍编码和几个概念。
1、编码：在使用遗传算法时，我们首先要对问题中的特征进行01编码。把特征特征映射成一个个只有01的字符串，这里的映射方式不唯一，含有灵感的艺术成分，有的人选取的01映射好，就可以很简单的解决问题，有的人选取的不恰当，就根本解决不了问题。具体映射要看具体问题。每个特征映射成一个只含有01的字符串后，将每个样本点的这些字符串拼接起来，就成了一个很长的字符串，这样原始的样本空间就映射成为了01编码空间里的点。这些01组成的小字符串就是基因，一个样本的长串就是染色体。如0001010111，这个就是一个染色体。
有了01编码后，我们就可以定义模式。
2、模式：我们将种群中的个体即染色体字符串中的相似样板称为“模式”，模式表示染色体串中某些特征位相同的结构。
如
3、模式阶： 模式H中确定位置的个数称为模式H的模式阶，记作O(H)。如
4、定义矩：模式中第一个确定位置和最后一个确定位置之间的距离，记作
5、模式匹配：如果一个模式
有了上面的定义，我们就可以引出模式定理的内容。如下：
设
其中:
证明：
对于模式H，在第t步中与模式H匹配的个体数量为：
这些个体在进行交叉和变异前被选作父体的数学期望为：
所以，若不进行交叉和变异，有：
考虑交叉，在交叉时，可以证明模式H被破坏的概率为：
所以模式保留的概率为：
所以上面t+1步时不考虑变异，有：
考虑变异，在交叉的基础上，模式H某个位置上不被变异的概率为：
所以模式H不被变异影响的概率为：
由泰勒展开得到：
所以将交叉和变异考虑在一起就有：
所以去掉右侧展开式中正的混合项，哟
当模式H具有低模式阶、短定义距以及平均适应度高于种群平均适应度。此时，有：
当模式阶很低，定义距很短，平均适应度高于种群平均适应度时，会有
所以
所以在子代中呈指数级增长。
所以模式定理得到证明。
可以看出，模式定理值证明了遗传算法的有效性，并没有给出算法最终会收敛到最优值。所以这里有个积木块假设。
积木块假设(building block hypothesis)： 遗传算法中，分别具有高适应值、短定义矩或低阶的模式（积木块）在遗传算子的作用下相互结合，形式同时具有高适应值、短定义矩和低阶的更优秀的模式，最终接近全局最优解。
积木块假设是基于以下两个基本前提的：
1、 表现型相近的个体具有相近的基因型；
2、 遗传算子间相对独立，相关性低。
其中：表现型的意思就是在没进行01编码映射前的样子，基因型就是编码后的样子。
有了上面的理论基础，我们给出遗传算法的流程如下。
遗传算法的流程图和流程如下：
1、对原始问题编码
原始问题描述的一般是一些现实问题的数学模型，这些模型大部分是一般的最优化模型，我们要对这个问题的可行域里面的数据进行01编码，做一个映射，将可行域里面的特征数据映射成含有01的字符串，具体怎么映射是不固定的，这里很具有艺术和灵感的成分，映射的好坏直接决定是否能解决问题。
2、产生初始种群
一般要解决的优化问题可行域都是连续的，所以，对于连续的可行域我们不可能进行遍历，只能先进行离散化，这里离散化的方法很多，一般采用随机的方式，也就是在连续的区间随机产生数据，这就是随机产生初始种群，当然这里也可以不随机。需要注意的是，这里产生的初始种群是在原始问题里面产生的，需要将其映射到编码空间。
3、计算适应性函数值
对于初始种群，我们要计算这些种群中所含染色体模式H的适应性函数值，这样用来判断，那些个体保留，那些个体被淘汰，所以这里的适应值函数是和目标函数挂钩的，但是，两者有一个区别，适应值函数一般都是正的，目标函数不一定为正。还有要注意的计算适应值函数值时，是在原始问题空间里面计算，不是在编码空间计算。
4、选择
设置一个适应值阈值，选择大于这个域字的个体作为接下来遗传的父体，淘汰小于这个域值的个体。
5、遗传运算
这里包括两部分，交叉和变异，交叉就是由上一步选择得到的父体，随机两两间进行k点交叉，这一步类似染色体的交叉遗传，得到新的个体染色体，对这个新的个体染色体，随机的对01字符串上的某一个节点进行变化，就是变异。这样就完成了遗传运算。
6、更新种群
进行完上面遗传运算后，原始的个体淘汰了一部分，且剩下的个体变成父体进行遗传运算，变成了新的种群。此时判读是否满足停止条件，若满足则停止，如不满足则进入第3步。
7、判断新种群是否满足停止条件
上面给出了整个遗传算法的步骤，但是遗传算法不能判断是否到达最优点附件，也就是说遗传算法不能判断是否收敛，那怎么才能停止运算呢，一般满足下面两个条件之一即可停止运算。
1、循环次数达到指定次数。
2、最大适应度值和平均适应度值变化不大、趋于稳定。
以上是对遗传算法的介绍，下面是遗传算法应用在特征选择上。
遗传算法进行特征选择
遗传算法特征选择的基本原理是用遗传算法寻找一个最优的二进制编码, 码中的每一位对应一个特征, 若第i位为“1”, 则表明对应特征被选取, 该特征将出现在估计器中, 为“0”, 则表明对应特征未被选取,该特征将不出现在分类器中。其基本步骤为:
1、 编码。采用二进制编码方法, 二进制码的每一位的值, “0”表示特征未被选中;“1”表示特征被选中。
2、初始群体的生成。随机产生N个初始串构成初始种群, 通常种群数确定为50～100。
3、 适应度函数。适应度函数表明个体或解的优劣性。针对特征选取问题, 适应度函数的构造非常重要, 它主要依据类别的可分性判据以及特征的分类能力。
4、 将适应度最大的个体, 即种群中最好的个体无条件地复制到下一代新种群中, 然后对父代种群进行选择、交叉和变异等遗传算子运算, 从而繁殖出下一代新种群其它n-1个基因串。通常采用转轮法作为选取方法, 适应度大的基因串选择的机会大, 从而被遗传到下一代的机会大, 相反, 适应度小的基因串选择的机会小, 从而被淘汰的机率大。交叉和变异是产生新个体的遗传算子, 交叉率太大, 将使高适应度的基因串结构很快被破坏掉, 太小则使搜索停止不前, 一般取为0.5～0.9。变异率太大, 将使遗传算法变为随机搜索, 太小则不会产生新个体, 一般取为0.01～0.1。
5、 如果达到设定的繁衍代数, 返回最好的基因串, 并将其作为特征选取的依据, 算法结束。否则, 回到4继续下一代的繁衍。
下面为具体的代码实现：
# -*- coding: utf-8 -*-
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris ##数据集
iris = load_iris()
x, y = iris.data, iris.target
data=pd.concat([pd.DataFrame(x),pd.DataFrame(y)],axis=1)
iris2 = np.mat(data)
X1 = iris2[0:50,0:4] #前50行为第一类别
X2 = iris2[50:100,0:4] #50到100行为第二类别
X3 = iris2[100:150,0:4] #100到150行为第三类别
m1 = np.mean(X1,axis = 0) # m1为第一类均值向量
m2 = np.mean(X2,axis = 0) # m2为第二类均值向量
m3 = np.mean(X3,axis = 0) # m3为第三类均值向量
m = np.mean(iris2,axis = 0) # m为总体均值向量
pc = 0.02 # pc为变异的概率
t = 100 #遗传算法迭代的次数
n = 50 #种群的个体数,要求大于20以保证具有随机性
#遗传算法
def GA(d):
population = np.zeros((n,4)) # 初始化种群
for i in range(n): # 定义种群的个体数为 n
a = np.zeros(4-d)
b = np.ones(d) # 将选择的d维特征定义为个体c中的1
c = np.append(a,b)
c = (np.random.permutation(c.T)).T # 随机生成一个d维的个体
population[i] = c # 初代的种群为 population，共有n个个体
# 遗传算法的迭代次数为t
fitness_change = np.zeros(t)
for i in range(t):
fitness = np.zeros(n) # fitness为每一个个体的适应度值
for j in range(n):
fitness[j] = Jd(population[j]) # 计算每一个体的适应度值
population = selection(population,fitness) # 通过概率选择产生新一代的种群
population = crossover(population) # 通过交叉产生新的个体
population = mutation(population) # 通过变异产生新个体
fitness_change[i] = max(fitness) #找出每一代的适应度最大的染色体的适应度值
# 随着迭代的进行，每个个体的适应度值应该会不断增加，所以总的适应度值fitness求平均应该会变大
best_fitness = max(fitness)
best_people = population[fitness.argmax()]
return best_people,best_fitness,fitness_change,population
#轮盘赌选择
def selection(population,fitness):
fitness_sum = np.zeros(n)
for i in range(n):
if i==0:
fitness_sum[i] = fitness[i]
else:
fitness_sum[i] = fitness[i] + fitness_sum[i-1]
for i in range(n):
fitness_sum[i] = fitness_sum[i] / sum(fitness)
#选择新的种群
population_new = np.zeros((n,4))
for i in range(n):
rand = np.random.uniform(0,1)
for j in range(n):
if j==0:
if rand<=fitness_sum[j]:
population_new[i] = population[j]
else:
if fitness_sum[j-1]<rand and rand<=fitness_sum[j]:
population_new[i] = population[j]
return population_new
#交叉操作
def crossover(population):
father = population[0:10,:]
mother = population[10:,:]
np.random.shuffle(father) # 将父代个体按行打乱以随机配对
np.random.shuffle(mother)
for i in range(10):
father_1 = father[i]
mother_1 = mother[i]
one_zero = []
zero_one = []
for j in range(4):
if father_1[j]==1 and mother_1[j]==0:
one_zero.append(j)
if father_1[j]==0 and mother_1[j]==1:
zero_one.append(j)
length1 = len(one_zero)
length2 = len(zero_one)
length = max(length1,length2)
half_length = int(length/2) #half_length为交叉的位数
for k in range(half_length): #进行交叉操作
p = one_zero[k]
q = zero_one[k]
father_1[p]=0
mother_1[p]=1
father_1[q]=1
mother_1[q]=0
father[i] = father_1 #将交叉后的个体替换原来的个体
mother[i] = mother_1
population = np.append(father,mother,axis=0)
return population
#变异操作
def mutation(population):
for i in range(n):
c = np.random.uniform(0,1)
if c<=pc:
mutation_s = population[i]
zero = [] # zero存的是变异个体中第几个数为0
one = [] # one存的是变异个体中第几个数为1
for j in range(4):
if mutation_s[j]==0:
zero.append(j)
else:
one.append(j)
a = np.random.randint(0,len(zero)) # e是随机选择由0变为1的位置
b = np.random.randint(0,len(one)) # f是随机选择由1变为0的位置
e = zero[a]
f = one[b]
mutation_s[e] = 1
mutation_s[f] = 0
population[i] = mutation_s
return population
#个体适应度函数 Jd(x)，x是d维特征向量(1*4维的行向量,1表示选择该特征)
def Jd(x):
#从特征向量x中提取出相应的特征
Feature = np.zeros(d) #数组Feature用来存 x选择的是哪d个特征
k = 0
for i in range(4):
if x[i] == 1:
Feature[k] = i
k+=1
#将4个特征从iris2数据集中取出重组成一个150*d的矩阵iris3
iris3 = np.zeros((150,1))
for i in range(d):
p = Feature[i]
p = p.astype(int)
q = iris2[:,p]
q = q.reshape(150,1)
iris3 = np.append(iris3,q,axis=1)
iris3 = np.delete(iris3,0,axis=1)
#求类间离散度矩阵Sb
iris3_1 = iris3[0:50,:] #iris数据集分为三类
iris3_2 = iris3[50:100,:]
iris3_3 = iris3[100:150,:]
m = np.mean(iris3,axis=0) #总体均值向量
m1 = np.mean(iris3_1,axis=0) #第一类的均值向量
m2 = np.mean(iris3_2,axis=0) #第二类的均值向量
m3 = np.mean(iris3_3,axis=0) #第二类的均值向量
m = m.reshape(d,1) #将均值向量转换为列向量以便于计算
m1 = m1.reshape(d,1)
m2 = m2.reshape(d,1)
m3 = m3.reshape(d,1)
Sb = ((m1 - m).dot((m1 - m).T) + (m2 - m).dot((m2 - m).T) + (m3 - m).dot((m3 - m).T))/3 #除以类别个数
#求类内离散度矩阵Sw
S1 = np.zeros((d,d))
S2 = np.zeros((d,d))
S3 = np.zeros((d,d))
for i in range(50):
S1 += (iris3_1[i].reshape(d,1)-m1).dot((iris3_1[i].reshape(d,1)-m1).T)
S1 = S1/50
for i in range(50):
S2 += (iris3_2[i].reshape(d,1)-m2).dot((iris3_2[i].reshape(d,1)-m2).T)
S2 = S2/50
for i in range(50):
S3 += (iris3_3[i].reshape(d,1)-m3).dot((iris3_3[i].reshape(d,1)-m3).T)
S3 = S3/50
Sw = (S1 + S2 + S3)/3
#计算个体适应度函数 Jd(x)
J1 = np.trace(Sb)
J2 = np.trace(Sw)
Jd = J1/J2
return Jd
if __name__ == '__main__':
# best_d = np.zeros(d) # judge存的是每一个维数的最优适应度
# fitness_change是遗传算法在迭代过程中适应度变化
# best是每一维数迭代到最后的最优的适应度，用于比较
for d in range(1,4):
print("\n")
best_people,best_fitness,fitness_change,best_population = GA(d)
choice = np.zeros(d)
k = 0
print("在取%d维的时候，通过遗传算法得出的最优适应度值为：%.6f"%(d,best_fitness))
print("选出的最优染色体为：")
print(best_people)
for j in range(4):
if best_people[j] == 1:
choice[k]=j+1
k+=1
print("选出的最优特征为：")
print(choice)
#画图
x = np.arange(0,t,1)
plt.xlabel('dimension')
plt.ylabel('fitness')
plt.ylim((0,50)) # y坐标的范围
plt.plot(x,fitness_change,'b')
在取1维的时候，通过遗传算法得出的最优适应度值为：16.041283
选出的最优染色体为：
[0. 0. 1. 0.]
选出的最优特征为：
[3.]
在取2维的时候，通过遗传算法得出的最优适应度值为：15.488503
选出的最优染色体为：
[0. 0. 1. 1.]
选出的最优特征为：
[3. 4.]
在取3维的时候，通过遗传算法得出的最优适应度值为：10.474305
选出的最优染色体为：
[0. 1. 1. 1.]
选出的最优特征为：
[2. 3. 4.]
2.3 拉斯维加斯方法（Las Vegas Wrapper）
拉斯维加斯算法和蒙特卡洛算法都是以赌场命名的随机算法，不同的点是，拉斯维加斯算法在给定的步数内，可能找到解，也可能找不到解，但蒙特卡洛算法在给定的步数内，一定会给出解，只不过有时解很精确，有时解很粗糙。若无步数限制，则两者都能给出精确解。
拉斯维加斯方法进行特征选择的步骤是，随机生成特征子集，对特征子集进行模型训练（通常会结合交叉验证），求的此时分类器的误差，比较误差和上一步误差的大小，若误差比上一步误差小，或者误差一样时，所选的特征数量比上一步小，则保留当前特征子集为备选，舍弃上一步特征子集，直到指定的步数结束。
# -*- coding: utf-8 -*-
"""
Created on 2020-09-29
@author:dingweijie
在建立SVM模型时以拉斯维加斯方法剔除特征
"""
from sklearn import svm
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import sklearn
import random
from sklearn.model_selection import train_test_split
from sklearn.datasets import load_iris ##数据集
from sklearn import cross_validation
from sklearn.model_selection import cross_val_score
iris = load_iris()
x, y = iris.data, iris.target
train_data,test_data,train_label,test_label =cross_validation.train_test_split(x,y, random_state=1, train_size=0.6,test_size=0.4) #sklearn.model_selection.
T=5000 #拉斯维加斯算法停止步数
t=0 #算法启始
E_scores=0 #启始分数为零
d=4#数据集的特征总维度
while t<T:
n_random=random.sample(range(1,4),1)
index = random.sample(range(0,4),n_random[0])
index.sort()
tmp_d=n_random
x_new=x[:,index]
classifier=svm.SVC(C=2,kernel='rbf',gamma=10,decision_function_shape='ovo') # ovr:一对多策略
scores = cross_val_score(classifier,x_new,y) #cv：默认是3折交叉验证，可以修改cv=5，变成5折交叉验证。
tmp_E_scores=scores.mean()
if tmp_E_scores>E_scores or (tmp_E_scores==E_scores and tmp_d < d ):
t=0
E_scores=tmp_E_scores
d=tmp_d
else:
t=t+1
print("特征所在列index：",index)
print("最优分数：",E_scores)
特征所在列index： [0, 3]
最优分数： 0.9869281045751634
Embedded（嵌入式选择）
嵌入式选择特征与前面的过滤式和包裹式不一样，没有特征选择和模型训练的明显界限，特征选择和模型训练同时进行且同时完成。
这种方式一般包括L1正则法和使用决策树的树模型。下面分别介绍。
L1正则法（LASSO ）
在建立模型时，我们一般会加上正则项来防止模型的过拟合问题，一般正则项时参数的
对于数据集
上面就是最简单的线性回归模型。为防止过拟合出现，我们引入正则项，有：
当p=2时，正则项就是
这里我们统一对岭回归和拉索回归进行分析。
引入正则化的优化目标：
可以看成是以下带约束条件的拉格朗日函数。
所以，其实就是在区域
为分析问题，我们考虑w在只有两维的情况下，因为两维可以在图片中画出来。
设
如下图：
由上图可以看出，使用
下面是
# -*- coding: utf-8 -*-
import pandas as pd
import numpy as np
from sklearn.cross_validation import train_test_split
from sklearn.preprocessing import MinMaxScaler
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
import matplotlib.pyplot as plt
#导入数据
data = load_iris()#读取iris数据集
X = data['data'][:100]#因为是二分类问题，所以我们只需要前一百个样本，每类50个样本
y = data['target'][:100].reshape((-1, 1))#读取他们对应的类别标签
X_train,X_test,y_train,y_test=train_test_split(X,y,test_size=0.3,random_state=0)
#特征值缩放
#归一化
mms=MinMaxScaler()
X_train_norm=mms.fit_transform(X_train)
X_test_norm=mms.fit_transform(X_test)
#标准化
stdsc=StandardScaler()
X_train_std=stdsc.fit_transform(X_train)
X_test_std=stdsc.fit_transform(X_test)
#L1正则化的逻辑斯蒂模型
lr=LogisticRegression(penalty='l1',C=0.1)#penalty='l2'
lr.fit(X_train_std,y_train)
print ('Training accuracy:',lr.score(X_train_std, y_train))
print ('Test accuracy:',lr.score(X_test_std, y_test))#比较训练集和测试集，观察是否出现过拟合
print (lr.intercept_)#查看截距，三个类别
print (lr.coef_)#查看权重系数，L1有稀疏化效果做特征选择
Training accuracy: 1.0
Test accuracy: 1.0
[0.]
[[0. 0. 1.41249054 0.33710232]]
决策树、随机森林、极限树
决策树和随机森林的基础知识和特征选择的方法我在《Machine Learning-决策树与集成学习》文章中已经介绍。这里只整理极限树的知识。
极限树（Extratree）
极限树和随机森林有很多相似的地方，算是随机森林的一种变种，修改随机森林的两个地方得到极限树：
1、样本不使用bagging，不进行随机抽取，而是使用全部样本。
2、在每一棵树的每一个节点进行分裂时，不进行任何数据运算，直接随机选取一个特征，然后随机给出这个特征的分裂阈值。
用极限树进行特征选择的python例子如下：
"""
Created on 2020-09-12
@author:dingweijie529
极限树模型筛选特征
"""
from sklearn import metrics
from sklearn.ensemble import ExtraTreesClassifier
import matplotlib.pyplot as plt
#引入数据集。以下三行替换为信评数据集。
iris = load_iris()
X_data = iris.data
y_data = iris.target
#训练集测试集模型划分。
model = ExtraTreesClassifier()
model.fit(X_data, y_data)
name_list = iris.feature_names
num_list = model.feature_importances_
fig, ax = plt.subplots(facecolor=(1, 1, 1),figsize=(16,8))
rects=plt.bar(range(len(num_list)), num_list)
# X轴标题
index=[0,1,2,3]
index=[float(c)+0.4 for c in index]
plt.ylim(ymax=1, ymin=0)
plt.xticks(index, name_list)
plt.ylabel("importances") #X轴标签
for rect in rects:
height = rect.get_height()
plt.text(rect.get_x() + rect.get_width() / 2, height, str(height), ha='center', va='bottom')
plt.show()
