Python: k-NN Feature Extraction について

k-NN Feature Extraction (k-近傍法を用いた特徴量抽出) という手法があるらしい。 これは、文字通り k-NN (k-Nearest Neighbor algorithm: k-近傍法) を特徴量の抽出に応用したもの。 興味深かったので、今回は自分でも Python を使って実装してみた。 手法について知ったのは、以下のブログを目にしたのがきっかけ。

upura.hatenablog.com

また、上記は以下のブログに記載のある R の実装を参考にしているとのことだった。

Feature Extraction with KNN • fastknn

ただ、先ほどの Python API では、特徴量を付与する対象のデータをパラメータとして指定できない点が気になった。 具体的には、以下のような交差検証を使った性能の計測が難しいのではと感じた。

• データセットを学習用と検証用に分割する
• 学習用データに対して k-NN Feature Extraction で特徴量を付与する
  • 付与するデータと付与に使うデータに分けながら順番に付与していく
• 特徴量が付与された学習用データを使ってモデルを学習させる
• 学習用データを使って検証用データに k-NN Feature Extraction で特徴量を付与する
• 特徴量が付与された検証用データを使ってモデルの性能を計測する

R の実装では、特徴量を付与するデータと付与に使うデータを別々に指定できる。 そこで、今回書いてみたものは R の実装にならって別々に指定できるようにしてみた。

使った環境は次の通り。

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

k-NN Feature Extraction について

ここからは k-NN Feature Extraction の具体的な手法について解説していく。 まず、k-NN Feature Extraction では、名前の通り k-NN (k-近傍法) を特徴量の抽出に応用している。 元々の k-NN は分類問題を解くためのアルゴリズムだった。

最初に k-NN では、分類したい対象のデータから最も近い k 点の学習データを探し出す。 k 点の k には 1 や 3 など、任意の数を用いる。 最も近い k 点の学習データが見つかったら、その学習データについているクラスを使って対象のデータを分類する。 k が 1 より大きいときは、多数決などを使って分類先のクラスを決める。

k-NN Feature Extraction では、上記のアルゴリズムのうち前半部分だけを利用する。 具体的には、最も近い k 点の学習データを探し出す部分。 このとき、最も近い k 点の学習データを探すために、データ間の距離を計算することになる。 k-NN Feature Extraction では、このデータ間の距離を特徴量として用いる。

少し k-NN と異なるのは、k 点の近傍データを探すのがクラスごとという点。 つまり、二値分類問題であれば、学習データの中のクラス 0 とクラス 1 ごとに近傍データをそれぞれ探す。 もちろん多値分類問題であれば、クラス数 c の数だけ近傍データを探すことになる。 さらに、k が 1 を越える場合には、1 ~ k までの全ての数で近傍データを探すようだ。 これによって k * c 個の特徴量が新たに生成されることになる。

文章だけだとなかなか分かりにくいと思うので図示してみる。 まず、以下のような感じで特徴量を付与したいデータとクラスが三つに分かれた学習データがあるとする。

特徴量を付与したいデータから、それぞれのクラスの近傍データまでの距離を計算して特徴量にする。

もし k > 1 であれば、複数の近傍点までの距離を足し合わせたものを特徴量にする。

いじょう。

下準備

実際に実装してみる前に、必要なパッケージを環境にインストールしておく。

$ pip install scikit-learn matplotlib

サンプルデータについて

特徴量の付与に使うサンプルデータについては、以下のような線形分離できないものを用いた。

上記を生成したのは次のようなサンプルコードになる。

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

import numpy as np
from matplotlib import pyplot as plt


def main():
    # サンプルデータを生成する
    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])

    # データを可視化する
    plt.scatter(X[y == 0, 0], X[y == 0, 1])
    plt.scatter(X[y == 1, 0], X[y == 1, 1])
    plt.show()


if __name__ == '__main__':
    main()

適当な名前をつけて実行すると、先ほどのグラフが出力される。

$ python visualize.py

k-NN Feature Extraction で特徴量を抽出してみる

それでは、実際に k-NN Feature Extraction を試してみよう。 以下のサンプルコードの knn_extract() 関数が特徴量を生成する関数になる。 基本的な使い方としては、学習データとしてxtr と、そのクラスデータとして ytr を渡す。 そして、特徴量を生成する対象のデータとして xte を渡す。

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

import functools
from concurrent.futures import ProcessPoolExecutor
import multiprocessing

import numpy as np
from sklearn.model_selection import KFold
from matplotlib import pyplot as plt


def _distance(a, b):
    """2 点間の距離を計算する"""
    distance = np.linalg.norm(a - b)
    return distance


def _knn_distance(xtr_c, target, k):
    """特定クラスの最近傍な教師データへの距離を計算する"""
    # 各要素への距離を計算する
    distances = np.array([_distance(target, x) for x in xtr_c])
    # 距離で昇順ソートする
    sorted_distances = np.sort(distances)
    # 先頭 k 件を取り出す (最近傍への距離 k 点)
    nearest_distances = sorted_distances[:k]
    # 距離を足し算して返す
    sum_distances = np.sum(nearest_distances)
    return sum_distances


def _knn_distance_fold(xtr_c, k, folds, target):
    """過学習を防ぐためにサンプリングしながら最近傍への距離の平均を計算する"""
    # 途中までの計算結果を代入するための入れ物
    distances = np.empty([folds])

    # K-分割することでサンプリングする
    kf = KFold(shuffle=True, n_splits=folds)
    for i, (train_index, _) in enumerate(kf.split(xtr_c)):
        # サンプリングした教師信号
        xtr_c_sampled = xtr_c[train_index]
        # ターゲットと教師信号との距離を計算する
        distance = _knn_distance(xtr_c_sampled, target, k)
        distances[i] = distance

    # 各サンプルでの距離の平均を返す
    average_distance = distances.mean()
    return average_distance


def knn_extract(xtr, ytr, xte, k=1, folds=5, nprocesses=-1):
    if nprocesses == -1:
        # デフォルトでは CPU のコア数で並列計算する
        nprocesses = multiprocessing.cpu_count()

    # 教師データに含まれているクラス (ラベル) の種類
    # NOTE: CV するときは含まれるクラスが偏らないようにすること
    classes = np.unique(ytr)

    # 生成する特徴量の入れ物を用意する
    features = np.zeros([len(xte), len(classes) * k])

    for i, class_ in enumerate(classes):
        # 処理対象のクラス
        xtr_c = xtr[ytr == class_]

        for j, k_n in enumerate(range(1, k + 1), 1):
            # クラス class_ で近傍数 k の距離を計算する関数に実引数を部分適用する
            f = functools.partial(_knn_distance_fold, xtr_c, k_n, folds)

            # クラス class_ で近傍数 k で計算した特徴ベクトル
            with ProcessPoolExecutor(max_workers=nprocesses) as executor:
                feature = executor.map(f, xte)
                features[:, i * j] = list(feature)

    return features


def main():
    # サンプルデータを生成する
    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])

    classes = np.unique(y)
    k = 1
    extracted_features = np.empty([0, len(classes) * k])

    # データを 5 分割して、それぞれで特徴量を作っていく
    skf = KFold(n_splits=5)
    for train_index, test_index in skf.split(X):
        # 学習用データと検証用データに分ける
        X_train, X_test = X[train_index], X[test_index]
        y_train, y_test = y[train_index], y[test_index]

        # 検証用データに対して、学習用データを使って特徴量を抽出する
        extracted_feature = knn_extract(X_train, y_train, X_test, k)
        extracted_features = np.append(extracted_features, extracted_feature, axis=0)

    # k-NN Feature Extraction の結果を可視化する
    plt.scatter(extracted_features[:, 0], extracted_features[:, 1], c=y)
    plt.show()


if __name__ == '__main__':
    main()

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

$ python knnfe.py

実行結果として得られるグラフは次の通り。 今回の問題は k = 1, c = 2 なので二次元の特徴量が新たに得られる。 それらを散布図にしてみると、ちゃんと線形分離できそうな形で特徴量が生成されている。

続いて Iris データセットに使ってみた場合についても紹介する。 まず、データセットの読み込み部分を以下のように変更する。

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

こちらも、四次元の特徴量から以下のような二次元の特徴量が得られる。 それなりに上手く分離できそうな雰囲気がある。

課題について

実際に試した上で気になった点としては、やはり計算量があげられる。 今回使ったようなサイズの小さいデータセットでも、実行してみると結構な時間がかかる。 これは最近傍の探索が素朴な実装になっているから。 そのため、実用上は ANN (Approximate Nearest Neighbor: 近似最近傍) を使う必要があると考えられる。

いじょう。

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

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 の計測方法になる。 ただ、文章だけだとちょっと分かりにくいかもと思ったので図示してみる。 最初のベースラインを測る部分は一般的な機械学習モデルの汎化性能を計測するやり方と同じ。

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

上記において勘違いしやすいポイントは 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

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

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

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

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

いじょう。

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

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)

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>

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

いじょう。

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

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/shadow は root ユーザか 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> なんかがそれ。

先ほど確認した自分自身のパスワードハッシュでは id が 6 になっていた。 つまり、ハッシュアルゴリズムとして 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>]

上記を見て分かる通り id が 2a に対応するハッシュアルゴリズムの 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.c の pw_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プログラミングの王道

詳解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: 地域インターネットレジストリ) を表している。

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

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

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

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

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

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

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

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

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

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

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

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

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

/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版

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

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によるデータ分析入門 第2版 ―NumPy、pandasを使ったデータ処理

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

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