Python: CatBoost を使ってみる

今回は CatBoost という、機械学習の勾配ブースティング決定木 (Gradient Boosting Decision Tree) というアルゴリズムを扱うためのフレームワークを試してみる。 CatBoost は、同じ勾配ブースティング決定木を扱うフレームワークの LightGBM や XGBoost と並んでよく用いられている。

CatBoost は学習にかかる時間が LightGBM や XGBoost に劣るものの、特にカテゴリカル変数を含むデータセットの扱いに定評がある。 ただし、今回使うデータセットはカテゴリカル変数を含まない点について先に断っておく。

使った環境は次の通り。

$ sw_vers                                     
ProductName:    Mac OS X
ProductVersion: 10.14.3
BuildVersion:   18D109
$ python -V            
Python 3.7.2

インストール

まずは pip を使って CatBoost をインストールする。 また、データセットの読み込みなどに使うため一緒に scikit-learn も入れておこう。

$ pip install catboost scikit-learn

基本的な使い方 (二値分類問題)

まずは Breast Cancer データセットを使った二値分類問題を CatBoost で処理してみよう。 以下のサンプルコードでは、CatBoost を学習させた上で汎化性能をホールドアウト検証で確認している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from catboost import CatBoost
from catboost import Pool

from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score


def main():
    # 乳がんデータセットを読み込む
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target
    # データセットを学習用と検証用に分割する
    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)
    # CatBoost が扱うデータセットの形式に直す
    train_pool = Pool(X_train, label=y_train)
    test_pool = Pool(X_test, label=y_test)
    # 学習用のパラメータ
    params = {
        # タスク設定と損失関数
        'loss_function': 'Logloss',
        # 学習ラウンド数
        'num_boost_round': 100,
    }
    # モデルを学習する
    model = CatBoost(params)
    model.fit(train_pool)
    # 検証用データを分類する
    # NOTE: 確率がほしいときは prediction_type='Probability' を使う
    y_pred = model.predict(test_pool, prediction_type='Class')
    # 精度 (Accuracy) を検証する
    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)


if __name__ == '__main__':
    main()

上記を実行してみよう。 各ラウンドでの学習データに対する損失と時間が表示される。 それぞれのラウンドでは決定木がひとつずつ作られている。

$ python helloworld.py 
Learning rate set to 0.100436
0: learn: 0.5477448   total: 108ms    remaining: 10.6s
1: learn: 0.4414628   total: 140ms    remaining: 6.85s
2: learn: 0.3561365   total: 172ms    remaining: 5.55s
...(snip)...
97:    learn: 0.0120305   total: 3.26s   remaining: 66.6ms
98:    learn: 0.0116963   total: 3.29s   remaining: 33.3ms
99:    learn: 0.0113142   total: 3.33s   remaining: 0us
Accuracy: 0.9532163742690059

最終的に精度 (Accuracy) で約 0.953 という結果が得られた。

ラウンド数を増やしてみる

先ほどの例ではラウンド数、つまり用いる決定木の数が最大で 100 本だった。 次は、試しにこの数を 1,000 まで増やしてみよう。

以下のサンプルコードではラウンド数を増やした上で、最も最適なラウンド数を用いて最終的な性能を確認している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from catboost import CatBoost
from catboost import Pool

from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score


def main():
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)

    train_pool = Pool(X_train, label=y_train)
    test_pool = Pool(X_test, label=y_test)

    params = {
        'loss_function': 'Logloss',
        'num_boost_round': 1000,  # ラウンド数を多めにしておく
    }

    model = CatBoost(params)
    # 検証用データに対する損失を使って学習課程を評価する
    # 成果物となるモデルには、それが最も良かったものを使う
    model.fit(train_pool, eval_set=[test_pool], use_best_model=True)

    y_pred = model.predict(test_pool, prediction_type='Class')

    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)


if __name__ == '__main__':
    main()

上記を実行してみる。

$ python valset.py 
Learning rate set to 0.070834
0: learn: 0.5872722   test: 0.5996097  best: 0.5996097 (0)   total: 109ms    remaining: 1m 49s
1: learn: 0.5025692   test: 0.5240925  best: 0.5240925 (1)   total: 140ms    remaining: 1m 9s
2: learn: 0.4279025   test: 0.4535490  best: 0.4535490 (2)   total: 172ms    remaining: 57s
...(snip)...
997:   learn: 0.0007559   test: 0.0891592  best: 0.0872851 (771) total: 37s  remaining: 74.1ms
998:   learn: 0.0007554   test: 0.0891578  best: 0.0872851 (771) total: 37s  remaining: 37ms
999:   learn: 0.0007543   test: 0.0891842  best: 0.0872851 (771) total: 37s  remaining: 0us

bestTest = 0.08728505181
bestIteration = 771

Shrink model to first 772 iterations.
Accuracy: 0.9590643274853801

最終的に、学習データの損失において最も優れていたのは 771 ラウンドで、そのときの精度 (Accuracy) は約 0.959 だった。 先ほどよりも精度が 0.6 ポイント改善していることが分かる。

学習の過程を可視化する

続いては、モデルが学習する過程を可視化してみよう。 次のサンプルコードでは 1,000 ラウンド回した場合の、学習データと検証データに対する損失の変化を可視化している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from catboost import CatBoost
from catboost import Pool

from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from matplotlib import pyplot as plt


def main():
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)

    train_pool = Pool(X_train, label=y_train)
    test_pool = Pool(X_test, label=y_test)

    params = {
        'loss_function': 'Logloss',
        'num_boost_round': 1000,
    }

    model = CatBoost(params)
    model.fit(train_pool, eval_set=[test_pool], use_best_model=True)

    y_pred = model.predict(test_pool, prediction_type='Class')

    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)

    # メトリックの推移を取得する
    history = model.get_evals_result()
    # グラフにプロットする
    train_metric = history['learn']['Logloss']
    plt.plot(train_metric, label='train metric')
    eval_metric = history['validation_0']['Logloss']
    plt.plot(eval_metric, label='eval metric')

    plt.legend()
    plt.grid()
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみる。

$ python visualize.py 
Learning rate set to 0.070834
0: learn: 0.5872722   test: 0.5996097  best: 0.5996097 (0)   total: 117ms    remaining: 1m 56s
1: learn: 0.5025692   test: 0.5240925  best: 0.5240925 (1)   total: 150ms    remaining: 1m 14s
2: learn: 0.4279025   test: 0.4535490  best: 0.4535490 (2)   total: 182ms    remaining: 1m
...(snip)...
997:   learn: 0.0007559   test: 0.0891592  best: 0.0872851 (771) total: 36.6s   remaining: 73.3ms
998:   learn: 0.0007554   test: 0.0891578  best: 0.0872851 (771) total: 36.6s   remaining: 36.7ms
999:   learn: 0.0007543   test: 0.0891842  best: 0.0872851 (771) total: 36.7s   remaining: 0us

bestTest = 0.08728505181
bestIteration = 771

Shrink model to first 772 iterations.
Accuracy: 0.9590643274853801

最終的な結果は先ほどと変わらない。

学習が完了すると、次のようなグラフが得られる。

上記を見ると 200 ラウンドを待たずに検証用データに対する損失は底を打っていることが分かる。

学習が進まなくなったら打ち切る

ここまで見てきた通り、CatBoost ではとりあえず多めのラウンドを回して、最終的に良かったものを使うことができる。 とはいえ、一旦学習が進まなくなると、そこからまた改善するということはそんなにない。 大抵の場合は大して改善しない状況が続くか、過学習が進みやすい。 そこで、続いては学習が進まなくなったらそこで打ち切る Early Stopping という機能を使ってみる。 この機能は LightGBM や XGBoost でも実装されており、よく用いられるものの一つとなっている。

次のサンプルコードでは Early Stopping を用いて 1,000 ラウンド回す中で、検証用データの損失が 10 ラウンド改善しなかったらそこで学習を打ち切るようになっている。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from catboost import CatBoost
from catboost import Pool

from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score


def main():
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)

    train_pool = Pool(X_train, label=y_train)
    test_pool = Pool(X_test, label=y_test)

    params = {
        'loss_function': 'Logloss',
        'num_boost_round': 1000,
        # 検証用データの損失が既定ラウンド数減らなかったら学習を打ち切る
        'early_stopping_rounds': 10,
    }

    model = CatBoost(params)
    model.fit(train_pool, eval_set=[test_pool])

    y_pred = model.predict(test_pool, prediction_type='Class')

    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)


if __name__ == '__main__':
    main()

上記を実行してみる。

$ python earlystop.py 
Learning rate set to 0.070834
0: learn: 0.5872722   test: 0.5996097  best: 0.5996097 (0)   total: 115ms    remaining: 1m 55s
1: learn: 0.5025692   test: 0.5240925  best: 0.5240925 (1)   total: 150ms    remaining: 1m 15s
2: learn: 0.4279025   test: 0.4535490  best: 0.4535490 (2)   total: 183ms    remaining: 1m
...(snip)...
136:   learn: 0.0139911   test: 0.0939635  best: 0.0932902 (128) total: 5.52s   remaining: 34.8s
137:   learn: 0.0139071   test: 0.0936872  best: 0.0932902 (128) total: 5.55s   remaining: 34.7s
138:   learn: 0.0138493   test: 0.0938802  best: 0.0932902 (128) total: 5.59s   remaining: 34.6s
Stopped by overfitting detector  (10 iterations wait)

bestTest = 0.09329017842
bestIteration = 128

Shrink model to first 129 iterations.
Accuracy: 0.9590643274853801

今度は 138 ラウンドで学習が止まったものの、最終的な汎化性能は先ほどと変わらないものが得られている。 このように Early Stopping は計算量を節約する上で重要な機能になっている。

scikit-learn インターフェースを使ってみる

CatBoost には scikit-learn のインターフェースもある。 続いてはそれを使ってみることにしよう。

以下のサンプルコードでは scikit-learn のインターフェースを備えた CatBoostClassifier を使っている。 これを用いることで使い慣れた scikit-learn のオブジェクトとして CatBoost を扱うことができる。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from catboost import CatBoostClassifier

from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score


def main():
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)

    # scikit-learn インターフェースを備えたラッパー
    model = CatBoostClassifier(num_boost_round=1000,
                               loss_function='Logloss',
                               early_stopping_rounds=10,
                               )
    model.fit(X_train, y_train, eval_set=[(X_test, y_test)])
    # 確率がほしいときは predict_proba() を使う
    y_pred = model.predict(X_test)

    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)


if __name__ == '__main__':
    main()

上記を実行してみる。

$ python skapi.py 
Learning rate set to 0.070834
0: learn: 0.5872722   test: 0.5996097  best: 0.5996097 (0)   total: 112ms    remaining: 1m 51s
1: learn: 0.5025692   test: 0.5240925  best: 0.5240925 (1)   total: 143ms    remaining: 1m 11s
2: learn: 0.4279025   test: 0.4535490  best: 0.4535490 (2)   total: 176ms    remaining: 58.4s
...(snip)...
136:   learn: 0.0139911   test: 0.0939635  best: 0.0932902 (128) total: 5.8s    remaining: 36.5s
137:   learn: 0.0139071   test: 0.0936872  best: 0.0932902 (128) total: 5.83s   remaining: 36.4s
138:   learn: 0.0138493   test: 0.0938802  best: 0.0932902 (128) total: 5.88s   remaining: 36.4s
Stopped by overfitting detector  (10 iterations wait)

bestTest = 0.09329017842
bestIteration = 128

Shrink model to first 129 iterations.
Accuracy: 0.9590643274853801

結果は先ほどと変わらない。

多値分類問題を処理してみる

ここまでは二値分類問題を扱ってきた。 続いては Iris データセットを用いて多値分類問題を解かせてみよう。

次のサンプルコードでは Iris データセットを CatBoost で分類している。 検証方法は先ほどと同じホールドアウト検証で、精度 (Accuracy) について評価している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from catboost import CatBoost
from catboost import Pool

from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score


def main():
    # Iris データセットを読み込む
    dataset = datasets.load_iris()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)

    train_pool = Pool(X_train, label=y_train)
    test_pool = Pool(X_test, label=y_test)

    params = {
        # 多値分類問題
        'loss_function': 'MultiClass',
        'num_boost_round': 1000,
        'early_stopping_rounds': 10,
    }

    model = CatBoost(params)
    model.fit(train_pool, eval_set=[test_pool])

    y_pred = model.predict(test_pool, prediction_type='Class')

    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)


if __name__ == '__main__':
    main()

上記の実行結果は次の通り。

$ python multiclass.py 
0: learn: -1.0663112  test: -1.0680680 best: -1.0680680 (0)  total: 65.3ms  remaining: 1m 5s
1: learn: -1.0369337  test: -1.0404371 best: -1.0404371 (1)  total: 74.7ms  remaining: 37.3s
2: learn: -1.0094257  test: -1.0172382 best: -1.0172382 (2)  total: 78.9ms  remaining: 26.2s
...(snip)...
268:   learn: -0.0561336  test: -0.1770632 best: -0.1768477 (260)    total: 986ms    remaining: 2.68s
269:   learn: -0.0557735  test: -0.1773530 best: -0.1768477 (260)    total: 991ms    remaining: 2.68s
270:   learn: -0.0556859  test: -0.1772400 best: -0.1768477 (260)    total: 994ms    remaining: 2.67s
Stopped by overfitting detector  (10 iterations wait)

bestTest = -0.1768476949
bestIteration = 260

Shrink model to first 261 iterations.
Accuracy: 0.9111111111111111

特徴量の重要度を可視化する

CatBoost が採用している勾配ブースティング決定木は名前の通り決定木の仲間なので特徴量の重要度を取得できる。 続いては特徴量の重要度を可視化してみよう。

以下のサンプルコードでは、先ほどと同じ条件で学習したモデルから特徴量の重要度を取得してグラフにプロットしている。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from catboost import CatBoost
from catboost import Pool

from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from matplotlib import pyplot as plt


def main():
    # Iris データセットを読み込む
    dataset = datasets.load_iris()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)

    train_pool = Pool(X_train, label=y_train)
    test_pool = Pool(X_test, label=y_test)

    params = {
        'loss_function': 'MultiClass',
        'num_boost_round': 1000,
        'early_stopping_rounds': 10,
    }

    model = CatBoost(params)
    model.fit(train_pool, eval_set=[test_pool])

    y_pred = model.predict(test_pool, prediction_type='Class')

    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)

    # 特徴量の重要度を取得する
    feature_importance = model.get_feature_importance()
    # 棒グラフとしてプロットする
    plt.figure(figsize=(12, 4))
    plt.barh(range(len(feature_importance)),
            feature_importance,
            tick_label=dataset.feature_names)

    plt.xlabel('importance')
    plt.ylabel('features')
    plt.grid()
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみる。

$ python featimp.py   
0: learn: -1.0663112  test: -1.0680680 best: -1.0680680 (0)  total: 61.7ms  remaining: 1m 1s
1: learn: -1.0369337  test: -1.0404371 best: -1.0404371 (1)  total: 67.2ms  remaining: 33.6s
2: learn: -1.0094257  test: -1.0172382 best: -1.0172382 (2)  total: 71ms remaining: 23.6s
...(snip)...
268:   learn: -0.0561336  test: -0.1770632 best: -0.1768477 (260)    total: 961ms    remaining: 2.61s
269:   learn: -0.0557735  test: -0.1773530 best: -0.1768477 (260)    total: 964ms    remaining: 2.61s
270:   learn: -0.0556859  test: -0.1772400 best: -0.1768477 (260)    total: 967ms    remaining: 2.6s
Stopped by overfitting detector  (10 iterations wait)

bestTest = -0.1768476949
bestIteration = 260

Shrink model to first 261 iterations.
Accuracy: 0.9111111111111111

得られたグラフは次の通り。

上記から Petal length と Petal width が分類において有効に作用していたことが確認できる。

回帰問題を処理してみる

続いては回帰問題を扱ってみる。 以下のサンプルコードでは Boston データセットを CatBoost で回帰している。 汎化性能の確認してはホールドアウト検証で RMSE (Root Mean Squared Error) を評価している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import math

from catboost import CatBoost
from catboost import Pool

from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error


def main():
    # Boston データセットを読み込む
    dataset = datasets.load_boston()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42)

    train_pool = Pool(X_train, label=y_train)
    test_pool = Pool(X_test, label=y_test)

    params = {
        # 損失関数に RMSE を使う
        'loss_function': 'RMSE',
        'num_boost_round': 1000,
        'early_stopping_rounds': 10,
    }

    model = CatBoost(params)
    model.fit(train_pool, eval_set=[test_pool])

    y_pred = model.predict(test_pool)

    # 最終的なモデルの RMSE を計算する
    mse = mean_squared_error(y_test, y_pred)
    print('RMSE:', math.sqrt(mse))


if __name__ == '__main__':
    main()

上記を実行してみる。

$ python regress.py 
0: learn: 24.2681654  test: 22.5038434 best: 22.5038434 (0)  total: 75.5ms  remaining: 1m 15s
1: learn: 23.6876829  test: 21.9483221 best: 21.9483221 (1)  total: 84.6ms  remaining: 42.2s
2: learn: 23.1173979  test: 21.3856082 best: 21.3856082 (2)  total: 92.1ms  remaining: 30.6s
...(snip)...
466:   learn: 2.4850147   test: 3.3697155  best: 3.3675366 (458) total: 4.11s   remaining: 4.69s
467:   learn: 2.4833511   test: 3.3700655  best: 3.3675366 (458) total: 4.13s   remaining: 4.69s
468:   learn: 2.4828831   test: 3.3701698  best: 3.3675366 (458) total: 4.14s   remaining: 4.69s
Stopped by overfitting detector  (10 iterations wait)

bestTest = 3.36753664
bestIteration = 458

Shrink model to first 459 iterations.
RMSE: 3.367536638941265

いじょう。

スマートPythonプログラミング: Pythonのより良い書き方を学ぶ

Python: k-NN Feature Extraction 用のライブラリ「gokinjo」を作った

表題の通り、k-NN Feature Extraction という特徴量抽出の手法に使う「gokinjo」という Python のライブラリを作った。 今回はライブラリの使い方について紹介してみる。

github.com

k-NN Feature Extraction で得られる特徴量は、Otto Group Product Classification Challenge という Kaggle のコンペで優勝者が使ったものの一つ。 手法自体の解説は以下のブログ記事を参照のこと。 なお、ブログ記事の中でも Python のコードを書いてるけど、よりナイーブな実装になっている。

blog.amedama.jp

以降、紹介に用いる環境は次の通り。 なお、ライブラリ自体にプラットフォームへの依存はない。 また、Python のバージョンは 3.6 以降をサポートしている。

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.14.2
BuildVersion:   18C54
$ python -V
Python 3.7.2

インストール

ライブラリは pip を使ってインストールできる。

$ pip install gokinjo

基本的な使い方

まずは簡単な例を使って使い方を解説する。

データを可視化するために matplotlib をインストールしておく。

$ pip install matplotlib 

Python のインタプリタを起動する。

$ python

まずは次のような二次元の特徴変数と二値の目的変数を持ったダミーデータを生成する。

>>> import numpy as np
>>> x0 = np.random.rand(500) - 0.5
>>> x1 = np.random.rand(500) - 0.5
>>> X = np.array(list(zip(x0, x1)))
>>> y = np.array([1 if i0 * i1 > 0 else 0 for i0, i1 in X])

このデータを可視化すると次のようになる。 データはラベルごとに XOR 的に配置されており、線形分離できないようになっている。

>>> from matplotlib import pyplot as plt
>>> plt.scatter(X[:, 0], X[:, 1], c=y)
<matplotlib.collections.PathCollection object at 0x11d47d048>
>>> plt.show()

上記のデータから、今回作ったライブラリを使って特徴量を抽出する。 抽出した特徴量は、近傍数 k とラベルの種類 C をかけた次元数を持っている。 今回であればデフォルトの近傍数 k=1 とラベル数 C=2 から k * C = 2 となって二次元のデータになる。

>>> from gokinjo import knn_kfold_extract
>>> X_knn = knn_kfold_extract(X, y)

抽出した上記のデータを可視化してみる。

plt.scatter(X_knn[:, 0], X_knn[:, 1], c=y)
plt.show()

元々のデータに比べると、もうちょっとで線形分離できそうなくらいの特徴量になっていることが分かる。

Iris データセットを使った例

先ほどと同様に Iris データセットでも特徴量を抽出して可視化してみることにしよう。 サンプルコードは次の通り。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from sklearn import datasets
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa

from gokinjo import knn_kfold_extract


def main():
    # データセットを読み込む
    dataset = datasets.load_iris()
    X = dataset.data
    y = dataset.target

    # k-NN 特徴量を取り出す
    extracted_features = knn_kfold_extract(X, y)

    # 取り出した特徴量を可視化する
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    ax.scatter(extracted_features[:, 0],
               extracted_features[:, 1],
               extracted_features[:, 2],
               c=y)
    ax.set_title('extracted features (3D)')

    axes = plt.gcf().get_axes()
    for axe in axes:
        axe.grid()

    plt.show()


if __name__ == '__main__':
    main()

上記に適当な名前をつけて実行する。

$ python knn_iris.py

抽出した特徴量を可視化したグラフは次の通り。 Iris データセットは多値分類問題なので k=1, C=3 より 3 次元のデータになっている。

こちらも良い感じに分離できそうな特徴量が抽出できている。

Digits データセットを使った例

可視化の次は実際に抽出した特徴量を分類問題に適用してみる。 以下のサンプルコードでは Digits データセットを使って特徴量をランダムフォレストで分類している。 Digits データセットというのは 0 ~ 9 までの手書き数字の画像が含まれたもので MNIST のミニチュアみたいな感じ。 パターンとしては「生データそのまま」「k-NN 特徴量」「t-SNE 特徴量」「t-SNE -> k-NN 特徴量」で試している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from sklearn import datasets
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_validate
from sklearn.manifold import TSNE

from gokinjo import knn_kfold_extract


def main():
    # データセットを読み込む
    dataset = datasets.load_digits()
    X = dataset.data
    y = dataset.target

    # 下準備
    clf = RandomForestClassifier(n_estimators=100, random_state=42)
    skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

    # 生データをそのまま使う場合
    score = cross_validate(clf, X, y, cv=skf)
    print('mean accuracy (raw):', score['test_score'].mean())

    # k-NN 特徴量に変換した場合
    X_knn_raw = knn_kfold_extract(X, y, random_state=42)
    score = cross_validate(clf, X_knn_raw, y, cv=skf)
    print('mean accuracy (raw -> k-NN):', score['test_score'].mean())

    # t-SNE に変換した場合
    tsne = TSNE(random_state=42)
    X_tsne = tsne.fit_transform(X)
    score = cross_validate(clf, X_tsne, y, cv=skf)
    print('mean accuracy (raw -> t-SNE):', score['test_score'].mean())

    # t-SNE をさらに k-NN 特徴量に変換した場合
    X_knn_tsne = knn_kfold_extract(X_tsne, y, random_state=42)
    score = cross_validate(clf, X_knn_tsne, y, cv=skf)
    print('mean accuracy (raw -> t-SNE -> k-NN):', score['test_score'].mean())


if __name__ == '__main__':
    main()

上記を実行した結果が以下の通り。 わずかな改善ではあるものの k-NN 特徴量が効いていることが分かる。

$ python knn_digits.py 
mean accuracy (raw): 0.975493504158074
mean accuracy (raw -> k-NN): 0.9777156040302326
mean accuracy (raw -> t-SNE): 0.9888661465635314
mean accuracy (raw -> t-SNE -> k-NN): 0.9894313077402828

バックエンドの切り替え

ライブラリの特徴として、k-NN アルゴリズムの実装が切り替えられる。 デフォルトではバックエンドが scikit-learn になっているけど、オプションとして Annoy にもできる。

Annoy をバックエンドとして使うときはインストールするときに環境として [annoy] を指定する。

$ pip install "gokinjo[annoy]"

そして、knn_kfold_extract() 関数の backend オプションに annoy を指定すれば良い。

knn_kfold_extract(..., backend='annoy')

バックエンドを切り替えることで、使うデータや環境によっては高速化につながるかもしれない。

分割しながらの生成が不要な場合

ちなみに k-NN Feature Extraction は目的変数にもとづいた特徴量抽出なので、生成するときに過学習を防ぐ工夫が必要になる。 具体的には、特徴量を生成する対象のデータを学習の対象に含めてはいけないという点。 ようするに Target Encoding するときと同じように k-Fold で分割しながら特徴量を生成する必要がある。 そのため、ここまでのサンプルコードでは、データを内部的に自動で分割しながら生成してくれる knn_kfold_extract() という関数を使ってきた。

もしデータ分割がいらない、もしくは自分でやるという場合には knn_extract() という関数が使える。 典型的なユースケースは Kaggle の test データに対して特徴量を生成する場合。 このときは train データを全て使って test データの特徴量を生成すれば良いため。

以下は knn_extract() を使って自前で分割しながら特徴量を生成する場合のサンプルコード。 生成した特徴量のソートとかでコードが分かりにくくなるので、ここでは Stratified でない k-Fold をシャッフルせずに使っている。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np
from sklearn import datasets
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import KFold
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_validate

from gokinjo import knn_extract


def main():
    # データセットを読み込む
    dataset = datasets.load_digits()
    X, y = dataset.data, dataset.target

    # 生成した特徴量を保存する場所
    number_of_classes = np.unique(y)
    k = 1
    X_knn = np.empty([0, len(number_of_classes) * k])

    # Leakage を防ぐために k-Fold で分割しながら生成する
    kf = KFold(n_splits=20, shuffle=False)  # 処理の単純化のためにシャッフルなし
    for train_index, test_index in kf.split(X, y):
        X_train, X_test = X[train_index], X[test_index]
        y_train, y_test = y[train_index], y[test_index]  # noqa

        # 学習用データで学習して、テストデータに対して特徴量を生成する
        X_test_knn = knn_extract(X_train, y_train, X_test, k=k)

        # 生成した特徴量を保存する
        X_knn = np.append(X_knn, X_test_knn, axis=0)

    # 生成した特徴量をランダムフォレストで学習してみる
    clf = RandomForestClassifier(n_estimators=100, random_state=42)
    skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    score = cross_validate(clf, X, y, cv=skf)
    print('mean accuracy (raw data)', score['test_score'].mean())
    score = cross_validate(clf, X_knn, y, cv=skf)
    print('mean accuracy (k-NN feature)', score['test_score'].mean())


if __name__ == '__main__':
    main()

上記の実行結果は次の通り。

$ python knn_nokfold.py 
mean accuracy (raw data) 0.975493504158074
mean accuracy (k-NN feature) 0.9782572815114076

scikit-learn インターフェース

ちなみに scikit-learn の Transformer を実装したインターフェースも用意してある。 こちらも、使うときは特徴量を生成するときに自前で k-Fold する必要がある点に注意が必要。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np
from sklearn import datasets
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import KFold
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_validate

from gokinjo.backend_sklearn import ScikitTransformer
# from gokinjo.backend_annoy import AnnoyTransformer  # Annoy 使うなら


def main():
    # データセットを読み込む
    dataset = datasets.load_digits()
    X, y = dataset.data, dataset.target

    # k-NN 特徴量を抽出する scikit-learn の Transformer として実装したやつ
    k = 1
    transformer = ScikitTransformer(n_neighbors=k)

    # 生成した特徴量を保存する場所
    number_of_classes = np.unique(y)
    X_knn = np.empty([0, len(number_of_classes) * k])

    # Leakage を防ぐために k-Fold で分割しながら生成する
    kf = KFold(n_splits=20, shuffle=False)  # 処理の単純化のためにシャッフルなし
    for train_index, test_index in kf.split(X, y):
        X_train, X_test = X[train_index], X[test_index]
        y_train, y_test = y[train_index], y[test_index]  # noqa

        # 学習用データで学習する
        transformer.fit(X_train, y_train)

        # テストデータに対して特徴量を生成する
        X_test_knn = transformer.transform(X_test)

        # 生成した特徴量を保存する
        X_knn = np.append(X_knn, X_test_knn, axis=0)

    # 生成した特徴量をランダムフォレストで学習してみる
    clf = RandomForestClassifier(n_estimators=100, random_state=42)
    skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    score = cross_validate(clf, X, y, cv=skf)
    print('mean accuracy (raw data)', score['test_score'].mean())
    score = cross_validate(clf, X_knn, y, cv=skf)
    print('mean accuracy (k-NN feature)', score['test_score'].mean())


if __name__ == '__main__':
    main()

実行結果は次の通り。

$ python knn_sklearn.py 
mean accuracy (raw data) 0.975493504158074
mean accuracy (k-NN feature) 0.9782572815114076

いじょう。

スマートPythonプログラミング: Pythonのより良い書き方を学ぶ

Linux Bridge でタグ VLAN を使う

今回は Linux Bridge (ネットワークブリッジ) でタグ VLAN (IEEE 802.1Q) を使う方法について。 設定が足らなくてだいぶ悩んだのでメモしておく。

使った環境は次の通り。

$ cat /etc/lsb-release 
DISTRIB_ID=Ubuntu
DISTRIB_RELEASE=18.04
DISTRIB_CODENAME=bionic
DISTRIB_DESCRIPTION="Ubuntu 18.04.1 LTS"
$ uname -r
4.15.0-29-generic
$ ip -V
ip utility, iproute2-ss180129
$ sudo apt-cache show iproute2 | grep -i version
Version: 4.15.0-2ubuntu1

下準備

最初に、今回使うパッケージをインストールしておく。

$ sudo apt-get -y install iproute2 tcpdump iputils-ping

Linux Bridge でタグ VLAN を有効にする

準備ができたところで、早速本題に入る。 Linux Bridge でタグ VLAN を扱うには vlan_filtering を有効にする必要がある。 vlan_filtering は ip(8) で有効にできる。

例えば ip(8) を使ってブリッジを作る場合には、次のように有効にできる。 これは br0 という名前でブリッジを作るときの例。

$ sudo ip link add dev br0 type bridge vlan_filtering 1

また、作ったあとからでも次のように有効にできる。 brctl(8) とか使って作った場合にも、これでいけるはず。

$ sudo ip link add dev br0 type bridge
$ sudo ip link set dev br0 type bridge vlan_filtering 1

ブリッジを作ったら状態を UP に設定する。

$ sudo ip link set br0 up

ここからは、上記の設定が正しく動いていることを確認していく。

アクセスポートを追加する

続いては、タグ VLAN を有効にした Linux Bridge に実際にアクセスポートを追加する。 追加するネットワークインターフェースには veth インターフェースを使う。 また、あとから ping(8) で疎通を確認できるように Network Namespace をつなげておく。

まずは Network Namespace を用意する。

$ sudo ip netns add ns1

続いて veth インターフェースのペアを作成する。

$ sudo ip link add ns1-veth0 type veth peer name br0-veth0

片方の veth インターフェースを、先ほど作った Network Namespace に所属させる。

$ sudo ip link set ns1-veth0 netns ns1

所属させたネットワークインターフェースを使える状態にして IP アドレスとして 192.0.2.1/24 を付与しておく。

$ sudo ip netns exec ns1 ip link set ns1-veth0 up
$ sudo ip netns exec ns1 ip addr add dev ns1-veth0 192.0.2.1/24

veth インターフェースのもう片方は、先ほど作った Linux Bridge に追加しておく。

$ sudo ip link set br0-veth0 master br0
$ sudo ip link set br0-veth0 up

準備ができたので、やっと VLAN の設定をしていく。 その前に、まずは bridge(8) で現状を確認しておく。 デフォルトでは、次のように VLAN ID として 1 が使われるように見えている。 これは、ようするに 802.1Q VLAN を使っていない状態。

$ bridge vlan show dev br0-veth0
port    vlan ids
br0-veth0    1 PVID Egress Untagged

そこで、bridge(8) を使ってポートに VLAN を設定する。 今回は Linux Bridge に追加したポートを使った通信で VLAN ID 100 が使われるように設定してみよう。 代わりに、既存のルールは削除しておく。

$ sudo bridge vlan add vid 100 dev br0-veth0 pvid untagged
$ sudo bridge vlan del dev br0-veth0 vid 1

設定は以下のようになった。

$ bridge vlan show dev br0-veth0
port    vlan ids
br0-veth0    100 PVID Egress Untagged

通信先 (対向) を用意する

ns1 に続いて、もう一つ ns2 という名前で Network Namespace を用意して、先ほどと同じようにセットアップしておく。 これは ns1 から ping(8) を打つときの相手になる。

$ sudo ip netns add ns2
$ sudo ip link add ns2-veth0 type veth peer name br0-veth1
$ sudo ip link set ns2-veth0 netns ns2
$ sudo ip netns exec ns2 ip link set ns2-veth0 up
$ sudo ip netns exec ns2 ip addr add dev ns2-veth0 192.0.2.2/24
$ sudo ip link set br0-veth1 master br0
$ sudo ip link set br0-veth1 up
$ sudo bridge vlan add vid 100 dev br0-veth1 pvid untagged
$ sudo bridge vlan del dev br0-veth1 vid 1

通信を観察する

準備ができたので、Linux Bridge の通信を tcpdump(1) でパケットキャプチャする。

$ sudo tcpdump -nel -i br0
tcpdump: verbose output suppressed, use -v or -vv for full protocol decode
listening on br0, link-type EN10MB (Ethernet), capture size 262144 bytes

パケットキャプチャの準備ができたら、別のターミナルから ns1 から ns2 に向かって ping を打ってみよう。

$ sudo ip netns exec ns1 ping -c 3 192.0.2.2
PING 192.0.2.2 (192.0.2.2) 56(84) bytes of data.
64 bytes from 192.0.2.2: icmp_seq=1 ttl=64 time=0.114 ms
64 bytes from 192.0.2.2: icmp_seq=2 ttl=64 time=0.060 ms
64 bytes from 192.0.2.2: icmp_seq=3 ttl=64 time=0.073 ms

--- 192.0.2.2 ping statistics ---
3 packets transmitted, 3 received, 0% packet loss, time 2036ms
rtt min/avg/max/mdev = 0.060/0.082/0.114/0.024 ms

打ち終わったら tcpdump(1) を実行していたターミナルに戻る。 すると、次のように 802.1Q のデータフレームが観測できているはず。 VLAN ID にも、ちゃんと 100 が使われていることが分かる。

$ sudo tcpdump -nel -i br0
...(略)...
18:20:18.031208 b6:05:b6:b6:c6:e3 > a2:bd:24:29:d8:76, ethertype 802.1Q (0x8100), length 102: vlan 100, p 0, ethertype IPv4, 192.0.2.1 > 192.0.2.2: ICMP echo request, id 1453, seq 1, length 64
18:20:18.031240 a2:bd:24:29:d8:76 > b6:05:b6:b6:c6:e3, ethertype 802.1Q (0x8100), length 102: vlan 100, p 0, ethertype IPv4, 192.0.2.2 > 192.0.2.1: ICMP echo reply, id 1453, seq 1, length 64
...(略)

めでたしめでたし。

スマートPythonプログラミング: Pythonのより良い書き方を学ぶ

Python: 自作ライブラリのパッケージングについて

今回は Python で自作したライブラリなどをパッケージングして、配布できる状態にする方法について書いてみる。

現在の Python では、パッケージングに setuptools というサードパーティ製のライブラリを使うのがデファクトスタンダードになっている。 この setuptools は、pip などを使ってパッケージをインストールするときにも必要になるので、実は気づかずに使っているという場合も多いかもしれない。 また、サードパーティ製といっても PyPA (Python Packaging Authority) というコミュニティが管理しているので準公式みたいな位置づけ。

今回使った環境は次の通り。

$ sw_vers                                      
ProductName:    Mac OS X
ProductVersion: 10.14.2
BuildVersion:   18C54
$ python -V              
Python 3.7.2
$ pip list | grep setuptools
setuptools 40.7.2

パッケージング用のサンプルコードを用意する

まずはパッケージングの対象となるサンプルコードを用意する必要がある。 そこで、今回は greet() という関数を一つだけ備えたパッケージとして example を用意した。 関数 greet() は、呼び出すと標準出力にメッセージをプリントする。

$ mkdir example
$ cat << 'EOF' > example/__init__.py
# -*- coding: utf-8 -*-

__version__ = '0.0.1'


def greet():
    """挨拶をする関数"""
    print('Hello, World!')
EOF

ちょっと紛らわしいんだけど Python においてパッケージというのは、パッケージングされたライブラリのことを意味していない。 ここでは、Python のコードが入ったディレクトリ、くらいの感覚で考えてもらえれば良い。

ちゃんとした Python におけるモジュールとパッケージという概念については以下を参照のこと。 簡単にいうとモジュールは Python ファイル (*.py) で、パッケージは __init__.py という名前のファイルが入ったディレクトリのことを指している。 blog.amedama.jp

上記で作ったパッケージは、今いるディレクトリからであれば Python の処理系からインポートして使うことができる。 これは、カレントワーキングディレクトリに Python のパスが通っているため。

$ python -c "import example; example.greet()"
Hello, World!

しかし、当然のことながら別の場所に移動してしまうと使えなくなる。 これは、先ほどのパッケージに Python のパスが通っていないため。

$ cd /tmp && python -c "import example; example.greet()" 
Traceback (most recent call last):
  File "<string>", line 1, in <module>
ImportError: No module named example

もちろん、明示的にパスを通してしまえば使うことはできる。 とはいえ、毎回こんなことをしたくはないはず。

$ PYTHONPATH=$HOME/workplace/packaging python -c "import example; example.greet()"
Hello, World!

そこで、作成したライブラリに必要なもの一式をまとめた上で、あらかじめパスの通っている場所に配置できるようにする仕組みが、今回扱うパッケージングというわけ。

サンプルコードをパッケージングする

それでは、先ほど作った example をパッケージングしてみよう。 setuptools を使ったパッケージングでは、まず setup.py という名前の Python ファイルを用意する。 このファイルは、セットアップスクリプトと呼ばれる。

$ cat << 'EOF' > setup.py 
#!/usr/bin/env python
# -*- coding: utf-8 -*-

from setuptools import setup


setup()
EOF

セットアップスクリプトの中では setuptools に含まれる関数 setup() を呼び出している。 引数には何も渡しておらず、これは比較的新しい書き方をしている。 古い書き方では setup() にパッケージに関する様々な情報を引数で渡すことになっていた。

では、パッケージに関する様々な情報は、代わりに何処で渡すのかというと setup.cfg という設定ファイルに記述する。 設定ファイルは ini フォーマットっぽい形式で書いていく。 パッケージの基本的な情報は [metadata] というセクションで扱う。

$ cat << 'EOF' > setup.cfg 
[metadata]
name = example
version = attr:example.__version__
author = Your Name
author_email = your-email@example.com
description = example is a example package
long_description = file:README.md
url = https://example.com/your/project/page
license = Apache License, Version 2.0
classifier =
    Development Status :: 1 - Planning
    Programming Language :: Python
    Intended Audience :: Developers
    Operating System :: POSIX
    Programming Language :: Python :: 2
    Programming Language :: Python :: 2.7
    Programming Language :: Python :: 3
    Programming Language :: Python :: 3.5
    Programming Language :: Python :: 3.6
    Programming Language :: Python :: 3.7
    
[options]
zip_safe = False
packages = find:
EOF

なお、この setup.cfg にメタデータなどを記述していく方式は setuptools のバージョン 30.0.3.0 (リリース日 2016年12月8日) から導入された。 この点は、インストール先となる環境の setuptools にもバージョン 30.0.3.0 以降が求められる点に注意が必要になる。 それより前のバージョンをサポートしたいときは、従来通り setup.py の中にメタデータなどを関数の引数として渡す。

setuptools.readthedocs.io

また、setup.cfg では、いくつかの項目において外部のファイルを読み込むこともできる。 例えば、先ほどの例において long_description では README.md の内容を読み込むよう指定している。 なので、ファイルを作っておこう。

$ cat << 'EOF' > README.md 
### What is this?

A long description for example project.
EOF

さて、これでパッケージングの準備ができた。 ディレクトリは、このような状態になっている。

$ ls
README.md   example     setup.cfg   setup.py

試しに pip を使ってパッケージをインストールしてみよう。 この操作は、ライブラリに必要なもの一式をまとめる作業と、インストールする作業を一気にやっているのに等しい。

$ pip install -U .

うまくいけば、次のように pip のインストール済みパッケージの一覧に表示される。

$ pip list | grep example
example    0.0.1  
$ pip show example
Name: example
Version: 0.0.1
Summary: example is a example package
Home-page: https://example.com/your/project/page
Author: Your Name
Author-email: your-email@example.com
License: Apache License, Version 2.0
Location: /Users/amedama/.virtualenvs/py37/lib/python3.7/site-packages
Requires: 
Required-by: 

また、全然別の場所でパッケージをインポートして使うこともできるようになる。 これは、pip によって作成したパッケージがパスの通った場所に配置されたため。

$ cd /tmp && python -c "import example; example.greet()"
Hello, World!

先ほど pip show したときの Location にあった通り、今回は以下のディレクトリにインストールされていた。

$ cat ~/.virtualenvs/py37/lib/python3.7/site-packages/example/__init__.py 
# -*- coding: utf-8 -*-

__version__ = '0.0.1'


def greet():
    """挨拶をする関数"""
    print('Hello, World!')

動作確認が終わったら、一旦パッケージをアンインストールしておこう。

$ pip uninstall -y example

ソースコード配布物 (sdist) にパッケージングする

先ほどはライブラリに必要なもの一式をまとめるのとインストールするのを pip コマンドで一気にやってしまった。 続いては、この工程を別々に分けてやってみよう。

まずは、基本となるソースコード配布物 (sdist) へのパッケージングから。 パッケージングは先ほど作ったセットアップスクリプトを使って sdist コマンドを実行する。

$ python setup.py sdist

すると dist というディレクトリ以下に tarball ができる。 これこそ、できあがったソースコード配布物のファイル。

$ ls dist 
example-0.0.1.tar.gz

このソースコード配布物からパッケージをインストールできる。

$ pip install dist/example-0.0.1.tar.gz
$ pip list | grep example              
example    0.0.1  

確認できたら、また一旦アンインストールしておこう。

$ pip uninstall -y example

Wheel フォーマットでパッケージングする

ソースコード配布物は基本ではあるものの、使うとちょっとめんどくさい場合もある。 具体的には Python/C API を使った拡張モジュールがライブラリに含まれているパターン。 ソースコード配布物では、拡張モジュールをビルドしていないソースコードの状態で同梱する。 すると、インストール先の環境にはそれをビルドするための開発ツール類一式が必要になる。

そこで、最近は次世代のパッケージング形式として Wheel フォーマットというものが使われ始めている。 このフォーマットでは拡張モジュールはあらかじめビルドされた状態で同梱される。 ただし、環境に依存する場合には各環境ごとに Wheel ファイルを用意する必要がある。

実際に Wheel ファイルをパッケージングしてみよう。 これにはセットアップスクリプトで bdist_wheel コマンドを実行する。

$ python setup.py bdist_wheel

すると、次の通り末尾が whl の Wheel ファイルが dist ディレクトリ以下に作られる。

$ ls dist | grep whl$
example-0.0.1-py3-none-any.whl

もちろんこの Wheel ファイルからもパッケージをインストールできる。

$ pip install dist/example-0.0.1-py3-none-any.whl
$ pip list | grep example                        
example    0.0.1

ちなみに、上記の Wheel ファイルは名前を見て分かる通り Python 3 専用としてビルドされている。 もし、Python 2 でも動く場合には setup.cfg の [wheel] セクションについて universal フラグを有効にする。

$ cat << 'EOF' > setup.cfg 
[metadata]
name = example
version = attr:example.__version__
author = Your Name
author_email = your-email@example.com
description = example is a example package
long_description = file:README.md
url = https://example.com/your/project/page
license = Apache License, Version 2.0
classifier =
    Development Status :: 1 - Planning
    Programming Language :: Python
    Intended Audience :: Developers
    Operating System :: POSIX
    Programming Language :: Python :: 2
    Programming Language :: Python :: 2.7
    Programming Language :: Python :: 3
    Programming Language :: Python :: 3.5
    Programming Language :: Python :: 3.6
    Programming Language :: Python :: 3.7
    
[options]
zip_safe = False
packages = find:

[wheel]
universal = 1
EOF

この状態でビルドすると Python 2/3 両対応の Wheel ファイルになる。

$ python setup.py bdist_wheel
$ ls dist | grep whl$
example-0.0.1-py2.py3-none-any.whl
example-0.0.1-py3-none-any.whl

動作が確認できたら、また一旦アンインストールしておこう。

$ pip uninstall -y example

パッケージングした成果物を PyPI に登録する

さて、パッケージングできるようになると、実は成果物を PyPI に登録して一般に公開できる。 この作業には PyPA 製の Twine というツールを使うのがおすすめ。 詳しくは、以下を参照のこと。

blog.amedama.jp

依存ライブラリを指定する

さて、話をちょっと戻して、ここからはパッケージングにおける色々なユースケースについて見ていく。

まずは、自作ライブラリが別のライブラリに依存している場合について。 この場合は [options] セクションに install_requires という項目で指定する。 例えば requests を依存ライブラリとして追加してみよう。

$ cat << 'EOF' > setup.cfg 
[metadata]
name = example
version = attr:example.__version__
author = Your Name
author_email = your-email@example.com
description = example is a example package
long_description = file:README.md
url = https://example.com/your/project/page
license = Apache License, Version 2.0
classifier =
    Development Status :: 1 - Planning
    Programming Language :: Python
    Intended Audience :: Developers
    Operating System :: POSIX
    Programming Language :: Python :: 2
    Programming Language :: Python :: 2.7
    Programming Language :: Python :: 3
    Programming Language :: Python :: 3.5
    Programming Language :: Python :: 3.6
    Programming Language :: Python :: 3.7
    
[options]
zip_safe = False
packages = find:
install_requires =
  requests

[wheel]
universal = 1
EOF

この状態でインストールしてみよう。

$ pip install -U .

すると、次のように requests が一緒にインストールされたことが分かる。

$ pip list | egrep "(example|requests)"
example    0.0.1     
requests   2.21.0    

環境に依存した依存ライブラリを指定する

続いては環境によって必要だったり不要だったりする依存ライブラリの指定について。 例えばリレーショナルデータベースを扱うアプリケーションだと、接続先が MySQL なのか PostgreSQL なのかで必要なドライバが変わる。

例として接続先が MySQL の環境なら mysqlclient をインストールする、という状況を想定してみよう。 この場合、[options.extras_require] セクションに環境名と必要なライブラリを指定していく。

$ cat << 'EOF' > setup.cfg 
[metadata]
name = example
version = attr:example.__version__
author = Your Name
author_email = your-email@example.com
description = example is a example package
long_description = file:README.md
url = https://example.com/your/project/page
license = Apache License, Version 2.0
classifier =
    Development Status :: 1 - Planning
    Programming Language :: Python
    Intended Audience :: Developers
    Operating System :: POSIX
    Programming Language :: Python :: 2
    Programming Language :: Python :: 2.7
    Programming Language :: Python :: 3
    Programming Language :: Python :: 3.5
    Programming Language :: Python :: 3.6
    Programming Language :: Python :: 3.7
    
[options]
zip_safe = False
packages = find:
install_requires =
  requests

[options.extras_require]
mysql =
  mysqlclient

[wheel]
universal = 1
EOF

この状態で、環境に mysql を指定してパッケージをアップデートしてみよう。

$ pip install -U ".[mysql]"

すると、mysqlclient がインストールされることが分かる。

$ pip list | grep mysqlclient
mysqlclient 1.4.1     

ソースコード以外のファイルを成果物に含める

続いては、ソースコード以外のファイルを成果物に含める方法について。 実は、デフォルトでは成果物にソースコード以外のファイルは含まれない。 試しに実験してみよう。

まずはパッケージに greet_from_file() という関数を追加する。 これはテキストファイルから読み込んだ内容を使ってメッセージをプリントする関数になっている。

$ cat << 'EOF' > example/__init__.py
# -*- coding: utf-8 -*-

from __future__ import print_function

import os

__version__ = '0.0.1'


def greet():
    """挨拶をする関数"""
    print('Hello, World!')


def greet_from_file():
    """テキストファイルの内容を使って挨拶する関数"""
    module_dir = os.path.dirname(__file__)
    filepath = os.path.join(module_dir, 'message.txt')
    with open(filepath, mode='r') as fp:
        print(fp.read(), end='')
EOF

続いて、表示するメッセージの元となるファイルを用意する。

$ cat << 'EOF' > example/message.txt 
Hello, World!
EOF

ローカルでは、これで動くようになる。

$ python -c "import example; example.greet_from_file()"
Hello, World!

では、インストールした上ではどうだろうか。

$ pip install -U .

試すと、そんなファイルないよと怒られる。 パッケージングした成果物の中にテキストファイルが含まれていないためだ。

$ cd /tmp && python -c "import example; example.greet_from_file()"
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "/Users/amedama/.virtualenvs/py37/lib/python3.7/site-packages/example/__init__.py", line 19, in greet_from_file
    with open(filepath, mode='r') as fp:
FileNotFoundError: [Errno 2] No such file or directory: '/Users/amedama/.virtualenvs/py37/lib/python3.7/site-packages/example/message.txt'

ソースコード以外のファイルを成果物に含めるには、まず setup.cfg の [options] セクションにおいて include_package_data の項目を有効にする。

$ cat << 'EOF' > setup.cfg 
[metadata]
name = example
version = attr:example.__version__
author = Your Name
author_email = your-email@example.com
description = example is a example package
long_description = file:README.md
url = https://example.com/your/project/page
license = Apache License, Version 2.0
classifier =
    Development Status :: 1 - Planning
    Programming Language :: Python
    Intended Audience :: Developers
    Operating System :: POSIX
    Programming Language :: Python :: 2
    Programming Language :: Python :: 2.7
    Programming Language :: Python :: 3
    Programming Language :: Python :: 3.5
    Programming Language :: Python :: 3.6
    Programming Language :: Python :: 3.7
    
[options]
zip_safe = False
include_package_data = True
packages = find:
install_requires =
  requests

[options.extras_require]
mysql =
  mysqlclient

[wheel]
universal = 1
EOF

その上で MANIFEST.in というファイルを用意する。 ここに、成果物に含めるソースコード以外のファイルを指定していく。

$ cat << 'EOF' > MANIFEST.in 
include example/message.txt
EOF

この状態で、もう一度試してみよう。

$  pip install -U .
$ cd /tmp && python -c "import example; example.greet_from_file()"
Hello, World!

今度はエラーにならず実行できた。

コマンドを追加する

続いてはパッケージを追加したときに新たにコマンドが使えるようにする方法について。

今回はコマンド用に別の設定ファイルを用意することにした。 まずは setup.cfg の [options] セクションに entry_points という項目を作る。

$ cat << 'EOF' > setup.cfg 
[metadata]
name = example
version = attr:example.__version__
author = Your Name
author_email = your-email@example.com
description = example is a example package
long_description = file:README.md
url = https://example.com/your/project/page
license = Apache License, Version 2.0
classifier =
    Development Status :: 1 - Planning
    Programming Language :: Python
    Intended Audience :: Developers
    Operating System :: POSIX
    Programming Language :: Python :: 2
    Programming Language :: Python :: 2.7
    Programming Language :: Python :: 3
    Programming Language :: Python :: 3.5
    Programming Language :: Python :: 3.6
    Programming Language :: Python :: 3.7
    
[options]
zip_safe = False
include_package_data = True
packages = find:
entry_points = file:entry_points.cfg
install_requires =
  requests

[options.extras_require]
mysql =
  mysqlclient

[wheel]
universal = 1
EOF

その上で、指定したのと同じ設定ファイルを用意する。 以下では example-greet というコマンドを実行したとき example パッケージの greet() 関数が呼ばれるようにしている。

$ cat << 'EOF' > entry_points.cfg
[console_scripts]
example-greet = example:greet
EOF

それでは、パッケージを更新してみよう。

$ pip install -U .

すると example-greet というコマンドが使えるようになっている。

$ example-greet
Hello, World!

テストコードを成果物に含めない

パッケージングしたくなるようなコードには、もちろんテストコードを書くことになる。 続いてはテストコードを成果物に含めない方法とともに、具体的なオペレーションを紹介してみる。

まずはテストをする環境向けの依存ライブラリを [options.extras_require] セクションで指定する。 テストフレームワークには pytest を使うことにした。 同時に、テストを成果物に含めないように [options.packages.find] セクションで tests というディレクトリを探索対象から除外する。

$ cat << 'EOF' > setup.cfg 
[metadata]
name = example
version = attr:example.__version__
author = Your Name
author_email = your-email@example.com
description = example is a example package
long_description = file:README.md
url = https://example.com/your/project/page
license = Apache License, Version 2.0
classifier =
    Development Status :: 1 - Planning
    Programming Language :: Python
    Intended Audience :: Developers
    Operating System :: POSIX
    Programming Language :: Python :: 2
    Programming Language :: Python :: 2.7
    Programming Language :: Python :: 3
    Programming Language :: Python :: 3.5
    Programming Language :: Python :: 3.6
    Programming Language :: Python :: 3.7
    
[options]
zip_safe = False
include_package_data = True
packages = find:
entry_points = file:entry_points.cfg
install_requires =
  requests

[options.extras_require]
mysql =
  mysqlclient
testing =
  pytest

[options.packages.find]
exclude =
  tests

[wheel]
universal = 1
EOF

続いてはテストコードを入れるパッケージを用意する。 中身でやっているのは、あまり本質的ではないことなので気にしなくて良いと思う。

$ mkdir tests
$ touch tests/__init__.py
$ cat << 'EOF' > tests/test_example.py   
#!/usr/bin/env python
# -*- coding: utf-8 -*-

import os
import sys

import pytest

from example import greet

if sys.version_info >= (3, 0, 0):
    # Python 3
    from io import StringIO
else:
    # Python 2
    from StringIO import StringIO


def test_greet():
    """挨拶は大事"""
    # 標準出力を入れ替える
    io = StringIO()
    sys.stdout = io

    # 挨拶する
    greet()

    # 挨拶した内容を確認する
    assert io.getvalue() == ('Hello, World!' + os.linesep)


if __name__ == '__main__':
    pytest.main([__file__])
EOF

準備ができたら、テスト環境向けにインストールする。

$ pip install -U ".[testing]"

修正したときの手間などを考えるとパッケージ本体は普段はアンインストールしておいた方が良いかも。

$ pip uninstall -y example

テストを走らせるにはテストランナーを実行する。

$ py.test
====================================== test session starts ======================================
platform darwin -- Python 3.7.2, pytest-4.2.0, py-1.7.0, pluggy-0.8.1
rootdir: /Users/amedama/Documents/temporary/packaging, inifile:
collected 1 item                                                                                

tests/test_example.py .                                                                   [100%]

=================================== 1 passed in 0.02 seconds ====================================

テストが実行できることが分かったので、続いては成果物の中身を確認してみよう。

$ rm -rf dist
$ python setup.py sdist
$ tar xvf dist/example-0.0.1.tar.gz -C dist
$ find dist/example-0.0.1 
dist/example-0.0.1
dist/example-0.0.1/PKG-INFO
dist/example-0.0.1/example.egg-info
dist/example-0.0.1/example.egg-info/PKG-INFO
dist/example-0.0.1/example.egg-info/not-zip-safe
dist/example-0.0.1/example.egg-info/SOURCES.txt
dist/example-0.0.1/example.egg-info/entry_points.txt
dist/example-0.0.1/example.egg-info/requires.txt
dist/example-0.0.1/example.egg-info/top_level.txt
dist/example-0.0.1/example.egg-info/dependency_links.txt
dist/example-0.0.1/example
dist/example-0.0.1/example/__init__.py
dist/example-0.0.1/example/message.txt
dist/example-0.0.1/MANIFEST.in
dist/example-0.0.1/README.md
dist/example-0.0.1/setup.py
dist/example-0.0.1/setup.cfg

成果物の中には tests というディレクトリが含まれていないことが分かる。

まとめ

• Python のパッケージングには setuptools というライブラリを使う
• パッケージングするときはセットアップスクリプト (setup.py) を用意する
• 最近 (>= 30.0.3.0) の書き方では、なるべく別の設定ファイル (setup.cfg) に内容を記述する
• パッケージングの成果物にはソースコード配布物 (sdist) と Wheel (whl) がある

スマートPythonプログラミング: Pythonのより良い書き方を学ぶ

Python: Unix における python* コマンドと処理系のバージョンについて (PEP394)

今回は Unix ライクなシステムにおける python* コマンドの振る舞いと、処理系のバージョンについて色々と書いてみる。

きっかけは、以下のブログ記事を目にしたため。

rheb.hatenablog.com

上記のブログでは Red Hat Enterprise Linux 8 に搭載される予定の Python と、それを扱うコマンド群について紹介されている。 そして、RHEL8 がそのような仕様に落ち着いた理由として以下の PEP (Python Enhancement Proposal) が挙げられている。

www.python.org

この PEP394 というドキュメントには、Unix ライクのシステムにおいて推奨される python* コマンドの振る舞いなどが記述されている。 そこで、今回は PEP394 に書かれている内容について紹介したい。 先のブログでは、話の切り口が「RHEL8 における仕様」なのでやむを得ないものの、やや説明が不足しており補足したい点もあった。

python* コマンドで実行される処理系のバージョンについて

まず、PEP394 には python* コマンドを実行したときに、どのバージョンの処理系が呼び出されるべきなのか書かれている。 例えば、次のようにコマンド名でバージョンを明示的に指定したときは、そのバージョンの処理系が呼ばれることになる。 もちろん「そのバージョンの処理系がインストールされていれば」という前提のもとで。

$ python2  # => Python 2.x の処理系が使われる
$ python3  # => Python 3.x の処理系が使われる

問題となるのは、バージョンが明示的に指定されていないとき。 一体、どのバージョンの処理系が呼び出されるべきなのだろうか。 その答えは「現時点では Python 2.x の処理系が使われることが推奨されている」になる。

$ python  # => 現時点では Python 2.x の処理系が使われることが推奨されている

現時点では？

太字にした通り、この挙動が推奨されているのは、あくまで現時点 (2019-01) の PEP394 ではという注釈つきになる。 なぜなら PEP394 には Future Changes to this Recommendation (この勧告に対する将来の変更) という項目があるため。

この項目には、次のような記述がある。(カッコ内は拙訳)

This recommendation will be periodically reviewed over the next few years, and updated when the core development team judges it appropriate. As a point of reference, regular maintenance releases for the Python 2.7 series will continue until at least 2020. (この勧告は将来の数年間に渡って定期的に見直されます。そして、コア開発チームが妥当と判断したときに更新されます。参考までに、Python 2.7 系の定常的なメンテナンスリリースは少なくとも 2020 年まで続きます。)

つまり、将来的にはバージョン指定なしの python コマンドで Python 3.x の処理系が使われることを推奨する内容に更新される余地がある。

例外について

ところで、上記の振る舞いは一部のディストリビューションと仮想環境では例外として扱われている。

一部のディストリビューションというのは、例えばこの PEP が勧告される以前から Python 3 移行を始めた、ある意味で優秀なコミュニティもあったりするため。

www.archlinux.org

もう一つの仮想環境というのは venv や virtualenv を使った環境を指している。 仮想環境の中では、環境を作るときに指定されたか、あるいはデフォルトの処理系が python コマンドで使われる。

$ python3 -m venv myvenv 
$ source myvenv/bin/activate
(myvenv) $ python -V
Python 3.7.2

シェバン (Shebang) について

PEP394 には、シェバン (Shebang) に関する記載もある。

シェバンにおいては、特定のバージョンの Python でしか動かないモジュールならバージョンを指定することが推奨されている。 例えば、Python 2 でしか動かないなら以下のようにする。

#!/usr/bin/env python2
...

Python 3 でしか動かないなら、こう。

#!/usr/bin/env python3
...

バージョン指定のないシェバンは、モジュールが両方のバージョンで動く場合に用いることが推奨されている。

#!/usr/bin/env python
...

インタプリタから別のインタプリタを起動する場合

また、subprocess なんかでインタプリタの中から別のインタプリタを呼ぶときは sys.executable の内容を使うのが良い。 ここには自身の処理系のパスが格納されている。

$ python3 -c "import sys; print(sys.executable)"
/usr/local/opt/python/bin/python3.7

特定のコマンド名 (python2, python3, python, etc) を使うと、誤って別の処理系を呼び出したりして思わぬバグが生じる恐れがある。

まとめ

• 前提: PEP394 の勧告は将来的に変更される余地がある
• 各コマンドで起動する処理系のバージョン
  • python2 # => Python 2.x
  • python3 # => Python 3.x
  • python # => 今のところ Python 2.x
• 上記の例外
  • 一部のディストリビューション (ArchLinux など)
  • 仮想環境 (venv, virtualenv)
• シェバン (Shebang) を書くとき
  • python2 # => Python 2.x のみで動作する
  • python3 # => Python 3.x のみで動作する
  • python # => 両方のバージョンで動作する
• インタプリタの中から別のインタプリタを起動するとき
  • sys.executable を使う

いじょう。

スマートPythonプログラミング: Pythonのより良い書き方を学ぶ

Python: XGBoost を使ってみる

XGBoost (eXtreme Gradient Boosting) は勾配ブースティング決定木 (Gradient Boosting Decision Tree) のアルゴリズムを実装したオープンソースのライブラリ。 最近は、同じ GBDT 系のライブラリである LightGBM にややお株を奪われつつあるものの、依然として機械学習コンペティションの一つである Kaggle でよく使われている。 今回は、そんな XGBoost の Python バインディングを使ってみることにする。

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.14.2
BuildVersion:   18C54
$ python -V
Python 3.7.2

もくじ

下準備

最初に、下準備として XGBoost と必要なパッケージを一通りインストールする。

そのために、まずは gcc をインストールしておく。 これは XGBoost が並列処理に OpenMP を使っている関係で必要になる。

$ brew install gcc@7

続いて XGBoost と、それ以外で今回使うパッケージをインストールする。

$ pip install xgboost scikit-learn matplotlib

乳がんデータセットを分類してみる

まずはハローワールド的な例として乳がんデータセットを使った二値分類 (Binary classification) から始める。

以下が XGBoost を使って乳がんデータセットを二値分類するサンプルコード。 なお、モデルの検証についてはここでの本題ではないことから、交差検証ではなくホールドアウト検証にとどめている。 最終的に検証用データで精度 (Accuracy) を確認している。 コードの説明についてはコメントを参照のこと。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import xgboost as xgb

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

"""XGBoost で二値分類するサンプルコード"""


def main():
    # 乳がんデータセットを読み込む
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target
    # データセットを学習用と検証用に分割する
    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)
    # XGBoost が扱うデータセットの形式に直す
    dtrain = xgb.DMatrix(X_train, label=y_train)
    dtest = xgb.DMatrix(X_test, label=y_test)
    # 学習用のパラメータ
    xgb_params = {
        # 二値分類問題
        'objective': 'binary:logistic',
        # 評価指標
        'eval_metric': 'logloss',
    }
    # モデルを学習する
    bst = xgb.train(xgb_params,
                    dtrain,
                    num_boost_round=100,  # 学習ラウンド数は適当
                    )
    # 検証用データが各クラスに分類される確率を計算する
    y_pred_proba = bst.predict(dtest)
    # しきい値 0.5 で 0, 1 に丸める
    y_pred = np.where(y_pred_proba > 0.5, 1, 0)
    # 精度 (Accuracy) を検証する
    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)


if __name__ == '__main__':
    main()

上記のサンプルコードに適当な名前をつけて保存した上で実行する。

$ python xgbhelloworld.py
[22:52:44] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 18 extra nodes, 0 pruned nodes, max_depth=5
[22:52:44] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 14 extra nodes, 0 pruned nodes, max_depth=4
[22:52:44] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 14 extra nodes, 0 pruned nodes, max_depth=4
...(snip)...
[22:52:44] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 2 extra nodes, 0 pruned nodes, max_depth=1
[22:52:44] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 2 extra nodes, 0 pruned nodes, max_depth=1
[22:52:44] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 2 extra nodes, 0 pruned nodes, max_depth=1
Accuracy: 0.9649122807017544

ホールドアウト検証の結果、精度として 0.965 前後の値が得られた。

学習過程を可視化する

先ほどの例では、いつの間にか学習が終わってモデルができたという感じだった。 そこで、続いては学習が進む過程をグラフで可視化してみる。

次のサンプルコードでは、学習用データと検証用データに対する損失を折れ線グラフで出力する。 そのためには、まず evals オプションで学習用データと検証用データを渡す。 その上で evals_result オプションに過程を記録するための辞書を渡す。 学習ラウンド数 (num_boost_round) は先ほどよりも多く 1000 まで増やした。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import xgboost as xgb

import numpy as np
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from matplotlib import pyplot as plt

"""XGBoost で学習の履歴を可視化するサンプルコード"""


def main():
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)

    dtrain = xgb.DMatrix(X_train, label=y_train)
    dtest = xgb.DMatrix(X_test, label=y_test)

    xgb_params = {
        'objective': 'binary:logistic',
        'eval_metric': 'logloss',
    }

    # 学習時に用いる検証用データ
    evals = [(dtrain, 'train'), (dtest, 'eval')]
    # 学習過程を記録するための辞書
    evals_result = {}
    bst = xgb.train(xgb_params,
                    dtrain,
                    num_boost_round=1000,  # ラウンド数を増やしておく
                    evals=evals,
                    evals_result=evals_result,
                    )

    y_pred_proba = bst.predict(dtest)
    y_pred = np.where(y_pred_proba > 0.5, 1, 0)
    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)

    # 学習の課程を折れ線グラフとしてプロットする
    train_metric = evals_result['train']['logloss']
    plt.plot(train_metric, label='train logloss')
    eval_metric = evals_result['eval']['logloss']
    plt.plot(eval_metric, label='eval logloss')
    plt.grid()
    plt.legend()
    plt.xlabel('rounds')
    plt.ylabel('logloss')
    plt.show()


if __name__ == '__main__':
    main()

先ほどと同じように、適当な名前をつけて実行する。

$ python xgbhistory.py
...(snip)...
Accuracy: 0.9649122807017544

すると、次のようなグラフが得られる。

学習が進むにつれて学習用データと検証用データの損失が減っていくことが分かる。 とはいえ、過学習はしていないようだけど 100 ラウンド前後で損失の減少は止まっているように見える。

損失が減らなくなったら学習を打ち切る

先ほどの例で分かる通り、損失が減らなくなったらそれ以上学習する必要はない。 それ以上回してしまうと、下手をすると損失が増える (過学習) することも考えられる。 そこで、続いては損失が減らなくなったタイミングで学習を打ち切るようにしてみる。

次のサンプルコードでは early_stopping_rounds オプションを指定することでそれを実現している。 このオプションを使うと、指定したラウンド数の間で損失が減らなかったときに学習を打ち切れる。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import xgboost as xgb

import numpy as np
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from matplotlib import pyplot as plt

"""XGBoost で early_stopping_rounds を使って学習ラウンド数を最適化するサンプルコード"""


def main():
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)

    dtrain = xgb.DMatrix(X_train, label=y_train)
    dtest = xgb.DMatrix(X_test, label=y_test)

    xgb_params = {
        'objective': 'binary:logistic',
        'eval_metric': 'logloss',
    }

    evals = [(dtrain, 'train'), (dtest, 'eval')]
    evals_result = {}
    bst = xgb.train(xgb_params,
                    dtrain,
                    num_boost_round=1000,
                    # 一定ラウンド回しても改善が見込めない場合は学習を打ち切る
                    early_stopping_rounds=10,
                    evals=evals,
                    evals_result=evals_result,
                    )

    y_pred_proba = bst.predict(dtest)
    y_pred = np.where(y_pred_proba > 0.5, 1, 0)
    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)

    train_metric = evals_result['train']['logloss']
    plt.plot(train_metric, label='train logloss')
    eval_metric = evals_result['eval']['logloss']
    plt.plot(eval_metric, label='eval logloss')
    plt.grid()
    plt.legend()
    plt.xlabel('rounds')
    plt.ylabel('logloss')
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。 ログ出力から、検証用データの損失に対して early stopping がかかるようになっていることが分かる。

$ python xgbearlystop.py
[23:05:43] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 18 extra nodes, 0 pruned nodes, max_depth=5
[0]    train-logloss:0.462396 eval-logloss:0.492899
Multiple eval metrics have been passed: 'eval-logloss' will be used for early stopping.

Will train until eval-logloss hasn't improved in 10 rounds.
[23:05:43] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 14 extra nodes, 0 pruned nodes, max_depth=4
...(snip)...
[23:05:44] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 2 extra nodes, 0 pruned nodes, max_depth=1
[100] train-logloss:0.006161  eval-logloss:0.09275
[23:05:44] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 2 extra nodes, 0 pruned nodes, max_depth=1
[101] train-logloss:0.00614   eval-logloss:0.093094
Stopping. Best iteration:
[91]  train-logloss:0.00637   eval-logloss:0.092265

Accuracy: 0.9649122807017544

今度は 101 ラウンドまでしか学習が進んでいない。 また、その中でも最も損失の少なかった 91 ラウンド目が結果として使われたことが分かる。

学習過程のグラフは次のようなものが得られた。

たしかに 100 前後の学習ラウンドまでしか表示されていない。

scikit-learn インターフェースを使ってみる

XGBoost には、ネイティブな API の他に scikit-learn 互換の API を持ったインターフェースもある。 続いては、これを使ってみよう。

次のサンプルコードでは scikit-learn 互換の API を備えた XGBClassifier を使っている。 ただし、実現していることは先ほどと全く変わらない。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import xgboost as xgb

from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from matplotlib import pyplot as plt

"""XGBoost の scikit-learn インターフェースを使ったサンプルコード (二値分類)"""


def main():
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)

    # scikit-learn API を備えた分類器
    clf = xgb.XGBClassifier(objective='binary:logistic',
                            # 'num_boost_round' の代わり
                            # adding 1 estimator per round
                            n_estimators=1000)
    # 学習する
    evals_result = {}
    clf.fit(X_train, y_train,
            # 学習に使う評価指標
            eval_metric='logloss',
            # 学習時に用いる検証用データ
            eval_set=[
                (X_train, y_train),
                (X_test, y_test),
            ],
            early_stopping_rounds=10,
            # 学習過程の記録はコールバック API で登録する
            callbacks=[
                xgb.callback.record_evaluation(evals_result)
            ],
            )

    y_pred = clf.predict(X_test)
    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)

    # 学習過程の名前は 'validation_{n}' になる
    train_metric = evals_result['validation_0']['logloss']
    plt.plot(train_metric, label='train logloss')
    eval_metric = evals_result['validation_1']['logloss']
    plt.plot(eval_metric, label='eval logloss')
    plt.grid()
    plt.legend()
    plt.xlabel('rounds')
    plt.ylabel('logloss')
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python xgbscikit.py
[0]    validation_0-logloss:0.607962  validation_1-logloss:0.615844
Multiple eval metrics have been passed: 'validation_1-logloss' will be used for early stopping.

Will train until validation_1-logloss hasn't improved in 10 rounds.
[1]   validation_0-logloss:0.538811   validation_1-logloss:0.555234
[2]   validation_0-logloss:0.480826   validation_1-logloss:0.5015
[3]   validation_0-logloss:0.430842   validation_1-logloss:0.458286
...(snip)...
[170] validation_0-logloss:0.008248   validation_1-logloss:0.094108
Stopping. Best iteration:
[160] validation_0-logloss:0.008465   validation_1-logloss:0.094088

Accuracy: 0.9649122807017544

学習ラウンド数は先ほどと違っているものの、最終的に得られた精度は変わっていないようだ。

学習過程のグラフは次のようなものが得られた。

多値分類問題を扱う

続いてはタスク設定として多値分類 (Multiclass classification) を試してみよう。 この場合、二値分類とは学習するときのパラメータが異なる。

次のサンプルコードでは XGBoost で Iris データセットを使った多値分類問題を扱っている。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import xgboost as xgb

from matplotlib import pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

"""XGBoost で多値分類するサンプルコード"""


def main():
    # Iris データセットを読み込む
    dataset = datasets.load_iris()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)

    dtrain = xgb.DMatrix(X_train, label=y_train)
    dtest = xgb.DMatrix(X_test, label=y_test)

    xgb_params = {
        # 多値分類問題
        'objective': 'multi:softmax',
        # クラス数
        'num_class': 3,
        # 学習用の指標 (Multiclass logloss)
        'eval_metric': 'mlogloss',
    }
    evals = [(dtrain, 'train'), (dtest, 'eval')]
    evals_result = {}
    bst = xgb.train(xgb_params,
                    dtrain,
                    num_boost_round=1000,
                    early_stopping_rounds=10,
                    evals=evals,
                    evals_result=evals_result,
                    )

    y_pred = bst.predict(dtest)
    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)

    train_metric = evals_result['train']['mlogloss']
    plt.plot(train_metric, label='train logloss')
    eval_metric = evals_result['eval']['mlogloss']
    plt.plot(eval_metric, label='eval logloss')
    plt.grid()
    plt.legend()
    plt.xlabel('rounds')
    plt.ylabel('logloss')
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python xgbmulticlass.py
[23:14:48] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 2 extra nodes, 0 pruned nodes, max_depth=1
[23:14:48] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 8 extra nodes, 0 pruned nodes, max_depth=3
[23:14:48] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 4 extra nodes, 0 pruned nodes, max_depth=2
[0]    train-mlogloss:0.742287    eval-mlogloss:0.765776
Multiple eval metrics have been passed: 'eval-mlogloss' will be used for early stopping.

Will train until eval-mlogloss hasn't improved in 10 rounds.
...
Stopping. Best iteration:
[16]  train-mlogloss:0.039852 eval-mlogloss:0.195911

Accuracy: 0.9333333333333333

乳がんデータセットに比べると、単純な分だけラウンド数がかなり少ないようだ。

得られた学習過程のグラフは次の通り。 学習用データの損失は減っているものの、検証用データの損失が減らない状況が生じていることから過学習の予兆が見られる。

特徴量の重要度を可視化する

XGBoost は決定木の仲間ということで特徴量の重要度 (Feature Importance) を可視化する機能を備えている。 次のサンプルコードでは、Iris データセットの分類にどの特徴量が有効だったのかを性能のゲインにもとづいて可視化している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import xgboost as xgb

from sklearn import datasets
from sklearn.model_selection import train_test_split
from matplotlib import pyplot as plt

"""XGBoost で特徴量の重要度を可視化するサンプルコード"""


def main():
    dataset = datasets.load_iris()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)

    # 可視化のために特徴量の名前を渡しておく
    dtrain = xgb.DMatrix(X_train, label=y_train,
                         feature_names=dataset.feature_names)
    dtest = xgb.DMatrix(X_test, label=y_test,
                        feature_names=dataset.feature_names)

    xgb_params = {
        'objective': 'multi:softmax',
        'num_class': 3,
        'eval_metric': 'mlogloss',
    }

    evals = [(dtrain, 'train'), (dtest, 'eval')]
    evals_result = {}
    bst = xgb.train(xgb_params,
                    dtrain,
                    num_boost_round=1000,
                    early_stopping_rounds=10,
                    evals=evals,
                    evals_result=evals_result,
                    )

    # 性能向上に寄与する度合いで重要度をプロットする
    _, ax = plt.subplots(figsize=(12, 4))
    xgb.plot_importance(bst,
                        ax=ax,
                        importance_type='gain',
                        show_values=False)
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python xgbfeatimp.py
[23:19:37] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 2 extra nodes, 0 pruned nodes, max_depth=1
[23:19:37] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 8 extra nodes, 0 pruned nodes, max_depth=3
[23:19:37] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 4 extra nodes, 0 pruned nodes, max_depth=2
[0]    train-mlogloss:0.742287    eval-mlogloss:0.765776
Multiple eval metrics have been passed: 'eval-mlogloss' will be used for early stopping.

Will train until eval-mlogloss hasn't improved in 10 rounds.
...(snip)...
Stopping. Best iteration:
[16]  train-mlogloss:0.039852 eval-mlogloss:0.195911

次のようなグラフが得られる。 どうやら petal length と petal width が分類する上で有効なようだ。

回帰問題を扱ってみる

続いては XGBoost で回帰問題を扱ってみる。 回帰問題を扱うときは学習時のパラメータとして渡す objective が reg から始まるようになる。

次のサンプルコードでは XGBoost で Boston データセットを回帰している。 学習と検証の評価指標には RMSE (Root Mean Squared Error) を用いた。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import math

import xgboost as xgb

from matplotlib import pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

"""XGBoost で回帰するサンプルコード"""


def main():
    # Boston データセットを読み込む
    dataset = datasets.load_boston()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        )

    dtrain = xgb.DMatrix(X_train, label=y_train)
    dtest = xgb.DMatrix(X_test, label=y_test)

    xgb_params = {
        # 回帰問題
        'objective': 'reg:linear',
        # 学習用の指標 (RMSE)
        'eval_metric': 'rmse',
    }
    evals = [(dtrain, 'train'), (dtest, 'eval')]
    evals_result = {}
    bst = xgb.train(xgb_params,
                    dtrain,
                    num_boost_round=1000,
                    early_stopping_rounds=10,
                    evals=evals,
                    evals_result=evals_result,
                    )

    y_pred = bst.predict(dtest)
    mse = mean_squared_error(y_test, y_pred)
    print('RMSE:', math.sqrt(mse))

    train_metric = evals_result['train']['rmse']
    plt.plot(train_metric, label='train rmse')
    eval_metric = evals_result['eval']['rmse']
    plt.plot(eval_metric, label='eval rmse')
    plt.grid()
    plt.legend()
    plt.xlabel('rounds')
    plt.ylabel('rmse')
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python xgbreg.py    
[23:25:09] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 22 extra nodes, 0 pruned nodes, max_depth=5
[0]    train-rmse:17.5096 eval-rmse:16.1546
Multiple eval metrics have been passed: 'eval-rmse' will be used for early stopping.

Will train until eval-rmse hasn't improved in 10 rounds.
...(snip)...
Stopping. Best iteration:
[33]  train-rmse:0.344815 eval-rmse:3.03055

RMSE: 3.0332711348600068

学習過程の損失をプロットしたグラフは次のようになった。

カスタムメトリックを扱う

これまでは XGBoost に組み込みで入っていた評価指標を用いて学習の進み具合を評価していた。 続いては自分で書いたカスタムメトリックを使って学習してみる。

次のサンプルコードでは、再び乳がんデータセットを使った二値分類を扱っている。 ただし、学習を評価するメトリックとして精度 (Accuracy) を計測するカスタムメトリックを使っている。 カスタムメトリックを扱うにはオプションとして feval に自分で書いた関数を渡す。 カスタムメトリックの関数は、評価指標の名前と値をタプルまたはリストで返すように作る。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import xgboost as xgb

import numpy as np
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from matplotlib import pyplot as plt

"""XGBoost で学習時にカスタムメトリックを使ったサンプルコード"""


def feval_accuracy(pred_proba, dtrain):
    """カスタムメトリックを計算する関数"""
    # 真のデータ
    y_true = dtrain.get_label().astype(int)
    # 予測
    y_pred = np.where(pred_proba > 0.5, 1, 0)
    # Accuracy を計算する
    acc = accuracy_score(y_true, y_pred)
    # メトリックの名前と数値を返す(最小化を目指すので 1 から引く)
    return 'accuracy', 1 - acc


def main():
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)

    dtrain = xgb.DMatrix(X_train, label=y_train)
    dtest = xgb.DMatrix(X_test, label=y_test)

    xgb_params = {
        'objective': 'binary:logistic',
        'eval_metric': 'logloss',
    }

    evals = [(dtrain, 'train'), (dtest, 'eval')]
    evals_result = {}
    bst = xgb.train(xgb_params,
                    dtrain,
                    num_boost_round=1000,
                    early_stopping_rounds=10,
                    evals=evals,
                    evals_result=evals_result,
                    # カスタマイズした評価関数を使う
                    feval=feval_accuracy,
                    )

    y_pred_proba = bst.predict(dtest)
    y_pred = np.where(y_pred_proba > 0.5, 1, 0)
    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)

    # デフォルトのメトリック
    _, ax1 = plt.subplots(figsize=(8, 4))
    train_metric = evals_result['train']['logloss']
    ax1.plot(train_metric, label='train logloss', c='r')
    eval_metric = evals_result['eval']['logloss']
    ax1.plot(eval_metric, label='eval logloss', c='g')
    ax1.set_ylabel('logloss')
    ax1.legend()
    ax1.set_xlabel('rounds')

    # カスタムメトリック
    ax2 = ax1.twinx()
    eval_custom_metric = evals_result['eval']['accuracy']
    ax2.plot(eval_custom_metric, label='eval accuracy', c='b')
    ax2.set_ylabel('accuracy')
    ax2.legend()

    plt.grid()
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみる。

$ python xgbcustmetric.py
[23:35:00] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 18 extra nodes, 0 pruned nodes, max_depth=5
[0]    train-logloss:0.462396 eval-logloss:0.492899  train-accuracy:0.017588    eval-accuracy:0.070175
Multiple eval metrics have been passed: 'eval-accuracy' will be used for early stopping.

Will train until eval-accuracy hasn't improved in 10 rounds.
...(snip)...
Stopping. Best iteration:
[14]  train-logloss:0.029256  eval-logloss:0.112603   train-accuracy:0.005025 eval-accuracy:0.046784

Accuracy: 0.9532163742690059

得られた学習過程のグラフは次の通り。

学習の評価指標として精度 (Accuracy) を使うと、早々に学習が打ち切られてしまっているようだ。 そのため、最終的に得られた精度も低いものになってしまっている。 これは、おそらく一つの要素が正しく分類できたか・できなかったに結果が大きく左右されてしまうため。 まあ、結果はそんなに良くないけど、こんな感じで書けるよという例として。

組み込みの交差検証の機能を使ってみる

XGBoost には、組み込みの交差検証 (Cross Validation) の機能がある。 正直、そんなに使いやすいものではないけど、とりあえずあるよということで紹介しておく。

次のサンプルコードでは、XGBoost に組み込みで入っている交差検証の機能を使っている。 機能は xgboost.cv() という関数が起点になる。 基本的なパラメータについては xgboost.train() で普通に学習するときとさほど変わらない。 関数から最終的に得られるのは、学習過程におけるメトリックの値の変化。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import xgboost as xgb

import numpy as np
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from matplotlib import pyplot as plt

"""XGBoost 組み込みの交差検証を使ったサンプルコード"""

def main():
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target

    # CV の中で分割するので丸ごと渡す
    dtrain = xgb.DMatrix(X, label=y)

    xgb_params = {
        'objective': 'binary:logistic',
        'eval_metric': 'logloss',
    }

    # 交差検証する
    history = xgb.cv(xgb_params,
                     dtrain,
                     num_boost_round=1000,
                     early_stopping_rounds=10,
                     nfold=10,
                     # 層化分割する
                     stratified=True,
                     # 検証の経過を出力する
                     verbose_eval=True,
                     )

    train_metric = history['train-logloss-mean']
    plt.plot(train_metric, label='train logloss')
    eval_metric = history['test-logloss-mean']
    plt.plot(eval_metric, label='eval logloss')
    plt.grid()
    plt.legend()
    plt.xlabel('rounds')
    plt.ylabel('logloss')
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみる。

$ python xgbcustmetric.py
[23:35:00] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 18 extra nodes, 0 pruned nodes, max_depth=5
[0]    train-logloss:0.462396 eval-logloss:0.492899  train-accuracy:0.017588    eval-accuracy:0.070175
Multiple eval metrics have been passed: 'eval-accuracy' will be used for early stopping.

Will train until eval-accuracy hasn't improved in 10 rounds.
...(snip)...
[23:39:47] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 2 extra nodes, 0 pruned nodes, max_depth=1

すると、次のような学習過程のグラフが得られる。

いじょう。

まとめ

今回は XGBoost の機能を色々と試してみた。

Kaggleで勝つデータ分析の技術

統計的学習の基礎 ―データマイニング・推論・予測―

Linuxで動かしながら学ぶTCP/IPネットワーク入門

Python: Hyperopt で機械学習モデルのハイパーパラメータを選ぶ

今回は、機械学習モデルのハイパーパラメータをチューニングするのに用いられる Python のフレームワークの一つとして Hyperopt を使ってみる。 このフレームワークは、機械学習コンペティションの一つである Kaggle でよく用いられるものとして知られている。

なお、このブログでは過去にハイパーパラメータのチューニングについて Bayesian Optimization や Optuna を使った例を扱ったことがある。

blog.amedama.jp

blog.amedama.jp

今回使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.14.2
BuildVersion:   18C54
$ python -V
Python 3.7.2

下準備

まずは下準備として、必要になるパッケージをインストールする。

$ pip install hyperopt matplotlib scikit-learn pandas

単純な例で動作を確認する

機械学習モデルのハイパーパラメータを扱う前に、まずはもっと単純な例から始める。 手始めに、次のような一つの変数 x を取る関数の最大値を探す問題について考えてみよう。 グラフを見ると、最大値は 2 付近にあるようだ。

上記の関数の最大値を Hyperopt で探すサンプルコードは次の通り。 Hyperopt では、最適化したい目的関数を用意して、それを hyperopt.fmin() に渡す。 以下のサンプルコードにおいて目的関数は objective() という名前で定義している。 チューニングするパラメータは、どんな分布でどういった値域を持つのかあらかじめ定義しておく。 以下では hp.uniform('x', -5, +15) がそれで、値域が -5 ~ +15 の一様分布を定義している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import math

from hyperopt import hp
from hyperopt import fmin
from hyperopt import tpe


def f(x):
    """最大値を見つけたい関数"""
    return math.exp(-(x - 2) ** 2) + math.exp(-(x - 6) ** 2 / 10) + 1 / (x ** 2 + 1)


def objective(args):
    """最小化したい目的関数"""
    return -1 * f(args)  # 元のタスクが最大化なので符号を反転する


def main():
    # 変数の値域を定義する
    space = hp.uniform('x', -5, +15)
    # 目的関数を最小化するパラメータを探索する
    best = fmin(fn=objective, space=space, algo=tpe.suggest, max_evals=100)
    # 結果を出力する
    print(best)


if __name__ == '__main__':
    main()

上記に適当な名前をつけて実行する。 すると、たしかに解析的に確認していた最大値として 2 前後の値が得られる。

$ python hellohp.py
{'x': 2.0069551331288276}

探索過程を記録する

先ほどの例では、最終的に最適化された値が得られておしまいだった。 とはいえ、実際に使う場面では、その値がどういった経緯を経て得られたのか重要になる場合もある。 そこで、Hyperopt では Trials というオブジェクトに過程を記録できる。

次のサンプルコードでは Trials のインスタンスに探索過程を記録した上で、最後にそれを出力している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import math

from hyperopt import hp
from hyperopt import fmin
from hyperopt import tpe
from hyperopt import Trials
from hyperopt import STATUS_OK


def f(x):
    """最大値を見つけたい関数"""
    return math.exp(-(x - 2) ** 2) + math.exp(-(x - 6) ** 2 / 10) + 1 / (x ** 2 + 1)


def objective(args):
    """最小化したい目的関数"""
    # 計算過程を追跡するときは目的関数の返り値を辞書にする
    return {
        'loss': -1 * f(args),  # 最小化する損失 (必須パラメータ)
        'status': STATUS_OK,  # 計算結果 (必須パラメータ)
    }


def main():
    # 変数の値域を定義する
    space = hp.uniform('x', -5, +15)
    # 探索過程を記録するオブジェクト
    trials = Trials()
    # 目的関数を最小化するパラメータを探索する
    fmin(fn=objective, space=space, algo=tpe.suggest, max_evals=100, trials=trials)
    # 探索過程を出力する
    for item in trials.trials:
        vals = item['misc']['vals']
        result = item['result']
        print('vals:', vals, 'result:', result)


if __name__ == '__main__':
    main()

実行すると、どのような値を探索して結果がどうだったのかが確認できる。

$ python trials.py
vals: {'x': [10.496518222653432]} result: {'loss': -0.14140262193324998, 'status': 'ok'}
vals: {'x': [13.094879827038248]} result: {'loss': -0.012312365591358516, 'status': 'ok'}
vals: {'x': [-2.0228004674552014]} result: {'loss': -0.19799926655007422, 'status': 'ok'}
...(snip)...
vals: {'x': [2.6501743089141496]} result: {'loss': -1.1054773169025645, 'status': 'ok'}
vals: {'x': [9.247128914308604]} result: {'loss': -0.35996620277329294, 'status': 'ok'}
vals: {'x': [7.18849940727293]} result: {'loss': -0.8872540482414317, 'status': 'ok'}

探索過程に色々な値を記録する

先ほど使った Trials オブジェクトには、実は自分で色々な値を記録することもできる。 次のサンプルコードでは、探索した時刻と関数の引数を記録している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import time
import math

from hyperopt import hp
from hyperopt import fmin
from hyperopt import tpe
from hyperopt import Trials
from hyperopt import STATUS_OK


def f(x):
    """最大値を見つけたい関数"""
    return math.exp(-(x - 2) ** 2) + math.exp(-(x - 6) ** 2 / 10) + 1 / (x ** 2 + 1)


def objective(args):
    """最小化したい目的関数"""
    return {
        'loss': -1 * f(args),  # 最小化する損失 (必須パラメータ)
        'status': STATUS_OK,  # 計算結果 (必須パラメータ)
        # 一定の制約の元に自由に入れて良い
        'estimate_time': time.time(),
        'vals': args,
    }


def main():
    # 変数の値域を定義する
    space = hp.uniform('x', -5, +15)
    # 探索過程を記録するオブジェクト
    trials = Trials()
    # 目的関数を最小化するパラメータを探索する
    fmin(fn=objective, space=space, algo=tpe.suggest, max_evals=100, trials=trials)
    # 探索過程を出力する
    for item in trials.trials:
        result = item['result']
        print('result:', result)


if __name__ == '__main__':
    main()

上記を実行すると、探索した時刻や引数が得られていることが分かる。

$ python params.py 
result: {'loss': -0.339479707912701, 'status': 'ok', 'estimate_time': 1548600261.9176042, 'vals': -1.407873925928973}
result: {'loss': -0.682422018336043, 'status': 'ok', 'estimate_time': 1548600261.919095, 'vals': 8.012065634262315}
result: {'loss': -0.055263914137346264, 'status': 'ok', 'estimate_time': 1548600261.920187, 'vals': -4.135976583356424}
...(snip)...
result: {'loss': -1.1031356060488282, 'status': 'ok', 'estimate_time': 1548600262.10813, 'vals': 2.65344333175562}
result: {'loss': -0.27738481344650867, 'status': 'ok', 'estimate_time': 1548600262.111007, 'vals': -1.6261239199291966}
result: {'loss': -0.14645158143456388, 'status': 'ok', 'estimate_time': 1548600262.1134, 'vals': -2.4222183825698775}

探索過程を可視化してみる

先ほどの例では、どんな風に探索が進められているのかがいまいちイメージしにくかった。 そこで、どんな風に探索しているのかをグラフ上にプロットしてみた。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import time
import math

from matplotlib import pyplot as plt

import numpy as np

from hyperopt import hp
from hyperopt import fmin
from hyperopt import tpe
from hyperopt import Trials
from hyperopt import STATUS_OK


def f(x):
    """最大値を見つけたい関数"""
    return math.exp(-(x - 2) ** 2) + math.exp(-(x - 6) ** 2 / 10) + 1 / (x ** 2 + 1)


def objective(args):
    """最小化したい目的関数"""
    return {
        'loss': -1 * f(args),  # 最小化する損失 (必須パラメータ)
        'status': STATUS_OK,  # 計算結果 (必須パラメータ)
        'vals': args,
    }


def main():
    # 変数の値域を定義する
    space = hp.uniform('x', -5, +15)
    # 探索過程を記録するオブジェクト
    trials = Trials()
    # 目的関数を最小化するパラメータを探索する
    fmin(fn=objective, space=space, algo=tpe.suggest, max_evals=100, trials=trials)

    # 検索過程をグラフで可視化してみる
    args = [result['vals'] for result in trials.results]
    values = [-result['loss'] for result in trials.results]
    # 元のグラフをプロットする
    X = [x for x in np.arange(-5, 15, 0.1)]
    y = [f(x) for x in X]
    plt.plot(X, y)
    # 探索した場所をプロットする
    plt.scatter(args, values)
    # グラフを描画する
    plt.xlabel('x')
    plt.ylabel('y')
    plt.grid()
    plt.show()


if __name__ == '__main__':
    main()

上記に適当な名前をつけて実行する。

$ python visualize.py

すると、次のようなグラフが得られる。

頂点付近以外にも、かなり色んな所を調べていることが伺える。

機械学習のモデルに適用してみる

続いては、ついに Hyperopt を機械学習のモデルに適用してみる。 モデルのアルゴリズムにはサポートベクターマシンを使った。 パラメータとしては kernel と C と gamma をチューニングする。

次のサンプルコードでは、目的関数でサポートベクターマシンを学習させて 5-Fold CV で Accuracy を評価している。 kernel は ['linear', 'rbf', 'poly'] の三種類から選ぶ。 C と gamma は常用対数で 0 ~ 100 と -10 ~ 1 の範囲で探索する。 2.303 を乗じているのは二進対数と常用対数に変換するため。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from functools import partial

from hyperopt import hp
from hyperopt import fmin
from hyperopt import tpe
from hyperopt import Trials
from hyperopt import space_eval

from sklearn import datasets
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_validate

from sklearn.svm import SVC


def objective(X, y, args):
    """最小化したい目的関数"""
    # モデルのアルゴリズムに SVM を使う
    model = SVC(**args)
    # Stratified 5 Fold Cross Validation
    kf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    scores = cross_validate(model, X=X, y=y, cv=kf)
    # 最小化なので符号を反転する
    return -1 * scores['test_score'].mean()


def main():
    # データセットを読み込む
    dataset = datasets.load_iris()
    X, y = dataset.data, dataset.target
    # 目的関数にデータセットを部分適用する
    f = partial(objective, X, y)
    # 変数の値域を定義する
    space = {
        # 底が自然対数なので 2.303 をかけて常用対数に変換する
        'C': hp.loguniform('C', 2.303 * 0, 2.303 * +2),
        'gamma': hp.loguniform('gamma', 2.303 * -2, 2.303 * +1),
        'kernel': hp.choice('kernel', ['linear', 'rbf', 'poly']),
    }
    # 探索過程を記録するオブジェクト
    trials = Trials()
    # 目的関数を最小化するパラメータを探索する
    best = fmin(fn=f, space=space, algo=tpe.suggest, max_evals=100, trials=trials)
    # 結果を出力する
    print(space_eval(space, best))


if __name__ == '__main__':
    main()

実行すると、次のような結果が得られる。 値は似通ったものになるはずだけど、結果は毎回異なるはず。

$ python svm.py 
{'C': 1.9877559946331484, 'gamma': 0.7522911704230325, 'kernel': 'linear'}

探索を可視化してみる

さきほどと同じように、どのような値を探索したのかグラフにプロットしてみよう。 次のサンプルコードでは、カーネルを linear に決め打ちした状態で、Hyperopt が探索した C と gamma を二次元でプロットしている。 標準化した Accuracy の高低は点の色で表現している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from functools import partial

from hyperopt import hp
from hyperopt import fmin
from hyperopt import tpe
from hyperopt import Trials
from hyperopt import STATUS_OK
from hyperopt import space_eval

from sklearn import datasets
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_validate
from sklearn.svm import SVC

import numpy as np

import pandas as pd

from matplotlib import pyplot as plt
from matplotlib import cm


def objective(X, y, args):
    """最小化したい目的関数"""
    # モデルのアルゴリズムに Linear SVM を使う
    model = SVC(kernel='linear', **args)
    # Stratified 5 Fold Cross Validation
    kf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    scores = cross_validate(model, X=X, y=y, cv=kf)
    # 最小化なので符号を反転する
    return {
        'loss': -1 * scores['test_score'].mean(),
        'status': STATUS_OK,
        'vals': args,
    }


def main():
    # データセットを読み込む
    dataset = datasets.load_iris()
    X, y = dataset.data, dataset.target
    # 目的関数にデータセットを部分適用する
    f = partial(objective, X, y)
    # 変数の値域を定義する
    space = {
        'C': hp.loguniform('C', 2.303 * 0, 2.303 * +2),
        'gamma': hp.loguniform('gamma', 2.303 * -2, 2.303 * +1),
    }
    # 探索過程を記録するオブジェクト
    trials = Trials()
    # 目的関数を最小化するパラメータを探索する
    best = fmin(fn=f, space=space, algo=tpe.suggest, max_evals=100, trials=trials)
    # 結果を出力する
    print('best parameters:', space_eval(space, best))

    # 探索過程を可視化する
    xs = [result['vals']['C'] for result in trials.results]
    ys = [result['vals']['gamma'] for result in trials.results]
    zs = np.array([-1 * result['loss'] for result in trials.results])
    s_zs = (zs - zs.min()) / (zs.max() - zs.min()) # 0 ~ 1 の範囲に正規化する
    # プロットする
    sc = plt.scatter(xs, ys, c=s_zs, s=20, zorder=10, cmap=cm.cool)
    plt.colorbar(sc)

    plt.xlabel('C')
    plt.xscale('log')
    plt.ylabel('gamma')
    plt.yscale('log')
    plt.grid()
    plt.show()

    print('-- accuracy summary --')
    print(pd.Series(np.array(zs).ravel()).describe())


if __name__ == '__main__':
    main()

上記を実行する。

$ python svmplot.py
best parameters: {'C': 2.3141971486483675, 'gamma': 0.014160000992695598}
-- accuracy summary --
count    100.000000
mean       0.975800
std        0.008350
min        0.960000
25%        0.966667
50%        0.980000
75%        0.980000
max        0.986667
dtype: float64

すると、次のようなグラフが得られる。

たしかに、場所によって性能に違いが出ているようだ。 見た感じ gamma よりも C の方が性能に違いを与えていそう。

複数のモデルから選ぶ

先ほどの例では、サポートベクターマシンの中でハイパーパラメータをチューニングしていた。 これだと、別のモデルを含めて評価したいときは別々に実行する必要がある。 そんなとき Hyperopt ではまとめて評価することもできる。

次のサンプルコードではランダムフォレスト、サポートベクターマシン、ロジスティック回帰を一度に評価している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from functools import partial

from hyperopt import hp
from hyperopt import fmin
from hyperopt import tpe
from hyperopt import Trials
from hyperopt import space_eval
from hyperopt.pyll.base import scope

from sklearn import datasets
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_validate

from sklearn.svm import SVC


def objective(X, y, args):
    """最小化したい目的関数"""
    classifiers = {
        'svm': SVC,
        'rf': RandomForestClassifier,
        'logit': LogisticRegression,
    }
    classifier = classifiers.get(args['model_type'])
    del args['model_type']
    model = classifier(**args)
    # Stratified 5 Fold Cross Validation
    kf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    scores = cross_validate(model, X=X, y=y, cv=kf)
    # 最小化なので符号を反転する
    return -1 * scores['test_score'].mean()


def main():
    # データセットを読み込む
    dataset = datasets.load_iris()
    X, y = dataset.data, dataset.target
    # 目的関数にデータセットを部分適用する
    f = partial(objective, X, y)
    # 変数の値域を定義する
    space = hp.choice('algorithms', [
        {
            'model_type': 'rf',
            'n_estimators': scope.int(hp.uniform('n_estimators', 1e+1, 1e+3)),
            'max_depth': scope.int(hp.uniform('max_depth', 1e+1, 1e+3)),
        },
        {
            'model_type': 'svm',
            'C': hp.uniform('C', 1e+0, 1e+2),
            'gamma': hp.lognormal('gamma', 1e-2, 1e+1),
        },
        {
            'model_type': 'logit',
            'solver': 'lbfgs',
            'multi_class': 'auto',
            'max_iter': 1000,
        }
    ])
    # 探索過程を記録するオブジェクト
    trials = Trials()
    # 目的関数を最小化するパラメータを探索する
    best = fmin(fn=f, space=space, algo=tpe.suggest, max_evals=100, trials=trials)
    # 結果を出力する
    print(space_eval(space, best))


if __name__ == '__main__':
    main()

上記を実行してみよう。 このパラメータの中ではサポートベクターマシンが最も良い結果が得られたようだ。

$ python multimodel.py
{'C': 98.52420177369001, 'gamma': 0.011789334572006299, 'model_type': 'svm'}

いじょう。

スマートPythonプログラミング: Pythonのより良い書き方を学ぶ