特征选择方法 (MIC互信息、DCOR距离相关系数、KPCA、随机森林贡献等)，含python代码

1.特征选择方法简介

在很多机器学习任务中，特征选择是一个至关重要的步骤，合适的特征选取不仅能够减少模型的复杂度与算力，同时还可以在一定程度上提高后续模型的性能，本文以二手船造价预测为具体事例，探讨如何通过更高效的特征选择法识别和保留最具影响力的特征，从而构建更高效、更准确的预测模型。

2.二手船造价预测任务描述

任务与相关数据介绍

在造船领域，对二手船的价格预测是一个较为重要的任务。其价格受到多重因素的影响，本文共列举了其中的11种因素：62%品位铁矿石期货价格、铜期货价格、锡期货价格、大豆期货价格、小麦期货价格、玉米期货价格、木材期货价格、WTI原油期货价格、美元指数、汇率Yuan/$、标准普尔500指数等，数据集如下图表格所示，共2011年一直收集到2023年，每隔一周收集一次数据：

本文的任务是，利用已知的11种因素，预测一周/一月之后的二手船价格。

3.特征选择方法介绍

本节介绍现有的特征选择法
我看了很多包括CSDN、知乎以及博客园等发现，当下介绍最多的特征选择方法主要分为三大类：过滤式方法、包裹式方法和嵌入式方法即：

对于过滤式方法：皮尔逊相关系数、卡方检验、ANOVA等主要用于线性数据多一些，对于复杂一点的非线性数据其效果并不是很好，因此本文对于这几种方法不予以介绍。只介绍能处理非线性数据的互信息方法(MIC)与距离相关系数法(DCOR)
对于包裹式的方法：其主要策略是首先建议一个基准模型(如SVM或LASSO回归等较为简单的模型)，用于评价特征选择的好坏。之后通过优化算法(如遗传算法)不断选择特征，找到最优的特征组合。这种方法虽然效果还行，但是可解释性太差了，比如你通过这种方法得到了一个最佳的特征组合如：选择小麦期货价格、玉米期货价格、木材期货价格、WTI原油期货价格、美元指数这几个特征效果最好。但是如果你导师问你，为啥他们几个好呢？那你只能弱弱的来一句：通过优化算法试出来的。对此本文不再介绍这类方法，感兴趣的同学可以看一下《基于传感器阵列的多组分 VOCs 检测及工业应用验证研究》这篇博士论文，有较为详细的介绍。

• (划重点)本文主要介绍亲测有效的：KPCA、互信息法、距离距离相关系数法、随机森林贡献度法。

4.相关算法与代码

本节介绍KPCA、MIC、DCOR和随机森林贡献度法的原理与代码

4.1 KPCA算法

首先是大家最熟悉的老当益壮KPCA法，核主成分分析(KPCA)是一种用于非线性数据降维的技术，它是经典PCA(主成分分析)的推广，能够处理数据集中非线性结构的情况。本文主要使用KPCA中的解释方差比概念来计算每种特则的相对重要程度。

(1). 核心思想

• 非线性映射: KPCA首先通过非线性映射$\varphi$将数据从原始空间$X$映射到一个高维特征空间$\mathcal{F}$，在这个空间中，数据可能表现出更好的线性可分性。
• 核技巧: 直接在高维空间中操作可能是不切实际的，因此KPCA使用核技巧，即通过核函数$k(x_i,x_j)=⟨\varphi (x_i),\varphi (x_j)⟩$来间接进行高维空间的内积计算。

(2).主要步骤

• 非线性映射：$x\to \varphi (x)$，将原始数据映射到高维特征空间

• 核函数选择：选择RBF（高斯）核函数：$k(x_i,x_j)=\exp (-\gamma \Vert x_i-x_j{\Vert }^2)$

• 数据中心化：

  • 在高维特征空间中对数据进行中心化处理
  • 通过调整核矩阵间接完成：$\tilde{K}=K-1_{N}K-K{1}_N+1_{N}K{1}_N$

• 构建核矩阵：计算核矩阵 $K$，其中 $K_{ij}=k(x_i,x_j)$

• 特征值分解：

  • 对中心化后的核矩阵 $\tilde{K}$ 进行特征值分解
  • 得到特征值 ${\lambda }_1\ge {\lambda }_2\ge ...\ge {\lambda }_n$ 和对应的特征向量

• 选择主成分：选取与最大特征值对应的特征向量作为主成分

• 投影数据：将原始数据投影到选定的主成分上，得到降维后的数据 $X_{kpca}$

• 计算方差解释比：

  • 对每个主成分 $i$：
    a. 创建仅包含第 $i$ 个主成分的零矩阵
    b. 将该主成分通过逆变换映射回原始空间：$X_{inverse}$
    c. 将 $X_{inverse}$ 重新投影到KPCA空间：$X_{proj}$
    d. 计算方差比：$\frac{\text{Var}(X_{proj,i})}{\text{Var}(X_{kpca,i})}$
  • 得到每个主成分的方差解释比列表

• 计算累积方差解释比：对方差解释比列表进行累积求和

(3).核心代码
首秀需要安装一下sklearn库，一般通过pip指令就可以安装：

# 对于非Anaconda用户
pip install -U scikit-learn

# 对于Anaconda用户
conda install scikit-learn

from sklearn.decomposition import KernelPCA
from sklearn.preprocessing import StandardScaler

# 标准化数据
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# 应用KPCA
n_components = X.shape[1]  # 使用所有可能的主成分
kpca = KernelPCA(n_components=n_components, kernel='rbf', fit_inverse_transform=True, alpha=0.1)
X_kpca = kpca.fit_transform(X_scaled)

# 计算每个主成分解释的方差比例
explained_variance_ratio = []
for i in range(n_components):
    X_inverse = kpca.inverse_transform(np.zeros((X.shape[0], n_components)))
    X_inverse[:, i] = X_kpca[:, i]
    X_proj = kpca.transform(X_inverse)
    explained_variance_ratio.append(np.var(X_proj[:, i]) / np.var(X_kpca[:, i]))

# 计算累积解释方差比例
cumulative_variance_ratio = np.cumsum(explained_variance_ratio)

(3).结果测试

4.2 基于MIC的互信息法

这个方法是经过我多次测试后最好使的方法。首先介绍一下互信息的概念，互信息（Mutual Information, MI）是衡量两个随机变量间依赖性的统计量，它描述了知道其中一个变量能减少多少关于另一个变量的不确定性，其公式在后面给出。在统计学和信息论中，互信息的计算通常是基于概率分布的，但直接从有限样本估计概率分布可能带来偏差，尤其是在高维度空间。Maximal Information Coefficient (MIC) 是一种旨在从有限样本中估计两个变量间互信息的非参数方法，它可以捕捉到线性和非线性关系。

(1). 核心思想

• 网格划分探索

  • 将数据点划分到不同大小和位置的网格中
  • 探索各种可能的网格划分方式

• 互信息最大化

  • 对每种网格划分，计算互信息
  • 寻找能使互信息最大化的网格划分方式

(2).主要步骤

• 数据离散化：首先，将连续变量划分为若干个区间，这一过程被称为离散化，目的是将连续的数据转化为可以计算互信息的形式。

• 互信息计算：对于每个可能的网格划分，计算互信息 $I(X;Y)$，这是衡量两个变量共享信息量的关键指标。互信息公式如下：
  $$I(X;Y)=\sum _{x\in X}\sum _{y\in Y}p(x,y)log\frac{p(x,y)}{p(x)p(y)}$$

  其中，$p(x,y)$为联合概率分布，而$p(x)$和$p(y)$分别为它们的边缘分布。

• 标准化处理：为了使 MIC 值具有可比性，需要将互信息值标准化。这一步骤是通过对互信息除以 log ⁡ ( m i n { n x , n y } ) \log(min\{n_x,n_y\}) log(min{nx​,ny​})实现的，其中$n_x$和$n_y$分别是$X$和$Y$的区间数。

• 最大化过程：在所有可能的网格划分中，选择使得标准化后的互信息值最大的划分，即为 MIC 值。MIC 的数学定义如下：
  M I C ( X , Y ) = max ⁡ n x ⋅ n y < B { I ∗ ( X , Y ) log ⁡ ( m i n { n x , n y } ) } MIC(X,Y) = \max_{n_x \cdot n_y < B} \left\{ \frac{I^*(X,Y)}{\log(min\{n_x,n_y\})} \right\} MIC(X,Y)=nx​⋅ny​<Bmax​{log(min{nx​,ny​})I∗(X,Y)​}

  其中，$B$是网格大小的上限，通常设定为数据点数量的函数，例如$B=n^{0.6}$，以平衡计算效率和精度。

(3).核心代码
首先需要安装一个库MINE，一般用pip指令就可以安装

# 对于非Anaconda用户
pip install minepy

# 对于Anaconda用户
conda install -c conda-forge minepy

通过MIC算法计算互信息的python实现：

# 导入MINE类，这是minepy库的一部分，专门用于计算MIC和其他相关统计量
from minepy import MINE

# 定义一个名为compute_mic的函数，接受两个参数：X（特征矩阵）和y（目标变量）
def compute_mic(X, y):
    # 初始化一个空列表，用于存储每一列与y之间的MIC得分
    mic_scores = []

    # 创建一个MINE实例，用于进行MIC的计算
    mine = MINE()

    # 遍历DataFrame X的所有列
    for column in X.columns:
        # 使用MINE的compute_score方法来计算当前列与目标变量y之间的MIC
        mine.compute_score(X[column], y)

        # 将计算得到的MIC得分追加到mic_scores列表中
        mic_scores.append(mine.mic())

    # 返回一个Series，其中索引为X的列名，值为相应的MIC得分
    return pd.Series(mic_scores, index=X.columns)

# 调用compute_mic函数，传入特征矩阵X和目标变量y，计算MIC得分
mic_scores = compute_mic(X, y)

(4).结果测试

4.3 距离相关系数 (Distance Correlation)

除了互信息外，另外一种比较好用的方法是距离相关系数(DCOR)。距离相关系数（Distance Correlation, 简称DCOR），由Gábor J. Székely等人在2007年提出，是一种革新性的统计量，用于衡量两个随机向量之间的依赖程度。与经典的Pearson相关系数相比，DCOR的最大优势在于它能够检测任何形式的依赖关系，而不仅仅是线性关系，从而提供了更为全面的依赖性评估。

(1). 核心思想

• 基于距离度量
  • 使用样本点间距离衡量变量关系
  • 可捕捉非线性依赖性
• 独立性检测
  • 仅当变量独立时为零
  • 克服Pearson相关系数局限
• 非参数特性
  • 不假设特定数据分布
  • 广泛适用于各类数据关系

(2).主要步骤

• 距离矩阵构建: 首先，对于两组样本 $X=\left(x_1,\dots ,x_n\right)和Y=\left(y_1,\dots ,y_n\right)$ ，分别计算它们之间的欧比距离矩阵：$$a_{ij}=\Vert x_i-x_j\Vert ,b_{ij}=\Vert y_i-y_j\Vert$$.
• 双重中心化：对距离矩阵进行双重中心化处理，即从每个元素中减去其所在行和列的平均值，再加回整个矩阵的平均值。这一过程旨在消除距离矩阵的偏移效应，公式如下:
  $$A_{ij}=a_{ij}-ai.-a.j+a..,Bij=bij-bi-bj+b$$.
  其中， $a_i,a_j,a_{ij}$分别表示 a_{ij} 的行平均、列平均和总平均
• 距离协方差与方差计算：基于双重中心化的距离矩阵，进一步计筫距离协方差和距离方差。距离协方差是衡量两个向量间距离相似性的关键指标，而距离方差则是单个向量内部距离变化的量化。计算公式如下:
  $$d{\text{Cov}}^2(X,Y)=\frac{1}{n^2}\sum _{i,j=1}^n{A}_{ij}{B}_{ij},d{\text{Var}}^2(X)=\frac{1}{n^2}\sum _{i,j=1}^n{A}_{ij}^2,d{\text{Var}}^2(Y)=\frac{1}{n^2}\sum _{i,j=1}^n{B}_{ij}^2$$
• 距离相关系数：最后，根据距离协方差和距离方差计算距离相关系数 dCor\funcapply(X,Y) ，它是衡量两个随机向量依赖性的标准化指标，取值范围为 [0,1] ，其中 0 表示完全独立， 1 表示存在确定性关系：
  $$d\text{Cor}(X,Y)=\frac{d\text{Cov}(X,Y)}{\sqrt{d\text{Var}(X)d\text{Var}(Y)}}$$

(3).核心代码

dcor库安装：

pip install dcor

# 距离相关系数法 (DCOR)
import pandas as pd
from dcor import distance_correlation

def compute_dcor(X, y):
    dcor_scores = []
    for column in X.columns:
        dcor_scores.append(distance_correlation(X[column], y))
    return pd.Series(dcor_scores, index=X.columns)

dcor_scores = compute_dcor(X, y)

(4).结果测试

4.4 随机森林贡献度法

随机森林模型的原理在很多博客里面都有讲到，我在这里就不在赘述，可以参考一下几个链接：
机器学习算法系列（十八）-随机森林算法（Random Forest Algorithm）
机器学习5—分类算法之随机森林（Random Forest）
我主要讲一下怎么利用该算法去计算不同变量的贡献度大小：

(1). 随机森林核心思想

• 随机森林是一种集成学习（Ensemble Learning）方法，其核心在于结合多个弱分类器（在随机森林中通常是决策树）的预测结果，以提高整体模型的准确性和稳定性。
• 在随机森林中，每棵树在构建时不仅使用数据的不同子集，而且在每个节点的分裂过程中只考虑特征的随机子集。算法不是在所有特征中寻找最优分割，而是在一个小的随机特征集合中寻找。通过限制每个节点可以考虑的特征数量，同时也促进了特征之间可能存在的冗余和相关性的处理。

(2).贡献度计算主要步骤

• Gini不纯度
  Gini不纯度是分类任务中常用的不纯度指标，定义如下：
  $$Gini(D)=1-\sum _{k=1}^K{p}_k^2$$
  其中，$D$是数据集，$K$是类别数量，$p_k$是属于第$k$类的概率。

• 熵
  熵是另一个常用于决策树的不纯度指标，主要用于分类任务，定义如下：
  $$Entropy(D)=-\sum _{k=1}^K{p}_k{\log }_2(p_k)$$

• 特征重要度计算

  假设我们有特征$X_j$，在某个节点上使用$X_j$分割数据集(D)成两个子集$D_l$和$D_r$。那么，使用特征$X_j$分割前后的不纯度减少可以表示为：
  $$\Delta I(X_j,D)=I(D)-\left(\frac{|D_l|}{|D|}I(D_l)+\frac{|D_r|}{|D|}I(D_r)\right)$$
  其中，$I(\cdot )$可以是Gini不纯度或者熵。

• 特征重要度

  在随机森林中，对于每个特征$X_j$，我们可以计算所有树中使用该特征进行分割时的不纯度减少总和，然后对所有特征的总和进行归一化。特征$X_j$的重要度$Im{p}_j$可以表示为：
  $$Im{p}_j=\frac{{\sum }_{t=1}^T{\sum }_{\text{all splits on }{X}_j}\Delta I(X_j,D_t)}{{\sum }_{j^{\prime }=1}^p{\sum }_{t=1}^T{\sum }_{\text{all splits on }{X}_{j^{\prime }}}\Delta I(X_{j^{\prime }},D_t)}$$
  其中，$T$是树的数量，$p$是特征的数量，$D_t$是树$t$的数据集。

(3). 核心代码

from sklearn.ensemble import RandomForestRegressor

# 随机森林特征重要性
rf = RandomForestRegressor(n_estimators=100, random_state=42)
rf.fit(X, y)
rf_importance = pd.Series(rf.feature_importances_, index=X.columns)

(4). 结果测试

5.结论

以上就是使用不同的算法对二手船造价进行分析的结果，共包含KPCA、MIC、DCOR以及随机森林贡献度四种方法。从每种方法的结果图中我们可以看出，不同算法的结果并不是一样的，甚至数值和排序上差了很多，就比如随机森林中标准普尔500贡献度异常的大。那么我们该怎么去选择算法呢？这里我提一个建议，可以去多尝试几种算法，选择其中结果较为相似的两种或几种作为最终结果。就比如在本文中，MIC和DCOR两种方法的结果很接近，那我们就有理由任务两个算法的结果都是正确的(就像你考完试之后和别人对数学题最后一大题答案，如果有两人或多人算出来的数是一致的就基本说明你们都做对了，除非你们心有灵犀~)。

---