単語分散表現(Skipgram, CBoW)のベクトル長について

Skipgram など単語の分散表現のベクトル長について

Skipgram, CBoW などの単語分散表現は、単語をベクトル空間に配置（埋め込む）するのですが、単語同士の類似度としてはコサイン類似度が用いられます。
言い換えると、単語ベクトル同士の角度が類似度となります（類似した単語ベクトルの角度は小さくなる）。

しかし、Gensim の Word2Vec model にはベクトル .syn0 と L2norm の .syn0norm が用意されています。
この使い分けの観点は何か？という問いが stackoverflow にあり、そこから「ベクトルの長さはどういう意味を持つのか」という考察が書かれています。

「non-normalized vector と normalized vector の使い分けがわからない問題について.」
stackoverflow.com
stats.stackexchange.com

端的に言うと

ベクトル長は異なるコンテキストに出現する頻度を反映している。 - 異なるコンテキストに何回も出現する単語は短いベクトルとなる - 少ない出現頻度の言葉は長いベクトルとなる

ということで、これは何を意味しているかと言うと、

単語ベクトルの方向が意味を捉えているだけでなく、長さも重要な情報を持っている。 特に、べクトルの長さは単語の出現頻度を内包しており、単語の重要度を測るのに有用である。

ということになります。
詳細はこの論文*1に詳しいですが、この単語出現数とベクトル長のプロットを見ると

単語出現数に対して、ベクトル長は短い間隔となり出現数が増えると共に徐々に上昇する、30 回以上の出現数でその傾向は減少傾向に転じる.

となっており、それに対して著者はこのような考察をしています。

単語出現数に対して、ベクトル長は短い間隔となり出現数が増えると共に徐々に上昇する、30 回以上の出現数でその傾向は減少傾向に転じる。 少ない出現頻度の言葉は首尾一貫した使われ方をする傾向があり、その結果ベクトルは長くなる。 高い出現頻度の言葉はたくさんの異なるコンテキストで使われる傾向があるため、平均処理が増えるだけ短いベクトルとなる。 この傾向は図が明確に示している。

単語が異なるコンテキストに出現したとき、ベクトルは異なる方向に動くように更新される。 従って、最終的なベクトルは異なるコンテキストでの重み付きの平均を表現したものとなる。 このように異なるコンテキストに出現する回数が増えれば、ベクトルは平均処理される回数が増えることから短くなる。 異なるコンテキストに何回も出現する単語は、少ない意味を内包するに違いない…という仮説になる。 一番良い例としては様々なコンテキストで無差別に出現するストップワードで、これらは出現頻度の割に、短いベクトル表現となる。

以上、Skipgram、CBoW など単語分散表現のベクトル長の意味についての考察でした。

*1:Measuring Word Significance using Distributed Representations of Words by Adriaan Schakel and Benjamin Wilson. https://arxiv.org/pdf/1508.02297.pdf

XGBoost における単調性制約

XGBoost の単調性制約

XGBoost などの Tree 系アルゴリズムは表現力が大きいため、容易に bias を小さくできます。

下の例では、増加するトレンドにある特徴量と目的変数がにノイズが乗ったデータに対して XGBoost でフィッティングすると、全体的には上昇トレンドを捉えますが、局所的に振動してしまいます。

これは、XGBoost が表現力（データにフィットする能力）が高いゆえにできることなのですが、一方で過度にデータにフィットしているとも言えます。

もちろん、それはモデル自身が判断することではなく、課題設定や事前情報、ビジネス上の考慮などで決まる場合が多いです。

ただ、諸理由によりこの「増加するトレンド」であったり「減少するトレンド」を制約として加え、制約の中で最大限にフィッティングさせたい場合に、この単調性の制約（単調に増加（減少）する制約）は有効です。

一般的に、このような制約を加える理由としては「精度向上のために」というよりも「解釈しやすくするため」という動機が多いように思います。

下図が単調性制約を与えてフィッティングしたモデルの結果です。 これを見ると、大局的な傾向は捉えていながら、局所的に増減する振動が消えています。

というようなことが下のページに書いてあるのですが、簡単に内容と実際のコードを書いておきます。

xgboost.readthedocs.io

例では次のデータを用いています（上図２つも同様です）

]

# jupyter notebook での実行
import numpy as np
import xgboost as xgb
import matplotlib.pylab as plt
from sklearn.cross_validation import train_test_split

%matplotlib inline

# データの作成
N = 1000
x1 = np.random.random(N)
x2 = np.random.random(N)
y = (5 * x1 + np.sin(10 * np.pi * x1)) - (5 * x2 + np.cos(10 * np.pi * x2)) + np.random.normal(0.0, 0.01)

# x1 の分布（単調増加）
plt.plot(x1, y, 'o', alpha=0.5)

# x2 の分布（単調減少）
plt.plot(x2, y, 'o', alpha=0.5)

# 学習データ
X = np.array([x1, x2]).T
X_train, X_test, y_train, y_test = train_test_split(X, y)

# 描画用
# 片方の変数を描画するときはもう片方を 0.5 で固定（影響を最小限にするため）
def plot_response_x1(model):
    
    # x1 (単調増加変数のレスポンス）
    _x1 = np.linspace(0, 1, 100)
    _x2 = np.zeros(100) + 0.5
    _X = np.array([_x1, _x2]).T
    
    _y = model.predict(_X)
    plt.plot(X_test[:,0], y_test, 'o', alpha=0.5)
    plt.plot(_x1, _y, color = 'black')

def plot_response_x2(model):

    # x2 (単調減少変数のレスポンス）
    _x1 = np.zeros(100) + 0.5
    _x2 = np.linspace(0, 1, 100)
    _X = np.array([_x1, _x2]).T
    
    _y = model.predict(_X)
    plt.plot(X_test[:, 1], y_test, 'o', alpha=0.5)
    plt.plot(_x2, _y, color = 'black')
    
# 単調性制約のない XGBoost で回帰モデルを構築します.
model = xgb.XGBRegressor()
model.fit(X_train, y_train)
plot_response_x1(model)
plot_response_x2(model)

# 単調性制約を課した XGBoost で回帰モデルを構築します.
model = xgb.XGBRegressor(monotone_constraints=(1,-1))
model.fit(X_train, y_train)
plot_response_x1(model)
plot_response_x2(model)

パラメータについて

monotone_constraints

1:単調増加制約、0:制約を加えない(デフォルト)、-1:単調減少制約、で説明変数毎にリストで設定します。例えば

• (1, 0): １つ目の説明変数に単調増加制約を与え、２つ目の説明変数には制約を与えない
• (0,-1): １つ目の説明変数に単調増加制約を与えず、２つ目の説明変数には単調減少制約を与える

となります。

tree_method

monotone_constraints を利用する場合は tree_method として次のうちのどれかを指定する必要があります。

• tree_method='exact','hist', gpu_hist

BayesSearchCV で XGBoost の early stopping 機能を使う

BayesSearchCV で XGBoost の early stopping 機能を使う

先日の BayesSearchCV イントロダクションの続きです。

XGBoost には early stopping という、

学習データとは別に、テストデータの Validation error が少なくとも <early_stopping_rounds> 回数連続して減少しないと、学習を停止する。

という機能があります。

early stopping は過学習を防ぐ学習方法として XGBoost だけでなく DeepNeuralNetwork でもよく使われる手法です。 端的にいうと、学習データだけではなくテストデータの Validation Error もモニタリングし、学習エラー(bias)と検証エラー(variance)のトレードオフが起きるところで学習を早期に終了させ、過学習を防ぐというものです。

注意すべきは、これは学習実行時のパラメータであり、いわゆるハイパーパラメータではないというところです。

というところで、sklearn では学習実行時のパラメータは fit_params という形で estimator.fit(X, y, **fit_params) の引数として学習実行時に指定されます。 （対してハイパーパラメータは GridSearch の search_param に引き渡されます）

BayesSearchCV でも学習実行時のパラメータを指定できますが、fit の引数ではなく、コンストラクタにて search_param と同様に引数で fit_params を与えます。

サンプルコードになります。

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import sys
import os
import numpy as np
import scipy as sp
import pandas as pd
from skopt import BayesSearchCV
from skopt.space import Real, Categorical, Integer
from sklearn.model_selection import train_test_split
import xgboost as xgb


# データセットを読み込み
from sklearn.datasets import load_breast_cancer
X, y = load_breast_cancer(return_X_y=True)

X_train, X_test, y_train, y_test = train_test_split(X, y, train_size=0.75, random_state=0)

estimator = xgb.XGBClassifier(
    booster='gbtree',
    silent=True,
    objective='binary:logistic',
    base_score=0.5,
    eval_metric='auc',
    n_jobs=-1
)

# Early Stopping のパラメータ設定
# eval_set に モニタリングするテストデータ・セットを設定
# early_stopping_rounds 回数メトリクスが減少（増加）しなければ学習を早期終了する
# eval_metric モニタリングする指標 - 最小化(RMSE, log loss, etc.) あるいは最大化(MAP, NDCG, AUC)は自動的に選ばれる
fit_params =  {
    'early_stopping_rounds': 5,
    'eval_metric':'error',
    'verbose':1,
    'eval_set':[[X_test, y_test]]
}

# パラメータ探索範囲設定
from skopt.space import Real, Categorical, Integer
param_grid = {
    'n_estimators': Integer(1, 4),
    'learning_rate': Real(0.01, 0.5, 'log-uniform'),
    'max_depth': Integer(0, 3)
}

clf = BayesSearchCV(
    estimator=estimator,
    search_spaces=param_grid,
    scoring='roc_auc',
    cv=3,
    n_jobs=-1,
    n_iter=10,
    verbose=1,
    refit=True,
    fit_params=fit_params # コンストラクタで設定すると estimator のフィッティング実行時に渡されます
)
def status_print(optim_result):
    '''
    ベイズ最適化のイテレーションで呼ばれるコールバック
    '''
    results    = pd.DataFrame(clf.cv_results_)
    new_params = results.tail(1).params.values[0]
    new_score  = results.tail(1).mean_test_score
    print('Model #%d/%d : ROC-AUC(newest/best) = %.4f / %.4f' % (
        len(results),
        clf.n_iter,
        new_score,
        clf.best_score_
    ))

# BayesSearchCV は sklearn の他の交差検証サーチ系クラス(RandomSearchCVなど)と同じAPIを持っています
clf.fit(X, y, callback=status_print)

print(pd.DataFrame(clf.cv_results_))

# 検証データのパフォーマンス
print(clf.score(X_test, y_test))

以上、Tips でした。

scikit-optimize の BayesSearchCV を用いたベイズ最適化によるハイパーパラメータ探索

ベイズ最適化でハイパーパラメータをチューニングする

ハイパーパラメータのチューニングではパラメータを様々に振りながら学習を繰り返し、一番パフォーマンスが良いところを探し出します。 学習時間はモデルやデータサイズなどによって異なりますが、効率よく探索するかしないかで良いパラメータを得るまでの時間が大きく変わります。 探索方法には様々あり、基本的なところではグリッドサーチ（全探索）、ランダムサーチ（ランダム探索）があります。

今回紹介するものは、ベイズ最適化の手法を用いて効率よく探索する枠組みです *1 。

BayesSearchCV

ここで紹介するのは CrossValidation で汎化精度を確かめながら、ハイパーパラメータを探索する scikit-optimize の BayesSearchCV です。

github.com

https://scikit-optimize.github.io/#skopt.BayesSearchCV

ハイパーパラメータの汎化精度は交差検証でわかります（*2）。 そのため、交差検証で汎化精度を見ながら最適なハイパーパラメータを探索します。

手順で言うと次のようになります。

1. ハイパーパラメータをランダムにピックアップする.
2. 交差検証で精度を見る
3. ベイズ最適化で次に探索するハイパーパラメータを定める（→2 へ）

これを全部やってくれるのが BayesSearchCV で、以下がサンプルコードとなります。

サンプル実装

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import sys
import os
import numpy as np
import scipy as sp
import pandas as pd
from skopt import BayesSearchCV
from skopt.space import Real, Categorical, Integer
from sklearn.model_selection import train_test_split
import xgboost as xgb


# データセットを読み込み
from sklearn.datasets import load_breast_cancer
X, y = load_breast_cancer(return_X_y=True)

X_train, X_test, y_train, y_test = train_test_split(X, y, train_size=0.75, random_state=0)

estimator = xgb.XGBClassifier(
    booster='gbtree',
    silent=True,
    objective='binary:logistic',
    base_score=0.5,
    eval_metric='auc',
    n_jobs=-1
)

# パラメータ探索範囲設定
# 探索範囲 min,max,distribution を space クラスでで指定(min, max,(distribution))
from skopt.space import Real, Categorical, Integer
param_grid = {
    'n_estimators': Integer(1, 4),
    'learning_rate': Real(0.01, 0.5, 'log-uniform'),
    'max_depth': Integer(0, 3)
}

# 参考
# https://www.kaggle.com/shaz13/santander-xgboost-bayesian-optimization/code
# https://www.kaggle.com/nanomathias/bayesian-optimization-of-xgboost-lb-0-9769
#

'''
GridSearchCV 系と比べた固有の API（一部）
- search_spaces : パラメータ探索範囲設定
- optimizer_kwargs : 最適化パラメータ。探索戦略やガウス過程パラメータなど。- 
- fit_params : clf.fit(X, y) に設定するパラメータ。学習実行時のパラメータ指定.
'''
clf = BayesSearchCV(
    estimator=estimator,
    search_spaces=param_grid,
    scoring='roc_auc',
    cv=3,
    n_jobs=-1,
    n_iter=10,
    verbose=1,
    refit=True
)
def status_print(optim_result):
    '''
    ベイズ最適化のイテレーションで呼ばれるコールバック
    '''
    results    = pd.DataFrame(clf.cv_results_)
    new_params = results.tail(1).params.values[0]
    new_score  = results.tail(1).mean_test_score
    print('Model #%d/%d : ROC-AUC(newest/best) = %.4f / %.4f' % (
        len(results),
        clf.n_iter,
        new_score,
        clf.best_score_
    ))

# BayesSearchCV は sklearn の他の交差検証サーチ系クラス(RandomSearchCVなど)と同じAPIを持っています
clf.fit(X, y, callback=status_print)

print(pd.DataFrame(clf.cv_results_))

# 検証データのパフォーマンス
print(clf.score(X_test, y_test))

探索の結果は clf.cv_resuts_ フィールドに格納されていて、pd.DataFrame にすると以下のようなデータとなっています。 この中には以下のデータが入っています（一部分列挙）。

• param_* は探索パラメータ
• mean_test_score は交差検証スコアの平均
• split{n}_test_score は交差検証の各分割でのスコア（これらを平均したものが mean_test_score）
• std_test_score は交差検証スコアの偏差

	mean_fit_time	mean_score_time	mean_test_score	split0_test_score	split1_test_score	split2_test_score	std_fit_time	std_score_time	std_test_score	param_colsample_bylevel	param_colsample_bytree
0	0.005673011	0.00219202	0.832636822	0.827494378	8.28E-01	0.842857143	0.000176986	0.00030934	0.007207978	0.432490483	0.101987568
1	0.005337397	0.002478361	0.5	0.5	0.5	0.5	0.000566812	0.000611297	0	0.45403928	0.928587396
2	0.005687555	0.002456427	0.5	5.00E-01	0.5	0.5	0.000361052	6.48E-05	0	0.471067626	0.880029631
3	0.005241315	0.00255092	0.5	5.00E-01	0.5	0.5	0.000284118	0.000263844	0	0.24114647	0.786112636
4	0.006174246	0.002198935	0.834788739	7.56E-01	0.832465381	0.915966387	0.000542449	0.000200357	0.065150251	0.699309626	0.752064985
5	0.007870754	0.002580643	0.5	0.5	0.5	0.5	0.000749173	0.000241813	0	0.932821363	0.339574308
6	0.004701535	0.00221777	0.5	0.5	0.5	0.5	0.000507002	0.000167713	0	0.818819988	0.79140138
7	0.004785299	0.002165318	0.5	0.5	5.00E-01	0.5	0.000287597	0.000235932	0	0.999741524	0.381359767
8	0.005346219	0.002468665	0.5	5.00E-01	0.5	0.5	0.000562654	0.000213206	0	0.55477075	0.31067584
9	0.005452951	0.00238506	0.5	5.00E-01	0.5	0.5	0.000326225	0.000221597	0	0.586112963	0.589381284

ベイズ最適化探索で得られた最適な結果は次のフィールドに格納されます。

• best_estimator_ : 最適なモデル
• best_params_ : 最適なハイパーパラメータ
• best_index : cv_results の最適なパラメータのインデックス
• best_score : 最適なモデルのスコア

まとめ

scikit-optimize の BayesSearchCV を用いて、ベイズ最適化によるハイパーパラメータ探索を試してみましたが、scikit-learn の RandomSearchCV や GridSearchCV と同じ使い方で、簡単にベイズ最適化を使えることがわかりました。

探索戦略のパラメータなどについては今後調べてみようと思います。

*1:東京大学の佐藤先生の講義が詳しいです https://www.youtube.com/watch?v=pQHWew4YYao

*2:≒交差検証精度の平均値 http://watanabe-www.math.dis.titech.ac.jp/users/swatanab/cross-val.html

jupyter notebook をコマンドラインで実行し HTML ファイルを生成する

はじめに

備忘録です。

最近は非エンジニアの分析官に一次分析内容を共有して定性分析を依頼したりする事が多いのですが、多数のモデルの場合 jupyter をある程度自動化したい場合があります。

私のよくやる流れとしては以下ですが、これを多数のモデルで PDCA を何回も回すことも珍しくありません。

1. jupyter で分析する.
2. html を出力して分析官に依頼する.
3. 分析官は自分の手慣れた方法（Excelなど）で分析、レポーティングする

今回はこういったことに役立つ、jupyter notebook をコマンドラインで実行し HTML ファイルを生成する方法を紹介します。

jupyter notebook をコマンドラインで実行する

コマンドラインで実行する

コマンドラインで実行することに関してはこちらが詳しいです。 qiita.com

コマンドラインで .pynb の HTML 化する

今回の HTML 化についてはこちらに殆ど書いています。

python - How export a Jupyter notebook to HTML from the command line? - Stack Overflow

# 以下を実行すると notebook.html ファイルが生成されます
$ jupyter nbconvert --execute --to html notebook.ipynb

私がよくやるのは、雛形となる jupyter notebook を作り、変更される値や設定を変数化し、コマンドラインで置換(perl や sed など)しながら様々な分析レポートを出力することです。

簡単な例

for model in model_A model_B .....;
do
    # 置換、或いは何らかのパーサー
    cat '雛形となる .ipynb' | perl -pe "s/置き換える文字列/${model}/s"  > ${model}.ipynb

    # html 生成
    jupyter nbconvert --execute --to html ${model}.ipynb
done

関数データ解析(FDA, Functional Data Analysis)とは

関数データ解析とは

関数データ解析とは

各個体や対象に対して、複数の離散点で時間や空間の変化に伴い観測・測定されたデータを、関数の集合として捉え、解析する方法
Ramsay and Silverman, 2005

とされています。

一般的にデータ解析では、データそのものの統計量をとったり、クラスタリングしたり、回帰分析したりします。 関数データ解析では、データを関数で表現して、その関数を改めてデータとして扱う手法です。

とはいえ、「データを関数で表現してそこで何か（分類や予測など）をする」という手法は、従来の機械学習にも使われていますし（例えばガウス線形回帰）、関数データ解析固有のものではありません。

関数データ解析では、特定の手法ではなくデータの関数表現をデータそのものの分析と同じように扱う方法論、ということになろうかと思います。

例えば、テキスト*1 では次のように書かれています。

Ramsay & Silverman (2005)は、関数データにおける平均、分散、共分散を定義した後、
主成分分析、線形モデル、正準相関分析、判別分析の関数データ対応版を扱っている。
さらに、関数の定義域を調整するRegistration (見当合わせ)、通常のデータを
関数データにする各種平滑化などについても詳細に検討している。

近年では、これらに加えクラスタリング（Kmeans の関数データ適用がよく使われるようです）、スパース解析（https://sites.google.com/site/hidetoshimatsui/field/fda）なども検討が進んでいます。

関数データ解析のながれ

主に次の二段階法が使われるようです。

1. 離散観測を連続の関数データで捉える（平滑化と呼ばれています）

   機械学習の手法で言うところのカーネル線形回帰などで、モデル関数を実データにフィッティングさせる。

2. 捉えた関数集合を分析する

   関数集合をデータ集合と同じように平均・分散・相関などを計算して分析する

所感

関数データ解析ではデータを関数で表現することで、データを集合で捉えより表現力のある構造（関数）で分析ができるようになります。 例えば従来のデータ解析でのデータの代表点（平均や中央値）や広がり（分散など）、関連（相関）などと同様の統計手法が関数データ解析においても使え、そこでは関数平均や関数分散、関数相関などの手法で分析ができます。

とは言え、私の中でもまだイメージの域を出ませんので、引き続き注目していきたいと思います。

*1:https://www.msi.co.jp/userconf/2007/pdf/07mizuta.pdf

Decision Tree - ID3(Iterative Dichotomiser 3 - 1979 John Ross Quinlan)

Decision Tree - ID3(Iterative Dichotomiser 3 - 1979 John Ross Quinlan)

(注意) これは私の勉強・備忘録のために記したものであり、間違いがあるやもしれません。どうぞご容赦ください。

はじめに

決定木のアルゴリズムに ID3 というものがあります。 ID3 は Quinlan が 1979 に発表したアルゴリズムで、この改良の C4.5 や、決定木を並列に多数構築した RandomForest、Boosting に活用した Gradient Boosted Trees など、発展型の元となったものです。

私の勉強のためにも、まずは ID3 から初めて、今 Kaggle でも最も定評・人気がある XGboost までたどってみようと思います。

ID3

Wikipediaja.wikipedia.org によると

その学習方法はオッカムの剃刀の原理に基づいている。すなわち最低限の仮説による事象の決定を行う。 「この方法は各独立変数に対し変数の値を決定した場合における平均情報量の期待値を求め、その中で最大のものを選びそれを木のノードにする操作を再帰的に行うことで実装される。」

とあります。これは ID3 に限ったことではなく決定木の基本的な考え方・原理です。

決定木は「最低限の仮説による事象の決定」とあるように、うまくコントロール*1すれば最小限の条件でその事象を表現できます。

特に、ヘテロジニアスデータ*2での適切に最適化された条件は人の感覚に合うため、人と機械のよいコンビネーションで PDCA を回せるなど、ビジネスの現場でもよく使われています。

実装

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import sys
import os
import math
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns

"""
    ID3(Iterative Dichotomiser 3 - 1979 John Ross Quinlan)
    その学習方法はオッカムの剃刀の原理に基づいている。すなわち最低限の仮説による事象の決定を行う。

    「この方法は各独立変数に対し変数の値を決定した場合における平均情報量の期待値を求め、
    その中で最大のものを選びそれを木のノードにする操作を再帰的に行うことで実装される。」
    
    https://ja.wikipedia.org/wiki/ID3
    制約: 離散でないといけない

    S: 訓練例の集合
    A: 属性の集合
    default: Yes / No の既定値
"""


# log の底は平均情報量の最大値が 1.0, 最小値が 0.0 になるように、分類数に合わせたほうが良い
# ここでは 3 とする.
def entropy(y, base=2):
    ''' 平均情報量の算出
    y=pd.Series
    target
    bird       0.4
    manmal     0.4
    reptile    0.2
    '''
    p = y.groupby(y).count() / y.shape[0]
    return -1 * sum(map(lambda x: x * math.log(x, base), p))

def mean_entropy(S, c, t='target'):
    ent = []
    for k, v in S.groupby(c).groups.items():
        # weighted average
        y = S.ix[v, t]
        ent.append(entropy(y, base=3) * y.shape[0] / S.shape[0])

    return np.sum(ent)


"""
 S : Tree
"""
class Node(object):

    def __init__(self, data, parent=None, col=None, feature=None, gain=0.):
        self.data   = data
        self.parent = parent
        self.col    = col
        self.feat   = feature

        # 子の分割基準
        self.gain   = gain
        self.childs = []

    def entropy(self, base=2):
        return entropy(self.data['target'], base=base)

    def shape(self):
        return self.data.shape

    def add_child(self, c):
        self.childs.append(c)

    def __repr__(self):
        return '%s(%s)\tH:%.3f Gain:%.3f' % (self.col, self.feat, self.entropy(), self.gain)


def ID3(S, A, t='target'):

    if S.data.shape[0] <= 0:
        # S が空集合
        return None
    
    elif len(A) <= 0:
        # Aが空集合
        return S

    elif S.entropy(base=3) <= 0.:
        return S

    else:
    
        # 1. 各説明変数で平均情報量を算出する <- このロジック mean_entropy の差し替えで C4.5 とかにできる
        Ms = map(lambda a: (a, mean_entropy(S.data, a)), A)

        # 2. 情報ゲインが最大な説明変数を選ぶ
        col, gain = max(map(lambda x: (x[0], S.entropy(3)-x[1]), Ms), key=lambda x:x[1])

        # 3. 木を分割する
        #_A = list(set(A) - set([col]))
        for f, x in S.data.groupby(col).groups.items():
            S.add_child(ID3(Node(S.data.ix[x], parent=S, col=col, feature=f, gain=gain), A))

        return S
    

def pretty_print_node(n, depth=0):
    if n is None:
        return 
    pref = '\t'
    print(pref*depth, n)
    for _n in n.childs:
        pretty_print_node(_n, depth+1)


if __name__ == '__main__':
    '''
    例題    食性（a1 ） 発生形態（a2 ）   体温（a3 ） 分類
    例題1（ペンギン） 肉食  卵生  恒温  鳥類
    例題2（ライオン） 肉食  胎生  恒温  哺乳類
    例題3（ウシ）　　 草食  胎生  恒温  哺乳類
    例題4（トカゲ）　 肉食  卵生  変温  爬虫類
    例題5（ブンチョウ）  草食  卵生  恒温  鳥類
    '''
    df = pd.DataFrame(
        [[ 'carnivorous', 'oviparity',  'homothermal',  'bird'],
         [ 'carnivorous', 'viviparity', 'homothermal',  'manmal'],
         [ 'herbivorous', 'viviparity', 'homothermal',  'manmal'],
         [ 'carnivorous', 'oviparity',  'cold_blooded', 'reptile'],
         [ 'herbivorous', 'oviparity',  'homothermal',  'bird']],
        columns = ['a1', 'a2', 'a3', 'target'],
        index = ['penguin', 'lion', 'caw', 'lizard','finch']
    )
    
    S = df
    A = ['a1', 'a2', 'a3'] # 属性集合
    X = df[['a1', 'a2', 'a3']]
    y = df.target
    
    Tree = ID3(Node(df, col='top'), A)
    pretty_print_node(Tree)
            

実行結果は

:decision_tree $ python id3.py
 top(None)      H:1.522 Gain:0.000
         a2(oviparity)  H:0.918 Gain:0.613
                 a3(cold_blooded)       H:-0.000 Gain:0.579
                 a3(homothermal)        H:-0.000 Gain:0.579
         a2(viviparity) H:-0.000 Gain:0.613

*1:実課題のデータはとても多様性があるため、ナイーブに決定木を最適化すると逆に条件が複雑になりオーバフィットしてしまいます。そのため、木の深さや単調性制約などを与えて機械的な最適化にある程度の歯止めをかけたりします

*2:ヘテロジニアスデータ:いわゆる異なる種類の説明変数で構成された表タイプのデータ、対するは画像ピクセルなど同種類のデータで構成されたホモジニアス・データ