CUBE SUGAR CONTAINER

技術系のこと書きます。

Python: 特徴量の重要度を Permutation Importance で計測する

学習させた機械学習モデルにおいて、どの特徴量がどれくらい性能に寄与しているのかを知りたい場合がある。 すごく効く特徴があれば、それについてもっと深掘りしたいし、あるいは全く効かないものがあるなら取り除くことも考えられる。 使うフレームワークやモデルによっては特徴量の重要度を確認するための API が用意されていることもあるけど、そんなに多くはない。 そこで、今回はモデルやフレームワークに依存しない特徴量の重要度を計測する手法として Permutation Importance という手法を試してみる。 略称として PIMP と呼ばれたりすることもあるようだ。

この手法を知ったのは、以下の Kaggle のノートブックを目にしたのがきっかけだった。

Permutation Importance | Kaggle

あんまりちゃんと読めてないけど、論文としては Altmann et al. (2010) になるのかな?

Permutation importance: a corrected feature importance measure | Bioinformatics | Oxford Academic

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.13.6
BuildVersion:   17G3025
$ python -V     
Python 3.7.1

Permutation Importance の考え方について

Permutation Importance の考え方はとてもシンプル。 ある特徴量が使い物にならないとき、モデルの性能がどれだけ落ちるかで重要度を計測する。 もし性能が大きく落ちるなら、その特徴量はモデルにおいて重要だと考えられるし、その逆もまたしかり。 下手をすると性能がむしろ上がって、ない方がマシということが分かったりするかもしれない。

この説明だけだと頭の中にクエスチョンマークがたくさん浮かぶと思うので、順を追って説明していく。 Permutation Importance では、まずデータセットを学習用のデータと検証用のデータに分割する。 その上で、学習用のデータを使ってモデルを学習させる。 そして、何も手を加えていない検証用のデータを使ってまずは性能を計測する。 性能の計測には、解決したい問題に対して任意の適切な評価指標を用いる。 ここまでの一連の流れは、一般的な機械学習モデルの汎化性能の計測方法と何も変わらない。 このときに測った性能を、後ほどベースラインとして用いる。

ベースラインが計測できたところで、続いては検証用データについて特定の特徴量だけを使い物にならない状態にする。 これには、手法の名前にもある通り "Permutation (交換・置換)" を用いる。 具体的には、特徴量をランダムにシャッフルすることで、目的変数に対して特徴量の相関が取れない状態にする。 特定の特徴量がシャッフルされた検証用データが完成したら、それを学習済みのモデルに与えて性能を計測する。 このとき、どれだけベースラインから性能が変化するかを確認する。

上記を全ての特徴量に対して実施して、ベースラインに対する性能の変化を特徴量ごとに比較する。 以上が Permutation Importance の計測方法になる。 ただ、文章だけだとちょっと分かりにくいかもと思ったので図示してみる。 最初のベースラインを測る部分は一般的な機械学習モデルの汎化性能を計測するやり方と同じ。

f:id:momijiame:20181118112928p:plain

ベースラインが計測できたら、続いて検証用データの特定の特徴量を Permutation で使い物にならない状態にする。 そして、その検証用データを使って性能を改めて検証する。 このとき、ベースラインからの性能の変化を観測する。

f:id:momijiame:20181118112944p:plain

上記において勘違いしやすいポイントは Permutation Importance でモデルの再学習は必要ないという点。 なんとなく大元のデータセットをシャッフルした上でモデルを再学習させて...と考えてしまうんだけど、実は必要ない。 モデルを再学習すると、それに計算コストも必要になるし、ベースラインで使ったモデルとの差異も大きくなってしまう。 そこで、学習済みのモデルと検証用データだけを使って性能を計測する。

実際に計測してみる

続いては実際に Permutation Importance を確認してみよう。

下準備として、使うパッケージをあらかじめインストールしておく。

$ pip install pandas scikit-learn matplotlib

今回はデータセットにみんな大好き Iris データセットを、モデルには Random Forest を選んだ。 結果がなるべく安定するように 5-fold CV でそれぞれ計算した結果を出している。 以下は手作業で計算処理を書いてるけど eli5 というパッケージに PIMP を計算するためのモジュールがあったりもする。

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

from collections import defaultdict

import numpy as np
import pandas as pd
from sklearn import datasets
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score
from sklearn.model_selection import StratifiedKFold
from matplotlib import pyplot as plt


def permuted(df):
    """特定のカラムをシャッフルしたデータフレームを返す"""
    for column_name in df.columns:
        permuted_df = df.copy()
        permuted_df[column_name] = np.random.permutation(permuted_df[column_name])
        yield column_name, permuted_df


def pimp(clf, X, y, cv=None, eval_func=accuracy_score):
    """PIMP (Permutation IMPortance) を計算する"""
    base_scores = []
    permuted_scores = defaultdict(list)

    if cv is None:
        cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

    for train_index, test_index in cv.split(X, y):
        # 学習用データと検証用データに分割する
        X_train, y_train = X.iloc[train_index], y.iloc[train_index]
        X_test, y_test = X.iloc[test_index], y.iloc[test_index]

        # 学習用データでモデルを学習する
        clf.fit(X_train, y_train)

        # まずは何もシャッフルしていないときのスコアを計算する
        y_pred_base = clf.predict(X_test)
        base_score = eval_func(y_test, y_pred_base)
        base_scores.append(base_score)

        # 特定のカラムをシャッフルした状態で推論したときのスコアを計算する
        permuted_X_test_gen = permuted(X_test)
        for column_name, permuted_X_test in permuted_X_test_gen:
            y_pred_permuted = clf.predict(permuted_X_test)
            permuted_score = eval_func(y_test, y_pred_permuted)
            permuted_scores[column_name].append(permuted_score)

    # 基本のスコアとシャッフルしたときのスコアを返す
    np_base_score = np.array(base_scores)
    dict_permuted_score = {name: np.array(scores) for name, scores in permuted_scores.items()}
    return np_base_score, dict_permuted_score


def score_difference_statistics(base, permuted):
    """シャッフルしたときのスコアに関する統計量 (平均・標準偏差) を返す"""
    mean_base_score = base.mean()
    for column_name, scores in permuted.items():
        score_differences = scores - mean_base_score
        yield column_name, score_differences.mean(), score_differences.std()


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

    # 計測に使うモデルを用意する
    clf = RandomForestClassifier(n_estimators=100)

    # Permutation Importance を計測する
    base_score, permuted_scores = pimp(clf, X, y)

    # 計測結果から統計量を計算する
    diff_stats = list(score_difference_statistics(base_score, permuted_scores))

    # カラム名、ベーススコアとの差、95% 信頼区間を取り出す
    sorted_diff_stats = sorted(diff_stats, key=lambda x: x[1])
    column_names = [name for name, _, _ in sorted_diff_stats]
    diff_means = [diff_mean for _, diff_mean, _ in sorted_diff_stats]
    diff_stds_95 = [diff_std * 1.96 for _, _, diff_std in sorted_diff_stats]

    # グラフにプロットする
    plt.plot(column_names, diff_means, marker='o', color='r')
    plt.errorbar(column_names, diff_means, yerr=diff_stds_95, ecolor='g', capsize=4)

    plt.title('Permutation Importance')
    plt.grid()
    plt.xlabel('column')
    plt.ylabel('difference')
    plt.show()


if __name__ == '__main__':
    main()

上記に適当な名前をつけて実行してみよう。

$ python pimp.py

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

f:id:momijiame:20181110230620p:plain

上記はベースラインに対する各特徴量の性能の平均的な落ち具合を可視化している。 Sepal length や Sepal width に比べると Petal length や Petal width の方が性能が大きく落ちていることが分かる。 ここから、このモデルでは Petal length や Petal width の方が識別に有効に使われていることが推測できる。

緑色のヒゲは、性能の落ち具合のバラつきが正規分布になるという仮定の元での 95% 信頼区間を示している。 実際に正規分布になるかは確認していないし、モデルにもよるだろうけど分散を確認するために描いてみた。

ちなみに、実際に色んなデータやモデルで計測してみると、意外と毎回結果が変わったりすることがある。 これは、そもそもモデルが決定論的に学習しないことだったり、試行回数が少ないことが原因として考えられる。 例えば計測する n-fold CV を何回も繰り返すなど、試行回数を増やすとより結果が安定しやすいかもしれない。

いじょう。

OpenSSH のコネクションが切れにくいように Keepalive を送る

たまに SSH のコネクションが頻繁に切れる環境があるので、定期的にデータを送受信することで切断されるのを防ぎたい。 これは OpenSSH のクライアントであれば ServerAliveInterval を設定することで実現できる。

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.13.6
BuildVersion:   17G3025
$ ssh -V
OpenSSH_7.6p1, LibreSSL 2.6.2

コマンドラインからであれば、次のように -o オプションで ServerAliveInterval に秒数を指定しながらリモートに接続する。

$ ssh -o ServerAliveInterval=60 <username>@<remotehost>

毎回コマンドラインで指定するのは面倒なので OpenSSH のコンフィグファイル (典型的には ~/.ssh/config) に記述しても良い。

Host *
  ServerAliveInterval 60

このオプションの詳細については ssh_config のセクション 5 に記載がある。

$ man 5 ssh_config

内容を以下に引用する。 この値が設定されていると、一定の期間にサーバからデータを受信しなかったときは SSH のコネクションを通してメッセージを送る、とある。

     ServerAliveInterval
             Sets a timeout interval in seconds after which if no data has been
             received from the server, ssh(1) will send a message through the
             encrypted channel to request a response from the server.  The default is
             0, indicating that these messages will not be sent to the server.

ちゃんとサーバまでメッセージが届けば、それに対する返答がくる。 それによってコネクションが維持されるという寸法。 ちなみに、メッセージを送信する回数に ServerAliveCountMax で上限を設けることもできるようだ。

いじょう。

OpenSSH[実践]入門 (Software Design plus)

OpenSSH[実践]入門 (Software Design plus)

Python: 実行環境が Jupyter Notebook か判定する

今回は Python の実行環境が Jupyter Notebook か否かを判定する方法について。 ベースのアイデアは以下の Stackoverflow からお借りした。

stackoverflow.com

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.13.6
BuildVersion:   17G3025
$ python -V
Python 3.6.7
$ pip list --format=columns | grep notebook
notebook          5.7.0

下準備

まずは今回使用するパッケージをインストールしておく。

$ pip install notebook ipywidgets tqdm

実行環境が Jupyter Notebook か否かをチェックする関数

早速だけど実行環境をチェックする関数は次の通り。

def is_env_notebook():
    """Determine wheather is the environment Jupyter Notebook"""
    if 'get_ipython' not in globals():
        # Python shell
        return False
    env_name = get_ipython().__class__.__name__
    if env_name == 'TerminalInteractiveShell':
        # IPython shell
        return False
    # Jupyter Notebook
    return True

この関数を使うと、例えば次のように環境ごとにインポートするパッケージや関数を切り替えることができる。

if is_env_notebook():
    # Jupyter Notebook environment
    from tqdm import tqdm_notebook as tqdm
else:
    # Other shells
    from tqdm import tdqm

関数の動作原理について

ここからは関数がどのように動作するかを確認していく。

まずは Jupyter Notebook を起動しておこう。

$ jupyter-notebook

新しい Notebook を作ったら、以下のコードを使ってグローバルスコープのオブジェクトを全て表示してみよう。

import pprint
pprint.pprint(globals())

すると、例えば次のような出力になる。 この中に get_ipython() という関数があることが確認できる。

{'In': ['', 'import pprint\npprint.pprint(globals())'],
 'Out': {},
 '_': '',
 '__': '',
 '___': '',
 '__builtin__': <module 'builtins' (built-in)>,
 '__builtins__': <module 'builtins' (built-in)>,
 '__doc__': 'Automatically created module for IPython interactive environment',
 '__loader__': None,
 '__name__': '__main__',
 '__package__': None,
 '__spec__': None,
 '_dh': ['/Users/amedama/Documents/temporary'],
 '_i': '',
 '_i1': 'import pprint\npprint.pprint(globals())',
 '_ih': ['', 'import pprint\npprint.pprint(globals())'],
 '_ii': '',
 '_iii': '',
 '_oh': {},
 'exit': <IPython.core.autocall.ZMQExitAutocall object at 0x104733c18>,
 'get_ipython': <bound method InteractiveShell.get_ipython of <ipykernel.zmqshell.ZMQInteractiveShell object at 0x104733dd8>>,
 'pprint': <module 'pprint' from '/Users/amedama/.pyenv/versions/3.6.7/lib/python3.6/pprint.py'>,
 'quit': <IPython.core.autocall.ZMQExitAutocall object at 0x104733c18>}

では、次に別のターミナルから通常の Python インタプリタを起動しよう。

$ python

先ほどと同じコードを実行してみる。 すると、こちらには get_ipython() 関数が見当たらない。 その他にも、いくつか Jupyter Notebook の環境ではあったオブジェクトが無いことも分かる。

>>> from pprint import pprint
>>> pprint(globals())
{'__annotations__': {},
 '__builtins__': <module 'builtins' (built-in)>,
 '__doc__': None,
 '__loader__': <class '_frozen_importlib.BuiltinImporter'>,
 '__name__': '__main__',
 '__package__': None,
 '__spec__': None,
 'pprint': <function pprint at 0x10b08de18>}

この特徴を元に、環境を一次切り分けできる。 少なくとも get_ipython() 関数がグローバルスコープになければ Jupyter Notebook の環境ではなさそう。

$ python

...

>>> 'get_ipython' in globals()
False

ただし、これだけだと上手くいかない場合がある。 それは IPython の環境で、IPython のインタプリタを起動するとグローバルスコープに get_ipython() 関数が見つかる。

$ ipython

...

In [1]: 'get_ipython' in globals()                                                                                                          
Out[1]: True

ただし、一つ救いがあって get_ipython() 関数の返り値が Jupyter Notebook とは異なっている。 IPython のインタプリタで返されるのは TerminalInteractiveShell というクラスのインスタンス。

In [2]: get_ipython()                                                                                                                       
Out[2]: <IPython.terminal.interactiveshell.TerminalInteractiveShell at 0x1080a7710>

それに対し Jupyter Notebook では ZMQInteractiveShell というクラスのインスタンスになっている。

$ jupyter-notebook

...

In [1]: 'get_ipython' in globals()
Out[1]: True

In [2]: get_ipython()
Out[2]: <ipykernel.zmqshell.ZMQInteractiveShell at 0x1081f0198>

そのため、クラスの違いを元に環境を識別できる。

いじょう。

Linux の UNIX パスワード認証について調べた

ふと Linux ディストリビューションのユーザ認証周りについて気になって、その中でも特に UNIX パスワード認証について調べてみた。 UNIX パスワード認証というのは、Linux に限らず Unix 系のディストリビューションで広く採用されているパスワードを使った認証の仕組み。 特に、ログイン用のパスワードを暗号化 (ハッシュ化) した上でパスワードファイル (/etc/passwd) やシャドウパスワードファイル (/etc/shadow) に保存するところが特徴となっている。 UNIX パスワード認証は Unix 系の色々なディストリビューションでそれぞれ実装されている。 その中でも、今回は Linux ディストリビューションの Ubuntu 18.04 LTS について見ていく。 先に断っておくと、かなり長い。

注: この記事では、やっていることがハッシュ化でも manpage の "encrypted" や "encryption" という表現を尊重して「暗号化」と記載している箇所がある

動作確認に使った環境は次の通り。

$ cat /etc/lsb-release
DISTRIB_ID=Ubuntu
DISTRIB_RELEASE=18.04
DISTRIB_CODENAME=bionic
DISTRIB_DESCRIPTION="Ubuntu 18.04.1 LTS"
$ uname -r
4.15.0-29-generic

パスワードファイルとシャドウパスワードファイルについて

太古の昔、ユーザのログインに使うパスワードはパスワードファイル (/etc/passwd) に記述されていたらしい。 しかし、今ではここにパスワードに関する情報が書かれることはほとんどない。 過去に暗号化されたパスワードが記載されていたフィールドには、代わりに x といった文字が載っている。

$ cat /etc/passwd | grep $(whoami)
vagrant:x:1000:1000:vagrant,,,:/home/vagrant:/bin/bash

なぜなら、このファイルの内容は誰でも読めるため。 仮にパスワードが暗号化 (ハッシュ化) されていたとしても、誰からでも読めるのはセキュリティ的には致命的といえる。 そこで、ログインに必要な情報はシャドウパスワードファイル (例えば /etc/shadow) に分離された。

$ ls -l /etc/passwd
-rw-r--r-- 1 root root 1662 Aug 24 09:07 /etc/passwd
$ ls -l /etc/shadow
-rw-r----- 1 root shadow 979 Aug 24 09:07 /etc/shadow

上記を見て分かる通り /etc/shadowroot ユーザか shadow グループのメンバー以外からは閲覧できない。

シャドウパスワードファイルにおいて、カンマで区切られた二つ目のカラムがユーザがログインに用いるパスワードを暗号化した文字列になる。 これは、ソルト付きのハッシュ化されたパスワードになっている。

$ sudo cat /etc/shadow | grep $(whoami)
vagrant:$6$xUMKs3vr$r/rGFiHf6.B0xsKDpustA2G9m1zGcWJnjZLMlyQeh.VDEpDEMEQ3AKSvrRv.NQEc2OBbban3BAHj4Nhwpa0t4.:17767:0:99999:7:::

ユーザがシステムにログインするときは、ユーザが入力した内容と上記の内容を突合して一致するかを確認することになる。

シャドウパスワードファイルの詳細については man コマンドで調べることができる。 man のセクション 5 には設定ファイルに関する情報が記載されている。

$ man 5 shadow

上記の内容を確認していくと、次のような記述がある。

       encrypted password
           Refer to crypt(3) for details on how this string is interpreted.

どうやらファイルに記載されている暗号化されたパスワードは crypt(3) に関連しているようだ。 man のセクション 3 は開発者向けのライブラリ関数に関する情報なので、それ用の manpage をインストールした上で確認する。

$ sudo apt-get -y install manpages-dev
$ man 3 crypt

上記を見ると暗号化されたパスワードは次のようなフォーマットになっていることが分かる。 id がハッシュアルゴリズムで salt がハッシュするときに使うソルト、そして最後の encrypted がハッシュ化されたパスワードになる。

$id$salt$encrypted

id については次のような記述がある。 使えるハッシュアルゴリズムは環境や glibc のバージョンによって異なるようだ。 glibc というのは C 言語の標準ライブラリである libc の、GNU Project による実装を指している。

       id identifies the encryption method used instead of DES and  this  then
       determines  how  the  rest  of the password string is interpreted.  The
       following values of id are supported:

              ID  | Method
              ─────────────────────────────────────────────────────────
              1   | MD5
              2a  | Blowfish (not in mainline glibc; added in some
                  | Linux distributions)
              5   | SHA-256 (since glibc 2.7)
              6   | SHA-512 (since glibc 2.7)

C 言語の標準ライブラリというのは、例えば多くの入門書で「おまじない」として紹介されている #include <stdio.h> なんかがそれ。

先ほど確認した自分自身のパスワードハッシュでは id6 になっていた。 つまり、ハッシュアルゴリズムとして SHA-512 が使われていることが分かる。

$ sudo cat /etc/shadow | grep $(whoami)
vagrant:$6$xUMKs3vr$r/rGFiHf6.B0xsKDpustA2G9m1zGcWJnjZLMlyQeh.VDEpDEMEQ3AKSvrRv.NQEc2OBbban3BAHj4Nhwpa0t4.:17767:0:99999:7:::

crypt(3) を使ってみよう

それでは、実際に crypt(3) を使ってみることにしよう。 ただ、C 言語で書くのではなく Python のラッパーである crypt モジュールを使うことにした。

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

$ python3

先ほど確認した自分自身の暗号化されたパスワードを再掲する。 ハッシュアルゴリズムは SHA-512 でソルトが xUMKs3vr だった。 そして、使われているパスワードは Vagrant の環境なので vagrant になっている。

$ sudo cat /etc/shadow | grep $(whoami)
vagrant:$6$xUMKs3vr$r/rGFiHf6.B0xsKDpustA2G9m1zGcWJnjZLMlyQeh.VDEpDEMEQ3AKSvrRv.NQEc2OBbban3BAHj4Nhwpa0t4.:17767:0:99999:7:::

上記の内容にもとづいて crypt モジュールの関数を呼び出してみよう。 ソルトについては暗号化されたパスワードの $ で区切られた部分をそのまま使って構わない。

>>> import crypt
>>> crypt.crypt('vagrant', '$6$xUMKs3vr$')
'$6$xUMKs3vr$r/rGFiHf6.B0xsKDpustA2G9m1zGcWJnjZLMlyQeh.VDEpDEMEQ3AKSvrRv.NQEc2OBbban3BAHj4Nhwpa0t4.'

見事に自分自身の暗号化されたパスワードと同じ文字列が生成できた。

ちなみに、システムで使えるハッシュアルゴリズムは crypt.methods で確認できる。

>>> crypt.methods
[<crypt.METHOD_SHA512>, <crypt.METHOD_SHA256>, <crypt.METHOD_MD5>, <crypt.METHOD_CRYPT>]

上記を見て分かる通り id2a に対応するハッシュアルゴリズムの Blowfish はこのシステムでは使えないらしい。

crypt(3) の実装を見てみよう

次は、念のため crypt(3) を実装している場所についても調べておくことにしよう。 特に興味がなければ、次のセクションまで飛ばしてもらって構わない。

先ほど見た通り、実装しているのは glibc と思われる。 そこで、まずは今使っているシステムの glibc のバージョンを調べよう。 検索パスに入っている共有ライブラリを ldconfig コマンドで調べたら、そこから libc を探す。

$ ldconfig -p | grep libc.so
    libc.so.6 (libc6,x86-64, OS ABI: Linux 3.2.0) => /lib/x86_64-linux-gnu/libc.so.6

あとは共有ライブラリをそのまま実行するとバージョンが確認できる。 どうやら、この環境ではバージョン 2.27 が使われているようだ。 先ほど crypt(3) のハッシュアルゴリズムとして SHA-512 が使えるのはバージョン 2.7 以降とあったので、その内容とも整合している。

$ /lib/x86_64-linux-gnu/libc.so.6
GNU C Library (Ubuntu GLIBC 2.27-3ubuntu1) stable release version 2.27.
Copyright (C) 2018 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.
There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE.
Compiled by GNU CC version 7.3.0.
libc ABIs: UNIQUE IFUNC
For bug reporting instructions, please see:
<https://bugs.launchpad.net/ubuntu/+source/glibc/+bugs>.

続いて公式サイトの記述に沿ってソースコードを取得する。

$ sudo apt-get -y install git
$ git clone git://sourceware.org/git/glibc.git
$ cd glibc/
$ git checkout release/2.27/master
Branch 'release/2.27/master' set up to track remote branch 'release/2.27/master' from 'origin'.
Switched to a new branch 'release/2.27/master'

とりあえず crypt(3) で検索してみると、それっぽいものが見つかる。

$ grep -r "crypt(3)" .
./crypt/ufc.c: * UFC-crypt: ultra fast crypt(3) implementation
./crypt/crypt-entry.c: * UFC-crypt: ultra fast crypt(3) implementation
./crypt/cert.c: * This crypt(3) validation program shipped with UFC-crypt
./crypt/crypt.c: * UFC-crypt: ultra fast crypt(3) implementation
./crypt/crypt_util.c: * UFC-crypt: ultra fast crypt(3) implementation
./crypt/README.ufc-crypt:- Crypt implementation plugin compatible with crypt(3)/fcrypt.
./crypt/README.ufc-crypt:- On most machines, UFC-crypt runs 30-60 times faster than crypt(3) when
./crypt/README.ufc-crypt:  that of crypt(3).
./crypt/README.ufc-crypt:such as the one said to be used in crypt(3).
./crypt/crypt-private.h: * UFC-crypt: ultra fast crypt(3) implementation
./crypt/crypt.h: * UFC-crypt: ultra fast crypt(3) implementation

ヘッダファイルに crypt() 関数らしき宣言があったり、各ハッシュアルゴリズムの実装も見受けられる。 関数のエントリポイントは crypt-entry.c にあって、そこから各実装に分岐しているみたい。

$ grep extern crypt/crypt.h | head -n 1
extern char *crypt (const char *__key, const char *__salt)
$ ls crypt | egrep "(md5|sha256|sha512)\.c"
md5.c
sha256.c
sha512.c

たしかに crypt(3) は glibc に実装があることが分かった。

crypt(3) は何処から使われているのか?

続いて気になるのは、暗号化されたパスワードを生成する crypt(3) が何処から使われているのかという点。 おそらくユーザがログインするタイミングで入力されたパスワードに対して crypt(3) を使っているはず。 そして、使った結果とシャドウパスワードファイルの内容を比較して、一致すればログイン成功という処理が何処かにあると考えられる。

一般的なシステムであれば、ユーザがログインする経路として代表的なものが二つ考えられる。 まず一つ目が、コンピュータに直接接続された物理的なディスプレイやキーボードを通してユーザ名とパスワードを入力するパターン。 この場合、端末デバイスには制御端末 (tty: teletypewriter) が割り当てられる。 そして二つ目が、SSH クライアントを通してネットワーク越しにユーザ名やパスワード、あるいは公開鍵の情報を渡すパターン。 この場合、端末デバイスには擬似端末 (pty: pseudo terminal) が割り当てられる。

制御端末 (tty: teletypewriter) からログインするパターン

制御端末 (tty) を通してログインするパターンから見ていく。 このパターンでは、まず Linux ディストリビューションの初期化プロセスから考える必要がある。 Ubuntu 18.04 LTS では init プロセス (pid に 1 を持つプロセス) として systemd が用いられる。

そして systemd は agetty(8) というプログラムを起動する。 このプログラムはログインプロンプトの表示部分を担当している。 システムのプロセスを確認すると、制御端末が割り当てられた状態で agetty(8) が起動していることが分かる。

$ ps auxww | grep [a]getty
root       678  0.0  0.1  16180  1656 tty1     Ss+  16:00   0:00 /sbin/agetty -o -p -- \u --noclear tty1 linux
$ pstree | head -n 3
systemd-+-VBoxService---7*[{VBoxService}]
        |-accounts-daemon---2*[{accounts-daemon}]
        |-agetty

このプログラムに関する詳細は man ページのセクション 8 の agetty に記載されている。

$ man 8 agetty

上記の内容を見ると、このプログラムは /bin/login を呼び出すとある。 どうやら、このプログラムが実際のログイン部分を担当しているようだ。

DESCRIPTION
       agetty  opens  a  tty port, prompts for a login name and invokes the /bin/login command.  It is normally invoked by
       init(8).

ということで、次はこの login(1) のソースコードが読みたい。 なので、まずはこのコマンドがどのパッケージに含まれているのかを dpkg コマンドで確認する。

$ dpkg -S /bin/login
login: /bin/login

どうやら、そのものズバリ login というパッケージに含まれているらしい。

パッケージ名が分かったので apt-get source コマンドを使ってソースコードを取り寄せる。 すると、どうやら login(1) の実体は shadow というプログラム群に含まれることが分かった。

$ sudo apt-get -y install dpkg-dev
$ sudo sed -i.back -e "s/^# deb-src/deb-src/g" /etc/apt/sources.list
$ sudo apt-get update
$ apt-get source login
Reading package lists... Done
Picking 'shadow' as source package instead of 'login'
NOTICE: 'shadow' packaging is maintained in the 'Git' version control system at:
https://anonscm.debian.org/git/pkg-shadow/shadow.git
Please use:
git clone https://anonscm.debian.org/git/pkg-shadow/shadow.git
to retrieve the latest (possibly unreleased) updates to the package.
...

ドキュメントを見ると、このプログラム群が Linux のシャドウパスワードを扱うためのものだということが分かる。 どうやら核心に近づいてきたようだ。

$ cd shadow-4.5/
$ head doc/HOWTO -n 6 | tail -n 2
  Linux Shadow Password HOWTO
  Michael H. Jackson, mhjack@tscnet.com

中身を確認すると src というディレクトリにコマンドラインから叩くプログラムがまとまっていることが分かった。 例えば login(1) もここにあるし、他には usermod(8)groupmod(8) なんかも見つかる。

$ ls src/ | grep login
login.c
login_nopam.c
nologin.c
sulogin.c
$ ls src/ | grep usermod
usermod.c
$ ls src/ | grep groupmod
groupmod.c

(2018/11/5 追記: 以降の説明はシステムで PAM を用いない場合の説明になっているため、一般的な環境の説明と異なっている。例えば、後述する関数が login コマンドにおいて呼び出されるフローは、ビルドする際に USE_PAM フラグを無効にした場合にしか通らない。)

改めて src/login.c を読み進めていくと、ログインパスワードを pw_auth() という関数に突っ込んでいる。 この関数は lib/pwauth.c で定義されていて、さらに lib/encrypt.cpw_encrypt() という関数を呼び出している。 そして、この pw_encrypt() 関数の中に、探し求めていた crypt(3) の呼び出しが見つかる! ついに見つけることができた。

$ grep "crypt (" lib/encrypt.c 
/*@exposed@*//*@null@*/char *pw_encrypt (const char *clear, const char *salt)
    cp = crypt (clear, salt);

さらに src/login.c 周辺を読み進めると crypt(3) で暗号化したパスワードは getpwnam(3)getspnam(3) と比較していることが分かった。 この関数は /etc/passwd/etc/shadow を読み取って、そこに記載されている暗号化されたパスワードを読み取るためのものだ。 関数を提供しているのは、まさにこの shadow で、ソースコードの実体は lib 以下に存在する。

それぞれの関数の情報は man コマンドから参照できる。

$ man 3 getspnam

その後、認証に成功するとパスワードファイルから読み取ったログインシェルを起動することもソースコードから確認できる。

$ grep -r shell login.c
    if (pwd->pw_shell[0] == '*') { /* subsystem root */
        pwd->pw_shell++;   /* skip the '*' */
        err = shell (tmp, pwd->pw_shell, newenvp); /* fake shell */
        /* exec the shell finally */
        err = shell (pwd->pw_shell, (char *) 0, newenvp);

これで Ubuntu 18.04 LTS におけるログインの流れをソースコードから確認できた。

  • ユーザからの入力パスワードを暗号化する
  • 上記をシャドウパスワードファイルに記載されている暗号化されたパスワードと比較する
  • 一致するときはログインシェルを起動する

実際に制御端末 (tty) からログインした状態で pstree を確認すると agetty(8) がいなくなって、代わりに login(1) 経由で bash(1) が起動している。

$ pstree
systemd─┬─VBoxService───7*[{VBoxService}]
...
        ├─login───bash

ちなみにシステムの /sbin/agetty/bin/login をリネームすると、なかなか面白い挙動になるので試してみるのもおすすめ。 例えば agetty が見つからない状態でシステムが起動するとログインプロンプトが表示されないことを確認できる。

擬似端末 (pty: pseudo terminal) からログインするパターン

続いては SSH クライアントを通してシステムにログインするパターン。 これは単純で、ようするに SSH デーモンに認証部分があるだろうことは容易に想像がつく。

なので、まずは sshd(8) の含まれるパッケージを確認する。

$ dpkg -S $(which sshd)
openssh-server: /usr/sbin/sshd

先ほどと同じ要領で OpenSSH のソースコードを取り寄せる。

$ sudo apt-get source openssh-server
Reading package lists... Done
Picking 'openssh' as source package instead of 'openssh-server'
NOTICE: 'openssh' packaging is maintained in the 'Git' version control system at:
https://salsa.debian.org/ssh-team/openssh
Please use:
git clone https://salsa.debian.org/ssh-team/openssh
to retrieve the latest (possibly unreleased) updates to the package.
...
$ cd openssh-7.6p1/

この中では、例えば crypt(3)xcrypt() 関数というプラットフォーム非依存な形でラップしていたりする。

$ grep " crypt(" */*.c
openbsd-compat/xcrypt.c:        crypted = crypt(password, salt);
openbsd-compat/xcrypt.c:        crypted = crypt(password, salt);
openbsd-compat/xcrypt.c:    crypted = crypt(password, salt);

あるいは各所で getspnam(3)getpwnam(3) が使われている。

$ grep "getspnam" *.c
auth.c:     spw = getspnam(pw->pw_name);
auth-shadow.c:  if ((spw = getspnam((char *)user)) == NULL) {
$ grep "getspnam" */*.c
openbsd-compat/xcrypt.c:    struct spwd *spw = getspnam(pw->pw_name);

ということで、こちらの経路に関しても crypt(3)getspnam(3) を駆使してユーザを認証しているらしいことが分かった。 プロセスを確認しても、ユーザがログインした状態では sshd(8) 経由で bash(1) が起動していることが分かる。

$ pstree | grep sshd
        |-sshd---sshd---sshd---bash-+-grep

まとめ

分かったことは次の通り。

  • シャドウパスワードファイル (例えば /etc/shadow) には暗号化されたユーザのログインパスワードが書き込まれている
  • 暗号化されたパスワードは crypt(3) で生成する
  • シャドウパスワードファイルに書き込まれている情報は getspnam(3) で読み取る
  • ユーザの入力内容から生成したパスワードハッシュと、シャドウパスワードファイルに書き込まれた内容を比較することで認証する

いじょう。

ふつうのLinuxプログラミング 第2版 Linuxの仕組みから学べるgccプログラミングの王道

ふつうのLinuxプログラミング 第2版 Linuxの仕組みから学べるgccプログラミングの王道

詳解UNIXプログラミング 第3版

詳解UNIXプログラミング 第3版

IPv4 アドレスの値段

今回は IPv4 アドレスの値段や、その売買に関する動向について調べてみた。

TL; DR

  • IPv4 アドレスは売買できる
  • オークションにおける IPv4 アドレスの売買単価は値上がり傾向にある
  • 2018 年 11 月 1 日現在、IPv4 アドレス一つあたり 17 ~ 20 USD ほどで取引されている
  • 近年は取引も活発になっている傾向にある

IPv4 アドレスの在庫枯渇と売買について

IPv4 アドレスの在庫が枯渇したというニュースが話題になってから、ずいぶんと時間が経ったように思う。 実際のところ IANA の中央在庫が枯渇したのは 2011 年 2 月 3 日なので、実に 7 年以上が経っている。

IANAにおけるIPv4アドレス在庫枯渇、およびJPNICの今後のアドレス分配について - JPNIC

IPv4 アドレスの在庫は、階層構造になったインターネットレジストリが管理している。 階層構造において下流に位置する組織は、上流の組織に対して必要に応じて割り当てを申請して、在庫の一部を受け取る。 しかし、在庫が枯渇してしまった上流の組織からは、新規の割り当てが制限される。 ちなみに、日本のインターネットレジストリである JPNIC の在庫は 2011 年 4 月 15 日に枯渇している。

IPv4アドレスの在庫枯渇に関して - JPNIC

新規の割り当てを受けられない状況において、ネットワーク事業者は既存の割り当て済みアドレスの利用効率を上げるなど、いくつかの対策が考えられる。 しかし、現実問題として既存の割り当て済みアドレスだけではどうにもならないようなケースも存在する。 例えば新たに IaaS のようなサービスを始めたり、設備を大幅に増強するような局面においては IPv4 アドレスが大量に必要となる恐れがある。

そのような場合の救済策として、他の事業者に割り当てられたアドレスを一定の条件下で譲り受ける (移転する) ことが可能になっている。 これは、過去に割り当てを受けたものの、既に不要となった IPv4 アドレスが企業などで死蔵されているケースなどがあるため。 もちろん、慈善事業ではないのでアドレスの移転を受ける場合には基本的に有償になると考えられる。 ただし、補足しておくとアドレスの移転ポリシーについては管轄のインターネットレジストリによって色々と事情が異なる。

そうした状況において、アドレスを買いたい事業者と売りたい事業者を仲介するブローカーも登場している。 前置きが長くなったけど、今回扱うのはそんなブローカーの一つである以下のサイト。 ここでは売却したい IPv4 アドレスのブロックがオークションにかけられている。

www.ipv4auctions.com

上記のサイトには、過去に売却された IPv4 アドレスのブロックに関する情報が記載されている。 そこで、今回はそれをスクレイピングして可視化してみることにした。 このサイトにおける売却は全世界のアドレス移転のごく一部を見ているにすぎないとはいえ、傾向については読み取れるんじゃないかと。

IPv4 アドレスの価格の推移

早速だけど、まずは IPv4 アドレスの価格の推移から。 横軸に売却された日付、縦軸に単一アドレスあたりの価格 (米ドル) をプロットした散布図にしてみた。 点の色は、中央在庫を管理している IANA の直下に位置する RIR (Regional Internet Registry: 地域インターネットレジストリ) を表している。

f:id:momijiame:20181031215956p:plain

2014 ~ 2016 年頃までは、だいたい 7 ~ 10 USD ほどで売却されていたアドレスが、最近では 17 ~ 20 USD 近くまで値上がりしていることが分かる。 仮にこの価格が下がるとすれば、もしかすると IPv4 アドレスが不要になりつつある兆しかもしれない。 他の要因ももちろんあるとはいえ IPv4 のインターネットが縮小期に入っている可能性はある。

また、このサイトでの売却は、北アメリカを管轄している ARIN (American Registry for Internet Numbers) のアドレスブロックで活発なようだ。 北アメリカはインターネット発祥の地とも言えるので、やはり初期に大きく割り当てられて死蔵されているアドレスブロックも多かったりするんだろうか?

売却されたアドレスの数・回数の推移

続いては、売却された IPv4 アドレスの数や回数の推移を時系列でプロットしてみた。

まず、以下は売却された IPv4 アドレスの数を /24 のブロック (232-24 = 256 個) 単位で時系列の棒グラフにしている。

f:id:momijiame:20181031220016p:plain

上記を見ると、売却される IPv4 アドレスの数は近年増加傾向にあることが分かる。

また、同様に売却の件数についても月ごとにプロットしてみた。

f:id:momijiame:20181031220023p:plain

こちらも近年は増加傾向にあることが見て取れる。

売却されたブロックサイズの推移

続いては、売却されたアドレスブロックのサイズについても月ごとに集計してプロットしてみた。 月単位での CIDR 長の最大値 (max)、最小値 (min)、平均値 (mean)、中央値 (median) を示している。

f:id:momijiame:20181031220006p:plain

これについては、時系列での傾向は特に見受けられなかった。 最も小さなブロックのサイズが期間を通じて /24 で固定なのは、それより小さいと BGP の経路フィルタで落とされる可能性が高いためだろうか?

ブロックサイズごとの価格の違い

ブロックサイズについて見たところで、次はブロックサイズと価格の関係をプロットしてみる。

f:id:momijiame:20181117192921p:plain

/17 だけやけに安いのは、そもそもオークションに出品された件数が少なく、それも初期に固まっているためだろう。 それ以外のブロックサイズについては、サイズごとに極端な傾向は見られない。 ただ、/20 ~ /24 の間に、ブロックサイズが大きいほどアドレス単価の中央値が安くなる傾向が見られるのは意外だった。 それというのも、ブロックサイズは大きい方が使い勝手が良いため、大きいほど高いのかと思っていた。 ブロックサイズが大きくなるほど取引にかかる総額も大きくなるため、価格が抑制されやすいという可能性はあるかもしれない。 あるいは、価格が上昇すると共に単に大きめのブロックサイズが出品されにくくなっているという疑似相関だろうか。

まとめ

今回は IPv4 アドレスのオークションサイトのデータをスクレイピングしてグラフとして可視化してみた。 最近は IPv4 アドレスの値段が値上がり傾向にあることが分かった。

備考

今回スクレイピングと可視化に使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.13.6
BuildVersion:   17G65
$ python -V          
Python 3.6.7

スクリプトを動作させるのに必要となるパッケージは次の通り。

$ pip install pandas lxml matplotlib tqdm percache

以下はスクレイピングと可視化に用いたスクリプト。

IPv4 アドレスの価格の推移

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

import datetime
import time

import percache
from tqdm import tqdm
import pandas as pd
from matplotlib import pyplot as plt


def _fetch(year):
    """オークションの履歴をスクレイピングする"""
    base_url = 'https://www.ipv4auctions.com/customer/account/previous?year={yyyy}'
    url = base_url.format(yyyy=year)
    dfs = pd.read_html(url)
    return dfs[1]


# スクレイピングした内容をローカルのディスクにキャッシュする
cache = percache.Cache('df-cache')


@cache
def auction_history(start_year, stop_year, fetch_interval=2):
    """各年の情報を結合して一枚のデータフレームにする"""
    dfs = []
    years = range(start_year, stop_year + 1)
    for year in tqdm(years):
        year_df = _fetch(year)
        dfs.append(year_df)
        # サービスに負荷をかけないようにインターバルをかける
        time.sleep(fetch_interval)

    # データフレームを結合する
    return pd.concat(dfs, axis=0)


def filter_outlier(df_):
    """データフレームから外れ値を取り除く"""
    # ブロックサイズが適切にパースできるものだけに絞る
    prefix_series = df_['BLOCK'].apply(lambda x: x[1:])
    return df_[prefix_series.str.isnumeric()]


def sold_region(df_):
    """取引された地域をパースする"""
    def parse(x):
        # 末尾にカンマの入ったデータが混ざっていたので取り除く
        return x[:-1] if x[-1] == ',' else x
    return df_.REGION.apply(parse)


def sold_date(df_):
    """取引された日付をパースする"""
    def parse(x):
        datetime_ = datetime.datetime.strptime(x, '%Y-%m-%d')
        return datetime.date(datetime_.year, datetime_.month, datetime_.day)
    return df_['SOLD DATE'].apply(parse)


def sold_yyyymm(df_):
    """取引された年・月をパースする"""
    def to_str(dt):
        return dt.strftime('%Y-%m')
    dt_series = pd.to_datetime(df_['SOLD DATE'])
    return dt_series.apply(to_str)


def price_per_addr(df_):
    """取引された価格をパースする"""
    def parse(x):
        ex_comma_x = x.replace(',', '')
        return float(ex_comma_x[1:])
    return df_['PRICE PER ADDRESS'].apply(parse)


def main():
    # データ取得・描画範囲
    YEAR_RANGE = (2014, 2018)

    # オークションの履歴を取得する
    df = auction_history(YEAR_RANGE[0], YEAR_RANGE[-1])

    # 外れ値を取り除く
    df = filter_outlier(df)

    # 可視化に使うデータを取り出す
    df = df.assign(sold_date=sold_date(df))
    df = df.assign(price_per_addr=price_per_addr(df))
    df = df.assign(sold_region=sold_region(df))
    df = df.assign(yyyymm=sold_yyyymm(df))

    # 地域ごとに集計・可視化する
    fig, ax = plt.subplots()
    regions = pd.unique(df.sold_region)
    for region in regions:
        region_df = df[df.sold_region == region]
        ax.plot_date(region_df.sold_date, region_df.price_per_addr, label=region)

    ax.set_title('IPv4 auction price graph')
    ax.legend()
    ax.set_ylabel('price per address (USD)')
    ax.set_xlabel('sold date')
    plt.show()


if __name__ == '__main__':
    main()

売却されたアドレスの数の推移

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

import time

import percache
from tqdm import tqdm
import pandas as pd
from matplotlib import pyplot as plt


def _fetch(year):
    """オークションの履歴をスクレイピングする"""
    base_url = 'https://www.ipv4auctions.com/customer/account/previous?year={yyyy}'
    url = base_url.format(yyyy=year)
    dfs = pd.read_html(url)
    return dfs[1]


# スクレイピングした内容をローカルのディスクにキャッシュする
cache = percache.Cache('df-cache')


@cache
def auction_history(start_year, stop_year, fetch_interval=2):
    """各年の情報を結合して一枚のデータフレームにする"""
    dfs = []
    years = range(start_year, stop_year + 1)
    for year in tqdm(years):
        year_df = _fetch(year)
        dfs.append(year_df)
        # サービスに負荷をかけないようにインターバルをかける
        time.sleep(fetch_interval)

    # データフレームを結合する
    return pd.concat(dfs, axis=0)


def class_c_blocks(df_):
    """取引された /24 ブロックの数を数える"""
    def parse(x):
        return 2 ** (24 - int(x[1:]))
    return df_['BLOCK'].apply(parse)


def filter_outlier(df_):
    """データフレームから外れ値を取り除く"""
    # ブロックサイズが適切にパースできるものだけに絞る
    prefix_series = df_['BLOCK'].apply(lambda x: x[1:])
    return df_[prefix_series.str.isnumeric()]


def sold_yyyymm(df_):
    """取引された年・月をパースする"""
    def to_str(dt):
        return dt.strftime('%Y-%m')
    dt_series = pd.to_datetime(df_['SOLD DATE'])
    return dt_series.apply(to_str)


def region(df_):
    """取引された地域をパースする"""
    def parse(x):
        # 末尾にカンマの入ったデータが混ざっていたので取り除く
        return x[:-1] if x[-1] == ',' else x
    return df_.REGION.apply(parse)


def main():
    # データ取得・描画範囲
    YEAR_RANGE = (2014, 2018)

    # オークションの履歴を取得する
    df = auction_history(YEAR_RANGE[0], YEAR_RANGE[-1])

    # 外れ値を取り除く
    df = filter_outlier(df)

    # 必要なデータを取り出す
    df = df.assign(blocks=class_c_blocks(df))
    df = df.assign(yyyymm=sold_yyyymm(df))
    df = df.assign(region=region(df))

    # 集計する
    pivot_df = df.pivot_table(index=['yyyymm'], columns=['region'], values=['blocks'], aggfunc='sum', fill_value=0)

    # 可視化する
    fig, ax = plt.subplots()
    for _, data in pivot_df.items():
        region_name = data.name[1]
        ax.bar(data.index, data, label=region_name)

    ax.set_title('Number of /24 blocks sold in the auction')
    ax.legend()
    ax.set_ylabel('sold /24 blocks')
    ax.set_xlabel('month')
    ax.xaxis.set_ticks(['{y}-06'.format(y=y) for y in range(YEAR_RANGE[0], YEAR_RANGE[-1] + 1)])
    plt.show()


if __name__ == '__main__':
    main()

売却された回数の推移

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

import time

import percache
from tqdm import tqdm
import pandas as pd
from matplotlib import pyplot as plt


def _fetch(year):
    """オークションの履歴をスクレイピングする"""
    base_url = 'https://www.ipv4auctions.com/customer/account/previous?year={yyyy}'
    url = base_url.format(yyyy=year)
    dfs = pd.read_html(url)
    return dfs[1]


# スクレイピングした内容をローカルのディスクにキャッシュする
cache = percache.Cache('df-cache')


@cache
def auction_history(start_year, stop_year, fetch_interval=2):
    """各年の情報を結合して一枚のデータフレームにする"""
    dfs = []
    years = range(start_year, stop_year + 1)
    for year in tqdm(years):
        year_df = _fetch(year)
        dfs.append(year_df)
        # サービスに負荷をかけないようにインターバルをかける
        time.sleep(fetch_interval)

    # データフレームを結合する
    return pd.concat(dfs, axis=0)


def class_c_blocks(df_):
    """取引された /24 ブロックの数を数える"""
    def parse(x):
        return 2 ** (24 - int(x[1:]))
    return df_['BLOCK'].apply(parse)


def filter_outlier(df_):
    """データフレームから外れ値を取り除く"""
    # ブロックサイズが適切にパースできるものだけに絞る
    prefix_series = df_['BLOCK'].apply(lambda x: x[1:])
    return df_[prefix_series.str.isnumeric()]


def sold_yyyymm(df_):
    """取引された年・月をパースする"""
    def to_str(dt):
        return dt.strftime('%Y-%m')
    dt_series = pd.to_datetime(df_['SOLD DATE'])
    return dt_series.apply(to_str)


def region(df_):
    """取引された地域をパースする"""
    def parse(x):
        # 末尾にカンマの入ったデータが混ざっていたので取り除く
        return x[:-1] if x[-1] == ',' else x
    return df_.REGION.apply(parse)


def main():
    # データ取得・描画範囲
    YEAR_RANGE = (2014, 2018)

    # オークションの履歴を取得する
    df = auction_history(YEAR_RANGE[0], YEAR_RANGE[-1])

    # 外れ値を取り除く
    df = filter_outlier(df)

    # 必要なデータを取り出す
    df = df.assign(blocks=class_c_blocks(df))
    df = df.assign(yyyymm=sold_yyyymm(df))
    df = df.assign(region=region(df))

    # 集計する
    pivot_df = df.pivot_table(index=['yyyymm'], columns=['region'], values=['blocks'], aggfunc='count', fill_value=0)

    # 可視化する
    fig, ax = plt.subplots()
    for _, data in pivot_df.items():
        region_name = data.name[1]
        ax.bar(data.index, data, label=region_name)

    ax.set_title('Number of sold counts in the auction')
    ax.legend()
    ax.set_ylabel('sold count')
    ax.set_xlabel('month')
    ax.xaxis.set_ticks(['{y}-06'.format(y=y) for y in range(YEAR_RANGE[0], YEAR_RANGE[-1] + 1)])
    plt.show()


if __name__ == '__main__':
    main()

売却されたブロックサイズの推移

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

import time

import percache
from tqdm import tqdm
import pandas as pd
import matplotlib.pyplot as plt


def _fetch(year):
    """オークションの履歴をスクレイピングする"""
    base_url = 'https://www.ipv4auctions.com/customer/account/previous?year={yyyy}'
    url = base_url.format(yyyy=year)
    dfs = pd.read_html(url)
    return dfs[1]


# スクレイピングした内容をローカルのディスクにキャッシュする
cache = percache.Cache('df-cache')


@cache
def auction_history(start_year, stop_year, fetch_interval=2):
    """各年の情報を結合して一枚のデータフレームにする"""
    dfs = []
    years = range(start_year, stop_year + 1)
    for year in tqdm(years):
        year_df = _fetch(year)
        dfs.append(year_df)
        # サービスに負荷をかけないようにインターバルをかける
        time.sleep(fetch_interval)

    # データフレームを結合する
    return pd.concat(dfs, axis=0)


def block_size(df_):
    """取引された IPv4 アドレスのブロックサイズをパースする"""
    def parse(x):
        return int(x[1:])
    return df_['BLOCK'].apply(parse)


def filter_outlier(df_):
    """データフレームから外れ値を取り除く"""
    # ブロックサイズが適切にパースできるものだけに絞る
    prefix_series = df_['BLOCK'].apply(lambda x: x[1:])
    return df_[prefix_series.str.isnumeric()]


def sold_yyyymm(df_):
    """取引された年・月をパースする"""
    def to_str(dt):
        return dt.strftime('%Y-%m')
    dt_series = pd.to_datetime(df_['SOLD DATE'])
    return dt_series.apply(to_str)


def region(df_):
    """取引された地域をパースする"""
    def parse(x):
        # 末尾にカンマの入ったデータが混ざっていたので取り除く
        return x[:-1] if x[-1] == ',' else x
    return df_.REGION.apply(parse)


def main():
    # データ取得・描画範囲
    YEAR_RANGE = (2014, 2018)

    # オークションの履歴を取得する
    df = auction_history(YEAR_RANGE[0], YEAR_RANGE[-1])

    # 外れ値を取り除く
    df = filter_outlier(df)

    # 取引された日付と価格を取得する
    df = df.assign(block_size=block_size(df))
    df = df.assign(yyyymm=sold_yyyymm(df))
    df = df.assign(region=region(df))

    # 集計する
    pivot_df = df.pivot_table(index=['yyyymm'], values=['block_size'], aggfunc=['max', 'mean', 'median', 'min'], fill_value=0)

    # 可視化する
    fig, ax = plt.subplots()
    for _, data in pivot_df.items():
        agg_type = data.name[0]
        ax.plot(data.index, data, label=agg_type)

    ax.legend()
    ax.set_ylabel('cidr')
    ax.set_xlabel('month')
    ax.xaxis.set_ticks(['{y}-06'.format(y=y) for y in range(YEAR_RANGE[0], YEAR_RANGE[-1] + 1)])
    plt.show()


if __name__ == '__main__':
    main()
#!/usr/bin/env python
# -*- coding: utf-8 -*-

import time

import percache
from tqdm import tqdm
import pandas as pd
from matplotlib import pyplot as plt


def _fetch(year):
    """オークションの履歴をスクレイピングする"""
    base_url = 'https://www.ipv4auctions.com/customer/account/previous?year={yyyy}'
    url = base_url.format(yyyy=year)
    dfs = pd.read_html(url)
    return dfs[1]


# スクレイピングした内容をローカルのディスクにキャッシュする
cache = percache.Cache('df-cache')


@cache
def auction_history(start_year, stop_year, fetch_interval=2):
    """各年の情報を結合して一枚のデータフレームにする"""
    dfs = []
    years = range(start_year, stop_year + 1)
    for year in tqdm(years):
        year_df = _fetch(year)
        dfs.append(year_df)
        # サービスに負荷をかけないようにインターバルをかける
        time.sleep(fetch_interval)

    # データフレームを結合する
    return pd.concat(dfs, axis=0)


def filter_outlier(df_):
    """データフレームから外れ値を取り除く"""
    # ブロックサイズが適切にパースできるものだけに絞る
    prefix_series = df_['BLOCK'].apply(lambda x: x[1:])
    return df_[prefix_series.str.isnumeric()]


def price_per_addr(df_):
    """取引された価格をパースする"""
    def parse(x):
        ex_comma_x = x.replace(',', '')
        return float(ex_comma_x[1:])
    return df_['PRICE PER ADDRESS'].apply(parse)


def sold_block_size(df_):
    """取引されたブロックのサイズ"""
    def parse(x):
        return int(x[1:])
    return df_['BLOCK'].apply(parse)


def main():
    # データ取得・描画範囲
    YEAR_RANGE = (2014, 2018)

    # オークションの履歴を取得する
    df = auction_history(YEAR_RANGE[0], YEAR_RANGE[-1])

    # 外れ値を取り除く
    df = filter_outlier(df)

    # 可視化に使うデータを取り出す
    df = df.assign(price_per_addr=price_per_addr(df))
    df = df.assign(sold_block_size=sold_block_size(df))

    # ブロックサイズごとに集計・可視化する
    fig, ax = plt.subplots()
    sizes = sorted(pd.unique(df.sold_block_size))
    block_size_series = [df[df.sold_block_size == size].price_per_addr for size in sizes]
    ax.boxplot(block_size_series, labels=sizes)

    ax.set_title('IPv4 auction price graph')
    ax.set_ylabel('price per address (USD)')
    ax.set_xlabel('sold block size')
    plt.show()


if __name__ == '__main__':
    main()

いじょう。

マスタリングTCP/IP 入門編 第5版

マスタリングTCP/IP 入門編 第5版

  • 作者: 竹下隆史,村山公保,荒井透,苅田幸雄
  • 出版社/メーカー: オーム社
  • 発売日: 2012/02/25
  • メディア: 単行本(ソフトカバー)
  • 購入: 4人 クリック: 34回
  • この商品を含むブログ (37件) を見る

Python: pandas-profiling でデータセットの概要を確認する

今回は pandas-profiling というパッケージを使ってみる。 このパッケージを使うと pandas の DataFrame に含まれる各次元の基本的な統計量や相関係数などを一度に確認できる。 最初にデータセットのサマリーを確認できると、その後の EDA (Exploratory Data Analysis: 探索的データ分析) の取っ掛かりにしやすいと思う。

使った環境は次の通り。

$ sw_vers 
ProductName:    macOS
ProductVersion: 12.4
BuildVersion:   21F79
$ python -V         
Python 3.9.13
$ pip3 list | grep pandas-profiling
pandas-profiling              3.2.0

下準備

まずは必要なパッケージをインストールしておく。

$ pip install pandas-profiling

サンプルのデータセットを用意する

サンプルとなるデータセットは Kaggle でおなじみの Titanic データセットにした。

Titanic - Machine Learning from Disaster | Kaggle

データセットは Web サイトか、あるいはコマンドラインツールからダウンロードしておく。

$ pip install kaggle
$ kaggle competitions download -c titanic
$ head train.csv 
PassengerId,Survived,Pclass,Name,Sex,Age,SibSp,Parch,Ticket,Fare,Cabin,Embarked
1,0,3,"Braund, Mr. Owen Harris",male,22,1,0,A/5 21171,7.25,,S
2,1,1,"Cumings, Mrs. John Bradley (Florence Briggs Thayer)",female,38,1,0,PC 17599,71.2833,C85,C
3,1,3,"Heikkinen, Miss. Laina",female,26,0,0,STON/O2. 3101282,7.925,,S
4,1,1,"Futrelle, Mrs. Jacques Heath (Lily May Peel)",female,35,1,0,113803,53.1,C123,S
5,0,3,"Allen, Mr. William Henry",male,35,0,0,373450,8.05,,S
6,0,3,"Moran, Mr. James",male,,0,0,330877,8.4583,,Q
7,0,1,"McCarthy, Mr. Timothy J",male,54,0,0,17463,51.8625,E46,S
8,0,3,"Palsson, Master. Gosta Leonard",male,2,3,1,349909,21.075,,S
9,1,3,"Johnson, Mrs. Oscar W (Elisabeth Vilhelmina Berg)",female,27,0,2,347742,11.1333,,S

pandas-profiling を使ってみる

それでは、実際に pandas-profiling を使ってみる。 まずは Python のインタプリタを起動する。

$ python

データセットの CSV を pandas で読み込む。

>>> import pandas as pd
>>> df = pd.read_csv('train.csv')

あとは読み込んだ DataFrame を pandas-profiling に渡すだけ。

>>> import pandas_profiling
>>> profile_report = pandas_profiling.ProfileReport(df)

解析結果は、次のように HTML で出力できる。

>>> profile_report.to_file('report.html')

一旦、インタプリタからは抜けておこう。

>>> exit()

あとは出力された HTML を確認する。

$ open report.html

このように、各次元の基本的な統計量や欠損値の有無、分布や相関係数などが一度に確認できる。

いじょう。

Python: scikit-learn の FeatureUnion を pandas の DataFrame と一緒に使う

今回は scikit-learnFeatureUnionpandasDataFrame を一緒に使うときの問題点とその解決策について。 scikit-learn の FeatureUnion は、典型的には Pipeline においてバラバラに作った複数の特徴量を一つにまとめるのに使われる機能。 この FeatureUnion を pandas の DataFrame と一緒に使おうとすると、ちょっとばかり予想外の挙動になる。 具体的には FeatureUnion の出力するデータが、本来なら DataFrame になってほしいところで numpy の ndarray 形式に変換されてしまう。 今回は、それをなんとか DataFrame に直す方法がないか調べたり模索してみた話。

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.13.6
BuildVersion:   17G65
$ python -V      
Python 3.7.0
$ pip list --format=columns | egrep -i "(scikit-learn|pandas)"
pandas          0.23.4 
scikit-learn    0.20.0 

下準備

まずは今回使うパッケージをインストールしておく。

$ pip install scikit-learn pandas

FeatureUnion について

まずは前提知識として FeatureUnion について紹介する。 これは Pipeline と一緒に使うことで真価を発揮する機能で、複数の特徴量を一つにまとめ上げるのに使われる。

次のサンプルコードでは Iris データセットを題材にして Pipeline と FeatureUnion を使っている。 この中では、元々 Iris データセットに含まれていた特徴量と、それらを標準化した特徴量を並列に準備して結合している。 動作確認のために Pipeline への入力データと出力データの情報を標準出力に書き出している。

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

from sklearn import datasets
from sklearn.pipeline import FeatureUnion
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import FunctionTransformer


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

    # パイプラインの入力データ
    print('input:')
    print(type(X), X.shape)
    print(X[:5])

    steps = [
        # 複数の特徴量を結合する
        ('union', FeatureUnion([
            # 生の特徴量をそのまま返す
            ('raw', FunctionTransformer(lambda x: x, validate=False)),
            # 標準化した特徴量を返す
            ('scaler', StandardScaler()),
        ])),
    ]
    pipeline = Pipeline(steps=steps)
    transformed_data = pipeline.fit_transform(X)

    # パイプラインの出力データ
    print('output:')
    print(type(transformed_data), transformed_data.shape)
    print(transformed_data[:5])


if __name__ == '__main__':
    main()

論より証拠ということで上記を実行してみよう。 すると、入力時点では 4 次元 150 行だったデータが、出力時点では 8 次元 150 行になっている。 これは、元々 Iris データセットに含まれていた 4 次元の特徴量に加えて、それらを標準化した特徴量が 4 次元加わっているため。

$ python featunion.py
input:
<class 'numpy.ndarray'> (150, 4)
[[5.1 3.5 1.4 0.2]
 [4.9 3.  1.4 0.2]
 [4.7 3.2 1.3 0.2]
 [4.6 3.1 1.5 0.2]
 [5.  3.6 1.4 0.2]]
output:
<class 'numpy.ndarray'> (150, 8)
[[ 5.1         3.5         1.4         0.2        -0.90068117  1.01900435
  -1.34022653 -1.3154443 ]
 [ 4.9         3.          1.4         0.2        -1.14301691 -0.13197948
  -1.34022653 -1.3154443 ]
 [ 4.7         3.2         1.3         0.2        -1.38535265  0.32841405
  -1.39706395 -1.3154443 ]
 [ 4.6         3.1         1.5         0.2        -1.50652052  0.09821729
  -1.2833891  -1.3154443 ]
 [ 5.          3.6         1.4         0.2        -1.02184904  1.24920112
  -1.34022653 -1.3154443 ]]

このように FeatureUnion を使うと特徴量エンジニアリングの作業が簡潔に表現できる。

pandas の DataFrame と一緒に使うときの問題点について

ただ、この FeatureUnion を pandas の DataFrame と一緒に使おうとすると問題が出てくる。 具体的には FeatureUnion を通すと、入力が DataFrame でも出力が numpy の ndarray の変換されてしまうのだ。

次のサンプルコードを見てほしい。 このコードでは Pipeline の入力データを pandas の DataFrame に変換している。

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

import pandas as pd

from sklearn import datasets
from sklearn.pipeline import Pipeline
from sklearn.pipeline import FeatureUnion
from sklearn.preprocessing import FunctionTransformer


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

    # データセットを pandas の DataFrame に変換する
    X_df = pd.DataFrame(X, columns=dataset.feature_names)

    # パイプラインの入力データ
    print('input:')
    print(type(X_df))
    print(X_df.head())

    steps = [
        ('union', FeatureUnion([
            ('raw', FunctionTransformer(lambda x: x, validate=False)),
        ])),
    ]
    pipeline = Pipeline(steps=steps)
    transformed_data = pipeline.fit_transform(X_df)

    # パイプラインの出力データ
    print('output:')
    print(type(transformed_data))
    print(transformed_data[:5])


if __name__ == '__main__':
    main()

上記を実行してみよう。 すると Pipeline の出力データが numpy の ndarray に変換されてしまっていることが分かる。

$ python dfissue.py
input:
<class 'pandas.core.frame.DataFrame'>
   sepal length (cm)        ...         petal width (cm)
0                5.1        ...                      0.2
1                4.9        ...                      0.2
2                4.7        ...                      0.2
3                4.6        ...                      0.2
4                5.0        ...                      0.2

[5 rows x 4 columns]
output:
<class 'numpy.ndarray'>
[[5.1 3.5 1.4 0.2]
 [4.9 3.  1.4 0.2]
 [4.7 3.2 1.3 0.2]
 [4.6 3.1 1.5 0.2]
 [5.  3.6 1.4 0.2]]

上記のような振る舞いになる原因は scikit-learn の以下のコードにある。 FeatureUnion では特徴量の結合に numpy.hstack() を使っているのだ。 この関数の返り値は numpy の ndarray になるので FeatureUnion の出力もそれになってしまう。

github.com

ColumnTransformer でも同じ問題が起こる

ちなみに FeatureUnion と似たような機能の ColumnTransformer を使っても同じ問題が起こる。 ColumnTransformer は処理対象のカラム (次元) を指定して特徴量が生成できる。

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

import pandas as pd

from sklearn import datasets
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import FunctionTransformer


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

    # データセットを pandas の DataFrame に変換する
    X_df = pd.DataFrame(X, columns=dataset.feature_names)

    # パイプラインの入力データ
    print('input:')
    print(type(X_df))
    print(X_df.head())

    steps = [
        ('union', ColumnTransformer([
            ('transparent', FunctionTransformer(lambda x: x, validate=False), [i for i, _ in enumerate(X_df)]),
        ])),
    ]
    pipeline = Pipeline(steps=steps)
    transformed_data = pipeline.fit_transform(X_df)

    # パイプラインの出力データ
    print('output:')
    print(type(transformed_data))
    print(transformed_data[:5])


if __name__ == '__main__':
    main()

上記を実行してみよう。 FeatureUnion のときと同じように出力データは numpy の ndarray 形式に変換されてしまった。

$ python columnissue.py 
input:
<class 'pandas.core.frame.DataFrame'>
   sepal length (cm)        ...         petal width (cm)
0                5.1        ...                      0.2
1                4.9        ...                      0.2
2                4.7        ...                      0.2
3                4.6        ...                      0.2
4                5.0        ...                      0.2

[5 rows x 4 columns]
output:
<class 'numpy.ndarray'>
[[5.1 3.5 1.4 0.2]
 [4.9 3.  1.4 0.2]
 [4.7 3.2 1.3 0.2]
 [4.6 3.1 1.5 0.2]
 [5.  3.6 1.4 0.2]]

解決方法について

この問題について調べると、同じように悩んでいる人たちがいる。 そして、例えば次のような解決策を提示している人がいた。

zablo.net

上記では、自分で Transformer と FeatureUnion を拡張することで対応している。 ただ、このやり方だと作業が大掛かりで scikit-learn の API 変更などに継続的に追随していく必要があると感じた。

そこで、別の解決策として「とてもダーティーなハック」を思いついてしまったので共有してみる。 この問題は、特徴量の結合に numpy の hstack() 関数を使っていることが原因だった。 だとすれば、この関数をモンキーパッチで一時的にすり替えてしまえば問題を回避できるのではないか?と考えた。

実際に試してみたのが次のサンプルコードになる。 このコードでは unittest.mock.patch() を使って、一時的に numpy.hstack()pd.concat() にすり替えている。 これなら出力データは pandas の DataFrame にできるのではないか?

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

from functools import partial
from unittest.mock import patch

import pandas as pd

from sklearn import datasets
from sklearn.pipeline import Pipeline
from sklearn.pipeline import FeatureUnion
from sklearn.preprocessing import FunctionTransformer


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

    # データセットを pandas の DataFrame に変換する
    X_df = pd.DataFrame(X, columns=dataset.feature_names)

    # パイプラインの入力データ
    print('input:')
    print(type(X_df))
    print(X_df.head())

    steps = [
        ('union', FeatureUnion([
            ('transparent', FunctionTransformer(lambda x: x, validate=False)),
        ])),
    ]
    pipeline = Pipeline(steps=steps)

    # numpy.hstack を一時的に pd.concat(axis=1) にパッチしてしまう
    horizontal_concat = partial(pd.concat, axis=1)
    with patch('numpy.hstack', side_effect=horizontal_concat):
        transformed_data = pipeline.fit_transform(X_df)

    print('output:')
    print(type(transformed_data))
    print(transformed_data.head())


if __name__ == '__main__':
    main()

上記を実行してみよう。 すると、見事に出力データが pandas の DataFrame になっていることが分かる。

$ python monkeypatch.py
input:
<class 'pandas.core.frame.DataFrame'>
   sepal length (cm)        ...         petal width (cm)
0                5.1        ...                      0.2
1                4.9        ...                      0.2
2                4.7        ...                      0.2
3                4.6        ...                      0.2
4                5.0        ...                      0.2

[5 rows x 4 columns]
output:
<class 'pandas.core.frame.DataFrame'>
   sepal length (cm)        ...         petal width (cm)
0                5.1        ...                      0.2
1                4.9        ...                      0.2
2                4.7        ...                      0.2
3                4.6        ...                      0.2
4                5.0        ...                      0.2

[5 rows x 4 columns]

この解決策の問題点としては、学習の過程において numpy.hstack() が使えなくなるところ。 pd.concat() は入力として pandas のオブジェクトしか渡すことができない。 もし numpy の ndarray を渡してしまうと、次のようにエラーになる。

>>> import numpy as n
>>> l1 = np.array([1, 2, 3])
>>> l2 = np.array([4, 5, 6])
>>> import pandas as pd
>>> pd.concat([l1, l2], axis=1)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/amedama/.virtualenvs/py37/lib/python3.7/site-packages/pandas/core/reshape/concat.py", line 225, in concat
    copy=copy, sort=sort)
  File "/Users/amedama/.virtualenvs/py37/lib/python3.7/site-packages/pandas/core/reshape/concat.py", line 286, in __init__
    raise TypeError(msg)
TypeError: cannot concatenate object of type "<class 'numpy.ndarray'>"; only pd.Series, pd.DataFrame, and pd.Panel (deprecated) objs are valid

巨大なパイプラインであれば numpy.hstack() をどこかで使っている可能性は十分に考えられる。 とはいえ、一旦試して駄目だったら前述した別の解決策に切り替える、というのもありかもしれない。

いじょう。

Pythonではじめる機械学習 ―scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎

Pythonではじめる機械学習 ―scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎

scikit-learnとTensorFlowによる実践機械学習

scikit-learnとTensorFlowによる実践機械学習

Pythonデータサイエンスハンドブック ―Jupyter、NumPy、pandas、Matplotlib、scikit-learnを使ったデータ分析、機械学習

Pythonデータサイエンスハンドブック ―Jupyter、NumPy、pandas、Matplotlib、scikit-learnを使ったデータ分析、機械学習