Ccmmutty logo
Commutty IT
・2 min read

０から決定木!

pythonMetrics
https://cdn.magicode.io/media/notebox/7720b971-b942-422c-9181-0d39dbf6ca2a.jpeg
決定木系モデルを一通り学ぶ。

1 決定木のアルゴリズム

決定木は、予測値の解釈性に優れている点で魅力的なモデルである。一連の質問に基づいて決断を下すという方法を繰り返すことでデータを分類する。
Classificationには分類木、Regressionには回帰木が対応する。
image.png

1.1 情報利得の最大化

質問によるデータの分割は、情報利得（information gain: 分割された集合の要素のばらつきの減少）が最大となる特徴量で分割する方法をとる。情報利得IGIGは以下のように定式化される。
IG(Dp,f)=I(Dp)−∑j=1mNjNpI(Dj)IG(D_p,f)=I(D_p)-\sum^m_{j=1}\frac{N_j}{N_p}I(D_j)
ffは分割を行う特徴量、DpD_pは親（分割前）のデータセット、DjD_jはj番目の子ノード（分割後）のデータセット、IIは不純度（impurity: 異なるクラスのサンプルがノードにどの程度でまざっているか定量化する指標）、NpN_pは親ノードのサンプル数、NjN_jはj番目の子ノードのサンプル数をそれぞれ指す。式の通り、IGIGは子ノードの不純度が低いほど大きくなる。**不純度指標をエントロピーとして、親ノードから3つ以上の子ノードを生み出すことを繰り返すアルゴリズムは"C4.5"と呼ばれることで有名である。**上式はm個の子ノードに分割する一般的な式だが、組み合わせ探索空間を減らすため、ほとんどのライブラリは二分決定木を実装している。つまり、親ノードDpD_pから、子ノードDleftD_{left}・DrightD_{right}に分割される。よって式は以下のようになる。
IG(Dp,f)=I(Dp)−NleftNpI(Dleft)−NrightNpI(Dright)IG(D_p,f)=I(D_p)-\frac{N_{left}}{N_p}I(D_{left})-\frac{N_{right}}{N_p}I(D_{right})
親ノードから左右２つの子ノードを生み出すことを繰り返すアルゴリズムは"CART"（Classfication And Regression Tree）と呼ばれることで有名である。
CARTでは分類タスクの時、twoing crietionという情報利得について別のメソッドをオプションに持っていて、
IG(Dp,f)=NleftNpNrightNp(∑i=1c∣p(i∣Dleft)−p(i∣Dright)∣)2IG(Dp,f)=\frac{N_{left}}{N_p}\frac{N_{right}}{N_p}(\sum^c_{i=1}|p(i|D_{left})-p(i|D_{right})|)^2
で規定される。
note: あらゆるライブラリで実装されているfeature_importance、つまり、特徴量ffの重要度は何を示しているのかというと、特徴量ffで分割することでどのくらい情報利得が生じたかを表している。
python
#collapse-show
import numpy as np

# 一般的な２分岐を想定する
def IG(impurity_func, target_p, target_left, target_right):
    N_left = len(target_left)
    N_right = len(target_right)
    N_p = len(target_p)
    I_p  = impurity_func(target_p)
    I_left = impurity_func(target_left)
    I_right = impurity_func(target_right)
    return I_p - N_left / N_p * I_left - N_right / N_p * I_right

def IG_twoing(target_left, target_right):
    N_left = len(target_left)
    N_right = len(target_right)
    N_p = N_left + N_right
    classes = np.unique(np.hstack([target_right,target_left]))
    I = 0
    for i in classes:
        I += np.abs(len(target_left[target_left==i])/N_left - len(target_right[target_right==i])/N_right)
    I = N_left * N_right / N_p**2 * I**2
    return I

1.2 不純度の指標

不純度（impurity）を表す指標は分類木では、ジニ不純度（Gini impurity）:IGI_G、、エントロピー（entropy）:IHI_H、分類誤差（classification error）:IEI_Eの３つ、回帰木では平均2乗和誤差（MSE）が存在する。
1.2.1 ジニ不純度
誤分類の確率を最小化する条件と解釈できる指標である。定義は以下。
IG(t)=∑i=1cp(i∣t)(1−p(i∣t))=1−∑i=1cp(i∣t)2I_G(t)=\sum^c_{i=1}p(i|t)(1-p(i|t))=1-\sum^c_{i=1}p(i|t)^2
最大になるのはクラスが完全に混合されている時である。
python
#collapse-show
def gini(target):
    N = len(target)
    _, counts = np.unique(target,return_counts=True)
    p = counts / N
    return 1 - np.sum(p**2)
1.2.2 エントロピー
すべての空でないクラスiiを対象にしたエントロピーの定義は以下のようになる。
IH(t)=−∑i=1cp(i∣t)log2p(i∣t)I_H(t)=-\sum^c_{i=1}p(i|t)log_2p(i|t)
p(i∣t)p(i|t)はノードttにおいてラベルがクラスiiのサンプルの割合を指す。エントロピーは平均情報量ともいわれる。
「情報量（じょうほうりょう）やエントロピー（英: entropy）は、情報理論の概念で、あるできごと（事象）が起きた際、それがどれほど起こりにくいかを表す尺度である。ありふれたできごと（たとえば「風の音」）が起こったことを知ってもそれはたいした「情報」にはならないが、逆に珍しいできごと（たとえば「曲の演奏」）が起これば、それはより多くの「情報」を含んでいると考えられる。情報量はそのできごとが本質的にどの程度の情報を持つかの尺度であるとみなすこともできる。」(by wiki)
ここで、 事象EEが起こる確率をP(E)P(E)とするとき、 事象 EE が起こったことを知らされたとき受け取る（選択）情報量I(E)I(E) を
I(E)=log⁡1P(E)=−log⁡P(E)I(E)=\log \frac {1}{P(E)}=-\log P(E)
と定義される。よってエントロピーは有限集合中、起きうる全事象から受け取る情報量をそれぞれの事象が起きる確率で重みづけ平均した情報量であり、平均情報量と呼ばれる所以である。
よって今回IHI_Hは有限集合であるノードttから受け取ることのできる平均情報量を表しているのである。
2値分類の場合、p(i=1∣t)=1or0p(i=1|t)=1or0でエントロピーは最小00になる。エントロピーが最大になるのはジニ不純度と同様にクラスが完全に混合されている場合、つまりp(i=1∣t)=0.5p(i=1|t)=0.5の時である。
python
#collapse-show
def entropy(target):
    N = len(target)
    _, counts = np.unique(target,return_counts=True)
    p = counts / N
    return -np.sum(p*np.log2(p))
1.2.3 分類誤差
以下のように定義される。
IE(t)=1−max(p(i∣t))I_E(t)=1-max(p(i|t))
多数派がどれだけ大多数か示すことで純度を測っている。
python
#collapse-show
def error(target):
    N = len(target)
    _, counts = np.unique(target,return_counts=True)
    p = counts / N
    return 1 - np.max(p)
以上3つの指標を比較してみよう。
ここでは2値分類を例として考える。クラス（0, 1）に対して親ノードのサンプル数を（40, 40）とする。以下の2通りの分岐を考えてみる。
A:(40,40)−−>(30,10)(10,30)A: (40,40)-->(30,10)(10,30)
B:(40,40)−−>(20,40)(20,0)B: (40,40)-->(20,40)(20,0)
それぞれのケースで情報利得はどのような値をとるのだろうか？
python
#collapse-show
target_p = np.hstack([np.zeros(40),np.ones(40)])
target_la = np.hstack([np.zeros(30),np.ones(10)])
target_ra = np.hstack([np.zeros(10),np.ones(30)])
target_lb = np.hstack([np.zeros(20),np.ones(40)])
target_rb = np.hstack([np.zeros(20),np.ones(0)])
ジニ不純度の場合、
python
#collapse-show
IG_a = IG(gini,target_p,target_la,target_ra)
IG_b = IG(gini,target_p,target_lb,target_rb)
print("A: {}, B: {}".format(IG_a,IG_b))
A: 0.125, B: 0.16666666666666669
となり、Bでの分割を優先することがわかる。実際問題そちらの方がより純粋である。
エントロピーの場合、
python
#collapse-show
IG_a = IG(entropy,target_p,target_la,target_ra)
IG_b = IG(entropy,target_p,target_lb,target_rb)
print("A: {}, B: {}".format(IG_a,IG_b))
A: 0.1887218755408671, B: 0.31127812445913283
となり、同様にBでの分割を優先することがわかる。このように、ジニ不純度とエントロピーは非常によく似た結果になることが知られていて、2択に時間をかける優先性はあまりない。
分類誤差の場合、
python
#collapse-show
IG_a = IG(error,target_p,target_la,target_ra)
IG_b = IG(error,target_p,target_lb,target_rb)
print("A: {}, B: {}".format(IG_a,IG_b))
A: 0.25, B: 0.25
となり、AとBを同等に評価していることがわかる。このように分類誤差はノードのクラス確率変化にあまり敏感ではないので決定木の成長に適していない。一方、過学習を防ぐため、決定木の分岐の深さに制限を設ける剪定（prune）には役立つ。
ちなみにtwoingはどうなるのか、
python
#collapse-show
IG_a = IG_twoing(target_la,target_ra)
IG_b = IG_twoing(target_lb,target_rb)
print("A: {}, B: {}".format(IG_a,IG_b))
A: 0.25, B: 0.33333333333333337
なるほど、一応Bをうまいこと優先しているようだ。
1.2.4 最小2乗和誤差
回帰に決定木を使用するには、連続値変数に適した不純度指標が必要である。そこで、ノードttの不純度指標として代わりにこの指標が使われ、
I(t)=MSE(t)=1Nt∑i∈Dt(y(i)−y^t)2I(t)=MSE(t)=\frac{1}{N_t}\sum_{i\in D_t}(y^{(i)}-\hat{y}_t)^2
y^t=1Nt∑i∈Dty(i)\hat{y}_t=\frac{1}{N_t}\sum_{i\in D_t}y^{(i)}
と定義される。ここで、NtN_tはノードttのサンプル数、DtD_tはノードttのサブセットインデックスの集合、y(i)y^{(i)}はラベル（真の目的値）、y^t\hat{y}_tはサンプルの平均（予測された目的値）として扱う。
決定木回帰の文脈で、この指標はよく「分割後のノード分散」と呼ばれる。分割条件はこれにならってよく「分散減少」（variance reduction）と呼ばれる。
python
#collapse-show
def mse(target):
    return np.mean((target - np.mean(target))**2)

2 決定木の実装

それでは単純な1本だけの決定木を実装してみよう。ここでは単純のため分類木に焦点をあてて実装してみる。（前に作った関数を流用したいのでクラス内メソッドに書き直しません、後決定木の図示はめんどいのでSklearnに任せます..）

2.1 最適な特徴量と閾値を選択する

親ノードから左右ノードへ分岐する際に情報利得が最も大きくなるように、分岐の基準とする特徴量とその閾値を求める。
python
#collapse-show
def search_best_split(impurity_func, data, target):   
    features = data.shape[1]
    best_thrs = None
    best_f = None
    IG_max = 0
 
    for feat_idx in range(features):
        values = data[:, feat_idx] 
        for val in values:
            target_left = target[values < val]
            target_right = target[values >= val]
            ig = IG(impurity_func, target, target_left, target_right)
            if IG_max < ig: # 情報利得の最大値を更新
                IG_max = ig
                best_thrs = val
                best_f = feat_idx
    return IG_max, best_thrs, best_f
今回はSklearnのIrisデータをお試しに使う。まずはデータを読み込む。
python
#collapse-show
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split

iris = load_iris()
features_name = iris.feature_names
data = iris.data
target = iris.target 
X_train, X_test, y_train, y_test = train_test_split(data, target, test_size=0.3, random_state=0)
先ほどの関数を使ってみると、
python
#collapse-show
IG_max, best_thrs, best_f = search_best_split(gini, X_train, y_train)
print('Information gain: {}, Best threshold: {}, Best feature: {}'.format(IG_max, best_thrs, features_name[best_f]))
Information gain: 0.3294995369039634, Best threshold: 3.0, Best feature: petal length (cm)
petal length (cm) 特徴量で閾値を3にした時に初めのベストな分割が行えるようだ。

2.2 再帰的にノード分割を行う

ノードクラスをつくって、クラス内でsplitメソッドにより子ノードを再帰的に繰り返す。これは左右分割後に、それぞれに対応するノードクラスをleftとright変数に格納してあげることで達成できる。
具体的に、変数left.method_nameはleftに格納されたノードクラスのmethod_nameを呼び出していて、そのメソッドでさらにleftノードクラスのleftとrightに分割されたノードクラスが格納され、さらに格納されたノードクラスのmethod_nameが呼び出されるというように繰り返される。制限をつけないと無限ループするので、ストップする基準としてmax_depthに達するか、左右子ノードの不純度が０（クラスが1種類しか含まない）という条件を使うことにする。
結果として学習後にこのノードクラスが成長後のすべてのノードを蓄えていることになる。また、IG_maxを特徴量ごとに足し合わせれば各特徴量の重要度を計算できる。
python
#collapse-show
feature_importances = {}

class Node:
    def __init__(self, data, target, max_depth, impurity):
        self.left = None
        self.right = None
        self.max_depth = max_depth
        self.depth = None
        self.data = data
        self.target = target
        self.threshold = None
        self.feature = None
        self.IG_max = None
        self.importance = {}
        self.label = np.argmax(np.bincount(target)) # ノード内で一番多いラベル（いわゆるそのノードでの予測値）
        self.impurity = impurity
        if self.impurity == 'gini':
            self.impurity_func = gini
        elif self.impurity == 'entropy':
            self.impurity_func = entropy
        elif self.impurity == 'error':
            self.impurity_func = error

    # 再帰的分割メソッド
    def split(self, depth):
        self.depth = depth
        self.IG_max, self.threshold, self.feature = search_best_split(self.impurity_func, self.data, self.target)
        print('Depth: {}, Sep at Feature: {},Threshold: {}, Label: {}'.format(self.depth, self.feature, self.threshold, self.label))
        
            
        # 各特徴量の重要度をグローバル変数に記録
        global feature_importances
        if self.IG_max != 0:
            if self.feature not in feature_importances:
                feature_importances[self.feature] = self.IG_max
            else:
                feature_importances[self.feature] += self.IG_max
        
        # 剪定（prune)して無限ループを防ぐ
        if self.depth == self.max_depth or self.IG_max == self.impurity_func(self.target):
            return       
        
        # 得られた特徴量と閾値にしたがてデータを2分割する
        idx_left = self.data[:, self.feature] >= self.threshold
        idx_right = self.data[:, self.feature] < self.threshold

        # 分割された左右ノードを対応する変数に格納し、再帰的に分割を行う
        self.left = Node(self.data[idx_left],  self.target[idx_left], self.max_depth, self.impurity)
        self.right = Node(self.data[idx_right], self.target[idx_right], self.max_depth, self.impurity)
        self.left.split(self.depth +1)
        self.right.split(self.depth +1)
    
    # 渡されたデータを末端（葉）ノードに達するまで分岐にそって流して、行きついた先の対応するラベルを予測値として返す
    def predict(self, data):
        if self.IG_max == self.impurity_func(self.target) or self.depth == self.max_depth:
            return self.label
        else:
            if data[self.feature] > self.threshold:
                return self.left.predict(data)
            else:
               return self.right.predict(data)

2.3 All in one

最後に決定木分類器クラスを作ってデータへの学習と予測メソッドを設定する
python
#collapse-show
class DesicionTreeClassifier:
    def __init__(self, max_depth, impurity):
        self.max_depth = max_depth
        self.impurity = impurity
        self.tree = None
    
    # ルートノードから分割を繰り返し決定木を成長させる
    def fit(self, data, target):
        initial_depth = 0
        self.tree = Node(data, target, self.max_depth, self.impurity)
        self.tree.split(initial_depth)
    
    # 成長した決定木でデータのラベルを予測する
    def predict(self, data):
        pred = []
        for s in data:
            pred.append(self.tree.predict(s))
        return np.array(pred)
最後にIrisデータセットで学習してみると、
python
#collapse-show
clf = DesicionTreeClassifier(3, 'gini')
clf.fit(X_train, y_train)
y_pred = clf.predict(X_test)
score = sum(y_pred == y_test)/float(len(y_test))
print('Classification accuracy: {}'.format(score))
Depth: 0, Sep at Feature: 2,Threshold: 3.0, Label: 2 Depth: 1, Sep at Feature: 2,Threshold: 5.0, Label: 2 Depth: 2, Sep at Feature: 2,Threshold: 5.1, Label: 2 Depth: 3, Sep at Feature: None,Threshold: None, Label: 2 Depth: 3, Sep at Feature: 0,Threshold: 6.7, Label: 2 Depth: 2, Sep at Feature: 3,Threshold: 1.7, Label: 1 Depth: 3, Sep at Feature: 1,Threshold: 3.2, Label: 2 Depth: 3, Sep at Feature: None,Threshold: None, Label: 1 Depth: 1, Sep at Feature: None,Threshold: None, Label: 0 Classification accuracy: 0.9777777777777777
のようになり、テストデータに対し98%程の正解率をだしており上手いこと成長しているのがわかる。また、左側のノードは深さが３まで、右側のノードは深さが１で止まっていることがわかる。各特徴量の重要度も求められていて、
python
#collapse-show
import matplotlib.pyplot as plt
%matplotlib inline

keys = [features_name[i] for i in feature_importances.keys()]
values = list(feature_importances.values())
plt.figure(figsize=(10,7))
plt.title('Feature importances', fontsize=15)
plt.bar(keys, values)
plt.show()
png
のように図示される。petal length (cm)、つまり花びらの長さが一番分類に寄与していることがわかる。
sklearnと比較してみると、
python
#collapse-show
from sklearn.tree import DecisionTreeClassifier as DecisionTreeClassifier2

clf2 = DecisionTreeClassifier2(max_depth=3)
clf2.fit(X_train,y_train)
y_pred = clf2.predict(X_test)
score = sum(y_pred == y_test)/float(len(y_test))
print('Classification accuracy: {}'.format(score))
Classification accuracy: 0.9777777777777777
となり、同様の精度まで決定木が成長していることがわかる。
またSkleanのライブラリでは決定木の可視化ツールも用意していて以下のように使える。
python
#collapse-show
import graphviz
from sklearn.tree import export_graphviz

dot_data = export_graphviz(
                        clf2,
                        class_names=iris.target_names,
                        feature_names=features_name,
                        filled=True,
                        rounded=True,
                        out_file=None
                    )
graph = graphviz.Source(dot_data)
graph.render("iris-tree", format="png")
'iris-tree.png'
python
#collapse-show
from IPython.display import Image,display_png
display_png(Image('iris-tree.png'))

Discussion

コメントにはログインが必要です。