実践機械学習勉強会 Part 5 決定木系モデル
決定木系モデルを一通り学ぶ。
決定木は、予測値の解釈性に優れている点で魅力的なモデルである。一連の質問に基づいて決断を下すという方法を繰り返すことでデータを分類する。
Classificationには分類木、Regressionには回帰木が対応する。
質問によるデータの分割は、情報利得（information gain: 分割された集合の要素のばらつきの減少）が最大となる特徴量で分割する方法をとる。情報利得$IG$は以下のように定式化される。
$$IG(D_p,f)=I(D_p)-\sum^m_{j=1}\frac{N_j}{N_p}I(D_j)$$
$f$は分割を行う特徴量、$D_p$は親（分割前）のデータセット、$D_j$はj番目の子ノード（分割後）のデータセット、$I$は不純度（impurity: 異なるクラスのサンプルがノードにどの程度でまざっているか定量化する指標）、$N_p$は親ノードのサンプル数、$N_j$はj番目の子ノードのサンプル数をそれぞれ指す。式の通り、$IG$は子ノードの不純度が低いほど大きくなる。不純度指標をエントロピーとして、親ノードから3つ以上の子ノードを生み出すことを繰り返すアルゴリズムは"C4.5"と呼ばれることで有名である。上式はm個の子ノードに分割する一般的な式だが、組み合わせ探索空間を減らすため、ほとんどのライブラリは二分決定木を実装している。つまり、親ノード$D_p$から、子ノード$D_{left}$・$D_{right}$に分割される。よって式は以下のようになる。 $$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)=\frac{N_{left}}{N_p}\frac{N_{right}}{N_p}(\sum^c_{i=1}|p(i|D_{left})-p(i|D_{right})|)^2$$
で規定される。
#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
不純度（impurity）を表す指標は分類木では、ジニ不純度（Gini impurity）:$I_G$、、エントロピー（entropy）:$I_H$、分類誤差（classification error）:$I_E$の３つ、回帰木では平均2乗和誤差（MSE）が存在する。
誤分類の確率を最小化する条件と解釈できる指標である。定義は以下。 $$I_G(t)=\sum^c_{i=1}p(i|t)(1-p(i|t))=1-\sum^c_{i=1}p(i|t)^2$$
最大になるのはクラスが完全に混合されている時である。
#collapse-show
def gini(target):
N = len(target)
_, counts = np.unique(target,return_counts=True)
p = counts / N
return 1 - np.sum(p**2)
すべての空でないクラス$i$を対象にしたエントロピーの定義は以下のようになる。 $$I_H(t)=-\sum^c_{i=1}p(i|t)log_2p(i|t)$$
$p(i|t)$はノード$t$においてラベルがクラス$i$のサンプルの割合を指す。エントロピーは平均情報量ともいわれる。
「情報量（じょうほうりょう）やエントロピー（英: entropy）は、情報理論の概念で、あるできごと（事象）が起きた際、それがどれほど起こりにくいかを表す尺度である。ありふれたできごと（たとえば「風の音」）が起こったことを知ってもそれはたいした「情報」にはならないが、逆に珍しいできごと（たとえば「曲の演奏」）が起これば、それはより多くの「情報」を含んでいると考えられる。情報量はそのできごとが本質的にどの程度の情報を持つかの尺度であるとみなすこともできる。」(by wiki)
ここで、 事象$E$が起こる確率を$P(E)$とするとき、 事象 $E$ が起こったことを知らされたとき受け取る（選択）情報量$I(E)$ を $$I(E)=\log \frac {1}{P(E)}=-\log P(E)$$ と定義される。よってエントロピーは有限集合中、起きうる全事象から受け取る情報量をそれぞれの事象が起きる確率で重みづけ平均した情報量であり、平均情報量と呼ばれる所以である。
よって今回$I_H$は有限集合であるノード$t$から受け取ることのできる平均情報量を表しているのである。
2値分類の場合、$p(i=1|t)=1or0$でエントロピーは最小$0$になる。エントロピーが最大になるのはジニ不純度と同様にクラスが完全に混合されている場合、つまり$p(i=1|t)=0.5$の時である。
#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))
#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)$$ $$B: (40,40)-->(20,40)(20,0)$$ それぞれのケースで情報利得はどのような値をとるのだろうか？
#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)])
ジニ不純度の場合、
#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))
となり、Bでの分割を優先することがわかる。実際問題そちらの方がより純粋である。
エントロピーの場合、
#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))
となり、同様にBでの分割を優先することがわかる。このように、ジニ不純度とエントロピーは非常によく似た結果になることが知られていて、2択に時間をかける優先性はあまりない。
分類誤差の場合、
#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とBを同等に評価していることがわかる。このように分類誤差はノードのクラス確率変化にあまり敏感ではないので決定木の成長に適していない。一方、過学習を防ぐため、決定木の分岐の深さに制限を設ける剪定（prune）には役立つ。
ちなみにtwoingはどうなるのか、
#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))
なるほど、一応Bをうまいこと優先しているようだ。
回帰に決定木を使用するには、連続値変数に適した不純度指標が必要である。そこで、ノード$t$の不純度指標として代わりにこの指標が使われ、 $$I(t)=MSE(t)=\frac{1}{N_t}\sum_{i\in D_t}(y^{(i)}-\hat{y}_t)^2$$ $$\hat{y}_t=\frac{1}{N_t}\sum_{i\in D_t}y^{(i)}$$ と定義される。ここで、$N_t$はノード$t$のサンプル数、$D_t$はノード$t$のサブセットインデックスの集合、$y^{(i)}$はラベル（真の目的値）、$\hat{y}_t$はサンプルの平均（予測された目的値）として扱う。
決定木回帰の文脈で、この指標はよく「分割後のノード分散」と呼ばれる。分割条件はこれにならってよく「分散減少」（variance reduction）と呼ばれる。
#collapse-show
def mse(target):
return np.mean((target - np.mean(target))**2)
それでは単純な1本だけの決定木を実装してみよう。ここでは単純のため分類木に焦点をあてて実装してみる。（前に作った関数を流用したいのでクラス内メソッドに書き直しません、後決定木の図示はめんどいのでSklearnに任せます..）
親ノードから左右ノードへ分岐する際に情報利得が最も大きくなるように、分岐の基準とする特徴量とその閾値を求める。
#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データをお試しに使う。まずはデータを読み込む。
#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)
先ほどの関数を使ってみると、
#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]))
petal length (cm) 特徴量で閾値を3にした時に初めのベストな分割が行えるようだ。
ノードクラスをつくって、クラス内でsplitメソッドにより子ノードを再帰的に繰り返す。これは左右分割後に、それぞれに対応するノードクラスをleftとright変数に格納してあげることで達成できる。
具体的に、変数left.method_nameはleftに格納されたノードクラスのmethod_nameを呼び出していて、そのメソッドでさらにleftノードクラスのleftとrightに分割されたノードクラスが格納され、さらに格納されたノードクラスのmethod_nameが呼び出されるというように繰り返される。制限をつけないと無限ループするので、ストップする基準としてmax_depthに達するか、左右子ノードの不純度が０（クラスが1種類しか含まない）という条件を使うことにする。
結果として学習後にこのノードクラスが成長後のすべてのノードを蓄えていることになる。また、IG_maxを特徴量ごとに足し合わせれば各特徴量の重要度を計算できる。
#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)
最後に決定木分類器クラスを作ってデータへの学習と予測メソッドを設定する
#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データセットで学習してみると、
#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))
のようになり、テストデータに対し98%程の正解率をだしており上手いこと成長しているのがわかる。また、左側のノードは深さが３まで、右側のノードは深さが１で止まっていることがわかる。各特徴量の重要度も求められていて、
#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()
のように図示される。petal length (cm)、つまり花びらの長さが一番分類に寄与していることがわかる。
sklearnと比較してみると、
#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))
となり、同様の精度まで決定木が成長していることがわかる。
またSkleanのライブラリでは決定木の可視化ツールも用意していて以下のように使える。
#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")
#collapse-show
from IPython.display import Image,display_png
display_png(Image('iris-tree.png'))
決定木は半エキスパートシステムみたいなもので、この図のように最適化された質問分岐を表示してくれたり、解釈性に非常に優れている。
Refrences
- Python Machine Learning, 2nd Edition
- http://darden.hatenablog.com/entry/2016/12/15/222447#%E5%86%8D%E5%B8%B0%E5%91%BC%E3%81%B3%E5%87%BA%E3%81%97
- http://hktech.hatenablog.com/entry/2018/10/05/004235
- https://techblog.nhn-techorus.com/archives/14801
- https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingClassifier.html
- https://en.wikipedia.org/wiki/Gradient_boosting
- https://qiita.com/Quasi-quant2010/items/10f7ad4ed2e11004990f
- https://zaburo-ch.github.io/post/xgboost/
- https://www.st-hakky-blog.com/entry/2017/08/07/163216
- https://ja.wikipedia.org/wiki/%E6%B1%BA%E5%AE%9A%E6%9C%A8
- https://ja.wikipedia.org/wiki/%E6%83%85%E5%A0%B1%E9%87%8F
- ftp://ftp.boulder.ibm.com/software/analytics/spss/support/Stats/Docs/Statistics/Algorithms/14.0/TREE-CART.pdf
- Bishop Pattern Recognition And Machine Learning 14.4 Tree-based Models p665