CUBE SUGAR CONTAINER

技術系のこと書きます。

Python: IsolationForest で教師なし学習の外れ値検知を試す

今回は教師なし学習で外れ値の検知に使える IsolationForest というアルゴリズムを試してみる。 このアルゴリズムの興味深いところは、教師データの中にある程度外れ値が含まれていても構わないという点。 つまり、アノテーションしていないデータをそのまま突っ込むことが許容されている。

IsolationForest のアルゴリズムでは、決定木を使った分類しやすさにもとづいてデータが正常か外れ値かを判断する。 外れ値は正常なデータに比べると数が少なく、特徴が大きく異なると仮定する。 だとすると、外れ値は正常なデータに比べて分類するのに木の深さがより多く必要と考える。

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.14.4
BuildVersion:   18E226
$ python -V            
Python 3.7.3

下準備

まずは scikit-learn と matplotlib をインストールしておく。

$ pip install scikit-learn matplotlib

ダミデータを用意する

IsolationForest を試すのに、以下のようなダミーデータを生成した。 特徴量は二次元で、ひとまず Blob でクラスタが 1 つだけある正常値の周囲に一様分布で外れ値を散らしている。

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

import numpy as np
from sklearn.datasets import make_blobs
from matplotlib import pyplot as plt


def main():
    # 正常値のデータを作る
    NORMAL_DATA_N = 1000
    args = {
        # 正常値のデータ点数
        'n_samples': NORMAL_DATA_N,
        # 特徴量の次元数
        'n_features': 2,
        # クラスタの数
        'centers': 1,
        # データを生成するための乱数
        'random_state': 42,
    }
    X_normal, y_normal = make_blobs(**args)

    # 外れ値のデータを作る
    OUTLIER_DATA_N = 200
    DISTRIBUTION_RANGE = 6
    X_outlier = np.random.uniform(low=-DISTRIBUTION_RANGE,
                                  high=DISTRIBUTION_RANGE,
                                  size=(OUTLIER_DATA_N, 2))
    y_outlier = np.ones((OUTLIER_DATA_N, ))

    # 正常値を中心に外れ値を分布させる
    X_outlier += X_normal.mean(axis=0)

    # くっつけて一つのデータセットにする
    X = np.concatenate([X_normal, X_outlier], axis=0)
    y = np.concatenate([y_normal, y_outlier], axis=0)

    plt.scatter(X[y == 0, 0],
                X[y == 0, 1],
                label='negative')
    plt.scatter(X[y == 1, 0],
                X[y == 1, 1],
                label='positive')
    plt.legend()
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみる。

$ python blob.py

すると、以下のような散布図が得られる。

f:id:momijiame:20190420122548p:plain

正常値の領域にも外れ値が含まれてしまっているけど、まあ完全に分離可能なデータなんてそうそうないので。

IsolationForest を適用してみる

それでは、先ほどのデータに IsolationForest を適用してみよう。

以下のサンプルコードでは、ホールドアウト検証でどのように分類されるか可視化すると共に評価指標を計算している。 IsolationForest のインスタンスで fit() メソッドを呼ぶときに目的変数 (y_train) に関する情報を渡していないのがポイント。 また、渡している説明変数 (X_train) の中には約 20% の外れ値が含まれている。

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

import numpy as np
from scipy import stats
from sklearn.datasets import make_blobs
from sklearn.ensemble import IsolationForest
from sklearn.model_selection import train_test_split
from sklearn.metrics import precision_score
from sklearn.metrics import recall_score
from sklearn.metrics import f1_score
from sklearn.metrics import matthews_corrcoef
from matplotlib import pyplot as plt


def main():
    # 正常値のデータを作る
    NORMAL_DATA_N = 1000
    args = {
        # 正常値のデータ点数
        'n_samples': NORMAL_DATA_N,
        # 特徴量の次元数
        'n_features': 2,
        # クラスタの数
        'centers': 1,
        # データを生成するための乱数
        'random_state': 42,
    }
    X_normal, y_normal = make_blobs(**args)

    # 外れ値のデータを作る
    OUTLIER_DATA_N = 200
    DIST_RANGE = 6
    X_outlier = np.random.uniform(low=-DIST_RANGE,
                                  high=DIST_RANGE,
                                  size=(OUTLIER_DATA_N, 2))
    y_outlier = np.ones((OUTLIER_DATA_N, ))

    # 正常値を中心に外れ値を分布させる
    X_outlier += X_normal.mean(axis=0)

    # くっつけて一つのデータセットにする
    X = np.concatenate([X_normal, X_outlier], axis=0)
    y = np.concatenate([y_normal, y_outlier], axis=0)

    # ホールドアウト検証するためにデータを分割する
    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)

    # IsolationForest
    clf = IsolationForest(n_estimators=100,
                          contamination='auto',
                          behaviour='new',
                          random_state=42)

    # 学習用データを学習させる
    clf.fit(X_train)

    # 検証用データを分類する
    y_pred = clf.predict(X_test)

    # IsolationForest は 正常=1 異常=-1 として結果を返す
    y_pred[y_pred == 1] = 0
    y_pred[y_pred == -1] = 1

    # 評価指標を計算する
    print('precision:', precision_score(y_test, y_pred))
    print('recall:', recall_score(y_test, y_pred))
    print('f1:', f1_score(y_test, y_pred))
    print('mcc:', matthews_corrcoef(y_test, y_pred))

    xx, yy = np.meshgrid(np.linspace(X_normal[:, 0].mean() - DIST_RANGE * 1.2,
                                     X_normal[:, 0].mean() + DIST_RANGE * 1.2,
                                     100),
                         np.linspace(X_normal[:, 1].mean() - DIST_RANGE * 1.2,
                                     X_normal[:, 1].mean() + DIST_RANGE * 1.2,
                                     100))
    Z = clf.decision_function(np.c_[xx.ravel(), yy.ravel()])
    Z = Z.reshape(xx.shape)

    plt.contourf(xx, yy, Z, cmap=plt.cm.Blues_r)
    threshold = stats.scoreatpercentile(y_pred,
                                        100 * (OUTLIER_DATA_N / NORMAL_DATA_N))
    plt.contour(xx, yy, Z, levels=[threshold],
                    linewidths=2, colors='red')

    plt.scatter(X_train[y_train == 0, 0],
                X_train[y_train == 0, 1],
                label='train negative',
                c='b')
    plt.scatter(X_train[y_train == 1, 0],
                X_train[y_train == 1, 1],
                label='train positive',
                c='r')

    plt.scatter(X_test[y_pred == 0, 0],
                X_test[y_pred == 0, 1],
                label='test negative (prediction)',
                c='g')
    plt.scatter(X_test[y_pred == 1, 0],
                X_test[y_pred == 1, 1],
                label='test positive (prediction)',
                c='y')

    plt.legend()
    plt.show()


if __name__ == '__main__':
    main()

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

$ python isoforest.py 
precision: 0.9215686274509803
recall: 0.7833333333333333
f1: 0.8468468468468469
mcc: 0.8229295359450786

同時に、以下のグラフが得られる。 正常値のクラスタの周囲に識別境界ができて、外れ値検出ができていることが分かる。

f:id:momijiame:20190420122902p:plain

正常値のクラスタを増やしてみる

正常値のクラスタが一つのときは、まあなんとなく上手くいきそうなことが分かった。 続いては二つに増やしてみることにしよう。

以下のサンプルコードでは、先ほどの内容を少し修正してダミーデータのクラスタを二つに増やしている。

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

import numpy as np
from scipy import stats
from sklearn.datasets import make_blobs
from sklearn.ensemble import IsolationForest
from sklearn.model_selection import train_test_split
from sklearn.metrics import precision_score
from sklearn.metrics import recall_score
from sklearn.metrics import f1_score
from sklearn.metrics import matthews_corrcoef
from matplotlib import pyplot as plt


def main():
    NORMAL_DATA_N = 1000
    args = {
        'n_samples': NORMAL_DATA_N,
        'n_features': 2,
        # クラスタの数を増やしてみる
        'centers': 2,
        'random_state': 42,
    }
    X_normal, _ = make_blobs(**args)
    y_normal = np.zeros((NORMAL_DATA_N, ))

    OUTLIER_DATA_N = 200
    DIST_RANGE = 6
    X_outlier = np.random.uniform(low=-DIST_RANGE,
                                  high=DIST_RANGE,
                                  size=(OUTLIER_DATA_N, 2))
    y_outlier = np.ones((OUTLIER_DATA_N, ))

    X_outlier += X_normal.mean(axis=0)

    X = np.concatenate([X_normal, X_outlier], axis=0)
    y = np.concatenate([y_normal, y_outlier], axis=0)

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

    clf = IsolationForest(n_estimators=100,
                          contamination='auto',
                          behaviour='new',
                          random_state=42)

    clf.fit(X_train)

    y_pred = clf.predict(X_test)

    y_pred[y_pred == 1] = 0
    y_pred[y_pred == -1] = 1

    print('precision:', precision_score(y_test, y_pred))
    print('recall:', recall_score(y_test, y_pred))
    print('f1:', f1_score(y_test, y_pred))
    print('mcc:', matthews_corrcoef(y_test, y_pred))

    xx, yy = np.meshgrid(np.linspace(X_normal[:, 0].mean() - DIST_RANGE * 1.2,
                                     X_normal[:, 0].mean() + DIST_RANGE * 1.2,
                                     100),
                         np.linspace(X_normal[:, 1].mean() - DIST_RANGE * 1.2,
                                     X_normal[:, 1].mean() + DIST_RANGE * 1.2,
                                     100))
    Z = clf.decision_function(np.c_[xx.ravel(), yy.ravel()])
    Z = Z.reshape(xx.shape)

    plt.contourf(xx, yy, Z, cmap=plt.cm.Blues_r)
    threshold = stats.scoreatpercentile(y_pred,
                                        100 * (OUTLIER_DATA_N / NORMAL_DATA_N))
    plt.contour(xx, yy, Z, levels=[threshold],
                    linewidths=2, colors='red')

    plt.scatter(X_train[y_train == 0, 0],
                X_train[y_train == 0, 1],
                label='train negative',
                c='b')
    plt.scatter(X_train[y_train == 1, 0],
                X_train[y_train == 1, 1],
                label='train positive',
                c='r')

    plt.scatter(X_test[y_pred == 0, 0],
                X_test[y_pred == 0, 1],
                label='test negative (prediction)',
                c='g')
    plt.scatter(X_test[y_pred == 1, 0],
                X_test[y_pred == 1, 1],
                label='test positive (prediction)',
                c='y')

    plt.legend()
    plt.show()


if __name__ == '__main__':
    main()

実行結果は次の通り。

$ python twoclusters.py 
precision: 0.5238095238095238
recall: 0.7333333333333333
f1: 0.611111111111111
mcc: 0.528680532637681

また、以下のグラフが得られる。

f:id:momijiame:20190420123453p:plain

評価指標は悪化しているものの、可視化したものを見るとちゃんと二つのクラスタの周囲に分離境界ができていることがわかる。

いじょう。

Python: scikit-learn の cross_val_predict() 関数で OOF な特徴量を作る

scikit-learn には cross_val_predict() という関数がある。 この関数は、教師データを k-Fold などで分割したときに OOF (Out of Fold) なデータの目的変数を推論する目的で使われることが多い。 なお、OOF (Out of Fold) というのは、k-Fold などでデータを分割した際に学習に使わなかったデータを指している。

scikit-learn.org

ちなみに、この関数は method オプションを指定すると fit() メソッドのあとに呼び出すメソッド名を自由に指定できる。 典型的には、分類問題において predict_proba を指定することで推論結果を確率として取得するために使うことが多い。

ただ、この method オプションには predict 系のメソッドに限らず、別に何を指定しても良いことに気づいた。 そこで、今回の記事では cross_val_predict() 関数を使って OOF な特徴量を作ってみることにした。 というのも、Target Encoding と呼ばれる目的変数を利用した特徴量抽出においては OOF することが Leakage を防ぐ上で必須になる。 もし、その作業を cross_val_predict() 関数と scikit-learn の Transformer を組み合わせることでキレイに書けるとしたら、なかなか便利なのではと思った。

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

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.14.4
BuildVersion:   18E226
$ python -V    
Python 3.7.3

下準備

まずは pandas と scikit-learn をインストールしておく。

$ pip install pandas scikit-learn

サンプルコード

以下に cross_val_predict() 関数と Transformer を組み合わせて Target Encoding した特徴量を抽出するサンプルコードを示す。 TargetMeanEncodingTransformer という名前で Target Encoding する Transformer クラスを定義している。 これをそのまま使ってしまうと Target Leakage が起こるため cross_val_predict() 関数と組み合わせることで OOF な特徴量を生成している。 ちなみに cross_val_predict() 関数の返り値は numpy の配列になるため、pandas と組み合わせて使うときはその点に考慮が必要になる。

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

import pandas as pd
from sklearn.base import BaseEstimator
from sklearn.base import TransformerMixin
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_val_predict


class TargetMeanEncodingTransformer(BaseEstimator, TransformerMixin):
    """Target mean encoding に使う scikit-learn の Transformer

    NOTE: Target leakage を防ぐために必ず OOF で特徴量を作ること"""

    def __init__(self, category, target=None):
        """
        :param category: エンコード対象のカラム名
        :param target: 目的変数のカラム名
        """
        self.category_ = category
        self.target_ = target
        # 学習したデータを保存する場所
        self.target_mean_ = None

    def fit(self, X, y=None):
        # エンコード対象の特徴量をカテゴリごとにグループ化する
        gdf = X.groupby(self.category_)
        # ターゲットの比率を計算する
        self.target_mean_ = gdf[self.target_].mean()
        # 自身を返す
        return self

    def transform(self, X, copy=None):
        # エンコード対象のカラムを取り出す
        target_series = X[self.category_]
        # エンコードする
        transformed = target_series.map(self.target_mean_.to_dict())
        # 特徴量を返す
        return transformed

    def get_params(self, deep=True):
        # NOTE: cross_val_predict() は Fold ごとに Estimator を作り直す
        #       そのため get_params() でインスタンスのパラメータを得る必要がある
        return {
            'category': self.category_,
            'target': self.target_,
        }


def main():
    # 元データ
    data = [
        ('Apple', 1),
        ('Apple', 1),
        ('Apple', 0),
        ('Banana', 1),
        ('Banana', 0),
        ('Banana', 0),
        ('Cherry', 0),
        ('Cherry', 0),
        ('Cherry', 0),
    ]
    columns = ['category', 'delicious']
    df = pd.DataFrame(data, columns=columns)

    transformer = TargetMeanEncodingTransformer(category='category',
                                                target='delicious')

    # k-Fold しないで適用した場合 (Target leak するのでやっちゃダメ!)
    encoded_without_kfold = transformer.fit_transform(df)
    encoded_without_kfold.name = 'target_mean_encoded'
    print('Without k-Fold')
    concat_df = pd.concat([df, encoded_without_kfold], axis=1, sort=False)
    print(concat_df)

    # k-Fold して適用した場合 (データ点数の関係で実質 LOO になってるので本来はこれもすべきではない)
    skf = StratifiedKFold(n_splits=3,
                          shuffle=True,
                          random_state=42)
    # NOTE: numpy の array として結果が返る点に注意する
    encoded_with_kfold = cross_val_predict(transformer,
                                           # NOTE: エンコード対象のカテゴリで層化抽出する
                                           df, df.category,
                                           cv=skf,
                                           # fit() した後に transform() を呼ぶ
                                           method='transform')
    print('With k-Fold')
    concat_df = pd.concat([df,
                           pd.Series(encoded_with_kfold,
                                     name='target_mean_encoded')
                           ], axis=1, sort=False)
    print(concat_df)


if __name__ == '__main__':
    main()

上記の実行結果は次の通り。 k-Fold ありのパターンでは、ちゃんと OOF で Target Encoding した特徴量が作れているようだ。

$ python oof.py 
Without k-Fold
  category  delicious  target_mean_encoded
0    Apple          1             0.666667
1    Apple          1             0.666667
2    Apple          0             0.666667
3   Banana          1             0.333333
4   Banana          0             0.333333
5   Banana          0             0.333333
6   Cherry          0             0.000000
7   Cherry          0             0.000000
8   Cherry          0             0.000000
With k-Fold
  category  delicious  target_mean_encoded
0    Apple          1                  0.5
1    Apple          1                  0.5
2    Apple          0                  1.0
3   Banana          1                  0.0
4   Banana          0                  0.5
5   Banana          0                  0.5
6   Cherry          0                  0.0
7   Cherry          0                  0.0
8   Cherry          0                  0.0

いじょう。

Python: LightGBM でカスタムメトリックを扱う

今回は LightGBM で、組み込みで用意されていない独自の評価指標 (カスタムメトリック) を扱う方法について。 ユースケースとしては、学習自体は別の評価指標を使って進めつつ、本来の目標としている評価指標を同時に確認するといったもの。 例えば、精度 (Accuracy) やマシューズ相関係数 (Matthews Correlation Coefficient) は、学習にそのまま用いることは難しい。 しかしながら、最終的な目標としている評価指標がそれらになっていることはよくある。

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.14.4
BuildVersion:   18E226
$ python -V                  
Python 3.7.3
$ pip list | grep -i lightgbm
lightgbm        2.2.3  

下準備

$ brew install libomp
$ pip install lightgbm scikit-learn matplotlib

独自の評価指標を用いる

以下のサンプルコードでは、学習には LogLoss を使いつつ、同時に Accuracy を計算している。 独自の評価指標を計算するときは train() 関数や cv() 関数で feval というオプションを用いる。 指定するのは評価指標を計算する関数で、引数はモデルが予測した値と学習に使ったデータの二つ。 データからは get_label() という関数で真のラベルを取得できる。 注意点として、モデルが予測した値は多値分類問題であっても一次元の配列になっているため reshape する必要がある。 評価指標を計算する関数では、返り値として評価指標の名前、スコア、そしてスコアが大きい方が優れているのか否かを表す真偽値を返す。

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

import lightgbm as lgb
from sklearn import datasets
from matplotlib import pyplot as plt
import numpy as np
from sklearn.model_selection import train_test_split


def accuracy(preds, data):
    """精度 (Accuracy) を計算する関数"""
    # 正解ラベル
    y_true = data.get_label()
    # 推論の結果が 1 次元の配列になっているので直す
    N_LABELS = 3  # ラベルの数
    reshaped_preds = preds.reshape(N_LABELS, len(preds) // N_LABELS)
    # 最尤と判断したクラスを選ぶ 
    y_pred = np.argmax(reshaped_preds, axis=0)
    # メトリックを計算する
    acc = np.mean(y_true == y_pred)
    # name, result, is_higher_better
    return 'accuracy', acc, True


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

    # 学習用データと検証用データに分割する
    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        shuffle=True,
                                                        random_state=42)

    # LightGBM のデータセットの表現に直す
    lgb_train = lgb.Dataset(X_train, y_train)
    lgb_eval = lgb.Dataset(X_test, y_test, reference=lgb_train)

    lgbm_params = {
        'objective': 'multiclass',
        'num_class': 3,
    }

    evals_result = {}
    lgb.train(lgbm_params,
              lgb_train,
              # メトリックを追跡する対象のデータセット
              valid_sets=[lgb_eval, lgb_train],
              # 上記の名前
              valid_names=['eval', 'train'],
              num_boost_round=1000,
              # メトリックの履歴を残すオブジェクト
              evals_result=evals_result,
              # 独自メトリックを計算する関数
              feval=accuracy,
              )

    # 組み込みのメトリック
    eval_metric_logloss = evals_result['eval']['multi_logloss']
    train_metric_logloss = evals_result['train']['multi_logloss']

    # カスタムメトリック
    eval_metric_acc = evals_result['eval']['accuracy']
    train_metric_acc = evals_result['train']['accuracy']

    # グラフにプロットする
    _, ax1 = plt.subplots(figsize=(8, 4))
    ax1.plot(eval_metric_logloss, label='eval logloss', c='r')
    ax1.plot(train_metric_logloss, label='train logloss', c='b')
    ax1.set_ylabel('logloss')
    ax1.set_xlabel('rounds')
    ax1.legend()

    ax2 = ax1.twinx()
    ax2.plot(eval_metric_acc, label='eval accuracy', c='g')
    ax2.plot(train_metric_acc, label='train accuracy', c='y')
    ax2.set_ylabel('accuracy')
    ax2.legend()

    plt.grid()
    plt.show()


if __name__ == '__main__':
    main()

サンプルコードでは、評価指標の推移を折れ線グラフとしてプロットしている。

上記を実行してみよう。 ログを確認すると、ちゃんと Accuracy に関する情報も出力される。

$ python lgbcm.py
...
[1000] train's multi_logloss: 8.23321e-05    train's accuracy: 1 eval's multi_logloss: 0.396726  eval's accuracy: 0.947368

上記を実行すると、以下のようなグラフが得られる。

f:id:momijiame:20190331004325p:plain

上記のグラフを見ると、学習に使った LogLoss に加えて Accuracy の推移も確認できる。 どうやら、イテレーション数が 100 を越えないあたりから検証データに対する LogLoss が増加に転じているようだ。 検証データに対する LogLoss は増加することなく現象し続けており、過学習を起こしていることがわかる。 検証データの LogLoss が増加するタイミングでは、学習データに対する Accuracy が増加すると共に検証データの Accuracy も低下している。

複数のカスタムメトリックを計測したい場合

先ほどの例ではカスタムメトリックが一つだったけど、複数確認したい場合もあるはず。 そのような場合には複数のメトリックをリストとして返す関数を作ることで対応できる。

以下のサンプルコードでは Breast Cancer データセットを用いて Accuracy, Precision, Recall をカスタムメトリックとして計測している。

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

import numpy as np
from sklearn import datasets
from sklearn.metrics import precision_score
from sklearn.metrics import recall_score
from sklearn.model_selection import train_test_split
from matplotlib import pyplot as plt
import lightgbm as lgb


def accuracy(preds, data):
    """精度 (Accuracy) を計算する関数"""
    y_true = data.get_label()
    y_pred = np.where(preds > 0.5, 1, 0)
    metric = np.mean(y_true == y_pred)
    return 'accuracy', metric, True


def precision(preds, data):
    """適合率 (Precision) を計算する関数"""
    y_true = data.get_label()
    y_pred = np.where(preds > 0.5, 1, 0)
    metric = precision_score(y_true, y_pred)
    return 'precision', metric, True


def recall(preds, data):
    """再現率 (Recall) を計算する関数"""
    y_true = data.get_label()
    y_pred = np.where(preds > 0.5, 1, 0)
    metric = recall_score(y_true, y_pred)
    return 'recall', metric, True


def metrics(preds, data):
    """複数の評価指標を計算するための関数"""
    # リストでまとめて返せば良い
    return [
        accuracy(preds, data),
        precision(preds, data),
        recall(preds, data),
    ]


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

    # 学習用データと検証用データに分割する
    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        shuffle=True,
                                                        random_state=42)

    # LightGBM のデータセットの表現に直す
    lgb_train = lgb.Dataset(X_train, y_train)
    lgb_eval = lgb.Dataset(X_test, y_test, reference=lgb_train)

    lgbm_params = {
        'objective': 'binary',
        'metric': 'binary_logloss',
    }

    evals_result = {}
    lgb.train(lgbm_params,
              lgb_train,
              num_boost_round=1000,
              early_stopping_rounds=100,
              # メトリックを追跡する対象のデータセット
              valid_sets=[lgb_eval, lgb_train],
              # 上記の名前
              valid_names=['eval', 'train'],
              # メトリックの履歴を残すオブジェクト
              evals_result=evals_result,
              # 独自メトリックを計算する関数
              feval=metrics,
              )

    # 組み込みのメトリック
    eval_metric_logloss = evals_result['eval']['binary_logloss']
    train_metric_logloss = evals_result['train']['binary_logloss']

    # カスタムメトリック (Accuracy)
    eval_metric_acc = evals_result['eval']['accuracy']
    train_metric_acc = evals_result['train']['accuracy']

    # カスタムメトリック (Precision)
    eval_metric_prec = evals_result['eval']['precision']
    train_metric_prec = evals_result['train']['precision']

    # カスタムメトリック (Precision)
    eval_metric_recall = evals_result['eval']['recall']
    train_metric_recall = evals_result['train']['recall']

    # グラフにプロットする
    _, ax1 = plt.subplots(figsize=(8, 4))
    ax1.plot(eval_metric_logloss, label='eval logloss')
    ax1.plot(train_metric_logloss, label='train logloss')
    ax1.set_ylabel('logloss')
    ax1.set_xlabel('rounds')
    ax1.legend()

    ax2 = ax1.twinx()
    ax2.plot(eval_metric_acc, label='eval accuracy')
    ax2.plot(train_metric_acc, label='train accuracy')
    ax2.plot(eval_metric_prec, label='eval precision')
    ax2.plot(train_metric_prec, label='train precision')
    ax2.plot(eval_metric_recall, label='eval recall')
    ax2.plot(train_metric_recall, label='train recall')
    ax2.set_ylabel('ratio')
    ax2.legend()

    plt.grid()
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python lgbmcm.py
...
[1]    train's binary_logloss: 0.58296   train's accuracy: 0.629108  train's precision: 0.629108   train's recall: 1   eval's binary_logloss: 0.585825 eval's accuracy: 0.622378   eval's precision: 0.622378  eval's recall: 1

以下のようなグラフが得られる。 たしかに、複数のメトリックが計測できていることが分かる。

f:id:momijiame:20190613060743p:plain

注意点

今回のように複数の評価指標を LightGBM に計算させるときは early_stopping_rounds との併用において注意が必要になる。 というのも early_stopping_rounds で early stopping の対象になるのが「いずれかの評価指標で条件を満たした場合」になっているため。 つまり、本来は LogLoss の値が増加に転じた場合に止めたいのに、まだ現象を続けている段階で別の評価指標が条件を満たすと学習が止まってしまう恐れがある。

Python: LightGBM の cv() 関数から学習済みモデルを得る

(2020-09-05 追記)

LightGBM v3.0.0 から cv() 関数に return_cvbooster オプションが追加されました。 これにより直接 CVBooster のインスタンスが取得できるため、下記のコールバックを使う必要はなくなりました。


勾配ブースティング決定木を扱うフレームワークの一つである LightGBM の Python API には cv() という関数がある。 この "cv" というのは Cross Validation の略で、その名の通り LightGBM のモデルを交差検証するための関数になっている。 具体的には、この関数にデータセットを渡すと、そのデータでモデルを学習させると共に、指定した評価指標について交差検証で評価できる。

今回は、この関数から交差検証の過程で学習させたモデルを手に入れる方法について書いてみる。 というのも、この関数が返すのは指定した評価指標を用いて計測された性能に関する情報に限られているため。 ようするに、交差検証の過程で学習させた学習済みのモデルを手に入れる方法が、標準では用意されていない。 しかしながら、交差検証で性能が確かめられたモデルを取得したい、というニーズはあるはず。

なお、結果から先に書くとコールバック関数を使うことで学習済みモデルを手に入れることができた。

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.14.4
BuildVersion:   18E226
$ python -V          
Python 3.7.3
$ pip list | grep -i lightgbm
lightgbm        2.2.3  

下準備

まずは LightGBM と scikit-learn をインストールしておく。

$ brew install libomp
$ pip install lightgbm scikit-learn

cv() 関数からコールバック関数を使って学習済みモデルを取り出す

論よりコードということで、以下に cv() 関数からコールバック関数を使って学習済みモデルを取り出すサンプルコードを示す。 LightGBM の train() 関数や cv() 関数には、ブースティングのイテレーションごとに呼ばれるコールバック関数が登録できる。 一般的には、コールバック関数は学習の過程などを記録するために用いられる。 しかしながら、コールバック関数にはモデルも渡されるため、今回のような用途にも応用が効く。 内容の細かい解説に関してはコメントの形でコードに含めた。

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

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


class ModelExtractionCallback(object):
    """lightgbm.cv() から学習済みモデルを取り出すためのコールバックに使うクラス

    NOTE: 非公開クラス '_CVBooster' に依存しているため将来的に動かなく恐れがある
    """

    def __init__(self):
        self._model = None

    def __call__(self, env):
        # _CVBooster の参照を保持する
        self._model = env.model

    def _assert_called_cb(self):
        if self._model is None:
            # コールバックが呼ばれていないときは例外にする
            raise RuntimeError('callback has not called yet')

    @property
    def boosters_proxy(self):
        self._assert_called_cb()
        # Booster へのプロキシオブジェクトを返す
        return self._model

    @property
    def raw_boosters(self):
        self._assert_called_cb()
        # Booster のリストを返す
        return self._model.boosters

    @property
    def best_iteration(self):
        self._assert_called_cb()
        # Early stop したときの boosting round を返す
        return self._model.best_iteration


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

    # デモ用にデータセットを分割する
    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        shuffle=True,
                                                        random_state=42)

    # LightGBM 用のデータセット表現に直す
    lgb_train = lgb.Dataset(X_train, y_train)

    # 学習済みモデルを取り出すためのコールバックを用意する
    extraction_cb = ModelExtractionCallback()
    callbacks = [
        extraction_cb,
    ]

    # データセットを 5-Fold CV で学習する
    lgbm_params = {
        'objective': 'multiclass',
        'num_class': 3,
    }
    # NOTE: 一般的には返り値の内容 (交差検証の結果) を確認する
    lgb.cv(lgbm_params,
           lgb_train,
           num_boost_round=1000,
           early_stopping_rounds=10,
           nfold=5,
           shuffle=True,
           stratified=True,
           seed=42,
           callbacks=callbacks,
           )

    # コールバックのオブジェクトから学習済みモデルを取り出す
    proxy = extraction_cb.boosters_proxy
    boosters = extraction_cb.raw_boosters
    best_iteration = extraction_cb.best_iteration

    # 各モデルの推論結果を Averaging する場合
    y_pred_proba_list = proxy.predict(X_test,
                                      num_iteration=best_iteration)
    y_pred_proba_avg = np.array(y_pred_proba_list).mean(axis=0)
    y_pred = np.argmax(y_pred_proba_avg, axis=1)
    accuracy = accuracy_score(y_test, y_pred)
    print('Averaging accuracy:', accuracy)

    # 各モデルで個別に推論する場合
    for i, booster in enumerate(boosters):
        y_pred_proba = booster.predict(X_test,
                                       num_iteration=best_iteration)
        y_pred = np.argmax(y_pred_proba, axis=1)
        accuracy = accuracy_score(y_test, y_pred)
        print('Model {0} accuracy: {1}'.format(i, accuracy))

if __name__ == '__main__':
    main()

上記の実行結果は次の通り。 5-Fold CV から得られた全モデルを Averaging する場合と、各モデルごとに推論させた場合の精度 (Accuracy) が示される。

$ python lgbcv.py
...(snip)...
Averaging accuracy: 1.0
Model 0 accuracy: 1.0
Model 1 accuracy: 1.0
Model 2 accuracy: 0.9473684210526315
Model 3 accuracy: 1.0
Model 4 accuracy: 1.0

なお、今回のコードで注意すべきなのは、非公開のクラスを利用しているところがある点。 具体的には ModelExtractionCallback#boosters_proxy から得られるオブジェクトが lightgbm.engine._CVBooster というクラスのインスタンスになっている。 このクラスは、名前の先頭にアンダースコア (_) が含まれることから、外部に API として公開しているものではないと考えられる。 そのため、このモデルのインターフェースが唐突に変わったとしても文句はいえない。

一応、今の _CVBooster がどういったコードになっているのか示しておく。 実装があるのは以下の場所。 このクラスは内部にモデルをリストの形で保持している。 そして、オブジェクトに何らかの呼び出しがあると、特殊メソッドの __getattr__() でトラップされる。 トラップされた呼び出しは、内部のモデルに対して順番に実行されて結果がリストとして返されることになる。

github.com

いじょう。

(2019-04-15 追記)

この記事に掲載したコードをベースに OOF Prediction の実装を追加した Kaggle カーネルをまつけんさん (Twitter: @Kenmatsu4) が公開してくださいました。 ありがとうございます。

www.kaggle.com

Google Cloud SDK の CLI でデフォルトのゾーンやリージョンを設定する

Google Cloud SDK の CLI を使っていると何度も --zone--region といったオプションを指定することになってめんどくさい。 いつも同じ場所を使うのであれば、デフォルトを指定しておくと便利になる。

使った環境は次の通り。

$ sw_vers                                                   
ProductName:    Mac OS X
ProductVersion: 10.14.3
BuildVersion:   18D109
$ gcloud version                             
Google Cloud SDK 239.0.0
bq 2.0.42
core 2019.03.17
gsutil 4.37

デフォルトのゾーンを指定する

デフォルトのゾーンは gcloud config set compute/zone で指定する。

$ gcloud config set compute/zone <ZONE>

例えば東京の C ゾーンを指定してみる。

$ gcloud config set compute/zone asia-northeast1-c

試しにゾーンの指定をしないで Compute Engine のインスタンスを立ち上げてみよう。

$ gcloud compute instances create gce-example \
  --preemptible \
  --machine-type f1-micro \
  --image-project ubuntu-os-cloud \
  --image-family ubuntu-1804-lts

いつもなら、プロンプトでどのゾーンに作る?って聞かれるけど、デフォルトがあると聞かれない。

ちゃんとデフォルトのゾーンでインスタンスが作られる。

$ gcloud compute instances list
NAME         ZONE               MACHINE_TYPE  PREEMPTIBLE  INTERNAL_IP  EXTERNAL_IP   STATUS
gce-example  asia-northeast1-c  f1-micro      true         10.146.0.4   34.85.89.254  RUNNING

デフォルトのリージョンを設定する

同じようにリージョンも設定できる。

$ gcloud config set compute/region <REGION>

例えば東京リージョンを指定するなら、こんな感じ。

$ gcloud config set compute/region asia-northeast1

デフォルトの設定を削除する

デフォルトの設定がいらなくなったときは、次のようにして削除できる。

$ gcloud config unset compute/region
$ gcloud config unset compute/zone

いじょう。

Google Compute Engine で SSH Port Forwarding する

今回は Google Compute Engine のインスタンスで SSH Port Forwarding する方法について。 SSH Port Forwarding を使うと、インスタンスのポートをインターネットに晒すことなく利用できる。

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.14.3
BuildVersion:   18D109
$ gcloud version          
Google Cloud SDK 239.0.0
bq 2.0.42
core 2019.03.17
gsutil 4.37

下準備

まずは下準備から。

Google Cloud SDK をインストールする

今回は Google Cloud SDK の CLI を使って Compute Engine を操作する。 なので、最初に Homebrew Cask を使って Google Cloud SDK をインストールしておく。

$ brew cask install google-cloud-sdk

認証する

インストールできたら gcloud auth login コマンドで Google Cloud Platform の認証をしておく。

$ gcloud auth login

プロジェクトを作る

今回の動作確認をするプロジェクトを用意する。

$ gcloud projects create gce-example-$(whoami)-$(date "+%Y%m%d")
$ gcloud config set project gce-example-$(whoami)-$(date "+%Y%m%d")

API と課金を有効にする

プロジェクトで API と課金を有効にする。 この操作だけは CLI からできないのでブラウザから行う。

$ open https://console.cloud.google.com/apis/dashboard

API が使えないと実行できないコマンドが上手くいけばおっけー。

$ gcloud compute instances list
Listed 0 items.

インスタンスを起動する。

適当にインスタンスを立ち上げる。

$ gcloud compute instances create gce-example \
  --preemptible \
  --zone asia-northeast1-a \
  --machine-type f1-micro \
  --image-project ubuntu-os-cloud \
  --image-family ubuntu-1804-lts

立ち上げたインスタンスにログインする。

$ gcloud compute ssh --zone asia-northeast1-a gce-example

試しに Jupyter Notebook をインストールして起動する。 これでリモートの TCP:8888 ポートで Jupyter Notebook のサービスが動く。

gce-example $ sudo apt-get update
gce-example $ sudo apt-get -y install jupyter-notebook
gce-example $ jupyter notebook

準備ができたら、続いて SSH Port Forwarding する。

gcloud compute ssh では -- 以降に OpenSSH のオプションをそのまま渡せる。 これを利用してリモートの TCP:8888 をローカルホストの TCP:8888 にマッピングする。 -N オプションを指定すると、リモートでコマンドを実行しないことを表す。

$ gcloud compute ssh --zone asia-northeast1-a gce-example \
  -- -N -L 8888:localhost:8888

あとはブラウザでローカルホストの TCP:8888 ポートにアクセスすれば、いつもの画面が見える。

$ open http://localhost:8888

f:id:momijiame:20190323151508p:plain

めでたしめでたし。

Google Compute Engine のカスタムイメージを作る

今回は Google Compute Engine でカスタムイメージを作る方法について。 カスタムイメージというのは、既存のベースイメージに何らかのカスタマイズを施したイメージを指す。 あらかじめカスタムイメージを作っておくことで環境構築の手間をはぶくことができる場合がある。

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.14.3
BuildVersion:   18D109
$ gcloud version
Google Cloud SDK 238.0.0
bq 2.0.42
core 2019.03.08
gsutil 4.37

下準備

あらかじめ、いくつか下準備をしておく。

Google Cloud SDK をインストールする

今回の操作は基本的に Google Cloud SDK の CLI を使う。

なので、まずは Homebrew で Google Cloud SDK をインストールしておく。

$ brew cask install google-cloud-sdk

Google Cloud Platform にログインする

Google Cloud SDK をインストールできたら、gcloud auth login コマンドでログインしておく。

$ gcloud auth login

プロジェクトを作る

今回の動作確認用にプロジェクトを用意する。

$ gcloud projects create gce-example-$(whoami)-$(date "+%Y%m%d")

プロジェクトを作ったら、デフォルトのプロジェクトに指定する。

$ gcloud config set project gce-example-$(whoami)-$(date "+%Y%m%d")

課金と API を有効にする

続いては課金と API を有効にする。

ここだけは Google Cloud SDK からはできないのでブラウザから。 Compute Engine API を選んで有効にする。

$ open https://console.cloud.google.com/apis/dashboard

動作を確認する

API が有効になっていないと呼べないコマンドを実行してみて、上手くいけば準備完了。

$ gcloud compute instances list
Listed 0 items.

カスタムイメージを作る

下準備が終わったので、ここからはカスタムイメージを作っていく。 あらかじめ、基本的な手順について解説する。 カスタムイメージを作るときは、まず既存の任意のイメージを使って VM インスタンスを立ち上げる。 このとき、インスタンスが使うディスクには永続ディスク (Persistent Disk) を指定する。 その上で、インスタンスにソフトウェアをインストールするなどのカスタマイズを施す。 カスタマイズが終わったら、インスタンスを停止して永続ディスクを切り離す。 この切り離した永続ディスクから、カスタムイメージを作れる。 あとはカスタムイメージから VM インスタンスを立ち上げられれば作業完了。

永続化ディスクは pd- というプレフィックスから始まる。 リージョンごとに使えるディスクの種類が異なるので、確認しておく。

$ gcloud compute disk-types list | grep asia-northeast1
local-ssd    asia-northeast1-b          375GB-375GB
pd-ssd       asia-northeast1-b          10GB-65536GB
pd-standard  asia-northeast1-b          10GB-65536GB
local-ssd    asia-northeast1-c          375GB-375GB
pd-ssd       asia-northeast1-c          10GB-65536GB
pd-standard  asia-northeast1-c          10GB-65536GB
local-ssd    asia-northeast1-a          375GB-375GB
pd-ssd       asia-northeast1-a          10GB-65536GB
pd-standard  asia-northeast1-a          10GB-65536GB

各ディスクの特性については以下のページに記載がある。

cloud.google.com

各ディスクを利用するときの料金は以下。

cloud.google.com

また、これから作るイメージも保管には料金がかかる。 詳細については、以下に記載がある。

cloud.google.com

さて、早速永続ディスクを使ってインスタンスを立ち上げてみよう。 以下では 10GBpd-standard を使って Ubuntu 18.04 LTS のインスタンスを起動している。 料金を抑えるために Preemptible なインスタンスにした。

$ gcloud compute instances create gce-example \
  --preemptible \
  --zone asia-northeast1-a \
  --machine-type f1-micro \
  --image-project ubuntu-os-cloud \
  --image-family ubuntu-1804-lts \
  --boot-disk-type pd-standard \
  --boot-disk-size 10GB

上手く立ち上がれば gcloud compute instances list に表示される。

$ gcloud compute instances list
NAME         ZONE               MACHINE_TYPE  PREEMPTIBLE  INTERNAL_IP  EXTERNAL_IP   STATUS
gce-example  asia-northeast1-a  f1-micro      true         10.146.0.3   34.85.54.239  RUNNING

立ち上げたインスタンスを操作するために SSH でログインする。

$ gcloud compute ssh gce-example

ログインして lsblk コマンドを実行すると sda という名前で永続ディスクが認識されているようだ。

$ lsblk
NAME    MAJ:MIN RM  SIZE RO TYPE MOUNTPOINT
loop0     7:0    0 91.1M  1 loop /snap/core/6531
loop1     7:1    0 56.7M  1 loop /snap/google-cloud-sdk/75
sda       8:0    0   10G  0 disk
├─sda1    8:1    0  9.9G  0 part /
├─sda14   8:14   0    4M  0 part
└─sda15   8:15   0  106M  0 part /boot/efi

さて、ここからはイメージをカスタマイズするフェーズに入る。 試しにインスタンスに sl コマンドをインストールしておこう。

$ sudo apt-get update
$ sudo apt-get -y install sl
$ which sl
/usr/games/sl

カスタマイズが終わったらログアウトする。

$ exit

インスタンスを削除する。 このとき、永続ディスクについては削除されないように --keep-disks オプションを指定する。

$ gcloud compute instances delete gce-example \
  --keep-disks boot

これでインスタンスは削除された。

$ gcloud compute instances list
Listed 0 items.

ただし永続ディスクについては残っている。

$ gcloud compute disks list
NAME         ZONE               SIZE_GB  TYPE         STATUS
gce-example  asia-northeast1-a  10       pd-standard  READY

上記の永続ディスクを元にカスタムイメージを作る。 以下では ubuntu1804lts-sl という名前でカスタムイメージを作っている。

$ gcloud compute images create ubuntu1804lts-sl \
  --source-disk gce-example \
  --source-disk-zone asia-northeast1-a \
  --family ubuntu-1804-lts

次のように、カスタムイメージができた。

$ gcloud compute images list | grep gce-example-$(whoami)-$(date "+%Y%m%d")
ubuntu1804lts-sl                                      gce-example-20190321  ubuntu-1804-lts                               READY

カスタムイメージができたら、元になった永続ディスクは削除してしまっても構わない。

$ gcloud compute disks delete gce-example \
  --zone asia-northeast1-a

カスタムイメージを使ってインスタンスを立ち上げる

続いてはカスタムイメージを使ってインスタンスを立ち上げてみよう。

$ gcloud compute instances create gce-example \
  --preemptible \
  --zone asia-northeast1-a \
  --machine-type f1-micro \
  --image ubuntu1804lts-sl

起動したらインスタンスにログインする。

$ gcloud compute ssh gce-example

確認すると、たしかに sl コマンドがインストールされている。 ちゃんとカスタマイズされた状態になっているようだ。

$ which sl
/usr/games/sl

確認が終わったら、インスタンスからログアウトする。

$ exit

カスタムイメージを削除する

不要になったカスタムイメージは gcloud compute images delete で削除できる。

$ gcloud compute images delete ubuntu1804lts-sl

前述した通りカスタムイメージの保管にもお金がかかるので注意する。

後片付け

あとは後片付けもお忘れなく。

$ gcloud compute instances delete gce-example
$ gcloud projects delete gce-example-$(whoami)-$(date "+%Y%m%d")

いじょう。