2019-02-09

このブログについて

はじめに

このブログでは，オープンソースを手本にしながら機械学習アルゴリズムを実装していきます。

決定木・ランダムフォレストからスタートして，ツリー系・アンサンブル学習系の手法を中心に扱っていこうと考えています。

f:id:tatamiyatamiyatatatamiya:20190209163341j:plain

2019-03-24

回帰する（ジコリューで行こう！#5）

ジコリューで行こう！決定木機械学習 Python

ヤナギの木の下に人影が…！

f:id:tatamiyatamiyatatatamiya:20190323175407p:plain

ということで，怪奇…もとい，回帰アルゴリズムを実装していきたいと思います。*1

tatamiya.hatenablog.com

前回までは決定木・ランダムフォレストの分類アルゴリズムを実装していきました。今回は連続的な値を予測する回帰を扱えるように拡張を行います。

指針

前回までのあらすじと今回やること

今まで扱ってきた分類問題では，不純度が小さくなるようなデータ分割方法を探索し，学習モデルを構築していました。

一方で，株価の値動きなどのように連続的な値をとるような場合でも，不純度の代わりに誤差を用いることで，今までと同様に学習モデルを作ることができます。

最小化する指標

分類問題では，データセット内に含まれる各クラスの構成率 $p_k (k \in \{1,2,3,...,K\})$ をもとに計算したGini係数

$\sum _ {k=1} ^{K} (1-p_k) p_k = 1 - \sum _{k=1} ^{K} p_k ^2$

の最小化を行いました。*2

一方回帰では，分割を行なった左右のデータ断片について，それぞれ平均値からの逸脱度を計算します。そして左右の逸脱度の和が最小になるように分割を行なっていきます。

逸脱度としては，ここでは二乗和を用いることにします。*3

分割の数式表現

N個の学習データセット $(X_i, y_i)_{i \in \{1,2,...,N\}}$ があるとします。 $X_i = (x_i^{(1)}, x_i^{(2)}, ..., x_i^{(D)})$ はD次元の特徴量ベクトルです。

この時，二乗和による評価指標は以下のようにかけます：

$E(j,s) = \sum_{i \in R_{\rm left}(j,s)} (y_{i} - c_{\rm left})^{2} + \sum_{i \in R_{\rm right}(j,s)} (y_{i} - c_{\rm right})^{2}$

ここで， $R_{\rm left/right}(j,s)$ は分割後左or右断片に入るデータインデックスの集合を， $c_{\rm left/right}$ は各断片内での平均値を表します：

$R_{\rm left}(j,s) = \{ i | x_i^{(j)} \le s \}, \quad R_{\rm right}(j,s) = \{ i | x_i^{(j)} > s \}$

$c_{\rm left/right} = \frac{1}{N_{\rm left/right}}\sum_{i \in R_{\rm left/right}(j,s)} y_{I}, \quad N_{\rm left/right} = \sum_{i \in R_{\rm left/right}(j,s)}1$

$j$ は分割に用いる特徴量のインデックス， $s$ は分割の閾値を表します。

最適な分割は，上記の評価指標 $E(j,s)$ が最小となる $(j,s)$ の組み合わせとして表されます。

実装

サンプルデータ

分類問題ではサンプルデータとしてirisを用いて，葉・萼片の長さと幅から品種を予測していました。

回帰問題を扱うにあたっては，同じくよく使われるデータセットとして， boston housingを用います。

データの説明については，こちらを参照願います。

scikit-learnから，以下のように読み込みます。

from sklearn import datasets

boston = datasets.load_boston()
X = boston.data
y = boston.target

13種類の特徴量からなるXと予測対象となる連続値yの組が，全部で506個含まれています。

最適な分割を探る

探索する特徴量を一つに固定した場合，二乗誤差和が最小になるような分割点は以下のような関数で求められます：

def find_optimal_division(x, y):
    list_criterion = []
    x_unique = np.unique(x)

    if len(x_unique) == 1:
        return None, func_criterion(y)

    for threshold in x_unique[:-1]:

        mask_divide = x > threshold
        y_left = y[mask_divide]
        y_right = y[~mask_divide]

        criterion_divide = criterion_lr(y_left, y_right, func_criterion)
        list_criterion.append(criterion_divide)
        
    array_criterion = np.array(list_criterion)
    i_div_opt = np.argmin(array_criterion)
    
    return x_unique[i_div_opt], array_criterion[i_div_opt]

基本的には分類の時と同じですが，指標の計算部分だけ異なります。

criterion_divide = criterion_lr(y_left, y_right, func_criterion)

ここの部分は，あらかじめ下記のように定義しておきます。

func_criterion = np.var

def criterion_lr(y_left, y_right, func_criterion):
    
    return (func_criterion(y_left) * len(y_left) + func_criterion(y_right) * len(y_right) ) / (len(y_left) + len(y_right))

誤差の二乗和を求めるにあたり，numpyのnp.varを使い分散を求めたうえでデータ数をかけるという方法をとりました。

無駄があるように見えますが，

numpyの最適化されている関数を有効活用できる
別の指標を扱う際にもコードの改変が少ない（後述）

というメリットがあります。

なお，求めた誤差二乗和をデータ数で割ることで，平均二乗誤差(Mean Squared Error, MSE)にしてあります。平均を取らなくても計算上問題ありませんが，こちらも別指標を使う際に共通化できるようにこのようなにしてあります。

これを例によってnumpyのapply_alogn_axis関数に適用することで，全特徴量から最適な分割を探索できます。

def divide(X, y):

    results = np.apply_along_axis(find_optimal_division, 0, X, y)

    arg_div = np.argmin(results[1])
    x_div = results[0, arg_div]
    criterion_opt = results[1, arg_div]

    return arg_div, x_div, criterion_opt

あとは今までと同一です。

使ってみる

以上をクラスにしてまとめ，実行してみました。コードはGitHubにあげてあります。

github.com

boston = datasets.load_boston()
X = boston.data
y = boston.target

from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y)

tree = MyTree(criterion='mse', max_depth=5, verbose=True)
tree.fit(X_train, y_train)

y_pred = tree.predict(X_test)
1 - np.sum((y_test - y_pred) ** 2) / len(y_test) / y_test.var()
# 0.7113242620957975

モデルの評価には，決定係数を用いました：

$R^{2} = 1 - \frac{\sum _ {i=1} ^{N} (y_i - y_i^{\rm (pred)})^{2}}{\sum _ {i=1} ^{N} (y_i - y_i^{\rm (mean)})^{2}}$

データ自身のばらつきと比べて，予測値と実際の値のズレが小さいほど，この値は1に近い値となります。

今回作成したモデルでは決定係数は0.71でした。モデルの誤差がデータ自身のばらつきの3割程度に抑えられているので，パラメータチューニングを行なっていないにも関わらずそれなりによい値のように思えます。

まとめ

回帰についても，評価指標の部分を変えるだけで実装することができました。

なお，評価指標の部分で，criterion_lr(y_left, y_right, func_criterion)という関数を導入しました。ここの部分の引数func_criterionに別の関数を入れてやることで，別の指標を使うことができます。

例えば，Gini係数を下記のように定義して代入すると，左右の分割断片断片でのGini係数の平均値（データ数で重み付き）を出すことができます。

def gini(x):
    _, counts = np.unique(x, return_counts=True)
    prob = counts / len(x)
    
    return 1 - (prob * prob).sum()

func_criterion = gini
criterion_divide = criterion_lr(y_left, y_right, func_criterion)

このようにして，異なる指標でも引数に入れる関数を変えるだけで共通のコードを使い回せるようにしてあります。

なお，MyTreeクラスにまとめ直す際に，引数にcriterionを加えました。これをginiもしくはmseにすることで，分類と回帰を切り替えられるようにしてあります。

今回は，決定木による回帰を，以前作成した分類との共通部分を意識しながら実装していきました。

しかしscikit-learnの実装では，分類はDecisionTreeClassifier，回帰はDecisionTreeRegressorと，異なるインターフェイスに分けています。とはいっても，共通する部分は継承をすることで，重複したコードが少なくなるように実装されているはずです。

自己流での実装もひと段落しましたので，次回からはいよいよscikit-learnの中身を覗いていきたいと思います。

*1:だんだん冒頭ネタもいらすとも苦しくなってきた…

*2:Gini係数の他に，情報エントロピー $-\sum_{k=1}^{K}p_k \log p_k$ を使う場合もあります。

*3:絶対値和を用いることもあります。

2019-03-17

木構造を整える（ジコリューで行こう！#番外編）

ジコリューで行こう！決定木機械学習 Python

こんにちは，木木木林（きききりん）です。

f:id:tatamiyatamiyatatatamiya:20190317194149p:plain

tatamiya.hatenablog.com

前回までの記事で一通りランダムフォレストの実装を行いました。

予告では今回は回帰を扱うとしていましたが，準備している最中に，分割の情報を記録するノードの構造をもう少し綺麗にできることに気がつきました。

そこで本記事では，今後の拡張性も意識して一般化つつ，木構造の実装整理をしていきたいと思います。

指針

前回までのあらすじと今回やること

前々回の記事では，決定木学習の際に各ノードで探した最適分割の情報は，ノードクラスインスタンスに保存していました。

tatamiya.hatenablog.com

その際に，中間ノードであればnode_internalクラス，葉ノードであればnode_leafクラスと２種類に分けて作成しました。それぞれのノードの役割は以下の通りです：

node_internal : 分割に用いる特徴量と閾値を保存
node_leaf : このノードにたどり着いたデータに割り当てるクラスラベルを保存

ノードにはそれぞれ固有の通し番号i_nodeを振っていました。各中間ノードには，次の行き先となる子ノードの番号も記録してあります。そして作成した全ノードは，keyをi_nodeとした辞書型オブジェクトdict_nodesとして保存していました。

予測を行う時には，i_node=0から順にデータを振り分けつつ，dict_nodeを参照しながら下流のノードへと進んでいきます。

最終的に葉ノードにたどり着けば予測値が得られますが，この際の葉ノード判定はクラスの名前を直接参照して行なっていました。

if node_next.__class__.__name__ == 'node_leaf':
    print(node_next.k_decided)

ノードクラス１種類だけでまとめられないか？

以上の方法では，以下の点で綺麗とは言い難い実装でした：

２種類のノードクラスを定義しなくてはいけない
中間ノード・葉ノードの区別はクラス名で無理やり行なっている
子ノードを参照する際にわざわざdict_nodesを解さなくてはいけない

これに対し，以下のようにノードクラスの構造を変えることで，よりスッキリした実装が実現しました：

ノードクラスを一つにまとめる
- 葉ノードかどうかの情報は，フラグ化してプロパティとして持っておく
子ノードの情報を取得する際は，ノード番号ではなくノードそのものを直接参照する

後者についてはイメージがつきにくいかもしれませんが，ノードクラスのインスタンスに別のノードクラスインスタンスを保存するという入れ子構造を取らせる形になります。

実装

実際に実装したものを見てみましょう。

class Node():
    
    def __init__(self, i_node, depth, criterion, value):

        self.i_node, self.depth = i_node, depth
        self.criterion, self.value = criterion, value
        
        self.is_leaf = False
        
        self.division = {}
        self.child = {0:None, 1:None}
        
    def set_division(self, i_feature, threshold, criterion_diff):
        
        self.division['i_feature'] = i_feature
        self.division['threshold'] = threshold
        self.division['criterion_diff'] = criterion_diff

以前中間ノード・葉ノード別々に作ったクラスと比べて，かなりスッキリしているかと思います。

大きな変更点は，

葉ノード判定用のis_leafプロパティを作成した
- デフォルトはFalseですが，葉ノードにすると決めた時点でTrueに書き換えます
分割の情報はdivisionプロパティに辞書型にしてまとめておく
- 学習の際には，set_divisionメソッドで一括して情報を入れられるようにしてあります。

そのほか，学習指標criterionも持っておくことにした。ここでいうcriterionは，今までいうと不純度の評価に使っていたGini係数です。今回このように一般化した表記に置き換えた箇所がいくつかありますが，Entropyなど別の指標への変更を見据えてこのようにしました。特に回帰だとGini係数が使えず平均二乗誤差(Mean Squared Error)や平均絶対誤差(Mean Absolute Error)などの指標に置き換えなくてはなりません。

また，今回の記事の内容では保存しておく必要はないものの，学習前後での指標値の差criterion_diffも記録します。これは特徴量の重要性を評価するfeature_importanceを計算する際に必要になります。 feature_importanceをはじめ特徴量の重要度や説明性については，また別の機会に解説したいと思います。

葉ノードor中間ノードの決定

今までのアルゴリズムは，葉ノードか中間ノードかはインスタンス生成時に決めていましたが，今回は後から決めます。

まず，今まで通りgo_on_division関数内でデータの最適な分割を探索します。そして，分割後の指標（今まででいうGini係数）の値をもとにフラグを書き換えます。

criterion_initial = node.criterion
arg_div, x_div, criterion_optimized = divide(X, y)
    
criterion_diff = criterion_initial - criterion_optimized
    
if criterion_diff < threshold_criterion:
    node.is_leaf = True

今までは分割後の値そのものに閾値を与えて分割を止めるか判断していましたが，より一般的な実装に近づけるために分割前後の指標値差で判断するように変更しました。これにより，指標自体が比較的大きい値で停留してしまう場合にも分割停止判断を下せるようになります。

子ノードの作成

上記の分割停止基準を満たなかった場合は，分割に用いた特徴量と閾値を保存し，子ノードを作成します。

node.set_division(arg_div, x_div, criterion_diff)

depth_next = node.depth + 1
list_divided = [(X_left, y_left), (X_right, y_right)]
for lr, divided in enumerate(list_divided):
    i_node += 1

    X_i, y_i = divided

    node_next = make_new_node(i_node, depth_next, y)
    node.child[lr] = node_next

ここで，あらかじめ定義しておいたmake_new_node関数を使って，子ノードの原型となる新たなNodeインスタンスnode_nextを作成しています（詳細は割愛します）。

注目すべきは，最後の行node.child[lr] = node_nextです。ここではNodeインスタンスnodeのchildプロパティに直接新しいインスタンスを入れています。

今までの実装では，

今いるnodeのプロパティに子ノードの番号を保存する
辞書型オブジェクトdict_nodesに子ノードの番号と新しく作成した子ノードインスタンスを保存する

の二段階必要でしたが，Nodeインスタンスに直接Nodeインスタンスを入れることで一回で済むようになりました。

あとは木の深さが閾値以下であれば再び分割を繰り返す，一定値に達したら先ほどと同様にis_leafフラグをTrueにして止める，という流れになります。

予測

Nodeクラスを整理したおかげで，予測の部分も比較的簡潔に書けるようになりました。

def pred_each_vector(x, node):

    node_current = node

    while node_current.is_leaf == False:
        division = node_current.division
        lr = int(x[division['i_feature']] > division['threshold'])
        node_current = node_current.child[lr]
    
    return node_current.value

一番最初のノード（根ノードともいう）から出発して，

分割の情報をnode_current.divisionにより取り出し，左右どちらの子ノードにデータを振り分けるか判断する
node_current.childから次のノードを取り出し，node_currentに置き換える

これをcurrent_node.is_leafがTrueになるまで続け，最後にnode_current.valueで葉ノードの出力値を取得すれば終了です。

まとめ

ノードインスタンスのプロパティに新しいノードインスタンスを入れて入れ子にしていくことで，比較的スッキリした木構造を作ることができました。

回帰の実装は次回行いたいと思います。

2019-03-10

木を見たので森を見る：ランダムフォレストを作る（ジコリューで行こう！#4）

ジコリューで行こう！ランダムフォレスト決定木機械学習 Python

御機嫌よろしくて，森林門郎（もりりんもんろう）ですの。*1

f:id:tatamiyatamiyatatatamiya:20190209162707p:plain

tatamiya.hatenablog.com

前回までの記事では，決定木を作りました。

今回は木をたくさん作って多数決をとるランダムフォレストを作っていきます。

指針

前回までのあらすじと今回やること

決定木の実装は，GitHubの方に上げてあります。*2

こちらを使って，以下のように学習・予測を行うことができます。

import numpy as np
from sklearn import datasets
from sklearn.model_selection import train_test_split

from tree_modules.tree_base import MyTree


iris = datasets.load_iris()
X = iris.data
y = iris.target

X_train, X_test, y_train, y_test = train_test_split(X, y)

tree = MyTree()
tree.fit(X_train, y_train)
'''
=== node 0 (depth 1): arg_div -> 2, x_div -> 1.9 ===
=== node 2 (depth 2): arg_div -> 2, x_div -> 4.8 ===
=== node 3 (depth 3): arg_div -> 3, x_div -> 1.6 ===
=== node 6 (depth 3): arg_div -> 3, x_div -> 1.7 ===
'''

y_pred = tree.predict(X_test)
(y_pred == y_test).sum() / len(y_test)
'''
0.9473684210526315
'''

今回は，ランダムに条件を変えて作った複数の決定木を集め，その中で多数決をとります。

一本一本の木をどう作るか？

ただ単に同じデータ・同じ学習方法で決定木を作ると，何回やっても同じ答えが返ってくる木しか作れないので，多数決をとる意味がありません。

したがって，少しずつ違う性質を持つ木を作る必要があります。

そのために，ランダムフォレストでは以下の方法を用います。

使うデータをランダムに選択する
木を分割する際の特徴量選択にランダムネスを加える

このうち1.については，N個のデータの中から重複を許してm個のデータをランダムにとってきて新しいデータセットを作成します。このような方法でデータセットを複製し統計量をとる手法を，ブートストラップサンプリングと呼びます。

２番目の分割特徴量についてですが，今までは各ノードで分割を行う際に，全特徴量ひとつひとつ調べて最も不純度が小さくなるものを探してきました。これに対してランダムフォレストでは，全部の特徴量は使わずにノードごとにランダムに選択し，ランダムサブセット（部分集合）の中から不純度が小さくなるものを探します。

これにより出来上がった一本一本の木の精度は，毎回全特徴量を使って探索した場合と比べて劣ってしまいます。しかし，ランダムよりはマシだけど少し精度の低いモデル（弱学習器）をたくさん集めて多数決をとることで全体の精度を高めるのが，アンサンブル学習の考え方です。

実装

データの復元抽出ランダムサンプリング（ブートストラップ）

numpyのrandom.choiceメソッドをreplace=Trueに設定して使えば，重複を許しつつランダムにデータをとってくることができます。

index_resample = np.random.choice(len(y), len(y), replace=True)
X_resample, y_resample = X[index_resample], y[index_resample]

特徴量の選出

numpyのrandom.choiceは使用特徴量の選択を各ノードで行う際にも使えます。

num_features = X.shape[1]
max_features = int(np.sqrt(num_features))

index_feat_chosen = np.random.choice(num_features, max_features, replace=False)

X_chosen = X_resample[:, index_feat_chosen]

今回は一度に同じ特徴量を複数選択しないように，日復元抽出replace=Falseとしました。ただ，replace=Trueにして選んでから重複を省くことで，指定した数より少ない特徴量から最適分割を探れるようにすることもできます。

なお，選択する特徴量の数ですが，一般には全特徴量数の平方根以下が好ましいとされているので，max_features = int(np.sqrt(num_features))としてあります。

決定木にランダム分割を組み込む

決定木の際は，各ノードでの分割をgo_on_dividing(X, y)という関数を再帰的に呼び出して行なっていました。

ランダムフォレスト用にこの部分を改造して，オプションに応じて全特徴量使うかランダム選択した中から分割するかを選べるようにします。

前回まで作った決定木はクラスに書き換えてあるので，go_on_dividing関数の代わりにMyTreeクラスの_go_on_dividingメソッドを書き換えます。中身は基本的に同じですが，global変数の代わりにself.から始まるクラス変数を使っています。

 def _go_on_dividing(self, X, y, depth=0):

     depth += 1
        
     if self.splitter == 'best':
         X_chosen = X
         index_feat_chosen = np.arange(X.shape[1])
        
     elif self.splitter == 'random':
            
         n_features = X.shape[1]
         num_feat_chosen = int(np.sqrt(n_features))
         index_feat_chosen = np.random.choice(n_features, num_feat_chosen, replace=False)
         X_chosen = X[:, index_feat_chosen]

     arg_div_tmp, x_div = self._divide(X_chosen, y)
     arg_div = index_feat_chosen[arg_div_tmp] # inevitable in case of RF

     # 以下省略

後述するMyTreeクラスインスタンスを作成する際にsplitterという引数をオプションに加え，全探索なら'best'，ランダム選択なら'random'を選べるようにします（defaultは'best'）。

学習と予測

学習の際には，こうしてランダム分割できるようにしたMyTreeを複数作成し，リストlist_treesに格納していきます。

n_estimatores=10

list_trees = []
length_data = len(y)
for i in range(0, n_estimators):
    
    index_resample = np.random.choice(length_data, length_data, replace=True)
    X_resample, y_resample = X[index_resample], y[index_resample]
    
    a_tree = MyTree(splitter='random', max_features=None)
    a_tree.fit(X_resample, y_resample)
    list_trees.append(a_tree)

予測の際は，list_treesから学習済みの木を一本ずつ取り出してデータを入れ，返ってきた予測ラベルをリストに入れていきます。

list_pred = []
for a_tree in list_trees:
    
    list_pred.append(a_tree.predict(X))

array_pred = np.array(list_pred)

あとで使いやすいように，numpy arrayとして持っておきます。

多数決

numpynのbincount関数とargmaxメソッドを使うことで，それぞれの木が出した予測ラベルの中で多数決を行います。

わかりやすくするために，まずは一件目の入力データに対する各木の予測ラベルを見てみます：

array_pred[:,0]
# array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

今は木の数n_estimatorsを10に設定しているので，要素数10の一次元arrayとなっています。これをnp.bincount関数に入れることで，各ラベルの投票数が返ってきます。

np.bincount(array_pred[:,0], minlength=3)
# array([10,  0,  0])

minlenghth=3としているのは，予測ラベルの数が３種類(setosa, versicolour, virginica)あるからです。今回のように投票数が0のラベルが存在する場合，この指定をしないと何番目のクラスの投票数なのかわからなくなります。 *3 この中から最も投票数の多かったクラスラベルを得るには，.argmax()をとります。

今の操作を全データでいっぺんに行うために，今までも何回か登場したnumpyのapply_along_axis関数を使います。

pole_result = np.apply_along_axis(func1d=np.bincount, axis=0, arr=array_pred, minlength=3)

axis=0とすることで，各列について行方向にnp.bincountを作用させます。

さらにもう一度，axis=0方向にargmax()をとることで，最も投票数の多かったラベルを全データで得ることができます：

pole_result.argmax(axis=0)

なお，今回は最多投票数のラベルを返すようにしましたが，各ラベルへの投票率で返すこともできます：

(pole_result / n_estimators).T

まとめ

以上をまとめて，MyForestクラスにしてみました。

X_train, X_test, y_train, y_test = train_test_split(X, y)

myrf = MyForest(n_estimators=100)
myrf.fit(X_train, y_train)

y_pred = myrf.predict(X_test)

y_pred
# array([1, 1, 2, 0, 0, 2, 0, 0, 0, 0, 0, 1, 0, 0, 2, 1, 1, 1, 1, 1, 0, 1,
#        2, 2, 2, 1, 1, 2, 2, 2, 2, 2, 2, 0, 2, 2, 2, 1])

y_test
# array([1, 1, 2, 0, 0, 2, 0, 0, 0, 0, 0, 2, 0, 0, 2, 1, 1, 1, 1, 1, 0, 1,
#        2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 0, 2, 2, 2, 1])

(y_pred == y_test).sum() / len(y_test)
# 0.9473684210526315

trainとtestの分け方にも依存しますが，比較的安定して精度が出るようです。

決定木さえできれば，ランダムフォレストは比較的簡単に作ることができました。

今まではクラス分類を行ってきましたが，決定技・ランダムフォレストは連続値を持つデータに対して回帰を行うこともできます。

次回は，回帰の実装を行いたいと思います。

*1:筆者が昔通っていた塾に，森先生という人がいた。ある日，廊下の黒板に「森林門郎ケツでかい！！！！」とデカデカと落書きがされていた。どちらかというと痩せ身の先生だったのだが…

*2:前回までの内容から，一部改変が含まれますがご容赦願います。

*3:今回の例でいうと，minlengthの指定がない場合array([10])のように出力される。しかし2番目のクラスのみ，もしくは３番目のみに票が集中した場合にも同様に出力されてしまう。minlength=3としておけば，それぞれarray([0,10,0]), array([0,0,10])と出力されるので，どのクラスに票が入っているのか判別することができる。

2019-03-03

木を育てるその3：正解ラベルなしのデータを入れて予測を行う（ジコリューで行こう！#3）

この木なんの木？気になる木。とっても不思議な木ですから。*1

f:id:tatamiyatamiyatatatamiya:20190303214003p:plain

tatamiya.hatenablog.com

前回・前々回に引き続き，決定木の実装を行なっていきます。

前回では既知のデータと正解ラベルを元に学習を行う部分ができましたので，今回は新しいデータを入れて予測する部分を作っていきたいと思います。

指針

前回までのあらすじと今回やること

関数go_on_dividingに特徴量ベクトルXと正解ラベルyを入れることで，木を深めながら最適な分割を探索していきます。

i_node = 0
go_on_dividing(X, y)
'''
=== node 0 (depth 1): arg_div -> 2, x_div -> 1.9 ===
[(1, 0, 2, 1.9, 0)] 0
=== node 2 (depth 2): arg_div -> 3, x_div -> 1.7 ===
=== node 3 (depth 3): arg_div -> 2, x_div -> 4.9 ===
[(1, 0, 2, 1.9, 1), (2, 2, 3, 1.7, 0), (3, 3, 2, 4.9, 0)] 1
[(1, 0, 2, 1.9, 1), (2, 2, 3, 1.7, 0), (3, 3, 2, 4.9, 1)] 2
[(1, 0, 2, 1.9, 1), (2, 2, 3, 1.7, 1)] 2
'''

これにより，

どの特徴量を，
どの値で分割し，
どこまで分割を繰り返せば良いか

が決まりましたので，今度はその情報を元に，正解ラベルなしの特徴量ベクトルを入れて，正解ラベルを予測したいと思います。

各ノードでの分割・予測の情報をどう持つか？

データの分割，もしくはクラス判定を行う部分を「ノード」と呼ぶことにしましたが，特に前者を「中間ノード」，後者を「葉ノード」と呼ぶことにします。

予測を行うにあたっては，各ノードにおいて以下の情報が必要になります。

中間ノードの場合

自身のノード番号
分割に用いる特徴量のindex
分割の閾値
True, Falseそれぞれの場合の次のノード番号

葉ノードの場合

自身のノード番号
クラスラベル（予測結果）

ここでは，ノードクラスを新たに作成してこれらの情報を格納していきたいと思います。

実装

ノードクラスを定義する

以下のようにして，中間ノードクラスnode_internalと葉ノードクラスnode_leafを作成しました。

両ノードに共通する部分はnode_basisとして定義し，継承させてあります。

class node_basis():
    
    def __init__(self, i_node, depth):
        
        self.i_node = i_node
        self.depth = depth


class node_internal(node_basis):
    
    def __init__(self, i_node, depth, i_feature, threshold):
        
        super().__init__(i_node, depth)
        self.i_feature = i_feature
        self.threshold = threshold
        
        self.node_child = {0:None, 1:None}
        

    def set_node_child(self, lr, i):
        self.node_child[lr] = i


class node_leaf(node_basis):
    
    def __init__(self, i_node, depth, k):
        
        super().__init__(i_node, depth)
        self.k_decided = k

上記では挙げませんでしたが，念の為木の深さdepthも持たせてあります。

中間ノードクラスnode_internalでは，分割する特徴量のindexをi_feature，分割の閾値をthresholdとしています。

そして，子ノードの番号をnode_childとして辞書型にして持たせています。 x[i_feature] <= thresholdならば次はnode_child[0]番目のノードに， x[i_feature] > thresholdならばnode_child[1]番目のノードに行きます。

これらは，学習の際にset_node_childメソッドを使って記録するようにします。

学習結果をノードクラスに記録していく

上述のように作成したノードクラスを，go_on_dividing関数に組み込んでいきます。

def go_on_dividing(X, y, depth=0,
                   threshold_gini=0.05, min_node_size=5, max_depth=3):
    
    global i_node, dict_nodes
    depth += 1

    arg_div, x_div = divide_tree(X, y)
    node_current = node_internal(i_node, depth, arg_div, x_div)
    dict_nodes[i_node] = node_current

    print("=== node {} (depth {}): arg_div -> {}, x_div -> {} ===".format(i_node, depth, arg_div, x_div))

    mask = X[:, arg_div] > x_div
    X_right, X_left = X[mask], X[~mask]
    y_right, y_left = y[mask], y[~mask]

    gini_left = gini(y_left)
    gini_right = gini(y_right)

    list_divided = [(X_left, y_left, gini_left), (X_right, y_right, gini_right)]

    for lr, divided in enumerate(list_divided):
        i_node +=1

        X_i, y_i, gini_i = divided
        if gini_i > threshold_gini and len(y_i)>min_node_size and depth+1 <= max_depth:

            node_current.set_node_child(lr, i_node)
            go_on_dividing(X_i, y_i, depth=depth)
        else:
            node_current.set_node_child(lr, i_node)
            feature_majority = np.bincount(np.array(y_i)).argmax()

            node_terminal = node_leaf(i_node, depth, feature_majority)
            dict_nodes[i_node] = node_terminal

大きく変わったのは，今までdiv_setにリストとして記録していた分割の過程を，辞書型のグローバル変数dict_nodesに入れるようにした部分です。 keyにはノード番号を，valueにはノードクラスのインスタンスを入れています。

このdict_nodesに新しいノードを加えていくのは，次のタイミングです。

1: 新しい分割を行った時

arg_div, x_div = divide_tree(X, y)
node_current = node_internal(i_node, depth, arg_div, x_div)
dict_nodes[i_node] = node_current

2 : 分割を止め，葉ノードを作成する時

feature_majority = np.bincount(np.array(y_i)).argmax()

node_terminal = node_leaf(i_node, depth, feature_majority)
dict_nodes[i_node] = node_terminal

予測を行う

これで学習結果を記録できるようになりました。実際にやってみますと，

i_node=0
dict_nodes = {}
go_on_dividing(X, y)
'''
=== node 0 (depth 1): arg_div -> 2, x_div -> 1.9 ===
=== node 2 (depth 2): arg_div -> 3, x_div -> 1.7 ===
=== node 3 (depth 3): arg_div -> 2, x_div -> 4.9 ===
'''

こうすると，グローバル変数として定義したdict_nodesは以下のようにkeyとvalueが追加されます。

dict_nodes
'''
{0: <__main__.node_internal at 0x1222dcfd0>,
 1: <__main__.node_leaf at 0x1222dccf8>,
 2: <__main__.node_internal at 0x1222dc9b0>,
 3: <__main__.node_internal at 0x1222dc908>,
 4: <__main__.node_leaf at 0x1222dcc88>,
 5: <__main__.node_leaf at 0x1222dc828>,
 6: <__main__.node_leaf at 0x1222dceb8>}
'''

中間ノードであれば，以下のようにして

分割する特徴量のindexと閾値
子ノードの番号

を確認できます：

dict_nodes[2].i_feature, dict_nodes[2].threshold
'''
(3, 1.7)
'''

dict_nodes[2].node_child
'''
{0: 3, 1: 6}
'''

葉ノードであればk_decidedで割り当てられたクラスを確認できます：

dict_nodes[1].k_decided
'''
0
'''

これを使って，新しいデータが入ってきた場合は以下のように予測を行います。現在いるノードnode_currentが中間ノードだとして，

x[node_current.i_feature] > node_current.thresholdを判定
- Trueなら，i_node_next = node_current.node_chile[1]番目のノードへ
- Falseなら，i_node_next = node_current.node_chile[0]番目のノードへ
node_next = dict_nodes[i_node_next]により次のノードを取得
- 中間ノードなら，node_currentをnode_nextに置き換えて上記を繰り返す
- 葉ノードならnode_next.k_decidedを予測結果として出力する

ここで問題になるのが，node_nextが中間ノードか葉ノードかをどうやって判定するかです。 pythonではクラスの名前を.__class__.__name__によって取得することができるので，これで判定してしまいます。

if node_next.__class__.__name__ == 'node_leaf':
    print(node_next.k_decided)

以上を実装すると，以下のようになります：

def pred_at_node(x, node):
    lr = int(x[node.i_feature] > node.threshold)
    i_node_next = node.node_child[lr]
    
    return i_node_next

def pred_each_vector(x, dict_nodes):
    node_current = dict_nodes[0]
    while True:
        node_next = dict_nodes[pred_at_node(x, node_current)]
        if node_next.__class__.__name__ == 'node_leaf':
            return node_next.k_decided
        else:
            node_current = node_next

試しに，特徴量ベクトルXから0番目のデータを取り出して適用してみます。

x_sample = X[0]
pred_each_vector(x_sample, dict_nodes)
'''
0
'''

さらに全データにこの関数を適用するために，前々回にも出てきたnumpyのapply_alogn_axis関数を使います。

y_pred = np.apply_along_axis(func1d=pred_each_vector, axis=1, arr=X, dict_nodes=dict_nodes)

ただし，ここでは学習に使ったのと同じデータで予測を行なっていることに注意してください。本来は学習データとは異なるデータで予測を行うべきであり，ここで精度が高かったとしてもあまり意味はありません。

まとめ

以上までの学習・予測の流れを整理して，MyTreeクラスにまとめました。詳細はGitHubにJupyter notebookをアップしましたのでそちらを参考にしてください。

これを使って，データを訓練用・テスト用に分けて学習し，精度を見てみたいと思います。

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, y)

tree = MyTree()
tree.fit(X_train, y_train)
'''
=== node 0 (depth 1): arg_div -> 2, x_div -> 1.9 ===
=== node 2 (depth 2): arg_div -> 2, x_div -> 4.8 ===
=== node 3 (depth 3): arg_div -> 3, x_div -> 1.6 ===
=== node 6 (depth 3): arg_div -> 3, x_div -> 1.7 ===
'''

y_pred = tree.predict(X_test)
(y_pred == y_test).sum() / len(y_test)
'''
0.9473684210526315
'''

なんだか機械学習やってる感じが出てきましたね。

以上で決定木による分類の実装は一通りできましたので，次回は複数の決定木を構築して多数決を行うランダムフォレストの実装を行いたいと思います。

*1:次の時代に新しい課税を。MISOJI Expire the Next!

2019-02-24

木を育てるその2：分割を繰り返して適当なところで止める（ジコリューで行こう！#2）

こんにちは，植木等 *1です。

わかっちゃいるけどやめられない。

f:id:tatamiyatamiyatatatamiya:20190224150754p:plain

この記事では前回に引き続き，決定木を自己流で実装していきたいと思います。

tatamiya.hatenablog.com

前回は不純度が最低になるような分割の探索を行いましたので，今回はこれを繰り返し，適当なところで止めてクラス分類を行いたいと思います。

指針

前回までのあらすじと今回やること

与えられた特徴量ベクトルXとtargetベクトルyを元に，分割する特徴量・各特徴量の値で全探索し，不純度が最低になる分割方法見つけ出すところまでを，以下の関数にまとめました。

途中の関数find_optimal_divisionは，各特徴量について最適な分割点を探索するものです。

def divide_tree(X, y):

    results = np.apply_along_axis(find_optimal_division, 0, X, y)

    arg_div = np.argmin(results[1])
    x_div = results[0, arg_div]

    return arg_div, x_div

返り値のうち，arg_divは分割する特徴量のindexを，x_divは分割の際の閾値です。

分割の指標となる不純度の計算には，Gini係数を用いました。

今度はこの最適な分割後の断片2つそれぞれに再度分割を施していき，より不純度を下げていきます。これには，先ほどの関数divide_treeを再帰的に呼び出していきます。

サンプルのデータとしては，引き続きirisを用います。

どこで分割を止めるか？

途中で不純度が十分下がりきった場合はそこで分割を止めます。

実際，一回分割した後の断片を見てみると，片方には0番目のクラス（品種名：Setosa）しか含まれていないことがわかります。

y[X[:, arg_div] <= x_div]
# array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

これ以上分割しても意味がないので，「不純度が一定値以下になったら分割をやめる」という条件を設けます。

また，データを細かく分けすぎてしまってもいけません。学習用にInputしたデータに過度に依存した分割ができてしまい，新しいデータを入れて予測を行う際に不適切な分割を行うことで精度が下がる恐れがあります。（過学習）

そのような事態を避けるためにも，

分割後の断片の大きさ
分割を行う回数（深さ）

に閾値を設け，適当なところで止めることにします。

実装

再帰呼び出しにより分割後の断片を再分割していく

まず，全体のイメージ的にはこんな実装になります。

def go_on_dividing_temp1(X, y):

    arg_div, x_div = divide_tree(X, y)
    
    mask = X[:, arg_div] > x_div
    X_right, X_left = X[mask], X[~mask]
    y_right, y_left = y[mask], y[~mask]
    
    go_on_dividing_tmp1(X_left, y_left)
    go_on_dividing_tmp1(X_right, y_right)

このgo_on_dividing_temp1関数では，最初に一回分割を行なってから，得られた特徴量が閾値より大きいか小さいかで右断片と左断片に分けます。

この際に，閾値以上のデータはTrue, それ以下のものはFalseとするbinaryのmaskを作成し，これをX, yの引数に撮りました。leftの方で~maskとしていますが，~はTrue, Falseを逆転させます。

そして，分けた断片を再度別々にgo_on_dividing_temp1関数に打ち込みます。

定義した関数の中でその関数自身を呼び出すこのやり方は，再帰呼び出しと呼ばれます。最初は何をやっているか理解が難しいかもしれませんが，何度か動かしていくうちに徐々に慣れてくると思います。

分割の停止条件を設ける

これで枝葉をどんどん伸ばしていくことができますが，このままでは分割できるデータがなくなるまで延々と分割を続けていってしまいます。

そのため，「指針」で触れた下記の３つの停止条件を付加します：

不純度が一定値以下
- 今回は各断片ごとのGini係数に閾値を設けます。
分割後の断片の大きさ
- 分割後の断片の長さlen(y_left), len(y_right)に閾値を設けます。
分割の回数
- 分割の深さをdepthとして変数に入れておき，こちらに閾値を設けます。

以上を実装すると，以下のようになります。

def go_on_dividing_temp2(X, y, depth=1,
                   threshold_gini=0.05, min_node_size=5, max_depth=3):        
    
    arg_div, x_div = divide_tree(X, y)
        
    mask = X[:, arg_div] > x_div
    X_right, X_left = X[mask], X[~mask]
    y_right, y_left = y[mask], y[~mask]
    
    gini_left = gini(y_left)
    gini_right = gini(y_right)
    
    list_divided = [(X_left, y_left, gini_left), (X_right, y_right, gini_right)]
    
    for divided in list_divided:
        
        X_i, y_i, gini_i = divided

        if gini_i > threshold_gini and len(y_i)>min_node_size and depth < max_depth:
            go_on_dividing_temp2(X_i, y_i, depth=depth+1)
        else:
            class_decided = np.bincount(y_i).argmax()

Gini係数，分割後断片のデータ長さ，分割の深さの閾値をそれぞれthreshold_gini, min_node_size, max_depthとしました。デフォルト引数でとりあえず適当な値を入れてあります。

分割の深さdepthはgo_on_dividing_temp2関数の引数として持っておき，デフォルトの値を1，呼び出す際に一つインクリメントするようにしました。

これらを合わせて，以下の条件が全部満たされている間だけ分割を続けるようにしてあります。

gini_i > threshold_gini
len(y_i)>min_node_size
depth < max_depth

y_i, gini_iには，分割後断片左右どちらかの正解クラスベクトルとGini係数が入ります。

もし３つの条件のうちどれか一つでも満たされない場合，分割をやめます。

そしてそれに合わせて，分割されたデータ断片がどのクラスに分類されるかを決めます。

ここでは，断片に含まれるクラスのうち最も多数派を占めるものを選び出します。

class_decided = np.bincount(y_i).argmax()

分割の様子を出力する

これで一通りの計算はできたのですが，このまま実行しても何も出力されません。

どのような分割を経て，
最終的にどのクラスに行き着いたか

を出力するにはどうすればいいでしょうか？

ここの部分は少し頭を使う必要があったのですが，最終的に次のようにしました：

def go_on_dividing(X, y, depth=1, div_set=None,
                   threshold_gini=0.05, min_node_size=5, max_depth=3):
    
    global i_node
    if div_set is None:
        div_set = []
    
    i_node_current = i_node
    
    arg_div, x_div = divide_tree(X, y)
    
    print("=== node {} (depth {}): arg_div -> {}, x_div -> {} ===".format(i_node, depth, arg_div, x_div))
    
    mask = X[:, arg_div] > x_div
    X_right, X_left = X[mask], X[~mask]
    y_right, y_left = y[mask], y[~mask]
    
    gini_left = gini(y_left)
    gini_right = gini(y_right)
    
    list_divided = [(X_left, y_left, gini_left), (X_right, y_right, gini_right)]
    
    for lr, divided in enumerate(list_divided):

        div_set_tmp = div_set.copy()
        div_set_tmp.append((depth, i_node_current, arg_div, x_div, lr))
        i_node +=1
        
        X_i, y_i, gini_i = divided
        if gini_i > threshold_gini and len(y_i)>min_node_size and depth < max_depth:
            go_on_dividing(X_i, y_i, depth=depth+1, div_set=div_set_tmp)
        else:
            class_decided = np.bincount(y_i).argmax()
            print(div_set_tmp, class_decided)

分割ごとに固有の通し番号i_nodeを振り，グローバル変数として持っておくことにしました。 Node（ノード）というのは，「節」や「結節点」という意味で，ここでは分割・もしくはクラス判定を行う部分を指します。

さらに，分割の経過をリストdiv_setに保存していきます。この中には，

何回目の分割(depth)の，
何番目のノード(i_node)で，
どの特徴量index(arg_div)を
閾値いくつ(x_div)で分割して，
左側・右側どちらの断片に入ったか？(lr)

の情報を入れていき，go_on_dividing関数の引数に入れて下流のノードに伝達していきます。

最終的に分割を止める際には，このdiv_setと，分類後のクラス番号を表示します。

これらを実行すると以下のようになります。

i_node = 0
go_on_dividing(X, y)
'''
=== node 0 (depth 1): arg_div -> 2, x_div -> 1.9 ===
[(1, 0, 2, 1.9, 0)] 0
=== node 2 (depth 2): arg_div -> 3, x_div -> 1.7 ===
=== node 3 (depth 3): arg_div -> 2, x_div -> 4.9 ===
[(1, 0, 2, 1.9, 1), (2, 2, 3, 1.7, 0), (3, 3, 2, 4.9, 0)] 1
[(1, 0, 2, 1.9, 1), (2, 2, 3, 1.7, 0), (3, 3, 2, 4.9, 1)] 2
[(1, 0, 2, 1.9, 1), (2, 2, 3, 1.7, 1)] 2
'''

ノード番号i_nodeはグローバル引数で定義したので，呼び出しの前に初期化しておく必要があります。

なお，リストdiv_setを更新する際に，

div_set_tmp = div_set.copy()
div_set_tmp.append((depth, i_node_current, arg_div, x_div, lr))

のように，いったんコピーdiv_set_tmpを作成してから更新・go_on_dividing関数に渡しています。

これをせずにdiv_setにそのままappendして引数に渡してしまうと，全ノードでこのリストが共有されてしまいます。そのため，最終的な分類結果を表示させるときに，自分がたどっていないノードでの分類過程まで出力されてしまいます。

ノード番号を振る際にグローバル変数を使っている点が気になりますが，クラスにしてまとめる際にうまく処理できればと思います。

なお，今回の実装ではノード番号は計算を行なった順に振っていますが，同じ深さのノード間では左から順に1ずつ増えていくように振る方がメジャーなようです（参考：みんな大好きはじめてのパターン認識第11章）*2。これに従ってノード番号を振るにはどうすれば良いかは，現在考え中です。

ともあれ以上で一通り決定木の学習のプロセスが実装できました。次回はこれをもとに新しいデータをいれてクラスラベル（品種）を予測できるようにしたいと思います。

*1:1960年代にヒットしたコメディーバンドハナ肇とクレイジーキャッツのボーカル。代表曲はスーダラ節，ハイそれまでョなど。2007年に死去。他のメンバーも気になって調べてみたところ，犬塚弘のみご存命(90歳)であった。

*2:個人的には数式の導入がそっけないのでそこまで好きではない。ただ，第11章の決定木〜アンサンブル学習はコンパクトにまとまっているので参考にしている。

2019-02-18

木を育てるその１：不純度を計算して最適な分割点を探す（ジコリューで行こう！#1）

f:id:tatamiyatamiyatatatamiya:20190217231702p:plain

どうも木多郎*1です。

今日は決定木を自己流で実装する第一回です。

基本的にこんな感じかな？でやっていくので，もしかしたらscikit-learnなどで用いられる一般的な実装とは，アルゴリズムからして異なっているかもしれません。

そこらへんはおいおい修正していくことにして，とりあえず実装してみます。

この記事では，Gini係数を計算して不純度を求め，最適なデータの分割方法を探る部分を実装します。

分割後のデータをさらに分割して枝葉を広げていくところは次回以降に行います。

アルゴリズムの全体像

ここでは以下のように学習を行なっていきます：

分割を行う特徴量を選ぶ
適当な位置で分割し，不純度を計算する
分割位置をずらしていって，不純度が最小になる点を見つける
全特徴量で同様に行い，不純度が最小になる特徴量・分割点の組み合わせを求める
分割後の断片で同じことを繰り返す

これを，適当な条件を満たすまで続けていきます。

停止の条件としては，とりあえず今は，

分割後の断片のデータ数が閾値より少ない
分割の深さが一定値以下
不純度が一定値以下

のいずれかを満たした場合，としようと思います。

今回扱うのは1~4の部分です。

使うデータ

無難にirisデータセットを使います。アヤメのがく片，花びらの長さと幅から品種を予測します。

読み込み

scikit-learnのサンプルのやつを読み込んできます。

from sklearn import datasets
import numpy as np

iris = datasets.load_iris()

X = iris.data
y = iris.target

データの説明

X : 各列に以下の値が入っています。

がく片長さ(cm)
がく片幅(cm)
花びら長さ(cm)
花びら幅(cm)

y: 3種類の品種がラベルづけされています。

0 -> Iris-Setosa
1 -> Iris-Versicolou
2 -> Iris-Virginica

実装

不純度の計算

一般的には不純度の計算にはentropyとGini係数が使われますが，とりあえず今回はGini係数の方でやってみます。

Gini係数の定義：

$\sum _ {k=1} ^{K} (1-p_k) p_k = 1 - \sum _{k=1} ^{K} p_k ^2$ 　*2

全K個あるクラスのうちk番目のクラスの確率を $p _ k$ と書きます。

この値は， $p_k$ の値が全てのkで等しい時（ $p_k=1/K$ ）に最大値を取り，特定のクラスしか存在しない場合0を取ります。

したがって，データセットを二分割した場合，Gini係数が小さいほど特定のクラスがマジョリティを占める綺麗な分割であることを意味します。

これを以下のように実装します。

def gini(y):
    _, counts = np.unique(y, return_counts=True)

    prob = counts / len(y)
    
    return 1 - (prob * prob).sum()

例えば，がく片幅が2.5cm以上かどうかで分割する場合，得られる不純度は0.637となります。

mask_divide = X[:, 1] > 2.5
y_upper = y[mask_divide]
y_lower = y[~mask_divide]

gini_divide = (gini(y_upper) * len(y_upper) + gini(y_lower) * len(y_lower)) / len(y)

二分割した２つの切片でそれぞれGini係数を計算し，データ数で重み付け平均を取っています。

分割しない場合の値0.666...と比べると，ほんの少ししか下がっていませんね。

分割点を動かす

上の例ではX[:, 1] > 2.5と適当に分割点を設けていましたが，今度は分割点をシステマティックに動かしながら同様に不純度を計算し，最小値をとる分割位置を探索します。

探索範囲をどう決めるかですが，適当なビン幅を決めて等分割していってもいいのですが，今回はデータの値そのものを使ってしまいます。

まずデータセットに入っているがく片幅の値X[:, 1]をユニーク化して，それぞれの値を分割の閾値にして不純度を計算していきます。

def find_optimal_division(x, y):
    list_gini = []
    x_unique = np.unique(x)

    for threshold in x_unique:

        mask_divide = x > threshold
        y_upper = y[mask_divide]
        y_lower = y[~mask_divide]

        gini_divide = (gini(y_upper) * len(y_upper) + gini(y_lower) * len(y_lower)) / len(y)

        list_gini.append(gini_divide)
        
    array_gini = np.array(list_gini)
    i_div_opt = np.argmin(array_gini)
    
    return x_unique[i_div_opt], array_gini[i_div_opt]

これにより，不純度が最小となる分割点とその時の不純度が得られました。

今回は，がく片幅3.3cmで分割するのが最適で，その際の不純度は0.540となりました。先ほどの適当な分割よりマシになりましたね。

find_optimal_division(X[:, 1], y)

全特徴量で調べる

上記ではがく片幅で分割してみましたが，今度は他の特徴量でもみていきましょう。

Pythonのリスト内包表記を使えば一行で書けます。

[find_optimal_division(X[:, i], y) for i in range(0, X.shape[1])]

これにより，下記のような結果が得られました：

[(5.4, 0.4389063317634746),
 (3.3, 0.539743283106115),
 (1.9, 0.3333333333333333),
 (0.6, 0.3333333333333333)]

花びらの長さか幅で分割すると最も不純度が下がるようです。

numpyのapply_along_axis関数を使って下記のように書き直し，不純度が最小になる特徴量（のindex）と分割の閾値を変数に格納しておきます。

results = np.apply_along_axis(find_optimal_division, 0, X, y)

arg_div = np.argmin(results[1])
x_div = results[0, arg_div]

出力：

arg_div, x_div
#(2, 1.9)

花びらの長さと幅が同率一位ですが，便宜上先にヒットする長さの方が選択されます。

今日のところはここまで。次回は分割を繰り返して枝を伸ばしていきます。

*1:北海道札幌市澄川にある老舗スープカレー屋の一つ。木造の店で野郎がカレーを作る店だからこの名前にした，とか。

*2:~~TeX表記がうまくできない。。。~~直せました！（参考）