Python: KMeans 法を実装してみる

KMeans 法は、機械学習における教師なし学習のクラスタリングという問題を解くためのアルゴリズム。 教師なし学習というのは、事前に教師データというヒントが与えられないことを指している。 その上で、クラスタリングというのは未知のデータに対していくつかのまとまりを作る問題をいう。

今回取り扱う KMeans 法は、比較的単純なアルゴリズムにも関わらず広く使われているものらしい。 実際に書いてみても、基本的な実装であればたしかにとてもシンプルだった。 ただし、データの初期化をするところで一点考慮すべき内容があることにも気づいたので、それについても書く。 KMeans 法の具体的なアルゴリズムについてはサンプルコードと共に後述する。

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

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.12.3
BuildVersion:   16D32
$ python --version
Python 3.5.3

依存パッケージをインストールする

あらかじめ、今回ソースコードで使う依存パッケージをインストールしておく。 グラフ描画ライブラリの matplotlib を Mac OS X で動かすには、pip でのインストール以外にもちょっとした設定が必要になる。

$ pip install scipy scikit-learn matplotlib
$ mkdir -p ~/.matplotlib
$ cat << 'EOF' > ~/.matplotlib/matplotlibrc
backend: TkAgg
EOF

まずは scikit-learn のお手本から

本来はお手本の例については後回しにしたいんだけど、今回は構成の都合で先に示しておく。 Python の機械学習系ライブラリである scikit-learn には KMeans 法の実装も用意されている。 自分で書いた実装をお披露目する前に、どんなものなのか見てもらいたい。

次のサンプルコードでは、ダミーのデータを生成した上でそれをクラスタリングしている。 scikit-learn にはデータセットを生成する機能も備わっている。 その中でも make_blobs() という関数は、いくつかのまとまりを持ったダミーデータを作ってくれる。 あとは、そのダミーデータを KMeans 法の実装である sklearn.cluster.KMeans で処理させている。 クラスタリングした処理結果は matplotlib を使って可視化した。

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

from matplotlib import pyplot as plt
from sklearn import datasets
from sklearn.cluster import KMeans


def main():
    # クラスタ数
    N_CLUSTERS = 5

    # Blob データを生成する
    dataset = datasets.make_blobs(centers=N_CLUSTERS)

    # 特徴データ
    features = dataset[0]
    # 正解ラベルは使わない
    # targets = dataset[1]

    # クラスタリングする
    cls = KMeans(n_clusters=N_CLUSTERS)
    pred = cls.fit_predict(features)

    # 各要素をラベルごとに色付けして表示する
    for i in range(N_CLUSTERS):
        labels = features[pred == i]
        plt.scatter(labels[:, 0], labels[:, 1])

    # クラスタのセントロイド (重心) を描く
    centers = cls.cluster_centers_
    plt.scatter(centers[:, 0], centers[:, 1], s=100,
                facecolors='none', edgecolors='black')

    plt.show()


if __name__ == '__main__':
    main()

上記のサンプルコードでは、クラスタ数として適当に 5 を選んだ。 ここで「クラスタ数を選ぶ」というのはポイントとなる。 KMeans 法では、あらかじめクラスタ数をハイパーパラメータとして指定するためだ。 つまり、データをいくつのまとまりとして捉えるかを人間が教えてやる必要がある。 クラスタリングの手法によっては、これをアルゴリズムがやってくれる場合もある。

それでは、上記のサンプルコードを実行してみよう。

$ python kmeans_scikit.py

すると、こんな感じのグラフが出力される。 それぞれの特徴ベクトルが所属するクラスタが、別々の色で表現されている。 また、それぞれのクラスタの真ん中らへんにある黒い円は、後述するセントロイドという概念を表している。 ちなみに、生成されるダミーデータは毎回異なるので、出力されるグラフも異なる。

自分で実装してみる

scikit-learn がクラスタリングしたお手本を見たところで、次は KMeans 法を自分で書いてみよう。

まず、KMeans 法ではセントロイド (重心) という概念が重要になる。 セントロイドというのは、文字通りクラスタの中心に置かれるものだ。 アルゴリズムの第一歩としては、このセントロイドを作ることになる。 セントロイドは、前述した通り各クラスタの中心に置かれる。 そのため、セントロイドはハイパーパラメータとして指定したクラスタ数だけ必要となる。 また、中心というのは、クラスタに属する特徴ベクトルの各次元でそれらの平均値にいる、ということ。

そして、データセットに含まれる特徴ベクトルは、必ず最寄りのセントロイドに属する。 ここでいう最寄りとは、特徴ベクトル同士のユークリッド距離が近いという意味だ。 属したセントロイドが、それぞれのクラスタを表している。

上記の基本的な概念を元に KMeans 法のアルゴリズムは次のようなステップに分けることができる

• 最初のセントロイドを何らかの方法で決める
• 特徴ベクトルを最寄りのセントロイドに所属させる
• 各クラスタごとにセントロイドを再計算する
• 上記 2, 3 番目の操作を繰り返す
• もしイテレーション前後で所属するセントロイドが変化しなければ、そこで処理を終了する

ここで、一番最初のステップである「最初のセントロイドを何らかの方法で決める」のがポイントだった。 例えば、次のページではあらかじめ各特徴ベクトルをランダムにクラスタリングしてから、セントロイドを決めるやり方を取っている。

tech.nitoyon.com

上記のやり方であれば、次のようなサンプルコードになる。 ここでは KMeans というクラスで KMeans 法を実装している。 インターフェースは、先ほどお手本として提示した scikit-learn と揃えた。 そのため main() 関数は基本的に変わっておらず、見るべきは KMeans クラスだけとなっている。

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

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


class KMeans(object):
    """KMeans 法でクラスタリングするクラス"""

    def __init__(self, n_clusters=2, max_iter=300):
        """コンストラクタ

        Args:
            n_clusters (int): クラスタ数
            max_iter (int): 最大イテレーション数
        """
        self.n_clusters = n_clusters
        self.max_iter = max_iter

        self.cluster_centers_ = None

    def fit_predict(self, features):
        """クラスタリングを実施する

        Args:
            features (numpy.ndarray): ラベル付けするデータ

        Returns:
            numpy.ndarray: ラベルデータ
        """
        # 初期データは各要素に対してランダムにラベルをつける
        pred = np.random.randint(0, self.n_clusters, len(features))

        # クラスタリングをアップデートする
        for _ in range(self.max_iter):

            # 各クラスタごとにセントロイド (重心) を計算する
            self.cluster_centers_ = np.array([features[pred == i].mean(axis=0)
                                              for i in range(self.n_clusters)])

            # 各特徴ベクトルから最短距離となるセントロイドを基準に新しいラベルをつける
            new_pred = np.array([
                np.array([
                    self._euclidean_distance(p, centroid)
                    for centroid in self.cluster_centers_
                ]).argmin()
                for p in features
            ])

            if np.all(new_pred == pred):
                # 更新前と内容を比較して、もし同じなら終了
                break

            pred = new_pred

        return pred

    def _euclidean_distance(self, p0, p1):
        return np.sum((p0 - p1) ** 2)


def main():
    # クラスタ数
    N_CLUSTERS = 5

    # Blob データを生成する
    dataset = datasets.make_blobs(centers=N_CLUSTERS)

    # 特徴データ
    features = dataset[0]
    # 正解ラベルは使わない
    # targets = dataset[1]

    # クラスタリングする
    cls = KMeans(n_clusters=N_CLUSTERS)
    pred = cls.fit_predict(features)

    # 各要素をラベルごとに色付けして表示する
    for i in range(N_CLUSTERS):
        labels = features[pred == i]
        plt.scatter(labels[:, 0], labels[:, 1])

    centers = cls.cluster_centers_
    plt.scatter(centers[:, 0], centers[:, 1], s=100,
                facecolors='none', edgecolors='black')

    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python kmeans_random.py

実行すると、次のようなグラフが表示される。 一見すると、上手くいっているようだ。

しかし、何度か実行すると、次のようなエラーが出ることがある。

$ python kmeans_random.py 
kmeans_random.py:41: RuntimeWarning: Mean of empty slice.
  for i in range(self.n_clusters)])
/Users/amedama/.virtualenvs/kmeans/lib/python3.5/site-packages/numpy/core/_methods.py:73: RuntimeWarning: invalid value encountered in true_divide
  ret, rcount, out=ret, casting='unsafe', subok=False)

表示されるグラフは、次のようなものになる。 何やら、ぜんぜん上手いことクラスタリングできていない。

上記が起こる原因について調べたところ、セントロイドの初期化に問題があることが分かった。 初期化時に特徴ベクトルに対してランダムにラベルをつけるやり方では、一つの特徴ベクトルも属さないセントロイドが生じる恐れがあるためだ。

ランダムにラベルをつけるやり方では、セントロイドはおのずとデータセット全体の平均あたりに生じやすくなる。 場合によっては、セントロイドが別のセントロイドに囲まれるなどして、一つの特徴ベクトルも属さないセントロイドがでてくる。 サンプルコードでは、そのような事態を想定していなかったのでセントロイドが計算できなくなって上手く動作しなかった。

先ほどの KMeans 法を図示しているページでも、特徴ベクトルの数 N に対してクラスタの数 k を増やして実行してみよう。 一つの特徴ベクトルも属さないクラスタが生じることが分かる。

初期化を改良してみる

一つの特徴ベクトルも属さないクラスタが生じることについて、問題がないと捉えることもできるかもしれない。 その場合には、一つも属さないクラスタのセントロイドを、前回の位置から動かなさないようにすることでケアできそうだ。

とはいえ、ここでは別のアプローチで問題を解決してみることにする。 それというのは、初期のセントロイドをデータセットに含まれるランダムな特徴ベクトルにする、というものだ。 つまり、最初に作るセントロイドの持つ特徴ベクトルがランダムに選んだ要素の特徴ベクトルと一致する。 こうすれば各セントロイドは、少なくとも一つ以上の特徴ベクトルが属することは確定できる。 ユークリッド距離がゼロになる特徴ベクトルが、必ず一つは存在するためだ。

次のサンプルコードでは、初期のセントロイドの作り方だけを先ほどと変更している。 ここでポイントとなるのは、最初にセントロイドとして採用する特徴ベクトルは重複しないように決める必要があるということだ。 そこで、初期のセントロイドは特徴ベクトルをシャッフルした上で先頭から k 個を取り出して決めている。

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

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


class KMeans(object):
    """KMeans 法でクラスタリングするクラス"""

    def __init__(self, n_clusters=2, max_iter=300):
        """コンストラクタ

        Args:
            n_clusters (int): クラスタ数
            max_iter (int): 最大イテレーション数
        """
        self.n_clusters = n_clusters
        self.max_iter = max_iter

        self.cluster_centers_ = None

    def fit_predict(self, features):
        """クラスタリングを実施する

        Args:
            features (numpy.ndarray): ラベル付けするデータ

        Returns:
            numpy.ndarray: ラベルデータ
        """
        # 要素の中からセントロイド (重心) の初期値となる候補をクラスタ数だけ選び出す
        feature_indexes = np.arange(len(features))
        np.random.shuffle(feature_indexes)
        initial_centroid_indexes = feature_indexes[:self.n_clusters]
        self.cluster_centers_ = features[initial_centroid_indexes]

        # ラベル付けした結果となる配列はゼロで初期化しておく
        pred = np.zeros(features.shape)

        # クラスタリングをアップデートする
        for _ in range(self.max_iter):
            # 各要素から最短距離のセントロイドを基準にラベルを更新する
            new_pred = np.array([
                np.array([
                    self._euclidean_distance(p, centroid)
                    for centroid in self.cluster_centers_
                ]).argmin()
                for p in features
            ])

            if np.all(new_pred == pred):
                # 更新前と内容が同じなら終了
                break

            pred = new_pred

            # 各クラスタごとにセントロイド (重心) を再計算する
            self.cluster_centers_ = np.array([features[pred == i].mean(axis=0)
                                              for i in range(self.n_clusters)])

        return pred

    def _euclidean_distance(self, p0, p1):
        return np.sum((p0 - p1) ** 2)


def main():
    # クラスタ数
    N_CLUSTERS = 5

    # Blob データを生成する
    dataset = datasets.make_blobs(centers=N_CLUSTERS)

    # 特徴データ
    features = dataset[0]
    # 正解ラベルは使わない
    # targets = dataset[1]

    # クラスタリングする
    cls = KMeans(n_clusters=N_CLUSTERS)
    pred = cls.fit_predict(features)

    # 各要素をラベルごとに色付けして表示する
    for i in range(N_CLUSTERS):
        labels = features[pred == i]
        plt.scatter(labels[:, 0], labels[:, 1])

    centers = cls.cluster_centers_
    plt.scatter(centers[:, 0], centers[:, 1], s=100,
                facecolors='none', edgecolors='black')

    plt.show()


if __name__ == '__main__':
    main()

それでは、上記のサンプルコードを実行してみよう。

$ python kmeans_select.py

これまでと変わらないけど、次のようなグラフが表示されるはず。

今度は N_CLUSTERS を増やしたり、何度やっても先ほどのようなエラーにはならない。

まとめ

今回は教師なし学習のクラスタリングという問題を解くアルゴリズムの KMeans 法を実装してみた。 KMeans 法はシンプルなアルゴリズムだけど、セントロイドをどう初期化するかは流派があるみたい。

はじめてのパターン認識

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

Python: データセットを標準化する効果を最近傍法で確かめる

データセットの標準化については、このブログでも何回か扱っている。 しかし、実際にデータセットを標準化したときの例については試していなかった。

blog.amedama.jp

blog.amedama.jp

そこで、今回は UCI の提供する小麦 (seeds) データセットを最近傍法で分類したときに、スコアが上がる様を見てみたいと思う。 あらかじめ、いくつかの説明変数が教師信号として与えられるので、そこから小麦の品種を分類 (Classification) する。

UCI Machine Learning Repository: seeds Data Set

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

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.12.3
BuildVersion:   16D32
$ python --version
Python 3.5.3

下準備

まずは今回使う Python のパッケージを pip でインストールする。 どれも科学計算系で定番のやつ。

$ pip install numpy pandas scipy scikit-learn

データセットを標準化しない場合

今回使う最近傍法というアルゴリズムでは、分類したい点から最も近くにあるデータの種別を使って分類する。 近さにはユークリッド距離を使うため、データセットの説明変数の大きさや単位に影響を受けやすい。 例えば説明変数の中に 1,000m ~ 10,000m を取る次元と 0.1cm ~ 1cm を取る次元があるとしよう。 当然ながら、ユークリッド距離を計算するとそのままでは前者の次元の影響が大きくなってしまう。 今回使うデータセットでも原理的には同じ問題が発生する。

次のサンプルコードでは、データセットを標準化しない状態で最近傍法を使った分類を実施している。 そして Leave-One-Out 法を使って、計測した汎化性能を出力する。

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

import pandas as pd
from sklearn.model_selection import LeaveOneOut
from sklearn.metrics import accuracy_score
from sklearn.neighbors import KNeighborsClassifier


def main():
    # 小麦データセットをダウンロードする
    dataset_url = 'https://archive.ics.uci.edu/ml/machine-learning-databases/00236/seeds_dataset.txt'  # noqa
    df = pd.read_csv(dataset_url, delim_whitespace=True, header=None)

    # データフレームから説明変数と目的変数を取り出す
    features = df.loc[:, :6].get_values()
    targets = df.loc[:, 7].get_values()

    # 予測した結果の入れ物
    predicted_labels = []

    # LOO で交差検証する
    loo = LeaveOneOut()
    for train, test in loo.split(features):
        train_data = features[train]
        target_data = targets[train]

        # k-NN 法を使う
        model = KNeighborsClassifier(n_neighbors=1)
        # 訓練データを学習させる
        model.fit(train_data, target_data)
        # テストデータを予測させる
        predicted_label = model.predict(features[test])
        # 予測した結果を追加する
        predicted_labels.append(predicted_label)

    # 正解率を出力する
    score = accuracy_score(targets, predicted_labels)
    print(score)


if __name__ == '__main__':
    main()

上記のサンプルコードの実行結果は次の通り。 データセットを標準化しない状態では約 90.5% の汎化性能が得られた。

$ python withoutnorm.py
0.904761904762

データセットを標準化する場合

それでは、次はデータセットを標準化する場合を試してみよう。 標準化されたデータセットでは、説明変数の各次元の値が平均 0 で標準偏差が 1 になる。 つまり、元々の単位や大きさは無くなってそれぞれの値の間隔の比率だけが残されることになる。

次のサンプルコードでは、先ほどとデータセットを標準化するところだけ変更している。 データセットの標準化には scikit-learn に用意されている StandardScaler を用いた。

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

import pandas as pd
from sklearn.model_selection import LeaveOneOut
from sklearn.metrics import accuracy_score
from sklearn.neighbors import KNeighborsClassifier
from sklearn.preprocessing import StandardScaler


def main():
    dataset_url = 'https://archive.ics.uci.edu/ml/machine-learning-databases/00236/seeds_dataset.txt'  # noqa
    df = pd.read_csv(dataset_url, delim_whitespace=True, header=None)

    features = df.loc[:, :6].get_values()
    targets = df.loc[:, 7].get_values()

    # データセットを Z-Score に標準化する
    sc = StandardScaler()
    sc.fit(features)
    normalized_features = sc.transform(features)

    predicted_labels = []

    loo = LeaveOneOut()
    for train, test in loo.split(normalized_features):
        train_data = normalized_features[train]
        target_data = targets[train]

        model = KNeighborsClassifier(n_neighbors=1)
        model.fit(train_data, target_data)

        predicted_label = model.predict(normalized_features[test])
        predicted_labels.append(predicted_label)

    score = accuracy_score(targets, predicted_labels)
    print(score)


if __name__ == '__main__':
    main()

上記の実行結果は次の通り。 今度は汎化性能が約 93.8% に上昇している。 データセットを標準化するだけで分類精度が 3.3% 上がった。

$ python withnorm.py
0.938095238095

まとめ

データセットを標準化することで、使う機械学習アルゴリズムによっては汎化性能を上げることができる。

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

統計: Python と R で重回帰分析してみる

今回は R と Python の両方を使って重回帰分析をしてみる。 モチベーションとしては、できるだけ手に慣れた Python を使って分析をしていきたいという気持ちがある。 ただ、計算結果が意図通りのものになっているのかを R の結果と見比べて確かめておきたい。

また、分析にはボストンデータセットを用いる。 このデータセットはボストンの各地区ごとの不動産の平均価格と、それに付随する情報が入っている。

今回使った環境は次の通り。 R と Python は、あらかじめインストール済みであることを想定している。

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.11.6
BuildVersion:   15G1212
$ python --version
Python 3.5.2
$ R --version
R version 3.3.2 (2016-10-31) -- "Sincere Pumpkin Patch"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin15.6.0 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under the terms of the
GNU General Public License versions 2 or 3.
For more information about these matters see
http://www.gnu.org/licenses/.

もし、まだ R をインストールしていないのであれば、こちらの記事が参考になるかも。

blog.amedama.jp

R

まずは R のやり方から。

最初に R を起動する。

$ R

続いてボストンデータセットの入っているライブラリ MASS を読み込む。

> library(MASS)

データセットの中身は次のようになっている。

> head(Boston)
     crim zn indus chas   nox    rm  age    dis rad tax ptratio  black lstat
1 0.00632 18  2.31    0 0.538 6.575 65.2 4.0900   1 296    15.3 396.90  4.98
2 0.02731  0  7.07    0 0.469 6.421 78.9 4.9671   2 242    17.8 396.90  9.14
3 0.02729  0  7.07    0 0.469 7.185 61.1 4.9671   2 242    17.8 392.83  4.03
4 0.03237  0  2.18    0 0.458 6.998 45.8 6.0622   3 222    18.7 394.63  2.94
5 0.06905  0  2.18    0 0.458 7.147 54.2 6.0622   3 222    18.7 396.90  5.33
6 0.02985  0  2.18    0 0.458 6.430 58.7 6.0622   3 222    18.7 394.12  5.21
  medv
1 24.0
2 21.6
3 34.7
4 33.4
5 36.2
6 28.7

各変数は次のような意味を持っている。 ただ、今回はボストンデータセットについて詳しく見ることが目的ではないので詳しくは立ち入らない。

1. CRIM: per capita crime rate by town
2. ZN: proportion of residential land zoned for lots over 25,000 sq.ft.
3. INDUS: proportion of non-retail business acres per town
4. CHAS: Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
5. NOX: nitric oxides concentration (parts per 10 million)
6. RM: average number of rooms per dwelling
7. AGE: proportion of owner-occupied units built prior to 1940
8. DIS: weighted distances to five Boston employment centres
9. RAD: index of accessibility to radial highways
10. TAX: full-value property-tax rate per $10,000
11. PTRATIO: pupil-teacher ratio by town
12. B: 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
13. LSTAT: % lower status of the population
14. MEDV: Median value of owner-occupied homes in $1000's
(https://archive.ics.uci.edu/ml/datasets/Housing より)

今回はシンプルに線形回帰モデルを使うので lm() 関数を使う。 medv は目的変数である「その地区の不動産の平均価格」になっている。 次の処理では、目的関数をその他の全ての変数を使って重回帰している。

> lmfit = lm(medv ~ ., data=Boston)

重回帰した結果の詳細は summary() 関数で取得できる。

> summary(lmfit)

Call:
lm(formula = medv ~ ., data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max
-15.595  -2.730  -0.518   1.777  26.199

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)  3.646e+01  5.103e+00   7.144 3.28e-12 ***
crim        -1.080e-01  3.286e-02  -3.287 0.001087 **
zn           4.642e-02  1.373e-02   3.382 0.000778 ***
indus        2.056e-02  6.150e-02   0.334 0.738288
chas         2.687e+00  8.616e-01   3.118 0.001925 **
nox         -1.777e+01  3.820e+00  -4.651 4.25e-06 ***
rm           3.810e+00  4.179e-01   9.116  < 2e-16 ***
age          6.922e-04  1.321e-02   0.052 0.958229
dis         -1.476e+00  1.995e-01  -7.398 6.01e-13 ***
rad          3.060e-01  6.635e-02   4.613 5.07e-06 ***
tax         -1.233e-02  3.760e-03  -3.280 0.001112 **
ptratio     -9.527e-01  1.308e-01  -7.283 1.31e-12 ***
black        9.312e-03  2.686e-03   3.467 0.000573 ***
lstat       -5.248e-01  5.072e-02 -10.347  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.745 on 492 degrees of freedom
Multiple R-squared:  0.7406, Adjusted R-squared:  0.7338
F-statistic: 108.1 on 13 and 492 DF,  p-value: < 2.2e-16

自由度調整済み決定係数が 0.7338 ということで、あまり上手く説明できていなさそうな感じ。

Python

続いて Python でのやり方について。 今回は statsmodels というライブラリを使うことにする。 ちなみに、線形回帰モデルは別のライブラリにも色々な実装がある。 ただ、フィッティングしたときに得られるパラメータが statsmodels は細かく取れて R に近いので、今回はこちらを使ってみた。

まずは必要なライブラリをインストールする。 今回は回帰に使わない scikit-learn を入れているのはボストンデータセットを読み込むため。

$ pip install statsmodels scikit-learn

最初に Python の REPL を起動しよう。

$ python

起動したら scikit-learn からボストンデータセットを読み込む。

>>> from sklearn import datasets
>>> dataset = datasets.load_boston()

このデータセットには次のようなメンバがある。

>>> dir(dataset)
['DESCR', 'data', 'feature_names', 'target']

target メンバが目的変数となる各地区ごとの不動産物件の平均価格となっている。

>>> dataset.target[:10]
array([ 24. ,  21.6,  34.7,  33.4,  36.2,  28.7,  22.9,  27.1,  16.5,  18.9])

そして説明変数が data に入っている。 ボストンを 506 地区に分割して 13 の変数を与えたデータになっている。

>>> dataset.data.shape
(506, 13)

変数名は feature_names から取得できる。

>>> dataset.feature_names
array(['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD',
       'TAX', 'PTRATIO', 'B', 'LSTAT'],
      dtype='<U7')

データセットを読み込めたところで statsmodels をインポートしよう。

>>> from statsmodels import api as sm

線形回帰モデルを使うときは OLS クラスを使う。 第一引数に目的変数を、第二変数に説明変数を渡す。

>>> model = sm.OLS(dataset.target, dataset.data)

準備ができたら fit() メソッドを使って分析する。

>>> result = model.fit()

分析結果からは summary() メソッドを使って R の summary() 関数に近い内容が取得できる。

>>> result.summary()
<class 'statsmodels.iolib.summary.Summary'>
"""
                            OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.959
Model:                            OLS   Adj. R-squared:                  0.958
Method:                 Least Squares   F-statistic:                     891.1
Date:                Thu, 22 Dec 2016   Prob (F-statistic):               0.00
Time:                        22:59:04   Log-Likelihood:                -1523.8
No. Observations:                 506   AIC:                             3074.
Df Residuals:                     493   BIC:                             3129.
Df Model:                          13
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
x1            -0.0916      0.034     -2.675      0.008        -0.159    -0.024
x2             0.0487      0.014      3.379      0.001         0.020     0.077
x3            -0.0038      0.064     -0.059      0.953        -0.130     0.123
x4             2.8564      0.904      3.160      0.002         1.080     4.633
x5            -2.8808      3.359     -0.858      0.392        -9.481     3.720
x6             5.9252      0.309     19.168      0.000         5.318     6.533
x7            -0.0072      0.014     -0.523      0.601        -0.034     0.020
x8            -0.9680      0.196     -4.947      0.000        -1.352    -0.584
x9             0.1704      0.067      2.554      0.011         0.039     0.302
x10           -0.0094      0.004     -2.393      0.017        -0.017    -0.002
x11           -0.3924      0.110     -3.571      0.000        -0.608    -0.177
x12            0.0150      0.003      5.561      0.000         0.010     0.020
x13           -0.4170      0.051     -8.214      0.000        -0.517    -0.317
==============================================================================
Omnibus:                      204.050   Durbin-Watson:                   0.999
Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1372.527
Skew:                           1.609   Prob(JB):                    9.11e-299
Kurtosis:                      10.399   Cond. No.                     8.50e+03
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 8.5e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
"""

めでたしめでたし、と行きたいところだけどちょっと待ってほしい。 先ほど取得した R の結果と内容が異なっている。

実は R の分析では説明変数にバイアス項という変数が自動で挿入されていた。 これは、不動産物件のベースとなる価格と捉えることができる。

statsmodels を使っているときにバイアス項を説明変数に挿入するときは statsmodels.add_constant() 関数を使う。 バイアス項を挿入した説明変数を X という名前で保存しておこう。 ついでに目的変数も Y という変数で保存しておく。

>>> X = sm.add_constant(dataset.data)
>>> Y = dataset.target

バイアス項を追加した説明変数を用いて、先ほどと同じように回帰してみよう。

>>> model = sm.OLS(Y, X)
>>> result = model.fit()
>>> result.summary()
<class 'statsmodels.iolib.summary.Summary'>
"""
                            OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.741
Model:                            OLS   Adj. R-squared:                  0.734
Method:                 Least Squares   F-statistic:                     108.1
Date:                Thu, 22 Dec 2016   Prob (F-statistic):          6.95e-135
Time:                        23:00:34   Log-Likelihood:                -1498.8
No. Observations:                 506   AIC:                             3026.
Df Residuals:                     492   BIC:                             3085.
Df Model:                          13
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const         36.4911      5.104      7.149      0.000        26.462    46.520
x1            -0.1072      0.033     -3.276      0.001        -0.171    -0.043
x2             0.0464      0.014      3.380      0.001         0.019     0.073
x3             0.0209      0.061      0.339      0.735        -0.100     0.142
x4             2.6886      0.862      3.120      0.002         0.996     4.381
x5           -17.7958      3.821     -4.658      0.000       -25.302   -10.289
x6             3.8048      0.418      9.102      0.000         2.983     4.626
x7             0.0008      0.013      0.057      0.955        -0.025     0.027
x8            -1.4758      0.199     -7.398      0.000        -1.868    -1.084
x9             0.3057      0.066      4.608      0.000         0.175     0.436
x10           -0.0123      0.004     -3.278      0.001        -0.020    -0.005
x11           -0.9535      0.131     -7.287      0.000        -1.211    -0.696
x12            0.0094      0.003      3.500      0.001         0.004     0.015
x13           -0.5255      0.051    -10.366      0.000        -0.625    -0.426
==============================================================================
Omnibus:                      178.029   Durbin-Watson:                   1.078
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              782.015
Skew:                           1.521   Prob(JB):                    1.54e-170
Kurtosis:                       8.276   Cond. No.                     1.51e+04
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.51e+04. This might indicate that there are
strong multicollinearity or other numerical problems.
"""

今度は R で得られた内容と一致している！

AIC

重回帰分析では、モデリングに用いる説明変数の選択基準として AIC という数値を小さくするやり方があるので、その計算方法について。

R

R では AIC() 関数が用意されているので、それに回帰結果を渡す。

> AIC(lmfit)
[1] 3027.609

Python

Python では statsmodels の回帰結果のインスタンスに aic というメンバがあって、そこに格納されている。 ちなみに AIC はサマリーの中にも表示されていた。

>>> result.aic
3025.6767200074596

標準化

標準化というのはデータセットを平均が 0 で標準偏差が 1 になるように加工すること。 説明変数が標準化されていると、それぞれの説明変数が目的変数に対してどれくらい影響を与えているかを見たりするのに分かりやすい。

R

R には標準化のための関数として scale() が用意されている。

> X <- scale(Boston[, -14])
> Y <- Boston[, 14]

説明変数を標準化して回帰してみよう。

> lmfit <- lm(Y ~ X)
> summary(lmfit)

Call:
lm(formula = Y ~ X)

Residuals:
    Min      1Q  Median      3Q     Max
-15.595  -2.730  -0.518   1.777  26.199

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.53281    0.21095 106.814  < 2e-16 ***
Xcrim       -0.92906    0.28269  -3.287 0.001087 **
Xzn          1.08264    0.32016   3.382 0.000778 ***
Xindus       0.14104    0.42188   0.334 0.738288
Xchas        0.68241    0.21884   3.118 0.001925 **
Xnox        -2.05875    0.44262  -4.651 4.25e-06 ***
Xrm          2.67688    0.29364   9.116  < 2e-16 ***
Xage         0.01949    0.37184   0.052 0.958229
Xdis        -3.10712    0.41999  -7.398 6.01e-13 ***
Xrad         2.66485    0.57770   4.613 5.07e-06 ***
Xtax        -2.07884    0.63379  -3.280 0.001112 **
Xptratio    -2.06265    0.28323  -7.283 1.31e-12 ***
Xblack       0.85011    0.24521   3.467 0.000573 ***
Xlstat      -3.74733    0.36216 -10.347  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.745 on 492 degrees of freedom
Multiple R-squared:  0.7406, Adjusted R-squared:  0.7338
F-statistic: 108.1 on 13 and 492 DF,  p-value: < 2.2e-16

説明変数が標準化されたことで、どの変数がどれくらい影響を与えているのかが分かりやすくなった。 例えば lstat などは影響が大きそうだ。

Python

Python の場合は scikit-learn の preprocessing() 関数を使うと簡単に標準化できる。

>>> X = dataset.data
>>> from sklearn import preprocessing
>>> X_scaled = preprocessing.scale(X)

標準化した上でバイアス項を追加して回帰してみよう。

>>> X_scaled_with_constant = sm.add_constant(X_scaled)
>>> model = sm.OLS(Y, X_scaled_with_constant)
>>> result = model.fit()
>>> result.summary()
<class 'statsmodels.iolib.summary.Summary'>
"""
                            OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.741
Model:                            OLS   Adj. R-squared:                  0.734
Method:                 Least Squares   F-statistic:                     108.1
Date:                Thu, 22 Dec 2016   Prob (F-statistic):          6.95e-135
Time:                        23:50:36   Log-Likelihood:                -1498.8
No. Observations:                 506   AIC:                             3026.
Df Residuals:                     492   BIC:                             3085.
Df Model:                          13
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const         22.5328      0.211    106.807      0.000        22.118    22.947
x1            -0.9204      0.281     -3.276      0.001        -1.472    -0.368
x2             1.0810      0.320      3.380      0.001         0.453     1.709
x3             0.1430      0.421      0.339      0.735        -0.685     0.971
x4             0.6822      0.219      3.120      0.002         0.253     1.112
x5            -2.0601      0.442     -4.658      0.000        -2.929    -1.191
x6             2.6706      0.293      9.102      0.000         2.094     3.247
x7             0.0211      0.371      0.057      0.955        -0.709     0.751
x8            -3.1044      0.420     -7.398      0.000        -3.929    -2.280
x9             2.6588      0.577      4.608      0.000         1.525     3.792
x10           -2.0759      0.633     -3.278      0.001        -3.320    -0.832
x11           -2.0622      0.283     -7.287      0.000        -2.618    -1.506
x12            0.8566      0.245      3.500      0.001         0.376     1.338
x13           -3.7487      0.362    -10.366      0.000        -4.459    -3.038
==============================================================================
Omnibus:                      178.029   Durbin-Watson:                   1.078
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              782.015
Skew:                           1.521   Prob(JB):                    1.54e-170
Kurtosis:                       8.276   Cond. No.                         9.82
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
"""

結果が R の内容と一致していることが分かる。

ちなみに標準化の操作自体は単純なので、自前でやることもできる。 数式に直すと、こんな感じ。 各要素から平均 () を引いて標準偏差 () で割る。

コードにするとこんな感じ。

>>> X_scaled = (X - X.mean(axis=0)) / X.std(axis=0)

交互作用モデル

交互作用モデルというのは、複数の変数をまとめて説明変数として扱うやり方。 ある説明変数と別の説明変数の相乗作用があるときに使うと良い。

R

R では lm() 関数を実行するときに説明変数を ^2 するだけで交互作用モデルを使える。

> lmfit = lm(medv ~ .^2, data=Boston)
> summary(lmfit)

Call:
lm(formula = medv ~ .^2, data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max
-7.9374 -1.5344 -0.1068  1.2973 17.8500

Coefficients:
                Estimate Std. Error t value Pr(>|t|)
(Intercept)   -1.579e+02  6.800e+01  -2.323 0.020683 *
crim          -1.707e+01  6.554e+00  -2.605 0.009526 **
zn            -7.529e-02  4.580e-01  -0.164 0.869508
indus         -2.819e+00  1.696e+00  -1.663 0.097111 .
chas           4.451e+01  1.952e+01   2.280 0.023123 *
nox            2.006e+01  7.516e+01   0.267 0.789717
rm             2.527e+01  5.699e+00   4.435 1.18e-05 ***
age            1.263e+00  2.728e-01   4.630 4.90e-06 ***
dis           -1.698e+00  4.604e+00  -0.369 0.712395
rad            1.861e+00  2.464e+00   0.755 0.450532
tax            3.670e-02  1.440e-01   0.255 0.798978
ptratio        2.725e+00  2.850e+00   0.956 0.339567
black          9.942e-02  7.468e-02   1.331 0.183833
lstat          1.656e+00  8.533e-01   1.940 0.053032 .
crim:zn        4.144e-01  1.804e-01   2.297 0.022128 *
crim:indus    -4.693e-02  4.480e-01  -0.105 0.916621
crim:chas      2.428e+00  5.710e-01   4.251 2.63e-05 ***
crim:nox      -1.108e+00  9.285e-01  -1.193 0.233425
crim:rm        2.163e-01  4.907e-02   4.409 1.33e-05 ***
crim:age      -3.083e-03  3.781e-03  -0.815 0.415315
crim:dis      -1.903e-01  1.060e-01  -1.795 0.073307 .
crim:rad      -6.584e-01  5.815e-01  -1.132 0.258198
crim:tax       3.479e-02  4.287e-02   0.812 0.417453
crim:ptratio   4.915e-01  3.328e-01   1.477 0.140476
crim:black    -4.612e-04  1.793e-04  -2.572 0.010451 *
crim:lstat     2.964e-02  6.544e-03   4.530 7.72e-06 ***
...(省略)...
rad:tax        3.131e-05  1.446e-03   0.022 0.982729
rad:ptratio   -4.379e-02  8.392e-02  -0.522 0.602121
rad:black     -4.362e-04  2.518e-03  -0.173 0.862561
rad:lstat     -2.529e-02  1.816e-02  -1.392 0.164530
tax:ptratio    7.854e-03  2.504e-03   3.137 0.001830 **
tax:black     -4.785e-07  1.999e-04  -0.002 0.998091
tax:lstat     -1.403e-03  1.208e-03  -1.162 0.245940
ptratio:black  1.203e-03  3.361e-03   0.358 0.720508
ptratio:lstat  3.901e-03  2.985e-02   0.131 0.896068
black:lstat   -6.118e-04  4.157e-04  -1.472 0.141837
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.852 on 414 degrees of freedom
Multiple R-squared:  0.9212, Adjusted R-squared:  0.9039
F-statistic: 53.18 on 91 and 414 DF,  p-value: < 2.2e-16

結果から説明変数の数が一気に増えてることが分かる。 これは 13 種類の説明変数から 2 種類を組み合わせで取り出して、それらを交互作用モデルの説明変数として追加しているため。

また、自由度調整済み決定係数が 0.9039 となっており、ただ 13 種類の変数を使うときよりも上手くモデリングできていることが分かる。 まあ、この値だけを見ていても過適合の恐れがあるけど。

Python

statsmodels には回帰するときのやり方を R のような書き方 (R-style らしい) で指定できる。 このやり方を使うとさくっと相互作用モデルも使えそうな感じ。

Fitting models using R-style formulas — statsmodels 0.7.0 documentation

ただまあ、これを使うよりも自前で用意する方が楽しそう。 なので、今回は自分で計算してみることにする。

Python で組み合わせを計算するときは itertools モジュールを使うと良い。 この中にある combinations() 関数を使えば簡単に組み合わせが取れる。 例えば説明変数の組み合わせを計算したいなら、こう。

>>> import itertools
>>> list(itertools.combinations(dataset.feature_names, 2))
[('CRIM', 'ZN'), ('CRIM', 'INDUS'), ('CRIM', 'CHAS'), ('CRIM', 'NOX'), ('CRIM', 'RM'), ('CRIM', 'AGE'), ('CRIM', 'DIS'), ('CRIM', 'RAD'), ('CRIM', 'TAX'), ('CRIM', 'PTRATIO'), ('CRIM', 'B'), ('CRIM', 'LSTAT'), ('ZN', 'INDUS'), ('ZN', 'CHAS'), ('ZN', 'NOX'), ('ZN', 'RM'), ('ZN', 'AGE'), ('ZN', 'DIS'), ('ZN', 'RAD'), ('ZN', 'TAX'), ('ZN', 'PTRATIO'), ('ZN', 'B'), ('ZN', 'LSTAT'), ('INDUS', 'CHAS'), ('INDUS', 'NOX'), ('INDUS', 'RM'), ('INDUS', 'AGE'), ('INDUS', 'DIS'), ('INDUS', 'RAD'), ('INDUS', 'TAX'), ('INDUS', 'PTRATIO'), ('INDUS', 'B'), ('INDUS', 'LSTAT'), ('CHAS', 'NOX'), ('CHAS', 'RM'), ('CHAS', 'AGE'), ('CHAS', 'DIS'), ('CHAS', 'RAD'), ('CHAS', 'TAX'), ('CHAS', 'PTRATIO'), ('CHAS', 'B'), ('CHAS', 'LSTAT'), ('NOX', 'RM'), ('NOX', 'AGE'), ('NOX', 'DIS'), ('NOX', 'RAD'), ('NOX', 'TAX'), ('NOX', 'PTRATIO'), ('NOX', 'B'), ('NOX', 'LSTAT'), ('RM', 'AGE'), ('RM', 'DIS'), ('RM', 'RAD'), ('RM', 'TAX'), ('RM', 'PTRATIO'), ('RM', 'B'), ('RM', 'LSTAT'), ('AGE', 'DIS'), ('AGE', 'RAD'), ('AGE', 'TAX'), ('AGE', 'PTRATIO'), ('AGE', 'B'), ('AGE', 'LSTAT'), ('DIS', 'RAD'), ('DIS', 'TAX'), ('DIS', 'PTRATIO'), ('DIS', 'B'), ('DIS', 'LSTAT'), ('RAD', 'TAX'), ('RAD', 'PTRATIO'), ('RAD', 'B'), ('RAD', 'LSTAT'), ('TAX', 'PTRATIO'), ('TAX', 'B'), ('TAX', 'LSTAT'), ('PTRATIO', 'B'), ('PTRATIO', 'LSTAT'), ('B', 'LSTAT')]

組み合わせは 通り。

>>> len(list(itertools.combinations(dataset.feature_names, 2)))
78

それぞれの説明変数を取り出すための添字にするとこんな感じ。

>>> list(itertools.combinations(range(dataset.data.shape[1]), 2))
[(0, 1), (0, 2), (0, 3), (0, 4), (0, 5), (0, 6), (0, 7), (0, 8), (0, 9), (0, 10), (0, 11), (0, 12), (1, 2), (1, 3), (1, 4), (1, 5), (1, 6), (1, 7), (1, 8), (1, 9), (1, 10), (1, 11), (1, 12), (2, 3), (2, 4), (2, 5), (2, 6), (2, 7), (2, 8), (2, 9), (2, 10), (2, 11), (2, 12), (3, 4), (3, 5), (3, 6), (3, 7), (3, 8), (3, 9), (3, 10), (3, 11), (3, 12), (4, 5), (4, 6), (4, 7), (4, 8), (4, 9), (4, 10), (4, 11), (4, 12), (5, 6), (5, 7), (5, 8), (5, 9), (5, 10), (5, 11), (5, 12), (6, 7), (6, 8), (6, 9), (6, 10), (6, 11), (6, 12), (7, 8), (7, 9), (7, 10), (7, 11), (7, 12), (8, 9), (8, 10), (8, 11), (8, 12), (9, 10), (9, 11), (9, 12), (10, 11), (10, 12), (11, 12)]

交互作用モデルで使うそれぞれの説明変数は、説明変数同士を乗算することで計算できる。

>>> b = np.array([X[:, 0] * X[:, 1]])

計算した説明変数を、既存の説明変数にカラムとして追加するには次のようにすれば良い。

>>> a = X
>>> np.concatenate((a, b.T), axis=1)
array([[  6.32000000e-03,   1.80000000e+01,   2.31000000e+00, ...,
          3.96900000e+02,   4.98000000e+00,   1.13760000e-01],
       [  2.73100000e-02,   0.00000000e+00,   7.07000000e+00, ...,
          3.96900000e+02,   9.14000000e+00,   0.00000000e+00],
       [  2.72900000e-02,   0.00000000e+00,   7.07000000e+00, ...,
          3.92830000e+02,   4.03000000e+00,   0.00000000e+00],
       ...,
       [  6.07600000e-02,   0.00000000e+00,   1.19300000e+01, ...,
          3.96900000e+02,   5.64000000e+00,   0.00000000e+00],
       [  1.09590000e-01,   0.00000000e+00,   1.19300000e+01, ...,
          3.93450000e+02,   6.48000000e+00,   0.00000000e+00],
       [  4.74100000e-02,   0.00000000e+00,   1.19300000e+01, ...,
          3.96900000e+02,   7.88000000e+00,   0.00000000e+00]])

これで説明変数が一つ増える。

>>> np.concatenate((a, b.T), axis=1).shape
(506, 14)

上記にもとづいて相互作用モデルで使う全ての変数を揃える。

>>> X = dataset.data
>>> for i, j in itertools.combinations(range(X.shape[1]), 2):
...     interaction_column = np.array([X[:, i] * X[:, j]])
...     X = np.concatenate((X, interaction_column.T), axis=1)
...
>>> X.shape
(506, 91)

そして、バイアス項を追加する。

>>> X = sm.add_constant(X)
>>> X.shape
(506, 92)

上記の説明変数を使って回帰する。

>>> model = sm.OLS(Y, X)
>>> result = model.fit()
>>> result.summary()
<class 'statsmodels.iolib.summary.Summary'>
"""
                            OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.921
Model:                            OLS   Adj. R-squared:                  0.904
Method:                 Least Squares   F-statistic:                     53.21
Date:                Thu, 22 Dec 2016   Prob (F-statistic):          5.85e-181
Time:                        19:16:55   Log-Likelihood:                -1197.3
No. Observations:                 506   AIC:                             2579.
Df Residuals:                     414   BIC:                             2967.
Df Model:                          91
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const       -159.4473     67.980     -2.346      0.019      -293.076   -25.819
x1           -16.9621      6.549     -2.590      0.010       -29.835    -4.089
x2            -0.0769      0.458     -0.168      0.867        -0.977     0.823
x3            -2.8163      1.695     -1.662      0.097        -6.148     0.516
x4            44.6397     19.519      2.287      0.023         6.271    83.008
x5            22.1561     75.150      0.295      0.768      -125.567   169.879
x6            25.4191      5.700      4.459      0.000        14.214    36.624
x7             1.2569      0.273      4.612      0.000         0.721     1.793
x8            -1.6398      4.604     -0.356      0.722       -10.690     7.410
x9             1.8056      2.462      0.733      0.464        -3.034     6.646
x10            0.0374      0.144      0.260      0.795        -0.246     0.320
x11            2.7520      2.849      0.966      0.335        -2.849     8.353
x12            0.1008      0.074      1.354      0.177        -0.046     0.247
x13            1.6587      0.853      1.944      0.053        -0.018     3.336
x14            0.4124      0.180      2.287      0.023         0.058     0.767
x15           -0.0467      0.448     -0.104      0.917        -0.927     0.834
x16            2.4377      0.571      4.271      0.000         1.316     3.560
x17           -1.2206      0.920     -1.327      0.185        -3.029     0.588
x18            0.2117      0.049      4.345      0.000         0.116     0.307
x19           -0.0031      0.004     -0.828      0.408        -0.011     0.004
x20           -0.1956      0.106     -1.848      0.065        -0.404     0.012
x21           -0.6599      0.581     -1.135      0.257        -1.803     0.483
x22            0.0349      0.043      0.815      0.416        -0.049     0.119
x23            0.4895      0.333      1.471      0.142        -0.164     1.143
x24           -0.0004      0.000     -2.513      0.012        -0.001 -9.71e-05
x25            0.0295      0.007      4.512      0.000         0.017     0.042
...(省略)...
x82         3.151e-05      0.001      0.022      0.983        -0.003     0.003
x83           -0.0439      0.084     -0.523      0.601        -0.209     0.121
x84           -0.0004      0.003     -0.165      0.869        -0.005     0.005
x85           -0.0251      0.018     -1.384      0.167        -0.061     0.011
x86            0.0078      0.003      3.135      0.002         0.003     0.013
x87        -2.303e-06      0.000     -0.012      0.991        -0.000     0.000
x88           -0.0014      0.001     -1.167      0.244        -0.004     0.001
x89            0.0012      0.003      0.349      0.727        -0.005     0.008
x90            0.0038      0.030      0.128      0.898        -0.055     0.062
x91           -0.0006      0.000     -1.510      0.132        -0.001     0.000
==============================================================================
Omnibus:                      121.344   Durbin-Watson:                   1.524
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              628.079
Skew:                           0.942   Prob(JB):                    4.11e-137
Kurtosis:                       8.123   Cond. No.                     1.18e+08
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.18e+08. This might indicate that there are
strong multicollinearity or other numerical problems.
"""

上記を見ると R で相互作用モデルを使って計算したときの内容と一致していることが分かる。

まとめ

今回は R と Python (statsmodels) で線形モデルを使った重回帰分析をしてみた。

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

Python: 多変数の関数から勾配法で最小値を探す

以前、このブログで一変数の関数から勾配法で最小値を探す記事を書いた。

blog.amedama.jp

このときは題材として という一変数の関数を扱った。 今回は、これを多変数の関数に拡張してみることにする。 ちなみに、この多変数というのは機械学習における多次元と同じ意味を持っている。 つまり、これができると多次元のデータセットで損失関数の出力から最小値を探すことができることを意味している。

使った環境は次の通り。

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.11.6
BuildVersion:   15G1108
$ python --version
Python 3.5.2

下準備

まずは、前回と同じように可視化に使う matplotlib と、数値計算用のライブラリである numpy をインストールしておこう。

$ pip install matplotlib numpy

多変数の関数を解析的に偏微分する

今回題材として扱う関数は にする。

勾配法では最初に関数を微分する必要があった。 これは、関数が多変数になったときも変わらない。 ただし、多変数になると名前だけは偏微分と呼ばれるようになる。 また、変数が複数あるので、それぞれの関数について微分することになる。

試しに、題材とする関数を について解析的に偏微分してみよう。 このとき、微分する対象ではない変数 は固定値と捉える。 そのため、結果は となる。

次のサンプルコードでは のときの を計算している。 とはいえ、上記の内容から結果は とすぐに分かる。

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


def f(x0, x1):
    """偏微分したい関数: f(x) = x_0^2 + x_1^2"""
    return x0 ** 2 + x1 ** 2


def df_x0(x0, x1):
    """関数 f の解析的な x0 についての偏微分: f'(x) = 2x_0"""
    return 2 * x0  # x_1 は固定値と捉えて微分することで消滅しており登場しない


def main():
    grad = df_x0(3.0, 4.0)  # ぶっちゃけ x_1 は影響を与えないので何を使っても良い
    print('解析的な x0 についての偏微分: {0}'.format(grad))


if __name__ == '__main__':
    main()

実行結果は次の通り。 まあ、当たり前だね。

解析的な x0 についての偏微分: 6.0

多変数の関数から数値解法で偏微分の近似値を得る

先ほどの例では人が数式を見て解析的に偏微分したけど、これだと機械学習アルゴリズムには適用が難しい。 近似値で良いから、先ほどの偏微分をコンピュータにやらせたい。 そこで、今度は前回のブログでも扱った中央差分を元にした数値微分を適用してみることにしよう。

次のサンプルコードを見てほしい。 ここでは、先ほどと同じように のときの を計算している。 ポイントとしては、偏微分というのは微分に使わない変数を固定値として扱うところだ。 なので を 4.0 に固定した関数 f_x0 を定義している。 そして、その関数を前回のブログ記事にも登場した数値微分の関数 numerical_diff で微分している。

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


def f(x0, x1):
    """偏微分したい関数: f(x) = x0^2 + x1^2"""
    return x0 ** 2 + x1 ** 2


def f_x0(x0):
    """関数 f を x0 について偏微分するために x1 を固定値にした関数"""
    return f(x0, 4.0)


def numerical_diff(f, x):
    """中央差分を元に数値微分する関数"""
    h = 1e-4
    return (f(x + h) - f(x - h)) / (2 * h)


def main():
    grad = numerical_diff(f_x0, 3.0)
    print('数値微分を使った x0 についての偏微分: {0}'.format(grad))


if __name__ == '__main__':
    main()

上記を実行してみよう。

数値微分を使った x0 についての偏微分: 6.00000000000378

おお！解析的な偏微分で得られた値とほとんど同じ値 (近似値) が得られている！

NumPy を使って汎用的にしてみる

先ほどのサンプルコードでは偏微分をするために、微分に関係しない変数を固定値にした関数を新たに用意していた。 しかし、これを毎回やるのは大変なので、今度は数値微分の関数 numerical_diffを偏微分用に拡張してみよう。

次のサンプルコードでは、偏微分する関数 f の引数を変更している。 具体的には、これまで別々の引数としていた と を NumPy の配列で渡すことにしている。 また、偏微分に拡張した関数 numerical_diff には引数 i を追加している。 これは、引数 x のうち何番目の要素で偏微分するかを表すインデックスになっている。 例えば なら について偏微分することを意味している。

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

import numpy as np


def f(x):
    """偏微分したい関数: f(x) = x0^2 + x1^2"""
    return x[0] ** 2 + x[1] ** 2


def numerical_diff(f, x, i):
    """中央差分を元に数値微分する関数 (偏微分)

    :param function f: 偏微分する関数
    :param numpy.ndarray x: 偏微分する引数
    :param int i: 偏微分する変数のインデックス
    """
    # 丸め誤差で無視されない程度に小さな値を用意する
    h = 1e-4
    # 偏微分する変数のインデックスにだけ上記の値を入れる
    h_vec = np.zeros_like(x)
    h_vec[i] = h
    # 数値微分を使って偏微分する
    return (f(x + h_vec) - f(x - h_vec)) / (2 * h)


def main():
    args = [3.0, 4.0]  # 偏微分に使う引数
    grad = numerical_diff(f, args, 0)  # 0 番目の要素 (x0) について偏微分する
    print('数値微分を使った x0 についての偏微分: {0}'.format(grad))


if __name__ == '__main__':
    main()

ポイントとしては偏微分する要素にだけ丸め誤差で無視されない程度に小さな値を加えて数値微分している点だ。

上記を実行してみよう。

数値微分を使った x0 についての偏微分: 6.00000000000378

先ほどと同じように上手く偏微分できていることが分かる。

偏微分を元に勾配を計算する

多変数の関数において、各変数について偏微分して得られた値をベクトルにまとめたものを勾配と呼ぶ。 次は、この勾配を計算してみよう。 とはいえ、ある変数について偏微分するやり方は先ほどの例で分かったので、あとはそれを繰り返すだけ。

次のサンプルコードでは numerical_gradient という関数で購買を計算している。 ポイントは特に無くて、各変数について淡々と偏微分するだけ。

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

import numpy as np


def f(x):
    """偏微分したい関数: f(x) = x0^2 + x1^2"""
    return x[0] ** 2 + x[1] ** 2


def numerical_diff(f, x, i):
    """中央差分を元に数値微分する関数 (偏微分)

    :param function f: 偏微分する関数
    :param numpy.ndarray x: 偏微分する引数
    :param int i: 偏微分する変数のインデックス
    """
    # 丸め誤差で無視されない程度に小さな値を用意する
    h = 1e-4
    # 偏微分する変数のインデックスにだけ上記の値を入れる
    h_vec = np.zeros_like(x)
    h_vec[i] = h
    # 数値微分を使って偏微分する
    return (f(x + h_vec) - f(x - h_vec)) / (2 * h)


def numerical_gradient(f, x):
    """勾配を計算する関数

    勾配というのは、全ての変数について偏微分した結果をベクトルとしてまとめたものを言う。
    """
    # 勾配を入れるベクトルをゼロで初期化する
    grad = np.zeros_like(x)

    for i, _ in enumerate(x):
        # i 番目の変数で偏微分する
        grad[i] = numerical_diff(f, x, i)

    # 計算した勾配を返す
    return grad


def main():
    grad = numerical_gradient(f, np.array([3.0, 4.0]))
    print('勾配: {0}'.format(grad))


if __name__ == '__main__':
    main()

実行すると のときの勾配が得られる。

勾配: [ 6.  8.]

題材となる関数を可視化してみる

さて、先ほどのサンプルコードで勾配を計算することができた。 とはいえ、何のために勾配を計算したのかよく分からないかもしれないので、ここで補足しておく。 先ほど計算した勾配には、それぞれの変数について偏微分したときの傾きが入っている。 これはつまり、各次元ごとのどちらに進めば値が小さくなるかという情報だ。

上記の説明だけでは分かりづらいかもしれないので、より原点に立ち返った説明もしてみよう。 例えば という三つの軸を持った三次元の空間をイメージしてほしい。 通常であれば変数の名前は になるけど、ここではあえて今回使った変数の名前にしている。 まず と が関数の入力で が出力だとしよう。 縦と横が と に対応していて高さが ということになる。 今回の問題は高さ が最も低くなる位置を と の平面の中から探したいということだ。

試しに、今回題材とする関数 を上記のように三次元で捉えて可視化してみることにしよう。

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

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import pyplot as plt
import numpy as np


def main():
    x0 = np.arange(-2, 2, 0.1)
    x1 = np.arange(-2, 2, 0.1)
    X0, X1 = np.meshgrid(x0, x1)
    Y = X0 ** 2 + X1 ** 2

    fig = plt.figure()
    ax = Axes3D(fig)
    ax.plot_wireframe(X0, X1, Y)

    plt.show()


if __name__ == '__main__':
    main()

上記を実行すると、次のようなグラフが得られる。 上記より、最小値は ということが分かる。

偏微分を元にした勾配の計算というのは、上記のグラフでいえば各点における傾きを という、それぞれの軸ごとに計算するという意味になる。

勾配法への組み込み

ここで、前回も扱った勾配法の数式を多変数に拡張して示す。

上記で現在の位置を としたとき、より最小値に近づいたものが になる。 は学習率と呼ばれるもので、どれくらいの勢いで最小値に近づいていくかを表している。 そして は を について偏微分した関数になる。

上記は、扱う変数が増えて多変数 (多次元) になったとしてもアルゴリズムの基本は全く変わらない。 単に、各次元ごとにそれぞれで勾配法を実行するだけだ。 次のサンプルコードでは初期位置を として、学習率 0.01 ステップ数 50 回で勾配法を実行している。

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

from matplotlib import pyplot as plt
import numpy as np


def f(x):
    """偏微分したい関数: f(x) = x0^2 + x1^2"""
    return x[0] ** 2 + x[1] ** 2


def numerical_diff(f, x, i):
    """中央差分を元に数値微分する関数 (偏微分)

    :param function f: 偏微分する関数
    :param numpy.ndarray x: 偏微分する引数
    :param int i: 偏微分する変数のインデックス
    """
    # 丸め誤差で無視されない程度に小さな値を用意する
    h = 1e-4
    # 偏微分する変数のインデックスにだけ上記の値を入れる
    h_vec = np.zeros_like(x)
    h_vec[i] = h
    # 数値微分を使って偏微分する
    return (f(x + h_vec) - f(x - h_vec)) / (2 * h)


def numerical_gradient(f, x):
    """勾配を計算する関数

    勾配というのは、全ての変数について偏微分した結果をベクトルとしてまとめたものを言う。
    """
    # 勾配を入れるベクトルをゼロで初期化する
    grad = np.zeros_like(x)

    for i, _ in enumerate(x):
        # i 番目の変数で偏微分する
        grad[i] = numerical_diff(f, x, i)

    # 計算した勾配を返す
    return grad


def gradient_descent(f, initial_position, learning_rate=0.1, steps=50):
    """勾配法で最小値を求める関数

    :param function f: 最小値を見つけたい関数
    :param numpy.ndarray initial_position: 関数の初期位置
    :param float learning_rate: 学習率
    :param int steps: 学習回数
    """
    # 現在地を示すカーソル
    x = initial_position

    # 学習を繰り返す
    for _ in range(steps):
        # 現在地の勾配 (どちらにどれだけ進むべきか) を得る
        grad = numerical_gradient(f, x)
        # 勾配を元にして現在地を移動する
        x -= learning_rate * grad

    # 最終的な位置を返す
    return x


def main():
    # 勾配法を使って関数 f() の最小値を探す (初期位置は 1, 2)
    min_pos = gradient_descent(f, [1, 2])
    print('勾配法が見つけた最小値: {0}'.format(min_pos))


if __name__ == '__main__':
    main()

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

勾配法が見つけた最小値: [  1.78405962e-05   3.56811923e-05]

最小値である に、どちらの次元についても近づけていることが分かる。

最小値に収束していく様子を可視化してみる

先ほどの結果から、最小値に収束させることができることは分かった。 次は、どのように収束していっているのかをグラフで可視化してみよう。

次のサンプルコードでは各ステップごとの と の位置を折れ線グラフで示している。

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

from matplotlib import pyplot as plt
import numpy as np


def f(x):
    """偏微分したい関数: f(x) = x0^2 + x1^2"""
    return x[0] ** 2 + x[1] ** 2


def numerical_diff(f, x, i):
    """中央差分を元に数値微分する関数 (偏微分)

    :param function f: 偏微分する関数
    :param numpy.ndarray x: 偏微分する引数
    :param int i: 偏微分する変数のインデックス
    """
    # 丸め誤差で無視されない程度に小さな値を用意する
    h = 1e-4
    # 偏微分する変数のインデックスにだけ上記の値を入れる
    h_vec = np.zeros_like(x)
    h_vec[i] = h
    # 数値微分を使って偏微分する
    return (f(x + h_vec) - f(x - h_vec)) / (2 * h)


def numerical_gradient(f, x):
    """勾配を計算する関数

    勾配というのは、全ての変数について偏微分した結果をベクトルとしてまとめたものを言う。
    """
    # 勾配を入れるベクトルをゼロで初期化する
    grad = np.zeros_like(x)

    for i, _ in enumerate(x):
        # i 番目の変数で偏微分する
        grad[i] = numerical_diff(f, x, i)

    # 計算した勾配を返す
    return grad


def gradient_descent(f, initial_position, learning_rate=0.1):
    """勾配法で最小値を求める関数

    :param function f: 最小値を見つけたい関数
    :param numpy.ndarray initial_position: 関数の初期位置
    :param float learning_rate: 学習率
    """
    # 現在地を示すカーソル
    x = initial_position

    # 学習を繰り返す
    while True:
        # 現在地の勾配 (どちらにどれだけ進むべきか) を得る
        grad = numerical_gradient(f, x)
        # 勾配を元にして現在地を移動する
        x -= learning_rate * grad
        # 現在の位置を返す
        yield np.copy(x)


def main():
    g = gradient_descent(f, [1, 2])

    x = range(50)
    y = np.array([next(g) for _ in x])

    plt.plot(x, y[:, 0], label='f_x0(x)')
    plt.plot(x, y[:, 1], label='f_x1(x)')
    plt.legend()
    plt.show()


if __name__ == '__main__':
    main()

上記の実行結果は次の通り。 それぞれの軸で最小値に収束していく様子が見て取れる。

まとめ

• 勾配法は多変数の関数でも使える
• ただし傾きを得るときの計算が偏微分になる
• 偏微分は解析的に計算する必要はなく数値解法での近似値が使える
• 各変数について偏微分した結果のベクトルを勾配という
• 勾配を元に各次元ごとに最小値を探すことができる

参考文献

ゼロから作るDeep Learning ―Pythonで学ぶディープラーニングの理論と実装

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

Python: データセットの標準化について

今回は機械学習とか統計で扱うデータセットの標準化について。

まずは、標準化されていない生のデータセットについて考えてみよう。 それらの多くは、次元によって数値の単位がバラバラだったり、あるいは大きさが極端に異なったりする。 これをそのまま扱ってしまうと、各次元を見比べたときにそれぞれの関係が分かりにくい。 また、機械学習においては特定の次元の影響が強く (または反対に弱く) 出てしまったりすることもあるらしい。 そこで、それぞれの次元のスケールを同じに揃えてやりたい。 これを標準化というようだ。

今回は「Zスコア」という標準化のやり方を扱う。 これは、一言で言ってしまえばデータセットの各要素から平均を引いて、標準偏差で割ったもの。 これをすると、データセットは平均が 0 で標準偏差・分散が 1 になる。

使った環境は次の通り。

$ python --version
Python 3.5.1

NumPy を使った標準化

まずは一番単純な NumPy の配列を直接使ったやり方から。

ひとまず pip を使って NumPy をインストールしておこう。

$ pip install numpy

Python の REPL を起動しよう。

$ python

NumPy をインポートしたら、次にZスコア化する対象となる配列を変数 l に代入する。 ようするに、これがデータセットを模したものということ。

>>> import numpy as np
>>> l = np.array([10, 20, 30, 40, 50])

この時点でデータセットの算術平均は 30 になっている。

>>> l.mean()
30.0

Zスコアによる標準化の第一段階では、まずデータセットの各要素から、データセットの平均値を引く。

>>> l2 = l - l.mean()

データセットの平均値は、これで 0 になった。

>>> l2.mean()
0.0

各要素から平均値を引いたので、中身はこんなことになる。

>>> l2
array([-20., -10.,   0.,  10.,  20.])

次にZスコア化の第二段階では標準偏差を扱う。 今のデータセットの標準偏差は約 14.14 だ。

>>> l2.std()
14.142135623730951

この標準偏差で、同じくデータセットの各要素を割る。

>>> l3 = l2 / l2.std()

これでデータベースの標準偏差は 1 になった。 きれいに 1 になっていないのは浮動小数点の計算誤差だろうね。

>>> l3.std()
0.99999999999999989

中身はこんなことになる。

>>> l3
array([-1.41421356, -0.70710678,  0.        ,  0.70710678,  1.41421356])

そして、これこそが標準化されたデータセットだ。 各要素の単位は、もはや元々のデータセットで扱われていたそれではなくなっている。 代わりに、単に各要素の相対的な位置関係を表したZスコアになっている。

SciPy を使った標準化

さっきは NumPy の配列を自分で標準化してみたけど、実際にはこの作業はライブラリが代わりにやってくれる。 ここでは SciPy を使った例を挙げておこう。

先ほどと同じように、今度は SciPy を pip コマンドでインストールする。

$ pip install scipy

インストールできたら Python の REPL を起動する。

$ python

さっきと同じようにデータセットを模した NumPy 配列を変数 l として用意しておく。

>>> import numpy as np
>>> l = np.array([10, 20, 30, 40, 50])

SciPy では、そのものずばり zscore() という名前の関数が用意されている。

>>> from scipy.stats import zscore

これに NumPy の配列を渡せば一発でZスコアに標準化できる。

>>> zscore(l)
array([-1.41421356, -0.70710678,  0.        ,  0.70710678,  1.41421356])

あっけない。

scikit-learn を使った標準化

同じように scikit-learn にも標準化のユーティリティが用意されている。

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

$ pip install scikit-learn

Python の REPL を起動したら…

$ python

データセットを模した NumPy 配列を用意して…

>>> import numpy as np
>>> l = np.array([10, 20, 30, 40, 50])

scikit-learn では StandardScaler というクラスを使って標準化する。

>>> from sklearn.preprocessing import StandardScaler

インスタンス化したら、データセットを模した配列を渡す。 これでデータセットの情報が StandardScaler のインスタンスにセットされる。

>>> sc = StandardScaler()
>>> sc.fit(l)

そして StandardScaler#transform() メソッドを使ってデータセットを標準化する。 返り値として標準化されたデータセットが返る。

>>> sc.transform(l)
array([-1.41421356, -0.70710678,  0.        ,  0.70710678,  1.41421356])

ばっちり。 ちなみに、上記を実行すると NumPy 配列の型が int から float に変換されたよ！っていう警告が出る。 けど、意図したものなので無視しておっけー。

まとめ

• データセットは標準化してから扱ったほうが良いらしい
• 標準化は「Zスコア」というものを使うのがメジャーっぽい
• 標準化のやり方には NumPy / SciPy / scikit-learn それぞれにやり方がある

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