專欄 | 基於 Jupyter 的特徵工程手冊：特徵選擇（二）

红色石头發表於2020-04-24
特徵工程

資料預處理後，我們生成了大量的新變數（比如獨熱編碼生成了大量僅包含0或1的變數）。但實際上，部分新生成的變數可能是多餘：一方面它們本身不一定包含有用的資訊，故無法提高模型效能；另一方面過這些多餘變數在構建模型時會消耗大量記憶體和計算能力。因此，我們應該進行特徵選擇並選擇特徵子集進行建模。

專案地址：

https://github.com/YC-Coder-Chen/feature-engineering-handbook

本文將介紹特徵工程第一種演算法：Filter Methods 過濾法（下）。

目錄：

1.1.1.5 Mutual Information (regression problem) 互資訊 (迴歸問題)

互資訊（Mutual Information）衡量變數間的相互依賴性。其本質為熵差，即 ( )− ( | )，即知道另一個變數資訊後混亂的降低程度 。當且僅當兩個隨機變數獨立時MI等於零。MI值越高，兩變數之間的相關性則越強。與Pearson相關和F統計量相比，它還捕獲了非線性關係。

公式:

  • 若兩個變數均為離散變數:

p( , )( , ) 為x和y的聯合概率質量函式 (PMF)， p ( )則為x的聯合概率質量函式 (PMF)。

  • 若兩個變數均為連續變數:

p( , )( , ) 為x和y的聯合概率密度函式 (PDF)，p ( )則為x的概率密度函式 (PDF)。連續變數情形下，在實際操作中，往往先對資料離散化分桶，然後逐個桶進行計算。

但是實際上，一種極有可能的情況是，x和y中的一個可能是離散變數，而另一個是連續變數。因此在sklearn中，它基於[1]和[2]中提出的基於k最臨近演算法的熵估計非引數方法。

[1] A. Kraskov, H. Stogbauer and P. Grassberger, “Estimating mutual information”. Phys. Rev. E 69, 2004.

[2] B. C. Ross “Mutual Information between Discrete and Continuous Data Sets”. PLoS ONE 9(2), 2014.

import numpy as np
from sklearn.feature_selection import mutual_info_regression
from sklearn.feature_selection import SelectKBest

# 直接載入資料集
from sklearn.datasets import fetch_california_housing
dataset = fetch_california_housing()
X, y = dataset.data, dataset.target # 利用 california_housing 資料集來演示
# 此資料集中，X，y均為連續變數，故此滿足使用MI的條件

# 選擇前15000個觀測點作為訓練集
# 剩下的作為測試集
train_set = X[0:15000,:].astype(float)
test_set = X[15000:,].astype(float)
train_y = y[0:15000].astype(float)

# KNN中的臨近數是一個非常重要的引數
# 故我們重寫了一個新的MI計算方程更好的來控制這一引數
def udf_MI(X, y):
    result = mutual_info_regression(X, y, n_neighbors = 5) # 使用者可以輸入想要的臨近數
    return result

# SelectKBest 將會基於一個判別方程自動選擇得分高的變數
# 這裡的判別方程為F統計量
selector = SelectKBest(udf_MI, k=2) # k => 我們想要選擇的變數數
selector.fit(train_set, train_y) # 在訓練集上訓練
transformed_train = selector.transform(train_set) # 轉換訓練集
transformed_train.shape #(15000, 2), 其選擇了第一個及第八個變數
assert np.array_equal(transformed_train, train_set[:,[0,7]])

transformed_test = selector.transform(test_set) # 轉換測試集
assert np.array_equal(transformed_test, test_set[:,[0,7]]);
# 可見對於測試集，其依然選擇了第一個及第八個變數

# 驗算上述結果
for idx in range(train_set.shape[1]):
    score = mutual_info_regression(train_set[:,idx].reshape(-1,1), train_y, n_neighbors = 5)
    print(f"第{idx + 1}個變數與因變數的互資訊為{round(score[0],2)}")
# 故應選擇第一個及第八個變數

第1個變數與因變數的互資訊為0.37

第2個變數與因變數的互資訊為0.03

第3個變數與因變數的互資訊為0.1

第4個變數與因變數的互資訊為0.03

第5個變數與因變數的互資訊為0.02

第6個變數與因變數的互資訊為0.09

第7個變數與因變數的互資訊為0.37

第8個變數與因變數的互資訊為0.46

1.1.1.6 Chi-squared Statistics (classification problem) 卡方統計量 (分類問題)

卡方統計量主要用於衡量兩個類別特徵之間的相關性。sklearn提供了chi2方程用於計算卡方統計量。其輸入的特徵變數必須為布林值或頻率（故對於類別變數應考慮獨熱編碼）。卡方統計量的零假設為兩個變數是獨立的，因為卡方統計量值越高，則兩個類別變數的相關性越強。因此，我們應該選擇具有較高卡方統計量的特徵。

公式:

其中， , 為在變數X上具有i-th類別值且在變數Y上具有j-th類別值的實際觀測點計數。 , 為利用概率估計的應在在變數X上具有i-th類別值且在變數Y上具有j-th類別值的觀測點數量。n為總觀測數， 為在變數X上具有i-th類別值的概率， 為在變數Y上具有j-th類別值的概率。

值得注意的是，通過解析原始碼，我們發現在sklearn中利用chi2計算出來的卡方統計量並不是統計意義上的卡方統計量。當輸入變數為布林變數時，chi2計算值為該布林變數為True時候的卡方統計量（我們將會在下文舉例說明）。這樣的優勢是，獨熱編碼生成的所有布林值變數的chi2值之和將等於原始變數統計意義上的卡方統計量。

舉個簡單的例子，假設一個變數I有0，1，2兩種可能的值，則獨特編碼後一共會產生3個新的布林值變數。這三個布林值變數的chi2計算出來的值之和，將等於變數I與因變數直接計算得出的統計意義上的卡方統計量。

解析sklearn中chi2的計算

# 首先，隨機生成一個資料集
import pandas as pd
sample_dict = {'Type': ['J','J','J',
                        'B','B','B',
                        'C','C','C','C','C'], 
               'Output': [0, 1, 0, 
                          2, 0, 1,  
                          0, 0, 1, 2, 2,]}
sample_raw = pd.DataFrame(sample_dict)
sample_raw #原始資料，Output是我們的目標變數，Type為類別變數


# 下面利用獨熱編碼生成布林變數，並利用sklearn計算每一個布林變數的chi2值
sample = pd.get_dummies(sample_raw)
from sklearn.feature_selection import chi2
chi2(sample.values[:,[1,2,3]],sample.values[:,[0]])
# 第一行為每一個布林變數的chi2值

(array([0.17777778, 0.42666667, 1.15555556]),
array([0.91494723, 0.8078868 , 0.56114397]))

# 下面直接計算原始變數Type與output統計學意義上的卡方統計量
# 首先，先統計每一個類別下出現的觀測數，用於建立列聯表
obs_df = sample_raw.groupby(['Type','Output']).size().reset_index()
obs_df.columns = ['Type','Output','Count']
obs_df

即列聯表（contingency table）為：

from scipy.stats import chi2_contingency
obs = np.array([[1, 1, 1], [2, 1, 2],[2, 1, 0]])
chi2_contingency(obs) # 第一個值即為變數Type與output統計學意義上的卡方統計量

(1.7600000000000002,

0.779791873961373,

4,

array([[1.36363636, 0.81818182, 0.81818182], [2.27272727, 1.36363636, 1.36363636], [1.36363636, 0.81818182, 0.81818182]]))

# 而chi2方程算出來的布林值之和為即為原始變數的統計意義上的卡方統計量
chi2(sample.values[:,[1,2,3]],sample.values[:,[0]])[0].sum() == chi2_contingency(obs)[0]

True

# 那麼sklearn中的chi2是如何計算的呢？
# 不妨以第一個生成的布林值為例，即Type為B
# chi2出來的值為0.17777778
# 而這與利用scipy以下程式碼計算出的計算一致
from scipy.stats import chisquare
f_exp = np.array([5/11, 3/11, 3/11]) * 3 # 預期頻數為 output的先驗概率 * Type為B 的樣本數
chisquare([1,1,1], f_exp=f_exp) # [1,1,1] 即Type為B 的樣本實際頻數
# 即sklearn 中的chi2 僅考慮了Type為B情形下的列連表

Power_divergenceResult(statistic=0.17777777777777778,pvalue=0.9149472287300311)

如何利用sklearn 來進行特徵選擇

import numpy as np
from sklearn.feature_selection import chi2
from sklearn.feature_selection import SelectKBest

# 直接載入資料集
from sklearn.datasets import load_iris # 利用iris資料作為演示資料集
iris = load_iris()
X, y = iris.data, iris.target
# 此資料集中，X為連續變數，y為類別變數
# 不滿足chi2的使用條件

# 將連續變數變為布林值變數以滿足chi2使用條件
# 不妨利用其是否大於均值來生成布林值（僅作為演示用）
X = X > X.mean(0)

# iris 資料集使用前需要被打亂順序
np.random.seed(1234)
idx = np.random.permutation(len(X))
X = X[idx]
y = y[idx]

# 選擇前100個觀測點作為訓練集
# 剩下的作為測試集
train_set = X[0:100,:]
test_set = X[100:,]
train_y = y[0:100]

# sklearn 中直接提供了方程用於計算卡方統計量
# SelectKBest 將會基於一個判別方程自動選擇得分高的變數
# 這裡的判別方程為F統計量
selector = SelectKBest(chi2, k=2) # k => 我們想要選擇的變數數
selector.fit(train_set, train_y) # 在訓練集上訓練
transformed_train = selector.transform(train_set) # 轉換訓練集
transformed_train.shape #(100, 2), 其選擇了第三個及第四個變數 
assert np.array_equal(transformed_train, train_set[:,[2,3]])

transformed_test = selector.transform(test_set) # 轉換測試集
assert np.array_equal(transformed_test, test_set[:,[2,3]]);
# 可見對於測試集，其依然選擇了第三個及第四個變數

# 驗證上述結果
for idx in range(train_set.shape[1]):
    score, p_value = chi2(train_set[:,idx].reshape(-1,1), train_y)
    print(f"第{idx + 1}個變數與因變數的卡方統計量為{round(score[0],2)}，p值為{round(p_value[0],3)}")
# 故應選擇第三個及第四個變數

第1個變數與因變數的卡方統計量為29.69，p值為0.0

第2個變數與因變數的卡方統計量為19.42，p值為0.0

第3個變數與因變數的卡方統計量為31.97，p值為0.0

第4個變數與因變數的卡方統計量為31.71，p值為0.0

1.1.1.7 F-Score (classification problem) F-統計量 (分類問題)

在分類機器學習問題中，若變數特徵為類別特徵，則我們可以使用獨熱編碼配合上述chi2方法選擇最重要的特徵。但若特徵為連續變數，則我們可以使用ANOVA-F值。ANOVA F統計量的零假設是若按目標變數（類別）分組，則連續變數的總體均值是相同的。故我們應選擇具有高ANOVA-F統計量的連續變數，因為這些連續變數與目標變數的關聯性強。

公式：

其中，SS(between)為組間的平方和，即組均值和總體均值之間的平方和。SS(error)是組內的平方和，即資料與組均值之間的平方和。m是目標變數的總類別數，n是觀測數。

import numpy as np
from sklearn.feature_selection import f_classif
from sklearn.feature_selection import SelectKBest

# 直接載入資料集
from sklearn.datasets import load_iris # 利用iris資料作為演示資料集
iris = load_iris()
X, y = iris.data, iris.target
# 此資料集中，X為連續變數，y為類別變數
# 滿足ANOVA-F的使用條件

# iris 資料集使用前需要被打亂順序
np.random.seed(1234)
idx = np.random.permutation(len(X))
X = X[idx]
y = y[idx]

# 選擇前100個觀測點作為訓練集
# 剩下的作為測試集
train_set = X[0:100,:]
test_set = X[100:,]
train_y = y[0:100]

# sklearn 中直接提供了方程用於計算ANOVA-F
# SelectKBest 將會基於一個判別方程自動選擇得分高的變數
# 這裡的判別方程為F統計量
selector = SelectKBest(f_classif, k=2) # k => 我們想要選擇的變數數
selector.fit(train_set, train_y) # 在訓練集上訓練
transformed_train = selector.transform(train_set) # 轉換訓練集
transformed_train.shape #(100, 2), 其選擇了第三個及第四個變數 
assert np.array_equal(transformed_train, train_set[:,[2,3]])

transformed_test = selector.transform(test_set) # 轉換測試集
assert np.array_equal(transformed_test, test_set[:,[2,3]]);
# 可見對於測試集，其依然選擇了第三個及第四個變數

# 驗證上述結果
for idx in range(train_set.shape[1]):
    score, p_value = f_classif(train_set[:,idx].reshape(-1,1), train_y)
    print(f"第{idx + 1}個變數與因變數的ANOVA-F統計量為{round(score[0],2)}，p值為{round(p_value[0],3)}")
# 故應選擇第三個及第四個變數

第1個變數與因變數的ANOVA-F統計量為91.39，p值為0.0

第2個變數與因變數的ANOVA-F統計量為33.18，p值為0.0

第3個變數與因變數的ANOVA-F統計量為733.94，p值為0.0

第4個變數與因變數的ANOVA-F統計量為608.95，p值為0.0

1.1.1.7 Mutual Information (classification problem) 互資訊 (分類問題)

【與1.1.1.5一樣】互資訊（Mutual Information）衡量變數間的相互依賴性。其本質為熵差，即 ( )− ( | )，即知道另一個變數資訊後混亂的降低程度 。當且僅當兩個隨機變數獨立時MI等於零。MI值越高，兩變數之間的相關性則越強。與Pearson相關和F統計量相比，它還捕獲了非線性關係。

公式:

  • 若兩個變數均為離散變數:

p( , )( , ) 為x和y的聯合概率質量函式 (PMF)， p ( )則為x的的聯合概率質量函式 (PMF)。

  • 若兩個變數均為連續變數:

p( , )( , ) 為x和y的聯合概率密度函式 (PDF)，p ( )則為x的的聯合概率密度函式 (PDF)。連續變數情形下，在實際操作中，往往先對資料離散化分桶，然後逐個桶進行計算。

但是實際上，一種極有可能的情況是，x和y中的一個可能是離散變數，而另一個是連續變數。因此在sklearn中，它基於[1]和[2]中提出的基於k最臨近演算法的熵估計非引數方法。

[1] A. Kraskov, H. Stogbauer and P. Grassberger, “Estimating mutual information”. Phys. Rev. E 69, 2004.

[2] B. C. Ross “Mutual Information between Discrete and Continuous Data Sets”. PLoS ONE 9(2), 2014.

import numpy as np
from sklearn.feature_selection import mutual_info_classif
from sklearn.feature_selection import SelectKBest

# 直接載入資料集
from sklearn.datasets import load_iris # 利用iris資料作為演示資料集
iris = load_iris()
X, y = iris.data, iris.target
# 此資料集中，X為連續變數，y為類別變數
# 滿足MI的使用條件

# iris 資料集使用前需要被打亂順序
np.random.seed(1234)
idx = np.random.permutation(len(X))
X = X[idx]
y = y[idx]

# 選擇前100個觀測點作為訓練集
# 剩下的作為測試集
train_set = X[0:100,:]
test_set = X[100:,]
train_y = y[0:100]

# KNN中的臨近數是一個非常重要的引數
# 故我們重寫了一個新的MI計算方程更好的來控制這一引數
def udf_MI(X, y):
    result = mutual_info_classif(X, y, n_neighbors = 5) # 使用者可以輸入想要的臨近數
    return result

# SelectKBest 將會基於一個判別方程自動選擇得分高的變數
# 這裡的判別方程為F統計量
selector = SelectKBest(udf_MI, k=2) # k => 我們想要選擇的變數數
selector.fit(train_set, train_y) # 在訓練集上訓練
transformed_train = selector.transform(train_set) # 轉換訓練集
transformed_train.shape #(100, 2), 其選擇了第三個及第四個變數 
assert np.array_equal(transformed_train, train_set[:,[2,3]])

transformed_test = selector.transform(test_set) # 轉換測試集
assert np.array_equal(transformed_test, test_set[:,[2,3]]);
# 可見對於測試集，其依然選擇了第三個及第四個變數

# 驗算上述結果
for idx in range(train_set.shape[1]):
    score = mutual_info_classif(train_set[:,idx].reshape(-1,1), train_y, n_neighbors = 5)
    print(f"第{idx + 1}個變數與因變數的互資訊為{round(score[0],2)}")
# 故應選擇第三個及第四個變數

第1個變數與因變數的互資訊為0.56
第2個變數與因變數的互資訊為0.28
第3個變數與因變數的互資訊為0.99
第4個變數與因變數的互資訊為1.02

本文首發於公眾號：AI有道（ID: redstonewill），歡迎關注！

相關文章