CUBE SUGAR CONTAINER

技術系のこと書きます。

Python: pandas でグループごとにデータをサンプリングする

取り扱うデータをサンプリングする機会は意外と多い。 ユースケースとしては、例えばデータが多すぎて扱いにくい場合や、グループごとに件数の偏りのある場合が挙げられる。 今回は pandas を使ってグループごとに特定の件数をサンプリングする方法について書いてみる。

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.14.3
BuildVersion:   18D109
$ python -V                 
Python 3.7.2
$ pip list | grep pandas    
pandas          0.24.2 

下準備

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

$ pip install pandas

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

$ python

そして、サンプルとなる DataFrame を用意する。 今回は野菜とくだものが入ったデータをサンプルとして使うことにした。

>>> import pandas as pd
>>> data = [
...   ('Apple', 'Fruit'),
...   ('Beetroot', 'Vegetable'),
...   ('Carrot', 'Vegetable'),
...   ('Date', 'Fruit'),
...   ('Eggplant', 'Vegetable'),
...   ('Fig', 'Fruit'),
... ]
>>> df = pd.DataFrame(data, columns=['name', 'category'])

グループを加味しないでサンプリングする

今回用意したデータは野菜とくだものが 3 つずつ入っている。

>>> df
       name   category
0     Apple      Fruit
1  Beetroot  Vegetable
2    Carrot  Vegetable
3      Date      Fruit
4  Eggplant  Vegetable
5       Fig      Fruit

これをそのままランダムサンプリングすると、野菜とくだものが均等に取り出される確率が最も高い。 とはいえ、それはあくまで期待値としての話なので実際にはそうならない場合も多い。 次のように、グループを加味しない状態で DataFrame からサンプリングすると偏ることもある。

>>> df.sample(n=4)
       name   category
0     Apple      Fruit
3      Date      Fruit
5       Fig      Fruit
1  Beetroot  Vegetable

こうなると都合が悪い。

グループを加味してサンプリングする

続いてはグループについて加味した上で特定の件数をサンプリングしてみる。

まずは DataFrame#groupby()category ごとにグルーピングする。

>>> gdf = df.groupby('category')

これで得られるのは DataFrameGroupBy というクラスのインスタンス。

>>> gdf
<pandas.core.groupby.generic.DataFrameGroupBy object at 0x117e6e1d0>
>>> gdf.size()
category
Fruit        3
Vegetable    3
dtype: int64

DataFrameGroupByapply() メソッドを使うと、グループ単位の DataFrame に対して処理を実行できる。 結果は、全ての DataFrame が一枚に結合された状態で得られる。 つまり、DataFrameGroupBy#apply() の中で DataFrame に対して sample() メソッドを呼び出せば良い。 例えば野菜とくだものを 2 件ずつ取り出したいときは次のようにする。

>>> gdf.apply(lambda x: x.sample(n=2))
                 name   category
category                        
Fruit     0     Apple      Fruit
          5       Fig      Fruit
Vegetable 1  Beetroot  Vegetable
          4  Eggplant  Vegetable

これで、野菜とくだものが必ず毎回 2 件ずつ取り出せるようになった。

なお、ランダムサンプリングの結果を安定させたいときは random_state オプションを指定すれば良い。

>>> gdf.apply(lambda x: x.sample(n=2, random_state=42))
                 name   category
category                        
Fruit     0     Apple      Fruit
          3      Date      Fruit
Vegetable 1  Beetroot  Vegetable
          2    Carrot  Vegetable

ちなみに、元のグループの比率を維持してサンプリングしたいときは、次のように DataFrame の長さを用いると良い。

>>> gdf.apply(lambda x: x.sample(n=round(len(x) * 0.5)))
                 name   category
category                        
Fruit     0     Apple      Fruit
          3      Date      Fruit
Vegetable 4  Eggplant  Vegetable
          2    Carrot  Vegetable

MultiIndex を解除する

ちなみに、このやり方で作った DataFrame はインデックスが MultiIndex になっている。

>>> sampled_df = gdf.apply(lambda x: x.sample(n=2))
>>> sampled_df.index
MultiIndex(levels=[['Fruit', 'Vegetable'], [0, 1, 4, 5]],
           codes=[[0, 0, 1, 1], [3, 0, 1, 2]],
           names=['category', None])

これはこれでグループごとに要素を取り出しやすくて良い。

>>> sampled_df.loc['Fruit']
    name category
5    Fig    Fruit
0  Apple    Fruit

もし、元の DataFrame のフォーマットと合わせたいときには DataFrame#reset_index() を使う。

>>> sampled_df.reset_index(level='category', drop=True)
       name   category
5       Fig      Fruit
0     Apple      Fruit
1  Beetroot  Vegetable
4  Eggplant  Vegetable

サンプリングしたい件数が実在する件数よりも多いとき

先ほどのデータはグループごとに件数が均衡だった。 しかし、現実世界のデータではこのようにキレイな分布をするものはあまりない。 次は分布に偏りがあるときに生じる悩みについて考えてみる。

次のデータは、先ほどのデータからくだものを 2 つほど取り除いたものを扱う。 これで、データに含まれるくだものは Fig (いちじく) だけになった。

>>> lack_df = df[~df.name.isin(['Apple', 'Date'])]
>>> lack_df
       name   category
1  Beetroot  Vegetable
2    Carrot  Vegetable
4  Eggplant  Vegetable
5       Fig      Fruit

先ほどと同じように category の列でグルーピングする。

>>> lack_gdf = lack_df.groupby('category')
>>> lack_gdf.size()
category
Fruit        1
Vegetable    3
dtype: int64

この状態で、野菜とくだものを 2 件ずつ取り出したいという場合について考えてみる。 しかし、データの中に存在するくだものは 1 件しかないので 2 件というサンプル数に満たない。 実際に実行してみると、次のようなエラーになってしまう。

>>> lack_gdf.apply(lambda x: x.sample(n=2))
Traceback (most recent call last):
...
ValueError: Cannot take a larger sample than population when 'replace=False'

こうしたときに、データがサンプル数に満たないものは全てそのまま取り出したいのであれば、次のようにする。 具体的には、組み込み関数 min() を使って、DataFrame の長さとサンプル数の少ない方に合わせてしまう。

>>> lack_gdf.apply(lambda x: x.sample(n=min(2, len(x))))
                 name
category             
Fruit     5       Fig
Vegetable 1  Beetroot
          2    Carrot

もし、重複 (繰り返し) を許してでもサンプル数にあわせてほしいときは、次のように replace オプションに True を指定する。 こういった、本来ある件数よりも多くサンプリングする処理はオーバーサンプリングと呼ばれる。

>>> lack_gdf.apply(lambda x: x.sample(n=2, replace=True))
                 name
category             
Fruit     5       Fig
          5       Fig
Vegetable 1  Beetroot
          4  Eggplant

いじょう。

SQLite3 のテーブルに CSV でデータを読み込む

メモリに乗り切らない程度のちょっとした集計をするのに SQLite3 のデータベースを使うのが意外と便利だなーと思っている今日このごろ。 サポートされている SQL の構文や関数が少ないするのはネックだけど、手軽さには変えられないという感じ。 今回は SQLite3 のテーブルに CSV のデータをインポートする方法について書いておく。

使った環境は以下の通り。

$ sw_vers          
ProductName:    macOS
ProductVersion: 12.1
BuildVersion:   21C52
$ uname -rm
21.2.0 arm64
$ sqlite3 --version
3.36.0 2021-06-18 18:58:49 d24547a13b6b119c43ca2ede05fecaa707068f18c7430d47fc95fb5a2232aapl

下準備

データとして Seaborn の Taxis データセットを使うので、あらかじめダウンロードしておく。 このとき、先頭行のヘッダは読み飛ばす。 これは、後ほど登場する SQLite3 の .import 命令が、ヘッダを読み飛ばすのに対応していないため。

$ brew install wget
$ wget -O - https://raw.githubusercontent.com/mwaskom/seaborn-data/master/taxis.csv | sed -e "1d" > /tmp/taxis.csv
$ head -n 5 /tmp/taxis.csv
2019-03-23 20:21:09,2019-03-23 20:27:24,1,1.6,7.0,2.15,0.0,12.95,yellow,credit card,Lenox Hill West,UN/Turtle Bay South,Manhattan,Manhattan
2019-03-04 16:11:55,2019-03-04 16:19:00,1,0.79,5.0,0.0,0.0,9.3,yellow,cash,Upper West Side South,Upper West Side South,Manhattan,Manhattan
2019-03-27 17:53:01,2019-03-27 18:00:25,1,1.37,7.5,2.36,0.0,14.16,yellow,credit card,Alphabet City,West Village,Manhattan,Manhattan
2019-03-10 01:23:59,2019-03-10 01:49:51,1,7.7,27.0,6.15,0.0,36.95,yellow,credit card,Hudson Sq,Yorkville West,Manhattan,Manhattan
2019-03-30 13:27:42,2019-03-30 13:37:14,3,2.16,9.0,1.1,0.0,13.4,yellow,credit card,Midtown East,Yorkville West,Manhattan,Manhattan

CSV ファイルにあわせてテーブルの定義を用意する。 データベースは taxis.db という名前で作成する。

$ cat << 'EOF' | sqlite3 taxis.db
CREATE TABLE IF NOT EXISTS taxis (
  pickup DATETIME,
  dropoff DATETIME,
  passengers INT,
  distance FLOAT,
  fare FLOAT,
  tip FLOAT,
  tolls FLOAT,
  total FLOAT,
  color TEXT,
  payment TEXT,
  pickup_zone TEXT,
  dropoff_zone TEXT,
  pickup_borough TEXT,
  dropoff_borough TEXT
);
EOF

データの区切り文字をカンマ (,) にしたら .import 命令を使って CSV ファイルをテーブルに読み込む。

$ cat << 'EOF' | sqlite3 taxis.db
.separator ,
.import /tmp/taxis.csv taxis
EOF

ちゃんと読み込まれているか確認しよう。

$ cat << 'EOF' | sqlite3 taxis.db
.headers on
SELECT * FROM taxis LIMIT 5
EOF
pickup|dropoff|passengers|distance|fare|tip|tolls|total|color|payment|pickup_zone|dropoff_zone|pickup_borough|dropoff_borough
2019-03-23 20:21:09|2019-03-23 20:27:24|1|1.6|7.0|2.15|0.0|12.95|yellow|credit card|Lenox Hill West|UN/Turtle Bay South|Manhattan|Manhattan
2019-03-04 16:11:55|2019-03-04 16:19:00|1|0.79|5.0|0.0|0.0|9.3|yellow|cash|Upper West Side South|Upper West Side South|Manhattan|Manhattan
2019-03-27 17:53:01|2019-03-27 18:00:25|1|1.37|7.5|2.36|0.0|14.16|yellow|credit card|Alphabet City|West Village|Manhattan|Manhattan
2019-03-10 01:23:59|2019-03-10 01:49:51|1|7.7|27.0|6.15|0.0|36.95|yellow|credit card|Hudson Sq|Yorkville West|Manhattan|Manhattan
2019-03-30 13:27:42|2019-03-30 13:37:14|3|2.16|9.0|1.1|0.0|13.4|yellow|credit card|Midtown East|Yorkville West|Manhattan|Manhattan

ばっちり。

Python: テストで SQLite3 のインメモリデータベースを使うときの問題点と解決策

今回は SQLite3 のインメモリデータベースをテストで使うときに生じる問題点と、その解決策について。

SQLite3 のインメモリデータベースを使うと、追加でソフトウェアをインストールしたり、データベースファイルを作ることなくリレーショナルデータベースを扱うことができる。 この点はリレーショナルデータベースを扱うソフトウェアを作るときに自動テストを組む上で望ましい特性といえる。 ただ、SQLite3 のインメモリデータベースには、制約が複数ある。 まず一つ目はコネクションを閉じると永続化した内容が消えてしまうところ。 そして、二つ目は異なるスレッドから内容を参照できないところ。

これらの解決策として、無理に制約のあるインメモリデータベースを使わずテンポラリファイルを使ってディスクに永続化することを提案する。

使った環境は次の通り。

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

インメモリデータベースの問題点について

まずはインメモリデータベースを使うときの問題点について書いていく。

Python のインタプリタを起動しておこう。

$ python

コネクションを閉じるとデータが消えてしまう問題

最初に、SQLite3 のインメモリデータベースへのコネクションを作成する。

>>> import sqlite3
>>> connection = sqlite3.connect(':memory:')

コネクションからカーソルオブジェクトを生成する。

>>> cursor = connection.cursor()

カーソルオブジェクトを使うことで SQL 文を実行できる。 まずはテーブルの作成から。

>>> create_query = """
... CREATE TABLE users (
...   id INTEGER,
...   age INTEGER NOT NULL,
...   name TEXT NOT NULL,
...   PRIMARY KEY (id)
... );
... """
>>> cursor.execute(create_query)
<sqlite3.Cursor object at 0x1037fa260>

続いて、作成したテーブルにレコードを追加する。

>>> insert_query = """
... INSERT INTO users
... VALUES
...   (1, 20, 'Alice'),
...   (2, 30, 'Bob'),
...   (3, 40, 'Carol');
... """
>>> cursor.execute(insert_query)
<sqlite3.Cursor object at 0x1037fa260>

この状態で、テーブルを SELECT する SQL 文を発行してみよう。 すると、ちゃんと永続化された内容が取得できる。

>>> select_query = """
... SELECT
...   *
... FROM users
... """
>>> cursor.execute(select_query)
<sqlite3.Cursor object at 0x1037fa260>
>>> cursor.fetchall()
[(1, 20, 'Alice'), (2, 30, 'Bob'), (3, 40, 'Carol')]

しかし、ここで一旦コネクションを閉じるとどうなるだろうか?

>>> connection.commit()
>>> connection.close()

改めてインメモリデータベースへのコネクションを開いてみよう。

>>> connection = sqlite3.connect(':memory:')
>>> cursor = connection.cursor()

そして SELECT 文を発行する。 すると、テーブルがないというエラーになってしまった。

>>> cursor.execute(select_query)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
sqlite3.OperationalError: no such table: users

このように、インメモリデータベースでは同じコネクションの中でしか永続化した内容が参照できない。

ちなみに、この問題は何もコネクションを閉じなければ問題ないかというと、そういう話でもない。 別のコネクションとして開いた場合にも問題になる。

まずはあるコネクションでテーブルを作ってレコードを追加しておく。

>>> conn1 = sqlite3.connect(':memory:')
>>> cur1 = conn1.cursor()
>>> create_query = """
... CREATE TABLE users (
...   id INTEGER,
...   age INTEGER NOT NULL,
...   name TEXT NOT NULL,
...   PRIMARY KEY (id)
... );
... """
>>> cur1.execute(create_query)
<sqlite3.Cursor object at 0x105c8b260>
>>> insert_query = """
... INSERT INTO users
... VALUES
...   (1, 20, 'Alice'),
...   (2, 30, 'Bob'),
...   (3, 40, 'Carol');
... """
>>> cur1.execute(insert_query)
<sqlite3.Cursor object at 0x105c8b260>

そして、別のコネクションを開いて SELECT 文を発行してみる。

>>> conn2 = sqlite3.connect(':memory:')
>>> cur2 = conn2.cursor()
>>> select_query = """
... SELECT
...   *
... FROM users
... """
>>> cur2.execute(select_query)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
sqlite3.OperationalError: no such table: users

すると、やはりテーブルがないとエラーになってしまう。

異なるスレッド間でデータベースを参照できない問題

続いては異なるスレッド間でデータベースを参照できない問題について。

まずは、先ほどと同じように動作確認用のデータベースを作成しておく。

>>> connection = sqlite3.connect(':memory:')
>>> cursor = connection.cursor()
>>> create_query = """
... CREATE TABLE users (
...   id INTEGER,
...   age INTEGER NOT NULL,
...   name TEXT NOT NULL,
...   PRIMARY KEY (id)
... );
... """
>>> cursor.execute(create_query)
<sqlite3.Cursor object at 0x1037fa260>
>>> insert_query = """
... INSERT INTO users
... VALUES
...   (1, 20, 'Alice'),
...   (2, 30, 'Bob'),
...   (3, 40, 'Carol');
... """
>>> cursor.execute(insert_query)
<sqlite3.Cursor object at 0x1037fa260>

続いて、次のようにデータベースのコネクションからレコードを SELECT する関数を定義する。

>>> def select(connection):
...     cursor = connection.cursor()
...     select_query = """
...     SELECT
...       *
...     FROM users
...     """
...     cursor.execute(select_query)
...     print(cursor.fetchall())
... 

上記の関数を別のスレッドから実行する。

>>> import threading
>>> t = threading.Thread(target=select, args=(connection,))

すると、次のようにエラーになる。 SQLite のオブジェクトは作成したスレッドでしか扱ってはいけない、という旨のメッセージが表示されている。 しかし、異なるコネクションを開いても内容が参照できないのは先ほど見た通り。

>>> t.start()
>>> Exception in thread Thread-1:
Traceback (most recent call last):
  File "/usr/local/Cellar/python/3.7.2_2/Frameworks/Python.framework/Versions/3.7/lib/python3.7/threading.py", line 917, in _bootstrap_inner
    self.run()
  File "/usr/local/Cellar/python/3.7.2_2/Frameworks/Python.framework/Versions/3.7/lib/python3.7/threading.py", line 865, in run
    self._target(*self._args, **self._kwargs)
  File "<stdin>", line 2, in select
sqlite3.ProgrammingError: SQLite objects created in a thread can only be used in that same thread. The object was created in thread id 4478109120 and this is thread id 123145389289472.

テンポラリファイルを使った解決策について

続いては、ここまで見てきた問題の解決策について提案する。 具体的には、無理に成約あるインメモリデータベースを使うのではなく、テンポラリファイルを使ってディスクに永続化するというもの。 Python には、使い終わった際に自動で削除されるファイルを扱うためのモジュールとして tempfile が用意されている。

まずは、このモジュールをインポートする。

>>> import tempfile

名前付きテンポラリファイルのオブジェクトをインスタンス化しよう。

>>> tfile = tempfile.NamedTemporaryFile()

このオブジェクトは、作成された時点でファイルシステムのどこかに対応するファイルが作成される。

>>> tfile.name
'/var/folders/1f/2k50hyvd2xq2p4yg668_ywsh0000gn/T/tmpc08x657m'

別のターミナルから確認すると、ちゃんとファイルができていることが分かる。

$ ls -alF /var/folders/1f/2k50hyvd2xq2p4yg668_ywsh0000gn/T/tmpc08x657m
-rw-------  1 amedama  staff  0  3  3 18:36 /var/folders/1f/2k50hyvd2xq2p4yg668_ywsh0000gn/T/tmpc08x657m

このファイルは、オブジェクトが GC されると自動的に削除される。

>>> del tfile

次のように、ファイルがなくなった。

$ ls -alF /var/folders/1f/2k50hyvd2xq2p4yg668_ywsh0000gn/T/tmpc08x657m
ls: /var/folders/1f/2k50hyvd2xq2p4yg668_ywsh0000gn/T/tmpc08x657m: No such file or directory

続いては、このテンポラリファイルを使って SQLite3 の動作を確認してみよう。 テンポラリファイルのパスを使って SQLite3 のオブジェクトを作成する。

>>> tfile = tempfile.NamedTemporaryFile()
>>> connection = sqlite3.connect(tfile.name)
>>> cursor = connection.cursor()

データベースにテーブルを作って、レコードを追加する。

>>> cursor.execute(create_query)
<sqlite3.Cursor object at 0x105d1f0a0>
>>> cursor.execute(insert_query)
<sqlite3.Cursor object at 0x105d1f0a0>

終わったらコネクションを閉じてしまおう。

>>> connection.commit()
>>> connection.close()

その上で、テンポラリファイルのパスを使って再度コネクションを開く。

>>> connection = sqlite3.connect(tfile.name)
>>> cursor = connection.cursor()

そして SELECT 文を発行すると、ちゃんと先ほど永続化した内容が取れることが分かる。

>>> cursor.execute(select_query)
<sqlite3.Cursor object at 0x105edd6c0>
>>> cursor.fetchall()
[(1, 20, 'Alice'), (2, 30, 'Bob'), (3, 40, 'Carol')]

テンポラリファイルを使い終わって削除すれば、永続化に使っていた SQLite3 のデータベースファイルは自動的に削除される。

>>> del tfile

実際にテストを書くときのイメージ

もう少し具体的にテストを書くときのイメージをつけるために pytest を使った例を示す。

まずは pytest をインストールしておく。

$ pip install pytest

以下は SQLite3 を使ってリレーショナルデータベースの動作を確認するテストのサンプルコード。

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

import tempfile
import sqlite3

import pytest


@pytest.fixture(scope='module', autouse=True)
def temp_dbfile():
    # テンポラリファイルを用意する
    dbfile = tempfile.NamedTemporaryFile()

    # SQLite3 のデータベースファイルとして使う
    conn = sqlite3.connect(dbfile.name)
    c = conn.cursor()

    # テーブルを作る
    create_query = """
      CREATE TABLE users (
        id INTEGER,
        age INTEGER NOT NULL,
        name TEXT NOT NULL,
        PRIMARY KEY (id)
      );
    """
    c.execute(create_query)

    # データを追加する
    insert_query = """
      INSERT INTO users VALUES
        (1, 20, 'Alice'),
        (2, 30, 'Bob'),
        (3, 40, 'Carol');
    """
    c.execute(insert_query)

    # 永続化する
    conn.commit()
    conn.close()

    return dbfile


def test_example(temp_dbfile):
    # テンポラリファイルの名前を使ってデータベースを開き直す
    conn = sqlite3.connect(temp_dbfile.name)
    c = conn.cursor()

    # データベースに永続化されているユーザを取り出す
    select_query = """
    SELECT
      *
    FROM users
    WHERE
      age > 35
    """
    c.execute(select_query)
    fetched_users = c.fetchall()

    # 取り出した内容を検証する
    assert [(3, 40, 'Carol')] == fetched_users


if __name__ == '__main__':
    pytest.main(['-v', __file__])

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

$ py.test -v test_example.py 
============================= test session starts ==============================
platform darwin -- Python 3.7.2, pytest-4.3.0, py-1.8.0, pluggy-0.9.0 -- /Users/amedama/.virtualenvs/py37/bin/python3.7
cachedir: .pytest_cache
rootdir: /Users/amedama/Documents/temporary, inifile:
collected 1 item                                                               

test_example.py::test_example PASSED                                     [100%]

=========================== 1 passed in 0.04 seconds ===========================

うまくテストがパスした。

注意事項

ここまで見てきた通り SQLite3 を使うとリレーショナルデータベースの自動テストが書きやすいのは、確かにそう。 ただ、リレーショナルデータベースは使うソフトウェアによってクセが強いので、結果をあまり信用しすぎない方が良い。 例えば、SQLite3 の自動テストで動作確認しながら開発していたけど、本番で使うデータベースと結合したときに動かないみたいなことが起こりうる。 これは、SQLAlchemy のような O/R マッパーを間に挟んでいても同じ。 そのため、本番で MySQL や PostgreSQL を使うのであれば、SQLite3 の自動テストはあくまで軽い動作確認程度にとどめた方が良い。 早めに Docker などを用いて本番で使うデータベースを使った E2E テストを書いておくと、だいぶ心理的には安心できる。

いじょう。

Python: Adversarial Validation について

最近、Kaggle などのデータ分析コンペで使われることの多い Adversarial Validation という手法について調べたり考えていたので書いてみる。

もくじ

背景

Adversarial Validation という手法は、データ分析コンペに存在する、ある課題を解決するために考案された。 その課題とは、提供される複数のデータセットの分布が異なる場合に、いかにして正しく予測するかというもの。

まず、データ分析コンペでは、一般に学習用データと検証用データが提供される。 学習用データには真のラベル情報が与えられるのに対して、検証用データにはそれがない。 参加者は学習用データでトレーニングしたモデルを用いて、検証用データを予測した結果をサブミットする。

サブミットした内容に対するフィードバックは、検証用データの一部を用いて評価した結果として得られる。 検証用データの一部を用いた評価は、Kaggle では Public Leader Board (Public LB) と呼ばれる。 ただし、これはあくまで一部を用いた評価に過ぎないため、検証用データの残りの部分を使ったときにどうなるかは分からない。 検証用データの残りの部分を使った評価は、コンペが終了したときに Private Leader Board (Private LB) として公表される。

Public LB にもとづいてモデルの良し悪しを評価すると、ときとして Private LB で大きく順位が下がることがある。 この現象は、Public LB の評価で用いられる一部のデータに対して、モデルが過剰適合した状態と考えられる。 これを防ぐには、Cross Validation (交差検証) などを用いたローカルでの汎化性能の評価 (Local CV) が重要になる。

しかし、コンペで提供される学習用データと検証用データは、ときとして分布が大きく異なる。 このとき何が起こるかというと、Local CV と Public LB / Private LB のスコアが相関しなくなる。 例えば Local CV では高いスコアが得られるのに、サブミットした結果では低いスコアしか得られない、といった事態が起こる。

Adversarial Validation

この問題を解決するために考案されたのが Adversarial Validation という手法。 Adversarial Validation では、まず二つのデータセットに本来の目的変数とは別の目的変数を設定する。 具体的には、どちらのデータセットに由来するデータなのかという情報を目的変数として付与する。 例えば学習用データが全て 0 で検証用データは全て 1 といった形で付与する。

そして、学習用データと検証用データを混ぜた状態で、両者がどちらに由来するものなのかを予測するモデルを作る。 学習用データと検証用データの分布が異なるほど、どちらに由来するかの予測は容易になる。 もし、全く同じ分布から得られたデータなのであれば、両者を区別することはそもそもできない。 この時点で、学習用データと検証用データの分布がどれくらい異なるのかが分かる。

もし分布が異なるのであれば、学習用データのサブセットを作ることになる。 サブセットを作る上での基準は、検証用データにより近いデータかどうか。 これは、先ほど確認した内容にもとづいて作成する。 つまり、学習用データの中から、検証用データと判断が難しかったものを選び出す。

具体的には、検証用データかどうかの確率を降順ソートして上から取り出せば良い。 その基準で取り出した学習用データのサブセットを使って Local CV することで、Public LB / Private LB との相関が改善する。 以上が Adversarial Validation という手法の説明になる。

試してみる

前置きがだいぶ長くなったけど、ここからは実際に Adversarial Validation を試してみよう。

使った環境は次の通り。

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

下準備

まずは、検証に使うパッケージをインストールする。

$ pip install scikit-learn matplotlib

二つのデータが同じ分布に由来するとき

まずは、二つのデータが全く同じ分布に由来するとき Adversarial Validation で何が起こるか確認する。

以下のサンプルコードでは二値分類問題に使うダミーのデータを生成している。 総数 5,000 点の二次元データで、割合は 1:1 の均衡データとなっている。 そのデータを 2,500 点ずつランダムにサンプリングして、これを学習用データと検証用データに模している。 これらのデータに、どちらに由来する要素なのかを表す目的変数 z_* を用意してランダムフォレストで二値分類として解いている。 均衡データなので、評価指標には単純に精度 (Accuracy) を用いた。

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

import numpy as np
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.model_selection import cross_validate
from sklearn.model_selection import StratifiedKFold
from matplotlib import pyplot as plt
from sklearn.ensemble import RandomForestClassifier


def main():
    # データ点数 5,000 のデータセットを作るパラメータ
    args = {
        'n_samples': 5000,
        'n_features': 2,
        'n_informative': 2,
        'n_redundant': 0,
        'n_repeated': 0,
        'n_classes': 2,
        'n_clusters_per_class': 1,
        'flip_y': 0,
        'class_sep': 1.5,
        'weights': [0.5, 0.5],
        'random_state': 42,
    }
    # 二次元のデータを生成する
    X, y = make_classification(**args)

    # データをランダムに二等分する
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5,
                                                        shuffle=True, random_state=42)

    # 分割したデータがどちらに由来したものか区別するためのラベル
    z_train = np.zeros(len(X_train))  # 片方は全て 0
    z_test = np.ones(len(X_test))  # もう片方は全て 1

    # 分割したデータを結合する
    X_concat = np.concatenate([X_train, X_test], axis=0)
    z_concat = np.concatenate([z_train, z_test], 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_concat, z_concat, cv=skf)
    # 5-Fold CV で評価した精度 (Accuracy) の平均
    print('Accuracy:', score['test_score'].mean())

    # 両データをプロットする
    plt.scatter(X_train[y_train == 0, 0],
                X_train[y_train == 0, 1],
                alpha=0.5,
                label='Train (Negative)')
    plt.scatter(X_train[y_train == 1, 0],
                X_train[y_train == 1, 1],
                alpha=0.5,
                label='Train (Positive)')
    plt.scatter(X_test[y_test == 0, 0],
                X_test[y_test == 0, 1],
                alpha=0.5,
                label='Test (Negative)')
    plt.scatter(X_test[y_test == 1, 0],
                X_test[y_test == 1, 1],
                alpha=0.5,
                label='Test (Positive)')
    plt.legend()
    plt.show()


if __name__ == '__main__':
    main()

上記の実行結果は次の通り。 精度 (Accuracy) は 0.5 を割っており、ランダムに答えたときの期待値よりも低い。

$ python samedist.py
Accuracy: 0.47759999999999997

二つのデータをプロットした散布図は次の通り。

f:id:momijiame:20190223142636p:plain

このように、全く同一の分布から得られたものであればデータの由来を予測することは困難といえる。

二つのデータが異なる分布に由来するとき

続いては二つのデータが異なる分布に由来するときに Adversarial Validation で何が起こるかを確認する。

次のサンプルコードでは、乱数を変えて 2,500 点のデータを二つ作っている。 二つのデータはそれなりに似ているものの結構違う、という絶妙なものを探してみた。 先ほどと同じようにどちらのデータに由来するものか目的変数を付与して、それをランダムフォレストで分類している。

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

import numpy as np
from sklearn.datasets import make_classification
from sklearn.model_selection import cross_validate
from sklearn.model_selection import StratifiedKFold
from matplotlib import pyplot as plt
from sklearn.ensemble import RandomForestClassifier


def main():
    # データ点数 2,500 のデータセットを作るパラメータ
    args = {
        'n_samples': 2500,
        'n_features': 2,
        'n_informative': 2,
        'n_redundant': 0,
        'n_repeated': 0,
        'n_classes': 2,
        'n_clusters_per_class': 1,
        'flip_y': 0,
        'class_sep': 1.5,
        'weights': [0.5, 0.5],
        'random_state': 42,
    }
    # 二次元のデータを生成する
    X_train, y_train = make_classification(**args)
    # 乱数を変えてもう一度作る
    args['random_state'] = 424242
    X_test, y_test = make_classification(**args)

    # 分割したデータがどちらに由来したものか区別するためのラベル
    z_train, z_test = np.zeros(len(X_train)), np.ones(len(X_test))

    # 結合する
    X = np.concatenate([X_train, X_test], axis=0)
    z = np.concatenate([z_train, z_test], 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, z, cv=skf)
    print('Accuracy:', score['test_score'].mean())

    # 両データをプロットする
    plt.scatter(X_train[y_train == 0, 0],
                X_train[y_train == 0, 1],
                alpha=0.5,
                label='Train (Negative)')
    plt.scatter(X_train[y_train == 1, 0],
                X_train[y_train == 1, 1],
                alpha=0.5,
                label='Train (Positive)')
    plt.scatter(X_test[y_test == 0, 0],
                X_test[y_test == 0, 1],
                alpha=0.5,
                label='Test (Negative)')
    plt.scatter(X_test[y_test == 1, 0],
                X_test[y_test == 1, 1],
                alpha=0.5,
                label='Test (Positive)')
    plt.legend()
    plt.show()


if __name__ == '__main__':
    main()

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

$ python diffdist.py 
Accuracy: 0.8498000000000001

今度は、それなりに分類できている。

データをプロットした散布図は次の通り。 Positive なデータの分布が、特に異なることが分かる。

f:id:momijiame:20190223143315p:plain

このように、分布が異なる場合には Adversarial Validation で、それを確認できる。

異なる分布でそのまま分類してみる

続いて、先ほど確認した分布の異なるデータを使って学習・予測してみよう。 これは Adversarial Validation を使わないアプローチで結果をサブミットしたときの Private LB のスコアに相当する。

次のサンプルコードでは、先ほど生成したのと同じデータセットを用意した。 そして、学習用データでランダムフォレストをトレーニングしてから、検証用データで予測した場合のスコアを確認している。

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

import numpy as np
from matplotlib import pyplot as plt
from sklearn.datasets import make_classification
from sklearn.metrics import accuracy_score
from sklearn.ensemble import RandomForestClassifier


def main():
    args = {
        'n_samples': 2500,
        'n_features': 2,
        'n_informative': 2,
        'n_redundant': 0,
        'n_repeated': 0,
        'n_classes': 2,
        'n_clusters_per_class': 1,
        'flip_y': 0,
        'class_sep': 1.5,
        'weights': [0.5, 0.5],
        'random_state': 42,
    }
    X_train, y_train = make_classification(**args)

    args['random_state'] = 424242
    X_test, y_test = make_classification(**args)

    clf = RandomForestClassifier(n_estimators=100,
                                 random_state=42)

    # 最初に作ったデータセットで学習する
    clf.fit(X_train, y_train)

    # 次に作ったデータセット (分布が異なる) でモデルを評価する
    y_pred = clf.predict(X_test)
    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)

    # モデルの識別境界を可視化する
    x_mesh, y_mesh = np.meshgrid(np.linspace(-5, 5, 50),
                                 np.linspace(-5, 5, 50))
    Z = clf.predict(np.c_[x_mesh.ravel(),
                          y_mesh.ravel()])
    Z = Z.reshape(x_mesh.shape)
    plt.contourf(x_mesh, y_mesh, Z, cmap=plt.cm.binary)

    # 元のデータをプロットする
    plt.scatter(X_train[y_train == 0, 0],
                X_train[y_train == 0, 1],
                alpha=0.5,
                label='Train (Negative)')
    plt.scatter(X_train[y_train == 1, 0],
                X_train[y_train == 1, 1],
                alpha=0.5,
                label='Train (Positive)')
    plt.scatter(X_test[y_test == 0, 0],
                X_test[y_test == 0, 1],
                alpha=0.5,
                label='Test (Negative)')
    plt.scatter(X_test[y_test == 1, 0],
                X_test[y_test == 1, 1],
                alpha=0.5,
                label='Test (Positive)')
    plt.legend()
    plt.show()

if __name__ == '__main__':
    main()

上記の実行結果は次の通り。 精度 (Accuracy) として 0.9516 という結果が得られた。

$ python nochange.py 
Accuracy: 0.9516

モデルの識別境界をプロットしたグラフが次の通り。

f:id:momijiame:20190223144159p:plain

Test (Positive) のデータの一部が誤って分類されていることが分かる。 とはいえモデルは Test に関する情報を知らないので妥当な結果といえる。

検証用データに似ているものを取り出す

続いては、Adversarial Validation を使って検証用データに似ているものを学習用データから取り出してみよう。

次のサンプルコードでは、Adversarial Validation で検証用データに近いと判断された学習用データの上位 1,000 件を取り出している。 ただし、取り出したデータを使った学習まではしておらずグラフへのプロットにとどまる。

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

import numpy as np
from sklearn.datasets import make_classification
from sklearn.model_selection import cross_val_predict
from matplotlib import pyplot as plt
from sklearn.ensemble import RandomForestClassifier


def main():
    args = {
        'n_samples': 2500,
        'n_features': 2,
        'n_informative': 2,
        'n_redundant': 0,
        'n_repeated': 0,
        'n_classes': 2,
        'n_clusters_per_class': 1,
        'flip_y': 0,
        'class_sep': 1.5,
        'weights': [0.5, 0.5],
        'random_state': 42,
    }
    X_train, y_train = make_classification(**args)

    args['random_state'] = 424242
    X_test, y_test = make_classification(**args)

    # 分割したデータがどちらに由来したものか区別するためのラベル
    z_train, z_test = np.zeros(len(X_train)), np.ones(len(X_test))

    X = np.concatenate([X_train, X_test], axis=0)
    z = np.concatenate([z_train, z_test], axis=0)

    # 要素が最初のデータセット由来なのか、次のデータセット由来なのか分類する
    clf = RandomForestClassifier(n_estimators=100,
                                 random_state=42)

    # 分類して確率を計算する
    z_pred_proba = cross_val_predict(clf, X, z,
                                     cv=5, method='predict_proba')

    # 先に生成したデータに対応する部分だけ取り出す
    z_train_pred_proba = z_pred_proba[:2500]
    # 後から作ったデータに近いと判断されたものを取り出す
    X_train_alike = X_train[np.argsort(z_train_pred_proba[:, 0])][:1000]
    # 上記の残り
    X_train_rest = X_train[np.argsort(z_train_pred_proba[:, 0])][1000:]

    # それぞれのデータをプロットする
    plt.scatter(X_test[:, 0],
                X_test[:, 1],
                alpha=0.5,
                label='Test')
    plt.scatter(X_train_rest[:, 0],
                X_train_rest[:, 1],
                alpha=0.5,
                label='NOT Resemblance')
    plt.scatter(X_train_alike[:, 0],
                X_train_alike[:, 1],
                alpha=0.5,
                label='Resemblance')

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


if __name__ == '__main__':
    main()

上記を実行してみる。

$ python resemblance.py

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

f:id:momijiame:20190223162721p:plain

Negative なデータにおいては中心部分が、Positive なデータにおいてはエッジ部分が似ていると判断されたようだ。

取り出したデータを使って学習してみる

続いては、先ほど取り出したデータを使って実際に学習してみよう。

以下のサンプルコードでは、先ほど取り出したのと同じデータを使ってモデルを学習している。 つまり、検証用データに近いと判断された学習用データの要素 1,000 点を使った。 そして、学習させたモデルを使って検証用データのスコアを計算している。 また、同時にモデルの識別境界も可視化した。

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

import numpy as np
from matplotlib import pyplot as plt
from sklearn.datasets import make_classification
from sklearn.metrics import accuracy_score
from sklearn.model_selection import cross_val_predict
from sklearn.ensemble import RandomForestClassifier


def main():
    args = {
        'n_samples': 2500,
        'n_features': 2,
        'n_informative': 2,
        'n_redundant': 0,
        'n_repeated': 0,
        'n_classes': 2,
        'n_clusters_per_class': 1,
        'flip_y': 0,
        'class_sep': 1.5,
        'weights': [0.5, 0.5],
        'random_state': 42,
    }
    X_train, y_train = make_classification(**args)

    args['random_state'] = 424242
    X_test, y_test = make_classification(**args)

    z_train, z_test = np.zeros(len(X_train)), np.ones(len(X_test))

    X = np.concatenate([X_train, X_test], axis=0)
    z = np.concatenate([z_train, z_test], axis=0)

    # 要素が最初のデータセット由来なのか、次のデータセット由来なのか分類する
    clf = RandomForestClassifier(n_estimators=100,
                                 random_state=42)
    z_pred_proba = cross_val_predict(clf, X, z,
                                     cv=5, method='predict_proba')

    # 先に生成したデータに対応する部分だけ取り出す
    z_train_pred_proba = z_pred_proba[:2500]
    # 後から作ったデータに近いと判断されたものを取り出す
    X_train_alike = X_train[np.argsort(z_train_pred_proba[:, 0])][:1000]
    y_train_alike = y_train[np.argsort(z_train_pred_proba[:, 0])][:1000]

    # 似ているデータで学習する
    clf.fit(X_train_alike, y_train_alike)
    # 後から作ったデータで評価する
    y_pred = clf.predict(X_test)
    # 精度 (Accuracy) を確認する
    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)

    # 識別境界を可視化する
    x_mesh, y_mesh = np.meshgrid(np.linspace(-5, 5, 50),
                                 np.linspace(-5, 5, 50))
    Z = clf.predict(np.c_[x_mesh.ravel(),
                          y_mesh.ravel()])
    Z = Z.reshape(x_mesh.shape)
    plt.contourf(x_mesh, y_mesh, Z, cmap=plt.cm.binary)

    plt.scatter(X_train[y_train == 0, 0],
                X_train[y_train == 0, 1],
                alpha=0.5,
                label='Train (Negative)')
    plt.scatter(X_train[y_train == 1, 0],
                X_train[y_train == 1, 1],
                alpha=0.5,
                label='Train (Positive)')
    plt.scatter(X_test[y_test == 0, 0],
                X_test[y_test == 0, 1],
                alpha=0.5,
                label='Test (Negative)')
    plt.scatter(X_test[y_test == 1, 0],
                X_test[y_test == 1, 1],
                alpha=0.5,
                label='Test (Positive)')
    plt.scatter(X_train_alike[:, 0],
                X_train_alike[:, 1],
                label='Model Learned')
    plt.legend()
    plt.show()


if __name__ == '__main__':
    main()

上記の実行結果は次の通り。 先ほど、学習用データをそのまま使ったパターンよりもスコアが改善している。

$ python advlearn.py 
Accuracy: 0.9788

識別境界を示すグラフは次の通り。

f:id:momijiame:20190223163622p:plain

ラベルごとに似ているデータを取り出す

先ほどの例ではラベルによらず検証用データに似たものを学習用データから 1,000 件取り出した。 ただ、可視化した内容を見ても分かる通り、似ているかどうかはラベルによって異なる場合がある。 例えば、先ほどの例では Negative なラベルに、より多くの似ているデータが含まれていた。 これでは、モデルに学習させるラベルの数が不均衡になる問題点がある。 なので、続いてはラベルごとに似ているデータを取り出してみることにしよう。

次のサンプルコードでは、ラベルごとに似ているデータを 200 件ずつ取り出していて学習に使っている。 取り出したデータの可視化と識別境界をプロットしている点については先ほどと変わらない。

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

import numpy as np
from matplotlib import pyplot as plt
from sklearn.datasets import make_classification
from sklearn.metrics import accuracy_score
from sklearn.model_selection import cross_val_predict
from sklearn.ensemble import RandomForestClassifier


def main():
    args = {
        'n_samples': 2500,
        'n_features': 2,
        'n_informative': 2,
        'n_redundant': 0,
        'n_repeated': 0,
        'n_classes': 2,
        'n_clusters_per_class': 1,
        'flip_y': 0,
        'class_sep': 1.5,
        'weights': [0.5, 0.5],
        'random_state': 42,
    }
    X_train, y_train = make_classification(**args)

    args['random_state'] = 424242
    X_test, y_test = make_classification(**args)

    z_train, z_test = np.zeros(len(X_train)), np.ones(len(X_test))

    X = np.concatenate([X_train, X_test], axis=0)
    z = np.concatenate([z_train, z_test], axis=0)

    # 要素が最初のデータセット由来なのか、次のデータセット由来なのか分類する
    clf = RandomForestClassifier(n_estimators=100,
                                 random_state=42)
    z_pred_proba = cross_val_predict(clf, X, z,
                                     cv=5, method='predict_proba')

    # 先に生成したデータに対応する部分だけ取り出す
    z_train_pred_proba = z_pred_proba[:2500]
    # 後から生成したデータへの近さにもとづいてソートする
    X_train_alike = X_train[np.argsort(z_train_pred_proba[:, 0])]
    y_train_alike = y_train[np.argsort(z_train_pred_proba[:, 0])]
    # 近いものをラベルごとに取り出す
    X_train_alike_neg = X_train_alike[y_train_alike == 0][:100]
    y_train_alike_neg = y_train_alike[y_train_alike == 0][:100]
    X_train_alike_pos = X_train_alike[y_train_alike == 1][:100]
    y_train_alike_pos = y_train_alike[y_train_alike == 1][:100]

    X_train_alike_concat = np.concatenate([X_train_alike_neg,
                                           X_train_alike_pos],
                                          axis=0)
    y_train_alike_concat = np.concatenate([y_train_alike_neg,
                                           y_train_alike_pos],
                                          axis=0)

    # 似ているデータで学習する
    clf.fit(X_train_alike_concat, y_train_alike_concat)
    # 後から作ったデータで評価する
    y_pred = clf.predict(X_test)
    # 精度 (Accuracy) を確認する
    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)

    # 識別境界を可視化する
    x_mesh, y_mesh = np.meshgrid(np.linspace(-5, 5, 50),
                                 np.linspace(-5, 5, 50))
    Z = clf.predict(np.c_[x_mesh.ravel(),
                          y_mesh.ravel()])
    Z = Z.reshape(x_mesh.shape)
    plt.contourf(x_mesh, y_mesh, Z, cmap=plt.cm.binary)

    plt.scatter(X_train[y_train == 0, 0],
                X_train[y_train == 0, 1],
                alpha=0.5,
                label='Train (Negative)')
    plt.scatter(X_train[y_train == 1, 0],
                X_train[y_train == 1, 1],
                alpha=0.5,
                label='Train (Positive)')
    plt.scatter(X_test[y_test == 0, 0],
                X_test[y_test == 0, 1],
                alpha=0.5,
                label='Test (Negative)')
    plt.scatter(X_test[y_test == 1, 0],
                X_test[y_test == 1, 1],
                alpha=0.5,
                label='Test (Positive)')
    plt.scatter(X_train_alike_concat[:, 0],
                X_train_alike_concat[:, 1],
                label='Model Learned')
    plt.legend()
    plt.show()


if __name__ == '__main__':
    main()

上記の実行結果は次の通り。 先ほどよりもスコアがさらに改善している。

$ python perlabel.py 
Accuracy: 0.9812

取り出したデータと識別境界をプロットしたグラフは次の通り。

f:id:momijiame:20190223164105p:plain

どれだけ取り出して学習させれば良いのか

ここまでで、検証用データに似たデータを学習用データから取り出して使うとスコアが上昇する可能性があることが分かった。 とはいえ、まだ疑問が残っている。 一体、どれだけ取り出して学習させると良い結果が得られるのだろうか。 先に断っておくと、まだ自分の中で結論は出ていない。 先ほど用いた上位 200 件というのは、特に根拠のない適当に選んだ数値だった。

試しに次のサンプルコードでは取り出す件数と精度の推移をプロットしている。 また、そのとき選ばれた末尾のデータは Adversarial Validation でどれくらい「検証用データっぽい」と判断されたのかも可視化した。

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

import numpy as np
from matplotlib import pyplot as plt
from sklearn.datasets import make_classification
from sklearn.metrics import accuracy_score
from sklearn.model_selection import cross_val_predict
from sklearn.ensemble import RandomForestClassifier


def main():
    args = {
        'n_samples': 2500,
        'n_features': 2,
        'n_informative': 2,
        'n_redundant': 0,
        'n_repeated': 0,
        'n_classes': 2,
        'n_clusters_per_class': 1,
        'flip_y': 0,
        'class_sep': 1.5,
        'weights': [0.5, 0.5],
        'random_state': 42,
    }
    X_train, y_train = make_classification(**args)

    args['random_state'] = 424242
    X_test, y_test = make_classification(**args)

    z_train, z_test = np.zeros(len(X_train)), np.ones(len(X_test))

    X = np.concatenate([X_train, X_test], axis=0)
    z = np.concatenate([z_train, z_test], axis=0)

    # 要素が最初のデータセット由来なのか、次のデータセット由来なのか分類する
    clf = RandomForestClassifier(n_estimators=100,
                                 random_state=42)
    z_pred_proba = cross_val_predict(clf, X, z,
                                     cv=5, method='predict_proba')

    # 先に生成したデータに対応する部分だけ取り出す
    z_train_pred_proba = z_pred_proba[:2500]
    # 後から生成したデータへの近さにもとづいてソートする
    sorted_train_index_test_similarity = np.argsort(z_train_pred_proba[:, 0])
    X_train_alike = X_train[sorted_train_index_test_similarity]
    y_train_alike = y_train[sorted_train_index_test_similarity]
    sorted_train_proba_test_similarity = np.sort(z_train_pred_proba[:, 0])

    accuracies = []
    probabilities = []
    picks = list(range(100, 2500, 100))
    for i in picks:
        # 近いものをラベルごとに取り出す
        X_train_alike_neg = X_train_alike[y_train_alike == 0][:i]
        y_train_alike_neg = y_train_alike[y_train_alike == 0][:i]
        X_train_alike_pos = X_train_alike[y_train_alike == 1][:i]
        y_train_alike_pos = y_train_alike[y_train_alike == 1][:i]

        X_train_alike_concat = np.concatenate([X_train_alike_neg,
                                               X_train_alike_pos],
                                              axis=0)
        y_train_alike_concat = np.concatenate([y_train_alike_neg,
                                               y_train_alike_pos],
                                              axis=0)

        # 似ているデータで学習する
        clf.fit(X_train_alike_concat, y_train_alike_concat)
        # 後から作ったデータで評価する
        y_pred = clf.predict(X_test)
        # 精度 (Accuracy) を確認する
        acc = accuracy_score(y_test, y_pred)
        accuracies.append(acc)
        # 取り出したデータの検証用データっぽさ (確率) を確認する
        proba = 1.0 - sorted_train_proba_test_similarity[i]
        probabilities.append(proba)

    # 結果をプロットする
    _, ax1 = plt.subplots(figsize=(8, 4))
    ax1.plot(picks, accuracies, c='b')
    ax1.set_xlabel('picks')
    ax1.set_ylabel('accuracy')

    ax2 = ax1.twinx()
    ax2.plot(picks, probabilities, c='r')
    ax2.set_xlabel('picks')
    ax2.set_ylabel('probability')

    plt.show()


if __name__ == '__main__':
    main()

なお、この精度の推移は Private LB の結果と同じ意味を持っているので、コンペ開催中にこのような探索をすることはできない。 代わりに、Public LB であれば Local CV との相関の推移を比較することはできると思う。

上記を実行してみる。

$ python pickacc.py

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

f:id:momijiame:20190223180541p:plain

どうやら、今回のデータでは上位 1,200 件を取り出すあたりから精度が大きく改善していることが分かる。 そして、結局上位 100 ~ 200 件を取り出すのが、最も良い結果が得られていた。 今回のケースでは、学習用データを大きく捨ててでも、より似ているデータを取り出す方が良い結果になるようだ。

ただし、上記は用いるデータによって結果は異なると考えられる。 おそらく Adversarial Validation でどれだけ高いスコアが得られたか、つまり分布が似ているか否かで取り出す量を変化させるのが良い気がしている。

Adversarial Validation の応用例

主にデータ分析コンペで用いられる Adversarial Validation だけど、おそらく現実世界でも応用が効くところはあると考えている。

例えば、運用中の機械学習モデルの再学習をするタイミングの判断に使えたりするんじゃないだろうか。 このアイデアでは、モデルの学習に使ったデータと、リアルタイムのデータを Adversarial Validation する。 モデルの学習に使ったデータとリアルタイムのデータの分布が近ければ、分類はできないはず。 しかし、分布が異なってくれば分類できるようになる。 この場合、運用中の機械学習モデルは未知の分布を処理していることになるので性能が低下している恐れがある。

このように、本来の目的変数とは異なる目的変数をデータに付与して分類させるという Adversarial Validation の考え方は、色々と応用が効きそうな気がしている。

いじょう。

Python: scikit-learn の cross_validate() 関数で独自の評価指標を計算する

今回は scikit-learncross_validate() 関数で、組み込みでは用意されていないような評価指標を計算する方法について書く。

使った環境は次の通り。

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

下準備

下準備として、事前に scikit-learn をインストールしておく。

$ pip install scikit-learn

cross_validate() 関数について

まずは cross_validate() 関数のおさらいから。 この関数は、その名の通り機械学習モデルの汎化性能を Cross Validation (交差検証) で評価するためにある。

次のサンプルコードでは、Breast Cancer データセットをランダムフォレストで学習させた場合の汎化性能を Stratified k-Fold CV で評価している。 cross_validate() 関数は、一度に複数の評価指標を計算できる点も特徴的。 以下では、組み込みで用意されている Accuracy, Precision, Recall, F1-Value について計算している。

#!/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


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

    # モデルを用意する
    clf = RandomForestClassifier(n_estimators=100, random_state=42)

    # Stratified k-Fold で汎化性能を評価する
    skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    # 評価する指標
    score_funcs = [
        'accuracy',
        'precision',
        'recall',
        'f1',
    ]
    # Cross Validation で検証する
    scores = cross_validate(clf, X, y, cv=skf, scoring=score_funcs)
    # 得られた指標を出力する
    print('accuracy:', scores['test_accuracy'].mean())
    print('precision:', scores['test_precision'].mean())
    print('recall:', scores['test_recall'].mean())
    print('f1:', scores['test_f1'].mean())


if __name__ == '__main__':
    main()

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

$ python cv.py 
accuracy: 0.9631088880338592
precision: 0.961883227291678
recall: 0.9805164319248826
f1: 0.9708395567654284

独自の評価指標を組み込む

続いて本題の独自の評価指標を組み込む方法について。 独自の評価指標を計算したいときは、真の値と推定した値を引数として受け取る関数を定義する。 そして、関数では引数を元に評価指標を計算して返す。

以下のサンプルコードでは、手始めに精度 (Accuracy) を計算する関数 accuracy_score を定義している。 精度 (Accuracy) の計算は組み込みでも用意されているけど、まずはそれと差異が出ないことを比較できるようにするため。 ポイントとして、定義した関数は make_scorer() 関数でラップする必要がある。

#!/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.metrics import make_scorer


def accuracy_score(y_true, y_pred):
    """オリジナルの評価関数 (ただし単なる精度の計算)"""
    accuracy = sum(y_true == y_pred) / len(y_true)
    return accuracy


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

    clf = RandomForestClassifier(n_estimators=100, random_state=42)

    skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

    score_funcs = {
        'builtin_accuracy': 'accuracy',
        # オリジナルの評価関数を登録する
        'custom_accuracy': make_scorer(accuracy_score),  # make_scorer() でラップする
    }

    scores = cross_validate(clf, X, y, cv=skf, scoring=score_funcs)
    print('built-in accuracy:', scores['test_builtin_accuracy'].mean())
    print('custom accuracy:', scores['test_custom_accuracy'].mean())


if __name__ == '__main__':
    main()

上記の実行結果は次の通り。 組み込みの精度 (Accuracy) 計算と同じ値が得られていることが分かる。

$ python acc.py 
built-in accuracy: 0.9631088880338592
custom accuracy: 0.9631088880338592

回帰問題でもやってみる

先ほどは分類問題だったので、念のため次は回帰問題でも確認しておこう。

以下のサンプルコードでは Boston データセットを ElasticNet 回帰で処理している。 独自の評価指標を計算する関数は rmse_score() という関数を用意した。 これは、文字通り RMSE (Root Mean Square Error: 平均二乗誤差平方根) を計算する関数になっている。

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

import math

from sklearn import datasets
from sklearn.linear_model import ElasticNet
from sklearn.model_selection import KFold
from sklearn.model_selection import cross_validate
from sklearn.metrics import mean_squared_error
from sklearn.metrics import make_scorer


def rmse_score(y_true, y_pred):
    """RMSE (Root Mean Square Error: 平均二乗誤差平方根) を計算する関数"""
    mse = mean_squared_error(y_true, y_pred)
    rmse = math.sqrt(mse)
    return rmse


def main():
    # Boston データセットを使った回帰問題
    dataset = datasets.load_boston()
    X, y = dataset.data, dataset.target

    # ElasticNet 回帰
    reg = ElasticNet(random_state=42)

    kf = KFold(n_splits=5, shuffle=True, random_state=42)

    score_funcs = {
        'rmse': make_scorer(rmse_score),
    }

    scores = cross_validate(reg, X, y, cv=kf, scoring=score_funcs)
    mean_rmse = scores['test_rmse'].mean()
    print('RMSE:', mean_rmse)


if __name__ == '__main__':
    main()

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

$ python rmse.py   
RMSE: 5.274746876673534

ちゃんと計算できているようだ。

めでたしめでたし。

機械学習のための特徴量エンジニアリング ―その原理とPythonによる実践

機械学習のための特徴量エンジニアリング ―その原理とPythonによる実践

  • 作者: Alice Zheng,Amanda Casari,株式会社ホクソエム
  • 出版社/メーカー: オライリージャパン
  • 発売日: 2019/02/23
  • メディア: 単行本(ソフトカバー)
  • この商品を含むブログを見る

Python: CatBoost を使ってみる

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

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

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

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

f:id:momijiame:20190216194828p:plain

上記を見ると 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

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

f:id:momijiame:20190216201841p:plain

上記から 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: 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()

f:id:momijiame:20190210021002p:plain

上記のデータから、今回作ったライブラリを使って特徴量を抽出する。 抽出した特徴量は、近傍数 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()

f:id:momijiame:20190210021013p:plain

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

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 次元のデータになっている。

f:id:momijiame:20190210021730p:plain

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

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

いじょう。