前言

在特徵的選擇過程中，如果學習器（基學習器）是樹模型的話，可以根據特徵的重要性來篩選有效的特徵。本文是對Random Forest、GBDT、XGBoost如何用在特徵選擇上做一個簡單的介紹。

各種模型的特徵重要性計算

Random Forests

• 袋外資料錯誤率評估
  RF的資料是boostrap的有放回取樣，形成了袋外資料。因此可以採用袋外資料（OOB）錯誤率進行特徵重要性的評估。
  袋外資料錯誤率定義為：袋外資料自變數值發生輕微擾動後的分類正確率與擾動前分類正確率的平均減少量。
  （1）對於每棵決策樹，利用袋外資料進行預測，將袋外資料的預測誤差記錄下來，其每棵樹的誤差為vote1,vote2,…,voteb
  （2）隨機變換每個預測變數，從而形成新的袋外資料，再利用袋外資料進行驗證，其每個變數的誤差是votel1,votel2,…votelb

• Gini係數評價指標 （和GBDT的方法相同）

GBDT

在sklearn中，GBDT和RF的特徵重要性計算方法是相同的，都是基於單棵樹計算每個特徵的重要性，探究每個特徵在每棵樹上做了多少的貢獻，再取個平均值。
在利用隨機森林對特徵重要性進行評估寫的比較清楚了，但是還是有一點小的問題，比較ensemble模型 零碎記錄中對原始碼的解析可以看出，前者計算中丟失了weighted_n_node_samples。

• 利用Gini計算特徵的重要性
  單棵樹上特徵的重要性定義為：特徵在所有非葉節在分裂時加權不純度的減少，減少的越多說明特徵越重要。
  沿用參考部落格裡的符號，我們將變數重要性評分（variable importance measures）用
  VIM$VIM$來表示，將Gini指數用GI$GI$來表示
  節點m的Gini指數的計算公式為：
  GIm=1−∑k=1|K|p2mk$$G{I}_m=1-\sum _{k=1}^{|K|}{p}_{mk}^2$$
  其中，K表示有K個類別，pmk$p_{mk}$表示節點m中類別k所佔的比例。直觀地說，就是隨便從節點m中隨機抽取兩個樣本，其類別標記不一致的概率。
  特徵Xj$X_j$在節點m$m$的重要性可以表示為加權不純度的減少
  VIMGinijm=Nm×GIm−Nl×GIl−Nr×GIr$$VI{M}_{jm}^{Gini}=N_m\times GIm-N_l\times GIl-N_r\times GIr$$
  其中，GIl$G{I}_l$和GIr$G{I}_r$分別表示分枝後兩個新節點的Gini指數。Nm$N_m$、Nl$N_l$、Nr$N_r$表示節點m、左孩子節點l和右孩子節點r的樣本數。
  如果，特徵
  Xj$X_j$在決策樹i中出現的節點在集合M中，那麼Xj$X_j$在第i顆樹的重要性為
  VIMij=∑m∈MVIMjm$$VI{M}_{ij}={\sum }_{m\in M}VI{M}_{jm}$$

~~如果這樣還不是很清晰的話，我們來舉個例子（李航統計學習方法表5.1）

import numpy as np
from sklearn.tree import DecisionTreeClassifier
from sklearn.externals.six import StringIO
from sklearn import tree
import pydotplus

clf = DecisionTreeClassifier()
x = [[1,1,1,1,1,2,2,2,2,2,3,3,3,3,3],
     [1,1,2,2,1,1,1,2,1,1,1,1,2,2,1],
     [1,1,1,2,1,1,1,2,2,2,2,2,1,1,1],
     [1,2,2,1,1,1,2,2,3,3,3,2,2,3,1]
     ]
y =  [1,1,2,2,1,1,1,2,2,2,2,2,2,2,1]
x = np.array(x)
x = np.transpose(x)
clf.fit(x,y)
print clf.feature_importances_
feature_name = ['A1','A2','A3','A4']
target_name = ['1','2']
dot_data = StringIO()
tree.export_graphviz(clf,out_file = dot_data,feature_names=feature_name,
                     class_names=target_name,filled=True,rounded=True,
                     special_characters=True)
graph = pydotplus.graph_from_dot_data(dot_data.getvalue())
graph.write_pdf("WineTree.pdf")
print('Visible tree plot saved as pdf.')

可以得到樹的劃分過程圖

特徵A3的重要性為 0.48×15−0.4444×9−0=3.2004$0.48\times 15-0.4444\times 9-0=3.2004$
特徵A2的重要性為 0.4444×9−0−0=3.9996$0.4444\times 9-0-0=3.9996$
特徵A1和A4的重要性都為0
所以該棵樹上所有節點總的加權不純度減少量為 3.2004+3.9996=7.3$3.2004+3.9996=7.3$
對其進行歸一化操作可以得到A1、A2、A3、A4的特徵重要性為

[ 0. 0.55555556 0.44444444 0. ]

這是單棵樹上特徵的計算方法，推廣到n棵樹

XGBoost

    def get_score(self, fmap='', importance_type='weight'):
        """Get feature importance of each feature.
        Importance type can be defined as:
            'weight' - the number of times a feature is used to split the data across all trees.
            'gain' - the average gain of the feature when it is used in trees
            'cover' - the average coverage of the feature when it is used in trees
        Parameters
        ----------
        fmap: str (optional)
           The name of feature map file
        """

        if importance_type not in ['weight', 'gain', 'cover']:
            msg = "importance_type mismatch, got '{}', expected 'weight', 'gain', or 'cover'"
            raise ValueError(msg.format(importance_type))

        # if it's weight, then omap stores the number of missing values
        if importance_type == 'weight':
            # do a simpler tree dump to save time
            trees = self.get_dump(fmap, with_stats=False)

            fmap = {}
            for tree in trees:
                for line in tree.split('\n'):
                    # look for the opening square bracket
                    arr = line.split('[')
                    # if no opening bracket (leaf node), ignore this line
                    if len(arr) == 1:
                        continue

                    # extract feature name from string between []
                    fid = arr[1].split(']')[0].split('<')[0]

                    if fid not in fmap:
                        # if the feature hasn't been seen yet
                        fmap[fid] = 1
                    else:
                        fmap[fid] += 1

            return fmap

        else:
            trees = self.get_dump(fmap, with_stats=True)

            importance_type += '='
            fmap = {}
            gmap = {}
            for tree in trees:
                for line in tree.split('\n'):
                    # look for the opening square bracket
                    arr = line.split('[')
                    # if no opening bracket (leaf node), ignore this line
                    if len(arr) == 1:
                        continue

                    # look for the closing bracket, extract only info within that bracket
                    fid = arr[1].split(']')

                    # extract gain or cover from string after closing bracket
                    g = float(fid[1].split(importance_type)[1].split(',')[0])

                    # extract feature name from string before closing bracket
                    fid = fid[0].split('<')[0]

                    if fid not in fmap:
                        # if the feature hasn't been seen yet
                        fmap[fid] = 1
                        gmap[fid] = g
                    else:
                        fmap[fid] += 1
                        gmap[fid] += g

            # calculate average value (gain/cover) for each feature
            for fid in gmap:
                gmap[fid] = gmap[fid] / fmap[fid]
            return gmap

在XGBoost中提供了三種特徵重要性的計算方法：

‘weight’ - the number of times a feature is used to split the data across all trees.
‘gain’ - the average gain of the feature when it is used in trees
‘cover’ - the average coverage of the feature when it is used in trees

簡單來說
weight就是在所有樹中特徵用來分割的節點個數總和；
gain就是特徵用於分割的平均增益
cover 的解釋有點晦澀，在[R-package/man/xgb.plot.tree.Rd]有比較詳盡的解釋：(https://github.com/dmlc/xgboost/blob/f5659e17d5200bd7471a2e735177a81cb8d3012b/R-package/man/xgb.plot.tree.Rd)：the sum of second order gradient of training data classified to the leaf, if it is square loss, this simply corresponds to the number of instances in that branch. Deeper in the tree a node is, lower this metric will be。實際上coverage可以理解為被分到該節點的樣本的二階導數之和，而特徵度量的標準就是平均的coverage值。

還是舉李航書上那個例子，我們用不同顏色來表示不同的特徵，繪製下圖

import xgboost as xgb
import numpy as np
x = [[1,1,1,1,1,2,2,2,2,2,3,