Python: scikit-learn のロジスティック回帰を使ってみる

最近、意外とロジスティック回帰が使われていることに気づいた。 もちろん世間にはもっと表現力のある分類器がたくさんあるけど、問題によってどれくらい複雑なモデルが適しているかは異なる。 それに、各特徴量がどのように働くか重みから確認したり、単純なモデルなのでスコアをベンチマークとして利用する、といった用途もあるらしい。 今回は、そんなロジスティック回帰を scikit-learn の実装で試してみる。

使った環境は次の通り。

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.13.6
BuildVersion:   17G65
$ python -V               
Python 3.6.6
$ pip list --format=columns | grep -i scikit-learn
scikit-learn    0.19.2 

下準備

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

$ pip install scikit-learn

乳がんデータセットをロジスティック回帰で分類してみる

以下にロジスティック回帰を使って乳がんデータセットを分類するサンプルコードを示す。 とはいえ scikit-learn は API が統一されているので、分類器がロジスティック回帰になってる以外に特筆すべき点はないかも。

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

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


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

    # ロジスティック回帰
    clf = LogisticRegression()
    # Stratified K-Fold CV で性能を評価する
    skf = StratifiedKFold(shuffle=True)
    scoring = {
        'acc': 'accuracy',
        'auc': 'roc_auc',
    }
    scores = cross_validate(clf, X, y, cv=skf, scoring=scoring)

    print('Accuracy (mean):', scores['test_acc'].mean())
    print('AUC (mean):', scores['test_auc'].mean())


if __name__ == '__main__':
    main()

実行してみよう。 だいたい精度 (Accuracy) の平均で 0.947 前後のスコアが得られた。

$ python logistic.py
Accuracy (mean): 0.9472570314675578
AUC (mean): 0.991659762496548

正直そんなに高くないけど、ロジスティック回帰くらい単純なモデルではこれくらいなんだなっていう指標にはなると思う。

ロジスティック回帰の利点

ロジスティック回帰の良いところは、モデルが単純で解釈も容易なところ。 例えば、基本的に線形モデルの眷属なので、各特徴量の重み (傾き) が確認できる。

実際に、それを体験してみよう。 今回例に挙げた乳がんデータセットは、腫瘍の特徴を記録した 30 次元の教師データだった。

>>> from sklearn import datasets
>>> dataset = datasets.load_breast_cancer()
>>> dataset.feature_names
array(['mean radius', 'mean texture', 'mean perimeter', 'mean area',
       'mean smoothness', 'mean compactness', 'mean concavity',
       'mean concave points', 'mean symmetry', 'mean fractal dimension',
       'radius error', 'texture error', 'perimeter error', 'area error',
       'smoothness error', 'compactness error', 'concavity error',
       'concave points error', 'symmetry error',
       'fractal dimension error', 'worst radius', 'worst texture',
       'worst perimeter', 'worst area', 'worst smoothness',
       'worst compactness', 'worst concavity', 'worst concave points',
       'worst symmetry', 'worst fractal dimension'], dtype='<U23')
>>> len(dataset.feature_names)
30

目的変数は 0 が良性で 1 が悪性となっている。

>>> import numpy as np
>>> np.unique(y)
array([0, 1])

ホールドアウト検証を使ってデータを分割したら、モデルを学習させよう。

>>> from sklearn.model_selection import train_test_split
>>> X, y = dataset.data, dataset.target
>>> X_train, X_test, y_train, y_test = train_test_split(X, y, shuffle=True, random_state=42)
>>> from sklearn.linear_model import LogisticRegression
>>> clf = LogisticRegression()
>>> clf.fit(X_train, y_train)
LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
          intercept_scaling=1, max_iter=100, multi_class='ovr', n_jobs=1,
          penalty='l2', random_state=None, solver='liblinear', tol=0.0001,
          verbose=0, warm_start=False)

すると、学習したモデルで切片 (LogisticRegression#intercept_) と重み (LogisticRegression#coef_) が確認できる。 重みは各特徴量ごとにある。

>>> clf.intercept_
array([0.40407439])
>>> clf.coef_
array([[ 2.18931343e+00,  1.51512837e-01, -1.57814199e-01,
        -1.03404299e-03, -1.29170075e-01, -4.23805008e-01,
        -6.47620520e-01, -3.37002545e-01, -1.97619418e-01,
        -3.23607668e-02, -6.88409834e-02,  1.48012177e+00,
         4.81243097e-02, -1.05177866e-01, -1.40690243e-02,
        -3.50323361e-02, -7.06715773e-02, -3.93587747e-02,
        -4.81468850e-02, -2.01238862e-03,  1.20675464e+00,
        -3.93262696e-01, -4.96613892e-02, -2.45385329e-02,
        -2.43248181e-01, -1.21314110e+00, -1.60969567e+00,
        -6.01906976e-01, -7.28573372e-01, -1.21974174e-01]])
>>> len(clf.coef_[0])
30

元のデータセットを標準化していないので、重みの値の大小については単純な比較が難しい。 ただ、特徴量がプラスであればその特徴量は悪性の方向に、反対にマイナスなら良性の方向に働く。

学習済みモデルの切片と重みから推論内容を確認する

ここからは学習済みモデルの切片と重みから計算した内容が推論と一致することを確認してみる。 ロジスティック回帰は線形回帰の式をシグモイド関数で 0 ~ 1 に変換したものになっている。

以下の式は線形回帰の式で、重みが w$w$ で切片が b$b$ に対応する。 ようするに特徴量と重みをかけて切片を足すだけ。

y=wX+b$y=wX+b$

上記をシグモイド関数に放り込むと値が 0 ~ 1 の範囲に収まる。

z=11+e−y$z=\frac{1}{1+e^{-y}}$

これがロジスティック回帰の出力となる。

実際に上記を学習済みモデルで確認してみよう。 例えば検証用データの先頭は悪性と判定されている。

>>> clf.predict([X_test[0]])
array([1])
>>> clf.predict_proba([X_test[0]])
array([[0.19165157, 0.80834843]])

悪性の確率は 0.80834843 になる。

正解を確認すると、たしかに悪性のようだ。

>>> y_test[0]
1

それでは学習済みモデルから計算した内容と上記が一致するかを確認してみよう。 まずは線形回帰の式を作る。 これが上記の y=wX+b$y=wX+b$ に対応する。

>>> import numpy as np
>>> y = np.sum(clf.coef_ * X_test[0]) + clf.intercept_
>>> y
array([1.4393142])

続いてシグモイド関数を定義しておく。

>>> def sigmoid(x):
...     return 1 / (1 + np.exp(-x))
... 

あとは、さきほど得られた結果をシグモイド関数に放り込むだけ。

>>> z = sigmoid(y)
>>> z
array([0.80834843])

結果は 0.80834843 となって、見事に先ほど得られた内容と一致している。

ばっちり。

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