主成分分析 - 静かなる名辞

SVMのsupport vectorを可視化してみた

2019-07-01T21:52:44+09:00

はじめに

　SVMはヒンジ関数を使ってマージン最大化を行い、境界付近のデータに基づいて分離超平面を決定する……ということはよく言われています。でも、実際のデータで確認している図はあまり見たことがありません。

　sklearnのSVMのドキュメントを読んでいたら、属性からサポートベクトル関連の情報が取れることがわかったので、いつものようにirisで見てみます。

見方

　ここに書いてあります。
sklearn.svm.SVC — scikit-learn 0.21.3 documentation

support_ : array-like, shape = [n_SV]
　　Indices of support vectors.
support_vectors_ : array-like, shape = [n_SV, n_features]
　　Support vectors.
n_support_ : array-like, dtype=int32, shape = [n_class]
　　Number of support vectors for each class.

　これを使うと便利です。support_とsupport_vectors_はどちらを使っても良いのですが、散布図でサポートベクトルとそれ以外の点を分けてプロットしたいという都合上、support_の方を使います。

　なお、小ネタとして、indiciesの配列がある場合、以下のようにすることでboolean maskに変換できます。今回は論理反転を使いたいので、これが使えると便利です。

>>> import numpy as np
>>> a = np.arange(10)
>>> idx = np.where(a > 4)
>>> idx
(array([5, 6, 7, 8, 9]),)
>>> mask = np.zeros(10, dtype="bool")
>>> mask[idx] = True
>>> mask
array([False, False, False, False, False,  True,  True,  True,  True,
        True])
>>> ~mask
array([ True,  True,  True,  True,  True, False, False, False, False,
       False])

　参考：python - How to invert numpy.where (np.where) function - Stack Overflow

コード

　素直にやります。以前やったRGB予測確率プロットも同時に出します。

sklearnとmatplotlibでiris（3クラス）の予測確率を可視化した話 - 静かなる名辞

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.svm import SVC

def main():
    iris = load_iris()

    # とりあえず二次元に
    X = PCA(n_components=2).fit_transform(iris.data)

    # 学習
    clf = SVC(gamma="scale", probability=True)
    clf.fit(X, iris.target)
    
    # support_をboolean maskにしておく
    support = np.zeros(X.shape[0], dtype="bool")
    support[clf.support_] = True
 
    # --散布図--
    ax = plt.subplot()
    cm = ListedColormap(["b", "g", "r"])

    # サポートベクトルとそれ以外を違うmarkerで
    ax.scatter(X[~support,0], X[~support,1], marker="2",
               c=iris.target[~support], cmap=cm, alpha=0.5)
    ax.scatter(X[support,0], X[support,1], marker=".", 
               c=iris.target[support], cmap=cm)

    # 確率のプロット
    XX, YY = np.meshgrid(np.arange(-5, 5, 0.025),
                         np.arange(-2, 2, 0.025))
    Z = clf.predict_proba(np.stack([XX.ravel(), YY.ravel()], axis=1))
    ZZ = np.flip(Z.reshape(XX.shape + (3, )), axis=1)
    ax.imshow(ZZ, alpha=0.1,
              aspect="auto", extent=(-5, 5, -2, 2))

    plt.savefig("result.png")
    
if __name__ == "__main__":
    main()

結果

　濃いめの点がサポートベクトル、薄い三菱マークがそれ以外です。なるほど、こうなっているのかという発見というか納得感が得られます。あと、よく見るとマージン最大化のため、緑色の領域がずいぶん青い側に寄っているみたいですね。

まとめ

　簡単に見れるので、2次元のデータでいろいろ突っ込んでみると面白いと思います。実際には全体の数割くらいのデータしか使わないで学習している様子がわかります。

sklearnとmatplotlibでiris（3クラス）の予測確率を可視化した話

2019-06-22T23:51:33+09:00

はじめに

　よく分類器の性質などを把握するために、2次元で可視化している図があります。

　特に予測確率なんかを平面的に出せるとかっこいいですよね。つまり、こういうのです。

Classifier comparison — scikit-learn 0.21.3 documentation

以前の記事より

君はKNN（k nearest neighbor）の本当のすごさを知らない - 静かなる名辞

　ただ、これが素直にできるのは2クラス分類までで、3クラス分類だと下のような図にしかなりません。

以前の記事より

【python】高次元の分離境界をなんとか2次元で見る - 静かなる名辞

　ということでずっと諦めていたのですが、ふと思いました。

「RGBに各クラスの予測確率あてればできるじゃん」

　簡単にできると思ったら思いの外手間取ったので、備忘録として書いておきます。

まずやる

　とりあえずirisを二次元でプロットします。この辺は定石どおりにやるだけです。

import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA

def main():
    iris = load_iris()
    pca = PCA(n_components=2)
    X = pca.fit_transform(iris.data)
    ax = plt.subplot()
    ax.scatter(X[:,0], X[:,1], c=iris.target, cmap="brg")
    plt.savefig("fig1.png")

if __name__ == "__main__":
    main()

fig1.png

　kNNを学習させて、まずは普通に分離境界を描きます。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.neighbors import KNeighborsClassifier

def main():
    iris = load_iris()
    pca = PCA(n_components=2)
    X = pca.fit_transform(iris.data)
    ax = plt.subplot()
    ax.scatter(X[:,0], X[:,1], c=iris.target, cmap="brg")
    
    clf = KNeighborsClassifier()
    clf.fit(X, iris.target)
    XX, YY = np.meshgrid(np.arange(-5, 5, 0.025),
                         np.arange(-2, 2, 0.025))
    Z = clf.predict(np.stack([XX.ravel(), YY.ravel()], axis=1))
    ZZ = Z.reshape(XX.shape)
    ax.pcolormesh(XX, YY, ZZ, alpha=0.05, cmap="brg", shading="gouraud")

    plt.savefig("fig2.png")

if __name__ == "__main__":
    main()

　参考：matplotlibのpcolormeshでalphaを小さくすると網目が出てくる対策 - 静かなる名辞

fig2.png

　さ、次はpredict_probaを呼ぶ訳ですが……pcolormeshとかこの辺の関数にはRGBのデータは渡せません。

matplotlib.pyplot.pcolormesh — Matplotlib 3.1.0 documentation

　しばし思案したあと、imshowならできると思いました。

　なにも考えずに書くと下のようなコードになります。

（これは動かないので注意してください）

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.neighbors import KNeighborsClassifier

def main():
    iris = load_iris()
    pca = PCA(n_components=2)
    X = pca.fit_transform(iris.data)
    ax = plt.subplot()
    ax.scatter(X[:,0], X[:,1], c=iris.target, cmap="brg")
    
    clf = KNeighborsClassifier()
    clf.fit(X, iris.target)
    XX, YY = np.meshgrid(np.arange(-5, 5, 0.025),
                         np.arange(-2, 2, 0.025))
    Z = clf.predict_proba(np.stack([XX.ravel(), YY.ravel()], axis=1))
    ZZ = Z.reshape(XX.shape + (3, ))
    ax.imshow(ZZ, alpha=0.2)

    plt.savefig("fig3.png")

if __name__ == "__main__":
    main()

fig3.png

　なにこれ？

　ああ、縮尺を合わせないといけないんですね。aspect, extentという引数でできそうです。

matplotlib.pyplot.imshow — Matplotlib 3.1.0 documentation

（これもちゃんと動かないので注意してください）

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.neighbors import KNeighborsClassifier

def main():
    iris = load_iris()
    pca = PCA(n_components=2)
    X = pca.fit_transform(iris.data)
    ax = plt.subplot()
    ax.set_xlim((-5, 5))
    ax.set_ylim((-2, 2))
    ax.scatter(X[:,0], X[:,1], c=iris.target, cmap="brg")
    
    clf = KNeighborsClassifier()
    clf.fit(X, iris.target)
    XX, YY = np.meshgrid(np.arange(-5, 5, 0.025),
                         np.arange(-2, 2, 0.025))
    Z = clf.predict_proba(np.stack([XX.ravel(), YY.ravel()], axis=1))
    ZZ = Z.reshape(XX.shape + (3, ))
    ax.imshow(ZZ, alpha=0.2,
              aspect="auto", extent=(-5, 5, -2, 2))

    plt.savefig("fig4.png")

if __name__ == "__main__":
    main()

　まず、先にax.set_xlimとax.set_ylimで図の範囲を指定し、そこにextentをあわせるようにしています。aspectはドキュメントを見た感じだとautoが無難そうに思います。

fig4.png

　どう見ても上下反転しているので、ZZを上下反転します。ついでに、マーカーの色を揃えることにします。

これが動くコードです

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.neighbors import KNeighborsClassifier

def main():
    iris = load_iris()
    pca = PCA(n_components=2)
    X = pca.fit_transform(iris.data)
    ax = plt.subplot()
    ax.set_xlim((-5, 5))
    ax.set_ylim((-2, 2))
    cm = ListedColormap(["b", "g", "r"])
    ax.scatter(X[:,0], X[:,1], c=iris.target, cmap=cm)
    
    clf = KNeighborsClassifier()
    clf.fit(X, iris.target)
    XX, YY = np.meshgrid(np.arange(-5, 5, 0.025),
                         np.arange(-2, 2, 0.025))
    Z = clf.predict_proba(np.stack([XX.ravel(), YY.ravel()], axis=1))
    ZZ = np.flip(Z.reshape(XX.shape + (3, )), axis=1)
    ax.imshow(ZZ, alpha=0.2,
              aspect="auto", extent=(-5, 5, -2, 2))

    plt.savefig("fig5.png")

if __name__ == "__main__":
    main()

fig5.png

　だいたい不満のない結果になりました。ここまで長かった。

他の分類器も試す

　せっかくなのでいろいろやってみます。SVM, ロジスティック回帰, ランダムフォレストを追加してやってみましょう。
　

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.neighbors import KNeighborsClassifier
from sklearn.svm import SVC
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier

def main():
    iris = load_iris()
    pca = PCA(n_components=2)
    X = pca.fit_transform(iris.data)

    fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(9, 9))
    knn = KNeighborsClassifier()
    svm = SVC(probability=True)
    lr = LogisticRegression()
    rfc = RandomForestClassifier(n_estimators=100)

    cm = ListedColormap(["b", "g", "r"])    
    XX, YY = np.meshgrid(np.arange(-5, 5, 0.025),
                         np.arange(-2, 2, 0.025))

    for ax, clf in zip(axes.ravel(), [knn, svm, lr, rfc]):
        ax.set_xlim((-5, 5))
        ax.set_ylim((-2, 2))
        ax.scatter(X[:,0], X[:,1], c=iris.target, cmap=cm)
    
        clf.fit(X, iris.target)
        Z = clf.predict_proba(np.stack([XX.ravel(), YY.ravel()], axis=1))
        ZZ = np.flip(Z.reshape(XX.shape + (3, )), axis=1)
        ax.imshow(ZZ, alpha=0.2,
                  aspect="auto", extent=(-5, 5, -2, 2))
        ax.set_title(clf.__class__.__name__)

    plt.tight_layout()
    plt.savefig("fig6.png")

if __name__ == "__main__":
    main()

fig6.png

　こんなもんか、という感じ。

まとめ

　やればできることはわかりました。もう少しかっこいい図にするには、さらなる工夫が要るのかもしれません。

【python】高次元の分離境界をなんとか2次元で見る

2019-05-24T04:01:31+09:00

はじめに

　分類器の特性を把握するために2次元データで分離境界を見るということが行われがちですが、高次元空間における分離器の特性を正確に表している訳ではありません。

　ということがずっと気になっていたので、なんとか高次元空間で分類させて2次元で見る方法を考えます。

方法

　PCAで2次元に落とせれば、線形変換で逆変換もできるので、それでやります。当然ながら情報は落ちますし、2次元でもなんとか見える程度のデータしか扱えませんが、妥協します。

　sklearnならinverse_transformという便利なメソッドがあるので、簡単です。

　というあたりまで考えた上で、こんなコードを書きました。

show_hyperplane.py

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

def show_hyperplane(dataset, clf, filename):
    pca = PCA(n_components=2)
    X = pca.fit_transform(dataset.data)
    plt.scatter(X[:,0], X[:,1], c=dataset.target)

    x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.01),
                         np.arange(y_min, y_max, 0.01))

    clf.fit(dataset.data, dataset.target)
    Z = clf.predict(
        pca.inverse_transform(np.c_[xx.ravel(), yy.ravel()]))
    plt.pcolormesh(xx, yy, Z.reshape(xx.shape),
                   alpha=0.03, shading="gouraud")
    plt.savefig(filename)

　汎用的に作ったので、これでいろいろなものを見てみようという算段です。

実験

　まずirisとSVM。

from show_hyperplane import show_hyperplane
from sklearn.datasets import load_iris
from sklearn.svm import SVC

iris = load_iris()
svm = SVC(C=50, gamma="scale")    
show_hyperplane(iris, svm, "iris_svm.png")

iris_svm.png

　特に興味深い知見は得られませんでした。

　次、irisとランダムフォレスト。

from show_hyperplane import show_hyperplane
from sklearn.datasets import load_iris
from sklearn.ensemble import RandomForestClassifier

iris = load_iris()
rfc = RandomForestClassifier(n_estimators=500, n_jobs=-1)    
show_hyperplane(iris, rfc, "iris_rf.png")

iris_rf.png

　ランダムフォレストで斜めの分離超平面の図を出したサイトはここくらいしかないのでは？　だからどうしたって話ですが。

　簡単なのでAdaBoostも試します。

from show_hyperplane import show_hyperplane
from sklearn.datasets import load_iris
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import AdaBoostClassifier

iris = load_iris()
ada = AdaBoostClassifier(
    base_estimator=DecisionTreeClassifier(max_depth=4),
    n_estimators=200)    
show_hyperplane(iris, ada, "iris_ada.png")

iris_ada.png

　面白いんですが、性能はいまいち悪そう。

　ちなみに、base_estimatorのパラメータでコロコロ結果が変わります。パラメータ設定については、以下の2記事を参照してください。

【python】sklearnのAdaBoostをデフォルトパラメータで使ってはいけない - 静かなる名辞
 AdaBoostとRandomForestの比較 - 静かなる名辞

　ただの決定木もやっておきましょう。

from show_hyperplane import show_hyperplane
from sklearn.datasets import load_iris
from sklearn.tree import DecisionTreeClassifier

iris = load_iris()
dtc = DecisionTreeClassifier()
show_hyperplane(iris, dtc, "iris_tree.png")

iris_tree.png

　つまらない。

　さて、irisは飽きてきたのでdigitsで同じことをやります。こちらは何しろ元が64次元で、2次元に落とすとかなり重なり合うので、カオスな結果になってくれそうです。

　が、その前にshow_hyperplane.pyをいじります。元のままだといろいろうまくいかなかったからです。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

def show_hyperplane(dataset, clf, filename):
    pca = PCA(n_components=2)
    X = pca.fit_transform(dataset.data)
    plt.scatter(X[:,0], X[:,1], s=5, c=dataset.target)

    x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.3),
                         np.arange(y_min, y_max, 0.3))

    clf.fit(dataset.data, dataset.target)
    Z = clf.predict(
        pca.inverse_transform(np.c_[xx.ravel(), yy.ravel()]))
    plt.pcolormesh(xx, yy, Z.reshape(xx.shape),
                   alpha=0.05, shading="gouraud")
    plt.savefig(filename)

　よし、やろう。

　まずSVM。今回からついでに学習データに対するスコアを見ます。コメントで記します。

from show_hyperplane import show_hyperplane
from sklearn.datasets import load_digits
from sklearn.svm import SVC

digits = load_digits()
svm = SVC(C=0.1, gamma="scale")    
score = svm.fit(
    digits.data, digits.target).score(
        digits.data, digits.target)
print(score) # => 0.9744017807456873
show_hyperplane(digits, svm, "digits_svm.png")

digits_svm.png

　あたりまえですが、64→2次元で情報落ちしているので、こんなふうにしか見えません。それでも、後々出てくるやつに比べればまともな方です。

　次。ランダムフォレスト。

from show_hyperplane import show_hyperplane
from sklearn.datasets import load_digits
from sklearn.ensemble import RandomForestClassifier

digits = load_digits()
rfc = RandomForestClassifier(n_estimators=500, n_jobs=-1)    
score = rfc.fit(
    digits.data, digits.target).score(
        digits.data, digits.target)
print(score) # => 1.0
show_hyperplane(digits, rfc, "digits_rfc.png")

digits_rfc.png

　これ面白いですね。ところどころ凹凸がありますが、それでもぱっと見SVMと同じくらい滑らかな分離超平面に見えます。高次元データほど強いというのもわかる気がします。

　アダブースト。

from show_hyperplane import show_hyperplane
from sklearn.datasets import load_digits
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import AdaBoostClassifier

digits = load_digits()
ada = AdaBoostClassifier(
    base_estimator=DecisionTreeClassifier(max_depth=3),
    n_estimators=200)    
score = ada.fit(
    digits.data, digits.target).score(
        digits.data, digits.target)
print(score) # 0.9660545353366722
show_hyperplane(digits, ada, "digits_ada.png")

digits_ada.png

　大丈夫なんかこれ。決定木のアダブーストはランダムフォレストと比べて個人的にいまいち信頼していないのですが、こういうの見るとその思いが強まります。

　決定木もやりましょう。irisではつまらなかったけど、こちらではどうなるでしょうか。
　

from show_hyperplane import show_hyperplane
from sklearn.datasets import load_digits
from sklearn.tree import DecisionTreeClassifier

digits = load_digits()
dtc = DecisionTreeClassifier()
score = dtc.fit(
    digits.data, digits.target).score(
        digits.data, digits.target)
print(score) # 1.0
show_hyperplane(digits, dtc, "digits_tree.png")

digits_tree.png

　あははー、なにこれカオス。デフォルトのまま高次元で使うな、ということですね。

まとめ

　元が64次元くらいでもだいぶ情報落ちするので、本当の高次元データでも使えるかというと微妙なのですが、それでもなんとなく傾向はつかめますし、面白かったです。

　SVMとランダムフォレストはどっちも優秀ですね。

【python】PCAと非負値行列因子分解のバイプロットを見比べる

2019-05-14T20:48:35+09:00

はじめに

　非負値行列因子分解は負の値が出現しないような行列に対して行える分解で、主成分分析とか因子分析に似ています。

　参考：
非負値行列因子分解（NMF）をふわっと理解する - Qiita

　上の記事によると、いいところとしては、

非負なので現実のデータに向く
非負なので解釈が楽
さらにスパースになる

　というあたりらしい。

　なので、PCAと比べます。sklearnを使います。

比較実験

　irisでやります。なんとかの一つ覚えです。

import numpy as np
import matplotlib.pyplot as plt

from sklearn.decomposition import PCA, NMF
from sklearn.datasets import load_iris

def main():
    iris = load_iris()

    pca = PCA(n_components=2)
    nmf = NMF(n_components=2)

    fig, axes = plt.subplots(nrows=1, ncols=2)

    for i, mname, method in zip([0,1], ["PCA", "NMF"], [pca, nmf]):
        X_2d = method.fit_transform(iris.data)

        # title
        axes[i].set_title("{} {}".format("iris", mname))

        # scatter
        axes[i].scatter(X_2d[:,0], X_2d[:,1], c=iris.target)

        # arrows
        pc0 = method.components_[0]
        pc1 = method.components_[1]

        pc0 = pc0 * (np.abs(X_2d[:,0]).max() / np.abs(pc0).max()) * 0.8 
        pc1 = pc1 * (np.abs(X_2d[:,1]).max() / np.abs(pc1).max()) * 0.8

        for j in range(pc0.shape[0]):
            axes[i].arrow(
                0, 0, pc0[j], pc1[j], color='r')
            axes[i].text(
                pc0[j]*1.1, pc1[j]*1.1, iris.feature_names[j], color='r')

    plt.show()

if __name__ == "__main__":
    main()

結果

　こうして見るとそんなに違いませんが、原点の右上だけが図に含まれるのが見ての通り特徴です。相違点としては

スケールの違いなどがわかるような気がする
PCAでは重なっているpetal width(cm)とpetal length(cm)の違いが出ている

　などがあるでしょうか。また、NMFではpetal width(cm)のy方向成分はゼロのようです。

メリット

　上の図を見る限りでは、別にどっちでも大差はなさそうだし、PCAの方が慣れているので意味的な解釈もしやすい気がします。

　非負であることが要請されるような特殊ケース以外は別にPCAでもこまらないという気もするのですが、実際のところどうなんでしょうね。

まとめ

　とにかく使えます。

本当は怖いSVMと交差検証

2019-01-26T01:21:13+09:00

概要

　SVMと交差検証を組み合わせて使うと、たとえ交差検証で高いスコアが出て汎化性能確保できた！　と思っても想像とかけ離れた分離超平面になっていることがままある。

　なのでこの組み合わせは少し怖いということを説明する。

コード

　irisを分類します。二次元で決定境界を可視化するために、irisを主成分分析を使って2次元に落としておきます。

　GridSearchCVを使って交差検証し、ベストパラメータを探ります。その後、ベストパラメータの分類器で、平面上で散布図と決定境界を可視化してみます。

　ちょっと長いけど、ざっくり読んでみてください。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.svm import SVC
from sklearn.model_selection import GridSearchCV
from sklearn.decomposition import PCA
from sklearn.datasets import load_iris

def main():
    # データの準備
    iris = load_iris()
    X, y = PCA(2).fit_transform(iris.data), iris.target

    # GridSearchCVの準備
    svm = SVC()
    params = {"C":10**np.arange(-3, 3, dtype=float),
              "gamma":10**np.arange(-3, 3, dtype=float)}
    gscv = GridSearchCV(svm, params, cv=8, iid=False, 
                        return_train_score=False,
                        verbose=1, n_jobs=-1)
    gscv.fit(X, y)

    # GridSearchCVの結果を表示する
    result_df = pd.DataFrame(gscv.cv_results_)
    print(result_df[["param_C", "param_gamma",
                     "rank_test_score", "mean_test_score"]
                ][result_df["rank_test_score"]==1])

    # 最良推定器をclfに代入
    clf = gscv.best_estimator_
    
    # 可視化の準備
    xmin, xmax, ymin, ymax = (X[:,0].min()-1, X[:,0].max()+1,
                              X[:,1].min()-1, X[:,1].max()+1)    
    x_ = np.arange(xmin, xmax, 0.01)
    y_ = np.arange(ymin, ymax, 0.01)
    xx, yy = np.meshgrid(x_, y_)

    # 予測
    zz = clf.predict(np.stack([xx.ravel(), yy.ravel()], axis=1)
                 ).reshape(xx.shape)

    # 可視化
    plt.pcolormesh(xx, yy, zz, cmap="winter", alpha=0.1, shading="gouraud")
    plt.scatter(X[:, 0], X[:, 1], c=y, edgecolors='k', cmap="winter")
    plt.savefig("result.png")

if __name__ == "__main__":
    main()

結果

　printされた結果は良さげなものでした。

   param_C param_gamma  rank_test_score  mean_test_score
22       1          10                1         0.972222

　0.972222なら悪くないaccuracyです。

　しかし、出力された画像は思っていたものとは違いました。

SVMの結果

　こわっ。どう考えてもこういうデータではないと思いますが、こういうことが現実に起こります。(-1.5, 0)あたりにデータが来たら、青かねずみ色のどちらかの色のグループに分類されてほしいところですが、実際は緑になってしまう訳です。

　ちょっと信頼ならないですね、SVM。

怖くない線形SVM

　ついでにいろいろな分類器を見てみます。

　LinearSVCをimportし、「# GridSearchCVの準備」以下「# 最良推定器をclfに代入」の上までを次のように書き換える。

    # GridSearchCVの準備
    svm = LinearSVC()
    params = {"C":10**np.arange(-3, 3, dtype=float)}
    gscv = GridSearchCV(svm, params, cv=8, iid=False, 
                        return_train_score=False,
                        verbose=1, n_jobs=-1)
    gscv.fit(X, y)

    # GridSearchCVの結果を表示する
    result_df = pd.DataFrame(gscv.cv_results_)
    print(result_df[["param_C",
                     "rank_test_score", "mean_test_score"]
                ][result_df["rank_test_score"]==1])

　結果。

  param_C  rank_test_score  mean_test_score
4      10                1         0.960317

　accuracyはわずかに下がるだけ。

線形SVMの結果

　はるかに納得感の高い結果になっています。こういうことがあるので、ほぼ線形分離可能なことがわかっているデータは、まずは線形なモデルで試すことがおすすめです。非線形なモデルで1%とかaccuracyをあげられるとしても、それで未知のデータに対して良い推定ができるかどうかは交差検証ではわからないのです*1。

怖いかどうか悩むランダムフォレスト

　LinearSVCと同様にRandomForestClassifierをimportし、同じ箇所を書き換えます。

    # GridSearchCVの準備
    clf = RandomForestClassifier(n_jobs=-1)
    params = {"n_estimators":10**np.arange(3),
              "min_samples_leaf":[1,2]}
    gscv = GridSearchCV(clf, params, cv=8, iid=False, 
                        return_train_score=False,
                        verbose=1, n_jobs=-1)
    gscv.fit(X, y)

    # GridSearchCVの結果を表示する
    result_df = pd.DataFrame(gscv.cv_results_)
    print(result_df[["param_n_estimators", "param_min_samples_leaf",
                     "rank_test_score", "mean_test_score"]
                ][result_df["rank_test_score"]==1])

　性能。SVMと比べると少し低下するかな。理由はよくわからないけど、木の本数10本で葉の最小サンプル数2のときが最高性能（本当になんで？　基本的に木が多い方が性能が高いはずなのだが・・・）。

  param_n_estimators param_min_samples_leaf  rank_test_score  mean_test_score
4                 10                      2                1         0.953373

　可視化。

ランダムフォレストの結果

　この状態ではrefitしているので、accuracyは1になるはずです。さて、ランダムフォレストの特徴は、決定木なので軸と90度の直線の組み合わせで決定境界が表現されることです。また、少ない分割でエントロピーが下がるような決定境界を追求していくので、データの全体の傾向はあまり見てくれません。

　それでもSVMよりはマシな感じでしょうか。

怖くない気がする多層パーセプトロン

　疲れてきたので最後。やり方は上と同じ。

    # GridSearchCVの準備
    clf = MLPClassifier(max_iter=3000)
    params = {"hidden_layer_sizes":[(5*x,) for x in range(1, 5)]}
    gscv = GridSearchCV(clf, params, cv=8, iid=False, 
                        return_train_score=False,
                        verbose=1, n_jobs=-1)
    gscv.fit(X, y)

    # GridSearchCVの結果を表示する
    result_df = pd.DataFrame(gscv.cv_results_)
    print(result_df[["param_hidden_layer_sizes",
                     "rank_test_score", "mean_test_score"]
                ][result_df["rank_test_score"]==1])

　同率一位が3つもありました。ま、どれでも良いか。

  param_hidden_layer_sizes  rank_test_score  mean_test_score
1                    (10,)                1          0.96627
2                    (15,)                1          0.96627
3                    (20,)                1          0.96627

多層パーセプトロンの結果

　まともな感じ。右側の決定境界はちょっと甘いかなー、という気はします。

　ただしパラメータ数の少ない隠れ層1層の多層パーセプトロンだからこれくらいの素直な結果に落ち着くのであって、深層学習は油断するとすぐ過学習するので注意が必要です。

SVMも怖くない！

　汎化重視のパラメータにすれば変なことにはならないので安心してください。基本的にCもgammaも低ければ低いほど過学習しづらくなります。

　今までと同じ箇所をこう書き換える。

    # SVMをfit
    clf = SVC(C=1, gamma=0.2)
    clf.fit(X, y)

汎化重視のパラメータのSVMの結果

　別に問題ないですよね。SVMそのものはこのようにパラメータで自由に分離超平面の複雑さを調整できる、優れた分類器です。

　ただし、交差検証で機械的にパラメータを決めてしまうとあまり良くない結果を招く可能性がある訳です。

現実的な話

　二次元で見ているから違和感があるのであって、高次元空間はサクサクメロンパン問題があるのでまた異なった挙動になります。

　↓サクサクメロンパン問題

次元の呪い、あるいは「サクサクメロンパン問題」 - 蛍光ペンの交差点[別館]

　直感的な理解は難しいと思いますが、こういうことが問題になるケースは少ないということはなんとなくわかります。

　とはいえ、「交差検証して最高性能だったSVMを投入したらトンチンカンな結果を出してくる」みたいなことが実際に起こらないとも限らないので、ある程度は注意した方が良いでしょう。

まとめ

　SVMはパラメータ設定を過学習するようにしたとき、他の手法と比べても群を抜いて変な結果になるのですが、何が駄目なんでしょうね。やっぱりカーネル使ってるからか。

　ちょっと注意が要るよなー、というお話でした。

*1:手持ちのデータの分布に依存する話ですが。未知のデータなんてないぜ！　といえるくらい豊富なデータがあればそれはそれで構いません

【python】sklearnのOneClassSVMを使って外れ値検知してみる

2018-12-14T11:11:26+09:00

はじめに

　OneClassSVMというものがあると知ったので使ってみます。

　「1クラスSVM？」と思われると思いますが、要するに異常検知・外れ値検出などで使う手法です。信頼区間を出すのに似ていますが、複雑な分布だったりそもそも分布が想定できないようなデータでも計算してくれるので、シチュエーションによっては役に立ちそうです。

　なお、わかりやすい記事があったので先に紹介しておきます。

異常検知のための One Class SVM - Qiita

実験

　異常検知・外れ値検出系で使える手法なので、センサデータの処理とか、為替や株価のアルゴリズム取引用の処理なんかをやると適当だと思いますが、私はそんなカッコいいデータは持っていません。

　なので、例によって例のごとく、irisをPCAで二次元に落としたデータを使います。

　使い方は簡単で、nuに異常値の割合を指定すれば良いようです。なんかドキュメントには意味深なことが書いてありますが、この理解で良さそうです。

　ちなみにデフォルトはnu=0.5なので、データの半数が異常値扱いになります。最初は、一体何事かと思いました。あと、predictすると正常値=1,異常値=-1という予測になります。

　ドキュメント

　sklearn.svm.OneClassSVM — scikit-learn 0.21.3 documentation

　コードは以下のとおりです。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.svm import OneClassSVM
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA

def main():
    iris = load_iris()
    pca = PCA(n_components=2)
    data = pca.fit_transform(iris.data)

    x = np.linspace(-5, 5, 500)
    y = np.linspace(-1.5, 1.5, 250)    
    X, Y = np.meshgrid(x, y)
    
    ocsvm = OneClassSVM(nu=0.1, gamma="auto")
    ocsvm.fit(data)
    df = ocsvm.decision_function(
        np.array([X.ravel(), Y.ravel()]).T).reshape(X.shape)
    preds = ocsvm.predict(data)

    plt.scatter(data[:,0], data[:,1], c=preds,
                cmap=plt.cm.RdBu, alpha=0.8)
    r = max([abs(df.min()), abs(df.max())])
    plt.contourf(X, Y, df, 10, vmin=-r, vmax=r,
                 cmap=plt.cm.RdBu, alpha=.5)
    plt.savefig("result.png")

if __name__ == "__main__":
    main()

　予測と決定関数を見るだけという手抜き。雰囲気はこれでわかると思うので、勘弁してください。

　なんかcontourfあたりでごちゃごちゃやっていますが、決定境界がcmapの中心と一致するように配慮しています。こうすることで、白色のあたり（というか青と赤の境界）が決定境界になります。

　余談ですが、このコードのためにlevelsをキーワード引数で指定しようとしたら、matplotlibのバグを踏みました。ひどい。

plt.contour levels parameter don't work as intended if receive a single int · Issue #11913 · matplotlib/matplotlib · GitHub

結果

　プロットされる図を示します。

result.png

　このように、お手軽に良さげな結果が得られます。分布の形状が複雑でもうまく推定できる訳です。良いですね。

まとめ

　SVMなので使いやすくて、うまく動くようです。手軽に良好な異常検知ができる手法としては、かなり便利だと思います。

【python】sklearnのFeatureAgglomerationを使ってみる

2018-12-10T03:56:41+09:00

はじめに

　FeatureAgglomerationは階層的クラスタリングを用いた教師なし次元削減のモデルです。特徴量に対して階層的クラスタリングを行い（つまり通常のサンプルに対するクラスタリングと縦横の向きが入れ替わる）、似ている特徴量同士をマージします。マージの方法はデフォルトでは平均のようです。

　使用例をあまり見かけませんが、直感的な次元削減方法なので何かしらの役に立つかもしれないと思って使ってみました。

sklearn.cluster.FeatureAgglomeration — scikit-learn 0.20.1 documentation

使い方

　パラメータは以下の通り。

class sklearn.cluster.FeatureAgglomeration(
    n_clusters=2, affinity=’euclidean’, memory=None, connectivity=None, 
    compute_full_tree=’auto’, linkage=’ward’, pooling_func=<function mean>)

　色々いじれるように見えますが、主要パラメータは2つだけです。

n_clusters

　PCAでいうところのn_componentsです。変換先の次元数を表します。

pooling_func

　似ている特徴量をマージする方法。callableが渡せます。何もしなければ平均が使われるので、平均より気の利いた方法を思いつく人以外はそのままで大丈夫です。

　あとは階層的クラスタリングのオプションが色々あります。それはそれで大切なものだと思いますが、今回は無視することにします。

実験

　もう何番煎じかわかりませんが、irisの2次元写像で試します。

import matplotlib.pyplot as plt

from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import FeatureAgglomeration

def main():
    iris = load_iris()

    pca = PCA(n_components=4)
    ss = StandardScaler()
    agg = FeatureAgglomeration(n_clusters=2)

    pca_X = pca.fit_transform(iris.data)
    agg_X = agg.fit_transform(
        ss.fit_transform(iris.data))

    print(pca.components_)
    print(agg.labels_)

    fig, axes = plt.subplots(nrows=1, ncols=2)
    axes[0].scatter(pca_X[:,0], pca_X[:,1], c=iris.target)
    axes[0].set_title("PCA")
    axes[1].scatter(agg_X[:,0], agg_X[:,1], c=iris.target)
    axes[1].set_title("FeatureAgglomeration\n{}".format(agg.labels_))
    plt.savefig("result.png")

if __name__ == "__main__":
    main()

　動作原理、目的と用途を考えると、事前にスケーリングしておいた方が恐らく無難です。

　printされた出力。

[[ 0.36138659 -0.08452251  0.85667061  0.3582892 ]
 [ 0.65658877  0.73016143 -0.17337266 -0.07548102]]
[0 1 0 0]

　FeatureAgglomerationは圧倒的に結果の解釈性が良いことがわかります。写像先の0次元目は元の0,2,3次元目の平均で*1、写像先の1次元目は元の1次元目ですね。こういうのはシチュエーション次第ですが、ちょっと嬉しいかもしれません。

　出力される画像。

プロットの結果

　概ねPCAと同等に使えています。うまく言葉で表現はできませんが、FeatureAgglomerationの方はなんとなくギザギザ感？みたいなものがあります。平均するとそうなる、というのがなんとなくわかる気もするし、わからない気もする。

考察

　結果の解釈性が良いのと、まがりなりにすべての特徴量の情報が結果に反映されるので、PCAより使いやすいシチュエーションはあると思います。分類前の次元削減とかで使ったときの性能とかは今回検討していませんが、たぶんそんなに良いということはないはず。

　あとドキュメントをあさっていたら、こんなページがあったので、

Feature agglomeration — scikit-learn 0.20.1 documentation

　真似してPCAでも同じものを出してみたら（コードはほとんど書き換えていないので省略。agglo = の行で代入するモデルをコメントアウトで切り替えて、あとlabels_の出力を外しただけです）、やっぱりFeatureAgglomerationはヘボかった（低次元で元の情報を保持することに関しては性能が低かった）です。

　10次元に落として元の情報をどこまで復元できるかという実験。

PCA

FeatureAgglomeration

　まあ、これは仕方ないか。

まとめ

　とにかく結果の解釈性の良さを活かしたい、とか、なにか特別な理由があって使う分には良いと思います。

*1:厳密にはどれか2つが先に平均されて、更に残りと平均されるはず。つまり3つの比重が違う順番はチェックしていないのでわかりませんが、children_属性をちゃんと読み取ればわかると思います

【python】sklearnのSparsePCAを使ってみる

2018-11-17T22:30:03+09:00

はじめに

　SparsePCAというものがあることを知ったので、使ってみようと思います。

SparsePCAとは？

　その名の通り、スパースな主成分分析です。スパースな主成分ベクトルを推定します。

Sparse PCA - Wikipedia

　原理などは理解しないで、カジュアルに使えるかどうか試してみるだけという趣旨です。なので「どうやって動いているの？」という質問には答えられません。許してください。

sklearnの実装

　きっちり存在しています（存在しなかったらこんな記事は書きませんが）。

sklearn.decomposition.SparsePCA — scikit-learn 0.20.1 documentation

　主要なパラメータとしては、以下のものがあります。

n_components

　PCAのと同じです。

alpha

　スパースPCAのキモで、L1正則化の強さを調整できます。

ridge_alpha

　こちらはtransformの際に使われるリッジ回帰（L2正則化）の正則化パラメータです。なんでリッジを使うのかは、実のところよくわかりません。

max_iter

　このパラメータがあるということは、最適化とか勾配法的なもので推定するのだな、というくらいに思っておきます。

normalize_components

　主成分ベクトルのノルムを1にするかどうか。Trueにしておくと良いと思います。

　結果に大きな影響を及ぼすのは上くらいだと思います。他のパラメータについてはドキュメントを参照してください。

実験

　今回はwineデータセットでやってみました。素のPCAでやった場合、alphaを0.5と5にした場合の結果をバイプロットで示します。

import matplotlib.pyplot as plt
from sklearn.datasets import load_wine
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA, SparsePCA

def biplot(X_2d, components, target, ax):
    r1 = 5
    r2 = 1.01
    for i, coef in enumerate(components.T):
        ax.arrow(0, 0, coef[0]*r1, coef[1]*r1, color='r')    
        ax.text(coef[0]*r1*r2, coef[1]*r1*r2, i, color='b', fontsize=8)

    ax.scatter(X_2d[:,0], X_2d[:,1], c=target, cmap="rainbow")

def main():
    wine = load_wine()
    ss = StandardScaler()
    X = ss.fit_transform(wine.data)

    pca = PCA(n_components=2)
    spca = SparsePCA(n_components=2,
                     max_iter=3000,
                     n_jobs=-1,
                     normalize_components=True)
    
    fig, axes = plt.subplots(figsize=(12, 6), nrows=1, ncols=3)

    X_2d = pca.fit_transform(X)
    biplot(X_2d, pca.components_, wine.target, axes[0])
    axes[0].set_title("PCA")

    for i,alpha in zip([1, 2], [0.5, 5]):
        spca.set_params(alpha=alpha)
        X_2d = spca.fit_transform(X)
        biplot(X_2d, spca.components_, wine.target, axes[i])
        axes[i].set_title("SPCA alpha={:.2f}".format(alpha))
    plt.savefig("result.png")

    # 図と突き合わせて確認するために特徴量の名前を出力しておく
    for i, name in enumerate(wine.feature_names):
        print(i, name)

if __name__ == "__main__":
    main()

　max_iterをきもち高めにしましたが、結果は数秒程度で出ました。

result.png

0 alcohol
1 malic_acid
2 ash
3 alcalinity_of_ash
4 magnesium
5 total_phenols
6 flavanoids
7 nonflavanoid_phenols
8 proanthocyanins
9 color_intensity
10 hue
11 od280/od315_of_diluted_wines
12 proline

　とりあえず、PCAの結果とSparsePCAの結果で左右が反転しているのに注意。

　あとは見ての通りで、alpha=0.5で一部の係数が主成分にべたっと張り付くようになり、alpha=5では大半の係数が主成分に張り付いています。これがSparsePCAの効果で、結果の解釈が容易になるということらしいです（この次元数だとあまり威力はありませんが、高次元では活躍しそうです）。

　ワインにはあまり詳しくないので、今回は結果を細かく解釈することはしませんが……。

まとめ

　使えることがわかりました。

GridSearchCV『の』パラメータ・チューニング高速化中心に

2018-06-18T22:11:44+09:00

はじめに

　機械学習でパラメータ・チューニングをしたい場合、グリッドサーチを行うのが定石とされています。sklearnではグリッドサーチはGridSearchCVで行うことができます。

sklearn.model_selection.GridSearchCV — scikit-learn 0.21.2 documentation

　それで何の問題もないかというと、さにあらず。

　グリッドサーチは計算コストの高い処理ですから*1、素直に書くとデータとアルゴリズム次第ではとんでもない処理時間がかかります。

　もちろん「寝ている間、出かけている間に回すから良い」と割り切るという方法もありますが、可能なら速くしたいですよね。

　そうすると、パラメータ・チューニングのために使うGridSearchCV『の』パラメータを弄り回すという本末転倒気味な目に遭います。そういうとき、どうしたら良いのかを、この記事で書きます。

　先に結論を言ってしまうと、本質的に計算コストの高い処理なので、劇的に速くすることは不可能です。それでも、ちょっとの工夫で2倍程度なら速くすることができます。その2倍で救われる人も結構いると思うので*2、単純なことですがまとめておきます。

はじめに
下準備とベースライン
cvを指定する（効果：大）
return_train_score=Falseする（効果：それなり）
まとめ
それでも時間がかかりすぎるときは

下準備とベースライン

　とりあえず、何も考えずに「GridSearchCVをデフォルトパラメタで使ってみた場合」の時間を測ります。

　そのためには適当なタスクを回してやる必要がありますが、今回はPCA+SVMでdigitsの分類でもやってみることにします。

　コードはこんな感じです。

import timeit
import pandas as pd
from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.svm import SVC
from sklearn.pipeline import Pipeline
from sklearn.model_selection import GridSearchCV

digits = load_digits()
svm = SVC()
pca = PCA(svd_solver="randomized")
pl = Pipeline([("pca", pca), ("svm", svm)])

params = {"pca__n_components":[10, 15, 30, 45],
          "svm__C":[1, 5, 10, 20], 
          "svm__gamma":[0.0001, 0.0005, 0.001, 0.01]}

def print_df(df):
    print(df[["param_pca__n_components",
              "param_svm__C", "param_svm__gamma", 
              "mean_score_time", 
              "mean_test_score"]])

def main1():
    clf = GridSearchCV(pl, params, n_jobs=-1)
    clf.fit(digits.data, digits.target)
    df = pd.DataFrame(clf.cv_results_)
    print_df(df)
    
if __name__ == "__main__":
    print(timeit.timeit(main1, number=1))

　色々なテクニックを使っているコードなので多少解説すると、とりあえずPipelineを使っています。また、GridSearchCV.cv_results_はそのままpandas.DataFrameに変換できる辞書として扱えることも利用しています。

www.haya-programming.com

　digits, svm, pca, pl, paramsの変数はmain関数の外でグローバル変数として作っていますが、これはあとでmain2とかmain3を作って使い回すための処置です。

　あと、速くするために必要と思われる常識的なこと（PCAのsvd_solver="randomized"とか、GridSearchCVのn_jobs=-1とか）はすでに実施しています。

　そんなことより本当にやりたいことは、この処理にどれだけ時間がかかるかを知ることです。そのために、timeitを使って時間を計測しています。

timeit --- 小さなコード断片の実行時間計測 — Python 3.7.4 ドキュメント

　さて、私の環境（しょぼいノートパソコン）ではこのプログラムの実行には42.2秒かかりました。

　これをベースラインにします。ここからどれだけ高速化できるかが今回のテーマです。

cvを指定する（効果：大）

　さて、GridSearchCVにはcvというパラメータがあります。default=3であり、この設定だと3分割交差検証になります。交差検証について理解していれば、特に不自然なところはないと思います。

　これを2にしてみます。交差検証できる最低の数字です。こうすると、

交差検証のループ回数が3回→2回になり、それだけで1.5倍速くなる
チューニング対象のモデルの計算量が学習データサイズnに対してO(n)以上なら、それ（nが小さくなること）によるご利益もある。なお予測データサイズmに対する予測時間は普通O(m)なので、影響はない

　この相乗効果で、高速化が期待できます。

　この方法のデメリットは学習データを減らしたことで性能が低めになることですが、チューニングのときはパラメータの良し悪し（スコアの大小関係）がわかれば良いので、あまり問題になりません。とにかくやってみます。

def main2():
    clf = GridSearchCV(pl, params, cv=2, n_jobs=-1)
    clf.fit(digits.data, digits.target)
    df = pd.DataFrame(clf.cv_results_)
    print_df(df)

if __name__ == "__main__":
    # print(timeit.timeit(main1, number=1))
    print(timeit.timeit(main2, number=1))

　上のコードと重複する部分は削ってあります。見比べると、ほとんど変わっていないことが、おわかりいただけるかと思います。

　この処置で、処理時間は28.0秒に改善しました。ちょっといじっただけで、2/3くらいに改善してしまった訳です。そして「mean_test_score」はやはり全体的に低くなりましたが、傾向は同じでした。よってパラメータチューニングには使えます。

return_train_score=Falseする（効果：それなり）

　さて、GridSearchCVはデフォルトの設定ではreturn_train_score='warn'になっています。「'warn'って何さ？」というと、こんな警告を出します。

FutureWarning: You are accessing a training score ('std_train_score'), which will not be available by default any more in 0.21. If you need training scores, please set return_train_score=True

　return_train_scoreは要するに学習データに対するスコアを計算するかどうかを指定できる引数です。この警告は割とくだらないことが書いてあるのですが、将来的にはこれがdefault=Falseにされるという警告です。

　基本的に、パラメータチューニングで見たいのはテストデータに対するスコアであるはずです。なのに、現在のデフォルト設定では学習データに対する評価指標も計算されてしまいます。

　これは無駄なので、return_train_score=Falseすると学習データに対する評価指標の計算分の計算コストをケチれます。予測時間なんてたかが知れていますが、それでも一応やってみましょう。

def main3():
    clf = GridSearchCV(pl, params, cv=2,
                       return_train_score=False,
                       n_jobs=-1)
    clf.fit(digits.data, digits.target)
    df = pd.DataFrame(clf.cv_results_)
    print_df(df)
    
if __name__ == "__main__":
    # print(timeit.timeit(main1, number=1))
    # print(timeit.timeit(main2, number=1))
    print(timeit.timeit(main3, number=1))

　この措置によって、処理時間は22.1秒まで短縮されました。ベースラインと比較すると1/2強の時間で済んでいる訳です。

まとめ

　この記事では

cv=2にする
return_train_score=Falseにする

　という方法で、パラメータチューニングの機能を損なわないまま2倍弱の速度の改善を実現しました。

工夫なし	cv=2	cv=2＆return_train_score=False
42.2秒	28.0秒	22.1秒

　このテクニックはきっと皆さんの役に立つと思います。

それでも時間がかかりすぎるときは

　そもそもグリッドサーチしようとしているパラメータ候補が多すぎる可能性があります。

　たとえば、3つのパラメータでそれぞれ10個の候補を調べるとなると、10*10*10=1000回の交差検証が行われます。いつまで経っても終わらない訳です。

　今回の記事では4*4*4=64回としましたが、これでもけっこう多い方だと思います。それでも解こうとしている問題が単純なので、デフォルトパラメータでも1分以内には処理できていますが、ちょっと重いモデルにちょっと多量のデータを突っ込んだりすると、もうダメです。何十分とか何時間もかかってしまいます。

　そういう場合、まずは粗いステップ（少ないパラメータ候補数）でざっくりパラメータチューニングしてしまい、どの辺りのパラメータが良いのかわかったら、その周辺に絞ってもう一回パラメータチューニングを行います。こういうのを二段グリッドサーチと言ったりします。

　あるいはベイズ最適化とか、他のアルゴリズムに走るのも一つの手かもしれません。

　粗いグリッドである程度チューニングしてから、RandomizedSearchCVを使うというのもいい手だと思います。

www.haya-programming.com

*1:なにせすべての組み合わせを計算する

*2:たとえば「今から回して、朝までにデータを出さないと教授への報告が間に合わないんだ！」みたいな状況を想定しています

【python】MeanShiftのbandwidthを変えるとどうなるか実験してみた

2018-05-28T22:18:19+09:00

　前回の記事ではMeanShiftクラスタリングを試してみました。

www.haya-programming.com

　このMeanShiftにはbandwidthというパラメータがあり、クラスタ数を決定する上で重要な役割を果たしているはずです。

　いまいち結果に納得がいかないというとき、bandwidthをいじって改善が見込めるのかどうか確認してみます。

プログラム

　例によってirisとwineで比較。簡単に書きました。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from sklearn.datasets import load_iris, load_wine
from sklearn.cluster import MeanShift, estimate_bandwidth
from sklearn.decomposition import PCA

def process(dataset, name):
    origin_bandwidth = estimate_bandwidth(dataset.data)
    rates = np.logspace(np.log10(0.2), np.log10(5), 11)
    fig, axes = plt.subplots(nrows=3, ncols=4, figsize=(24,18))

    PCA_X = PCA().fit_transform(dataset.data)
    for target in range(3):
        axes[0,0].scatter(PCA_X[dataset.target==target, 0],
                        PCA_X[dataset.target==target, 1],
                        c=cm.Paired(target/3))
    axes[0,0].set_title("original label", fontsize=28)

    for r, ax in zip(rates, axes.ravel()[1:]):
        ms = MeanShift(bandwidth=r*origin_bandwidth, n_jobs=-1)
        y = ms.fit_predict(dataset.data)
        n_cluster = ms.cluster_centers_.shape[0]
        for target in range(n_cluster):
            ax.scatter(PCA_X[y==target, 0],
                       PCA_X[y==target, 1],
                       c=cm.Paired(target/n_cluster))
        ax.set_title("r:{0:.3f} b:{1:.3f}".format(
            r, origin_bandwidth), fontsize=28)
    fig.savefig(name+".png")

def main():
    iris = load_iris()
    wine = load_wine()

    process(iris, "iris")
    process(wine, "wine")

if __name__ == "__main__":
    main()

　bandwidthをsklearn.cluster.estimate_bandwidthの推定値（デフォルトで用いられる値）の1/5倍から5倍まで変化させ、結果をプロットします。

結果

　プロットされた結果を示します。

　結果の図の見方は、まずタイトルが

　sklearn.cluster.estimate_bandwidthによる推定値

　かけた比率

　という風に対応しており、あとは便宜的に2次元上に主成分分析で写像した散布図が、クラスタごとに色分けされて出ています。一枚目が本来のクラスに基づく色分け、r=1の図が推定値による色分けです。

　まずiris。

iris.png

　きれいに元通りになるrは今回見た中にはありませんでした。クラスタ数的にはr=0.525とr=0.725の間くらいで3クラスタになりそうですが、この図を見るとそれでうまく元通りまとまるかは疑問です。

　次にwine。

wine.png

　こちらもうまく元通りにはならないようです。そもそもデータが悪いという話はあると思います。

結論

　確かにクラスタ数は変わるが、クラスタリングの良し悪しが改善するかはなんともいえないですね。

　データをスケーリングしたり、もっと色々頑張ると改善は見込めるかもしれません。

【python】sklearnのMeanShiftクラスタリングを試してみる

2018-05-27T18:23:20+09:00

はじめに

　MeanShiftはクラスタリングアルゴリズム。クラスタ数を自動で決定してくれるという長所がある。

　理論的には最急降下法で各クラスタの極大点を探していく感じらしいです。わかりやすい解説があったので、リンクを張っておきます（ただし私自身はすべては読み込めていない）。

Mean Shift Clustering

　このMeanShiftはsklearnに実装されているので、簡単に試してみることができます。

　sklearn.cluster.MeanShift — scikit-learn 0.20.1 documentation

　sklearnのトイデータで遊んでみましょう。

はじめに
使い方
実験
- プログラム
- 結果
結論

使い方

　いつものsklearnのモデルです。fitしてpredictするだけ。

　いつだったかFuzzy C-Meansをやったときは苦労しましたが、とりあえずそんな心配は要りません。

　となると気になるのはパラメータですが、指定しなくても

bandwidth

　勝手に推定される

seeds

　乱数のシードなので勝手に決まる。指定するときは[n_samples, n_features]が必要。

bin_seeding

　よくわからないけど、初期値の選び方みたいな。上のと関係がありそう。Trueにすると速くなるらしい。デフォルトのFalseの方がアルゴリズムとしては厳密なはず（よくわからんけど）。

min_bin_freq

　これも上のと関係がありそう。わかるようなわからないような感じ。

cluster_all

　すべての点をクラスタリングするかどうか。default=Trueなので敢えてFalseにする理由は・・・（高速化なんだろうな）。

n_jobs : int

　いつもの並列化数

　本質的な挙動に関わるのはbandwidthで、あとは高速化のために計算を端折るための引数がいっぱいあるだけっぽいです。

　そしてbandwidthも勝手に推定してくれるので、敢えて指定する必要性を感じません（推定の良し悪しがどうかという議論はありますが）。

　今回はn_jobs以外すべてデフォルトでやってみます。

実験

　iris, wineで見てみる。せっかくなのでK-Meansと比較してみます。ついでに、入力をスケーリングすると結果が変わるかも見ます。

プログラム

　いろいろ手抜きをしています。が、とにかく結果は出ます。

import matplotlib.pyplot as plt
from sklearn.datasets import load_iris, load_wine
from sklearn.cluster import MeanShift, KMeans
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

def main():
    iris = load_iris()
    wine = load_wine()

    kmeans = KMeans(n_clusters=3, n_jobs=-1)
    mean_shift = MeanShift(n_jobs=-1)
    s_kmeans = Pipeline([("scaler", StandardScaler()), 
                         ("kmeans", KMeans(n_clusters=3, n_jobs=-1))])
    s_mean_shift = Pipeline([("scaler", StandardScaler()), 
                             ("meanshift", MeanShift(n_jobs=-1))])


    # iris and wine
    pca = PCA(n_components=2)
    for dataset, dataset_name in zip([iris, wine], ["iris", "wine"]):
        fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(20,8))
        axes = axes.ravel()

        PCA_X = pca.fit_transform(dataset.data)
        origin_y = dataset.target
        km_y = kmeans.fit_predict(dataset.data)
        ms_y = mean_shift.fit_predict(dataset.data)
        s_km_y = s_kmeans.fit_predict(dataset.data)
        s_ms_y = s_mean_shift.fit_predict(dataset.data)
        n_clusters = [3, 3, mean_shift.cluster_centers_.shape[0],
                      3, s_mean_shift.named_steps.meanshift.cluster_centers_.shape[0]]

        for i, (y, name, n_cluster) in enumerate(
                zip([origin_y, km_y, ms_y, s_km_y, s_ms_y], 
                    ["original", "k-means", "mean-shift",
                     "scaling+k-means", "scaling+mean-shift"],
                    n_clusters)):

            for target in range(n_cluster):
                axes[i].scatter(PCA_X[y==target, 0],
                                PCA_X[y==target, 1],
                                c="rgb"[target])
            axes[i].set_title("{0} n_cluster:{1}".format(name, n_cluster))
        plt.savefig("{0}.png".format(dataset_name))

if __name__ == "__main__":
    main()

結果

　まずirisの結果から。

iris.png

　MeanShiftは2クラスタと解釈しているようです。本来のデータとは異なりますが、敢えて人の目で見ると妥当な結果な気もします。この場合、スケーリングによる変化は微々たるものです。

　次に、wineの方。

wine.png

　一見するとoriginalが悪すぎるように見えますが、PCAでむりやり二次元に写像しているためです。クラスタリング自体は写像前のオリジナルの空間で行っているため、影響はありません。

　この場合、合格と言って良いのは入力をスケーリングしたKMeansだけで、あとはダメダメです。特徴量の次元数が大きいから、うまく動いていないのでしょうか。ちょっと謎。

結論

　良いか悪いかの判断はつきかねますが、できることはわかりました。

　たぶんbandwidthを変えるとコロコロ結果が変わるのでしょう。どんな感じで変わるのかは、今後気が向いたときに検証しようと思います。

　→やりました。
www.haya-programming.com

【python】複数の特徴をまとめるFeatureUnion

2018-05-15T14:41:24+09:00

　単一の入力データから、複数の処理方法で幾つもの異なる特徴量が得られる・・・というシチュエーションがある。

　この場合、「どれが最善か」という観点でどれか一つを選ぶこともできるけど、そうすると他の特徴量の情報は捨ててしまうことになる。総合的な性能では他に一歩譲るが、有用な情報が含まれている特徴量がある・・・というような場合は、ちょっと困る。

　こういう状況で役に立つのがFeatureUnion。特徴抽出や次元削減などのモデルを複数まとめることができる。

　結果はConcatenateされる。Concatenateというのがわかりづらい人もいると思うけど、たとえば手法1で10次元、手法2で20次元の特徴量ベクトルが得られたら、これをそのまま横に繋げて30次元のベクトルとして扱うということ。

sklearn.pipeline.FeatureUnion — scikit-learn 0.20.1 documentation

　ちなみに、こいつはsklearn.pipeline以下に存在する。Pipelineの兄弟みたいな扱い。引数の渡し方とかもほとんど同じである。

　簡単に試してみよう。digitsの分類を行うことにする。PCA+GaussianNB, LDA+GNB, FeatureUnion(PCA, LDA)+GNBの3パターンでスコアを見比べる。

import warnings
warnings.filterwarnings('ignore')

from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn.naive_bayes import GaussianNB
from sklearn.pipeline import Pipeline, FeatureUnion
from sklearn.model_selection import cross_validate, StratifiedKFold

def main():
    digits = load_digits()
    
    pca = PCA(n_components=30)
    lda = LDA()
    gnb = GaussianNB()
    
    pca_gnb = Pipeline([("pca", pca), ("gnb", gnb)])
    lda_gnb = Pipeline([("lda", lda), ("gnb", gnb)])
    pca_lda_gnb = Pipeline([("reduction", FeatureUnion([("pca", pca),
                                                        ("lda", lda)])),
                            ("gnb", gnb)])

    scoring = {"p": "precision_macro",
               "r": "recall_macro",
               "f":"f1_macro"}

    for name, model in zip(["pca_gnb", "lda_gnb", "pca_lda_gnb"], 
                           [pca_gnb, lda_gnb, pca_lda_gnb]):

        skf = StratifiedKFold(shuffle=True, random_state=0)
        scores = cross_validate(model, digits.data, digits.target,
                                cv=skf, scoring=scoring)
        
        p = scores["test_p"].mean()
        r = scores["test_r"].mean()
        f = scores["test_f"].mean()
        print(name)
        print("precision:{0:.3f} recall:{1:.3f} f1:{2:.3f}".format(p,r,f))

if __name__ == "__main__":
    main()

　結果は、

pca_gnb
precision:0.947 recall:0.944 f1:0.945
lda_gnb
precision:0.955 recall:0.953 f1:0.953
pca_lda_gnb
precision:0.959 recall:0.957 f1:0.957

　ちょっと微妙だけど、誤差ではないみたい。このように比較的手軽に性能を改善できることがわかる（効くかどうかはケースバイケースだけど）。

複数の目的変数で回帰を行う方法

2018-05-07T11:51:44+09:00

はじめに

　回帰分析を行う際、複数の目的変数に対して回帰をしたい場合があります。普通のモデルではできないのでちょっと面食らいますが、やり方は色々あるようです。

はじめに
目的変数の数だけ回帰モデルを作る方法
複数の目的変数に対応したモデルを使う
まとめ

目的変数の数だけ回帰モデルを作る方法

　単純に考えると、一つの目的変数を出力する回帰モデルを目的変数の数だけ用意してやれば、所要を達しそうです。

　python+sklearnを使えば、これに対応したモデルが最初から用意されています。

sklearn.multioutput.MultiOutputRegressor — scikit-learn 0.20.2 documentation

　コンストラクタには好きな回帰モデルを渡してあげることができます。それが目的変数の数だけコピーされ、内部で束ねられて回帰に使われます*1。

複数の目的変数に対応したモデルを使う

　上の方法は単純ですが、回帰モデルの中には自然に複数の出力に対応しているものもあります。

　そういったモデルを使うことにどんなメリットがあるのか？　というと、まず目的変数の数だけ回帰モデルを作るのと比べて無駄が減るので、計算コストがケチれる可能性があります（あくまでも「可能性」の話）。

　また、複数存在する目的変数の間に何らかの相関性があれば、それも踏まえて上手く学習することでモデルの性能が上がる可能性があります（こちらもあくまでも「可能性」）。

　そういった複数の目的変数に対応したモデルを幾つか紹介します。すべては網羅しきれないので、その点はご承知ください。

正準相関分析

　正準相関分析はこの手の話で出てくる代表的なモデルです。単純な手法ですが、けっこう奥深いといえば奥深いです。

　参考（過去に書いた記事）：【python】正準相関分析（Canonical Correlation Analysis）を試してみる - 静かなる名辞

　これの良いところは、説明変数と目的変数*2のそれぞれでPCAみたく新たな軸を張り、次元削減を行ってくれることです。説明変数ン百次元、目的変数20次元みたいなケースだったとしても、次元削減の効果で「わかりやすい」結果が得られる可能性があります。つまり、現象を説明するモデルとしては非常に適しています。

　欠点は、回帰モデルとして考えると性能が高いと言えるかは微妙なこと、非線形への対応は基本的にはないことです。カーネルPCAみたくカーネル法で非線形対応させたモデルもありますが、良さげなライブラリが見当たらないのと、そこまでするなら他の手法を使いたいという気持ちがあるので紹介しません。

　sklearnのモデルはこれです。上に書いた通りカーネル正準相関の実装はありません。

sklearn.cross_decomposition.CCA — scikit-learn 0.20.2 documentation

　predictメソッドでXからYを予測できるので、普通に回帰に使えます。　

　入出力が割と線形なデータで、「説明」を重視したいときは使えると思います。

ランダムフォレスト回帰

　なぜかランダムフォレスト回帰は複数出力に対応しています。解説論文を見つけたので貼っておきます。興味のある方はどうぞ（私は読んでいません）。

　とにかく使いたければsklearnのRandomForestRegressorはそのまま使えます。目的変数も説明変数と同様に配列で入れてあげてください。

3.2.4.3.2. sklearn.ensemble.RandomForestRegressor — scikit-learn 0.20.2 documentation

多層パーセプトロン（ニューラルネットワーク回帰）

　ニューラルネットですからできて当然。複数出力にするためにやることといったら出力層ユニット数を増やすだけですから、一番シンプルかもしれません。これもsklearnのがそのまま使えます。

sklearn.neural_network.MLPRegressor — scikit-learn 0.20.2 documentation

まとめ

　複数の目的変数に対して回帰を行う場合について、2種類の方法を説明しました。

単純に目的変数の数だけ回帰モデルを用意する方法
複数の目的変数を出力できるモデルを用いる方法

　どちらが良いかは一概には言えません。データや目的に応じて、あるいは実際に走らせてみて評価指標や計算コストを勘案して考える必要があります。複数の目的変数に最初から対応したモデルの方が良いような気もしますが、そうとも言えないんじゃという話もあったりします。

　でもまあ、色々な選択肢があることは良いことです。いろいろ勘案して選べば良いでしょう。適当ですがこんな感じでシメます。

*1:ということだと思う・・実装読んでいないので断言しかねます

*2:回帰の記事なのでそう呼ぶが、ぶっちゃけ妥当ではない。CCAの枠組みではどっちがどっちでも大して構わないのだし

【python】SOMのライブラリSomocluはかなりおすすめ

2018-04-07T16:12:49+09:00

　SOM(Self-organizing maps：自己組織化写像)は割と古めの、データの可視化手法です（それ以外にも使えると思いますが）。

　今回はpythonのSOMライブラリSomocluを使ってみたら、けっこう良かったというネタです。

SOMの概要
ライブラリがない
それでも頑張ってググった
使ってみた
今どきSOMなんか使うの？（蛇足パート）
まとめ

SOMの概要

　昨今は深層学習が流行りですが、SOM、自己組織化写像は敢えて言えば単層学習とでも言うべきでしょうか。平面上だったり立体状（まあ理屈の上では何次元でも定義できる）に並べたニューロンにデータをマッピングします。それ以上の説明はwikipediaとか、ググれば色々出てくるページを読んでください。

wikipedia

自己組織化写像 - Wikipedia

九州工業大学大学院の先生が書いた読みやすかったページ

http://www.brain.kyutech.ac.jp/~furukawa/data/SOMtext.pdf

わかりやすい解説

子供でもわかる「自己組織化マップ」

ライブラリがない

　SOM、けっこう面白い性質があるみたいなのて使ってみたいのですが、ググってみるとpythonで使えそうなライブラリがとにかくあまり出てきません。

SOMPY

　申し訳ないけど、ちょっと使いづらかった。というかインストールしても挙動が変な感じだった。
GitHub - sevamoo/SOMPY: A Python Library for Self Organizing Map (SOM)

sompy

　日本人の方が実装されたようです。率直に言って「作ってみた」レベルで、実用にはどうかという感じ
自己組織化マップ(SOM)のPythonライブラリsompyを公開しました - 俺とプログラミング

PyMVPA

　多変量解析のためのそれなりに大きいライブラリで、SOMも実装されている。これが使えればよかったのだと思うが、python2系のサポートしかないので没・・・。
Self-organizing Maps — PyMVPA 2.6.1.dev1 documentation

　他にも色々あったのですが、割愛。古い手法なので、敢えて作ろうという人がいないのかな・・・。

　というか、SOMでググると「実装してみた」系の記事はたくさん出てくるのに、まともに使えるライブラリは出てこないというの、かなり異常というか残念というか・・・。

それでも頑張ってググった

　Somocluというのを見つけました。

Introduction — Somoclu 1.7.5 documentation

　ウリの部分を適当に訳したり訳さなかったりしつつ抜粋

OpenMPとCUDAがサポートされていてGPUでも計算できる
当然マルチプラットフォームでLinux, macOS, and Windowsでサポートされている
「Planar and toroid maps」平面とドーナツみたいな形のSOM両方が作れる
「Rectangular and hexagonal grids」四角と六角形がいける
「Gaussian or bubble neighborhood functions」近傍の計算を効率化する系のがある
「Visualization of maps, including those that were trained outside of Python.」
マップの初期化にはPCAが使える

　すごく良さそう。あと、pythonに依存しないツールでコマンドラインから直接コマンドで叩けます。pythonバインディングもあるよ、という位置づけ。真剣に開発されてる感じです。

使ってみた

　とりあえず使ってみました。SOMの可視化結果でよく見るU-matrixという奴を出します。以下のコードで動きました。

# coding: UTF-8
import numpy as np

from somoclu import Somoclu
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA

def main():
    # データを読み込む
    dataset = load_iris()
    X = dataset.data
    y = dataset.target
   
    # SOMに入れる前にPCAして計算コスト削減を測る（iris程度では無駄） 
    pca = PCA(n_components=0.95) 
    X = pca.fit_transform(X)

    # SOMの定義
    n_rows = 16
    n_cols = 24
    som = Somoclu(n_rows=n_rows, n_columns=n_cols,
                  initialization="pca", verbose=2)

    # 学習
    som.train(data=X, epochs=1000)

    # U-matrixをファイル出力
    som.view_umatrix(labels=y, bestmatches=True,
                     filename="umatrix.png")

if __name__ == "__main__":
    main()

　説明不要な感じ。コードも直感的だし、特に不満がないです。

　こんな画像が出てきます。

U-matrix

　この画像の見方は色の濃淡が重要で、色の明るい部分は相対的に縮尺が縮んでおり、逆に暗い部分は縮尺が相対的に大きい訳です。PCAで可視化した結果を参考に貼っておきます。

PCAによるirisの可視化結果

　紫がラベル0に、緑と黄色が1と2に対応している訳です。SOMを使うと、このようにデータの構造を捉えることができます。

　使いやすいし動作もまともだし、Somocluは素晴らしいライブラリです。SOMが必要になったら積極的に使っていきたいところ。

今どきSOMなんか使うの？（蛇足パート）

　t-SNEみたいなよくできた手法があるのに今更SOM？　と思う方もおられるかと思いますが、SOMはSOMでメリットがあると感じています。

　というのは、t-SNEはけっきょくパラメタに依存するし、ミクロな構造を捉えるのは得意でもマクロな構造はどこまで正しいのか？　という問題があるからです。

　例として、digitsを可視化してみます。

# coding: UTF-8
import numpy as np

from sklearn.datasets import load_digits
from sklearn.manifold import TSNE
from sklearn.decomposition import PCA
from somoclu import Somoclu
import matplotlib.pyplot as plt

def main():
    print("loading data")
    digits = load_digits()
    pca = PCA(n_components=0.95)
    pca_data = pca.fit_transform(digits.data)

    # tsneで可視化
    print("tsne")
    tsne = TSNE()
    X = tsne.fit_transform(pca_data)
    fig, ax = plt.subplots()
    plt.scatter(X[:,0], X[:,1], c=digits.target/10)
    
    i = 0
    for xy, l in zip(X, digits.target):
        if i%8 == 0: # 描画されるtextが多いと汚いので省く
            ax.annotate(l, xy=xy)
        i += 1
    plt.savefig("tsne_digits.png")

    # somで可視化
    print("som")
    # データを適当に省く
    sample_index = np.random.choice(X.shape[0], 400, replace=False)
    sample_X = pca_data[sample_index]
    sample_y = digits.target[sample_index]

    # som
    som = Somoclu(n_rows=30, n_columns=40,
                  initialization="pca")
    som.train(data=sample_X, epochs=1000)
    som.view_umatrix(labels=sample_y, bestmatches=True,
                     filename="som_digits.png")

if __name__ == "__main__":
    main()

t-SNEで可視化したdigits

SOMで可視化したdigits

　一見するとt-SNEは同じラベルごとにまとまっていて綺麗なんですが、形の似ている数字が近くに来るのはむしろSOMの方という気もします。0の周りに5,6,9が来るというのは（数字の形を考えると）妥当そうですね。主観的になってしまいますが、SOMも捨てたものではないという気がします。

まとめ

　SOMとSomocluは良いのでみんな使おう。

【python】sklearnのPCAでloading（主成分負荷量）を計算する

2018-03-31T01:24:28+09:00

　PCA（主成分分析）のloading*1がほしいときがあります。

　sklearnでは一発では出ません。

　ドキュメントはここ。
sklearn.decomposition.PCA — scikit-learn 0.21.2 documentation

PCA.components_は確かにあるけど・・・
loadingを計算しよう
罠だった
共分散行列のときはどうするのか
loadingを使うと何が良いのか

PCA.components_は確かにあるけど・・・

　結論から先に言うと、PCA.components_はノルム1の固有ベクトルです。ノルムを測ってみましょう。

>>> import numpy as np
>>> from sklearn.datasets import load_iris
>>> from sklearn.decomposition import PCA
>>> iris = load_iris()
>>> pca = PCA(n_components=2)
>>> pca.fit(iris.data)
PCA(copy=True, iterated_power='auto', n_components=2, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)
>>> pca.components_
array([[ 0.36158968, -0.08226889,  0.85657211,  0.35884393],
       [ 0.65653988,  0.72971237, -0.1757674 , -0.07470647]])
>>> np.linalg.norm(pca.components_, axis=1)
array([1., 1.])

　まあ、loadingも固有ベクトルには違いないのですが、ノルムを整えてやる必要があります。

loadingを計算しよう

　教科書などによく書いてあることですが、第主成分に対応する元の変数のloadingは次の式で出せます。

　は固有値。 eigenvectorは固有ベクトルで、元の変数の数だけ次元がありますから、これで良いわけです（雑な説明ですが・・・）。

　ということで、pythonで同様にやってみましょう。固有値はexplained_varianceに入っています。

>>> pca.components_*np.c_[np.sqrt(pca.explained_variance_)] # 縦ベクトルに変換する必要がある
array([[ 0.74322652, -0.16909891,  1.76063406,  0.73758279],
       [ 0.32313741,  0.35915163, -0.08650963, -0.03676921]])

　できました。これがloadingです。・・・と思ったけど、1を超えちゃってますね。なんてこった。

罠だった

　固有値は分散なので、データのスケールに依存します。

　とりあえず入力をスケーリングしてみよう。上の式は相関行列から行くときのものでした。なのでこれで平気なはず。

>>> from sklearn.preprocessing import StandardScaler as SS
>>> ss = SS()
>>> data = ss.fit_transform(iris.data)
>>> pca.fit(data)
PCA(copy=True, iterated_power='auto', n_components=2, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)
>>> pca.components_*np.c_[np.sqrt(pca.explained_variance_)]
array([[ 0.89421016, -0.45081822,  0.99500666,  0.96822861],
       [ 0.35854928,  0.89132754,  0.02031465,  0.06299656]])

　1を超えなくなくてめでたし、ということよりも、数字が変わったことの方が問題で、これで本当に正しいのかという疑念が生じてきました。

　確認のために元の特徴と主成分の相関係数を直接測ってみます。

>>> X = pca.fit_transform(data)
>>> np.corrcoef(np.hstack([iris.data, X]), rowvar=False)
array([[ 1.00000000e+00, -1.09369250e-01,  8.71754157e-01,  8.17953633e-01,  8.91224479e-01,  3.57352114e-01],
       [-1.09369250e-01,  1.00000000e+00, -4.20516096e-01,  -3.56544090e-01, -4.49312976e-01,  8.88351481e-01],
       [ 8.71754157e-01, -4.20516096e-01,  1.00000000e+00,  9.62757097e-01,  9.91684422e-01,  2.02468206e-02],
       [ 8.17953633e-01, -3.56544090e-01,  9.62757097e-01,  1.00000000e+00,  9.64995787e-01,  6.27862218e-02],
       [ 8.91224479e-01, -4.49312976e-01,  9.91684422e-01,  9.64995787e-01,  1.00000000e+00,  2.08904471e-17],
       [ 3.57352114e-01,  8.88351481e-01,  2.02468206e-02,  6.27862218e-02,  2.08904471e-17,  1.00000000e+00]])

　下の二行の4列目までを見てください。微妙に誤差があるような気はしますが（小数点以下3桁でずれてきてるので微妙ってほど微妙でもないが）、たぶん同じ数字になっています。

　微妙な誤差については、丸め誤差などが累積した、実は計算間違ってる、といった理由が考えられます。前者ならまだ許せるけど、後者はやだな・・・。

共分散行列のときはどうするのか

　どうするんだろうね・・・。

　2019/06/22追記
　手順は増えますが、スケールを考慮すれば同様に行えるようです。

出典：http://manabukano.brilliant-future.net/document/text-PCA.pdf　p.10

loadingを使うと何が良いのか

　相関係数なので、「どれくらい効いてるか」がよくわかります。よくある0.3以下なら～とか0.7以上なら～という論法ができます。それだけといえばそれだけ。

　このように取扱が大変なので、固有ベクトルのまま解釈した方が楽かもという気がしてきました。各主成分の寄与率はexplained_variance_ratio_で得られる訳だし、寄与率の大きい軸の固有ベクトルの大きい次元を見ればどんな感じかはわかるし・・・。

　でもまあ、一応（入力をスケーリングすれば）大体出るということはわかったので、これでよしとします。

　共分散行列でやるときのやり方は、どなたか詳しい方に教えて頂けると助かります。

*1:主成分負荷量、あるいは因子負荷量とも（なぜか）言われますが、この記事ではloadingで通します。けっきょくヘタに和訳しないのがいちばんわかりやすい

【python】sklearnで因子分析を試す

2018-03-31T00:22:11+09:00

　pythonで因子分析をやる人はあまりいないようだが、sklearnにはしっかりモデルが存在している。ついさっき気づいた。

sklearn.decomposition.FactorAnalysis — scikit-learn 0.20.1 documentation

　因子分析自体は前からどんなものなのか興味があり、かといってググるとRだったりSPSSだったりばっかり出てきて辟易していたのだが、sklearnにあると都合が良い。さっそく使ってみよう。

とりあえずirisをプロットする
とりあえずcomponentsを見る
使えることはわかった

とりあえずirisをプロットする

　私だけでも何十回もやってきた、世界中では何万回とやられてきたirisの二次元可視化をやってみる。

　次のようなコードを書いた。

# coding: UTF-8

from copy import deepcopy
from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA, FactorAnalysis as FA
from sklearn.pipeline import Pipeline
import matplotlib.pyplot as plt

def decomp_and_plot(dataset, model, file_name):
    X = model.fit_transform(dataset.data)
    plt.figure()
    plt.scatter(X[:,0], X[:,1], c=dataset.target/len(dataset.target_names))
    plt.savefig(file_name)
    
def main():
    iris = load_iris()

    ss = StandardScaler()
    pca = PCA(n_components=2)
    pl = Pipeline([("scaler", ss), ("pca", deepcopy(pca))])
    fa = FA(n_components=2, max_iter=5000)

    decomp_and_plot(iris, pca, "pca_plt.png")
    decomp_and_plot(iris, pl, "spca_plt.png")
    decomp_and_plot(iris, fa, "fa_plt.png")

if __name__ == "__main__":
    main()

　PCA、変数をスケーリングしたPCA（相関行列を使うことと等価）、因子分析でそれぞれplotしてみる。

　結果はこれ。

PCAの結果

PCA(相関行列)の結果

　相関行列はぱっと見いまいち（この絵一枚でダメかどうかは判断できないが）。

因子分析の結果

　うーん、相関行列のとも違うし、なんとも言い難いというか、素人目にはぶっちゃけあんまり良くないように見えるのだが、確率モデルなのでノイズの存在を仮定して見るとこうなるということだろう。

とりあえずcomponentsを見る

　次のようなmain2を作り、実行した。

def main2():
    iris = load_iris()

    print(iris.feature_names)
    print("pca")
    pca = PCA(n_components=2)
    pca.fit(iris.data)
    print(pca.components_)

    print("fa")
    fa = FA(n_components=2, max_iter=5000)
    fa.fit(iris.data)
    print(fa.components_)

　結果

['sepal length (cm)', 'sepal width (cm)', 'petal length (cm)', 'petal width (cm)']
pca
[[ 0.36158968 -0.08226889  0.85657211  0.35884393]
 [ 0.65653988  0.72971237 -0.1757674  -0.07470647]]
fa
[[ 0.72577591 -0.17754023  1.75733754  0.73196365]
 [-0.37036948 -0.24060118  0.02793388  0.04121372]]

　プロット結果から予想される通り、両者のcomponentsはよく似通っている。

　これがloadingなのかどうかはぶっちゃけよくわからないのだが（というか1を超えてくる時点でたぶん違うのだろうが）、とりあえずloadingだと思って解釈する。

　第一因子は花弁の長さと幅、がく片の長さに対応しているので花の大きさに対応しているっぽい。花の大きさとがく片の幅はなぜか若干反比例する。

　第二因子は花弁に関する係数が小さいので、がく片の大きさを表す因子と言って良さそうである。

　こんなところか。

使えることはわかった

　だから何？　って言われると、正直答えに窮しますが・・・とにかく使えます。主成分分析で良いじゃんと言われたら何も言い返せません。
　

【python】pythonで主成分分析のバイプロット

2018-03-28T23:13:05+09:00

　バイプロット（Biplot）という主成分分析（PCA）の結果の可視化方法があります。

　すごく大雑把に言うと、PCAによる写像の前の空間の各特徴（軸）が写像先の空間のどこに向いているかを可視化する方法です。

　具体的には、主成分ベクトル（因子負荷量などを使う場合もあります）と散布図を同じ図にplotします。これらを組み合わせることで、元の空間の性質が二次元（もしかしたら3次元）で手に取るようにわかります*1。

　バイプロットはR言語だと簡単に描けるらしいのですが、我らがpythonには（少なくとも一般的なライブラリには）そんな便利なものはありません。ちょっと困るのですが、幸い英語圏にはちらほらやりかたの情報があります。しかし、それはそれでページごとにやってることが違ったりして、（申し訳ないのですが）微妙に信用できなかったりします。

　で、けっきょく自分で書いてみることにしました。なお、参考にしたのはこの辺です。

方針

　まずsklearnの公式ドキュメントをできるだけ良く読み込みます。

sklearn.decomposition.PCA — scikit-learn 0.22.1 documentation

　PCA.components_が固有ベクトルであり、データをセンタリングしてこれと掛けるとPCAの出力が出てくることは前回の記事で確認しました。

　固有ベクトル行列が主成分*元のデータの特徴という形になっているとして、横に見ると負荷量（みたいなもの。本当は対応する固有値のsqrtを掛け算してやらないといけない）に、縦に見ると元の写像先で表現された特徴の軸になります。

　つまり、その軸をプロットするだけです。

　なお、この辺は微妙に議論があるようです。私もこれがどこまで正しい方法なのかは自信が持てません。

　参考：
色々と考えてみる: 文系のための「主成分分析の可視化」（２）

　だけど今回は、データをセンタリングしてPCAを学習させた上で、各軸に対応するone-hot vectorを渡してtransformしたら確かに上に書いた方法通りで上手く行きました（biplotの線の上に載った）。なので、「これで良いんだろう」と勝手に判断しました。どこまで妥当かはよくわからないんですけど。

実装

　こんな感じで書きました。

# coding: UTF-8

from sklearn.datasets import load_iris, load_wine
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

import matplotlib.pyplot as plt

def biplot(dataset, scale=False, arrow_mul=1, text_mul=1.1):
    if scale:
        ss = StandardScaler()
        X = ss.fit_transform(dataset.data)
    else:
        X = dataset.data

    if hasattr(dataset, "feature_names"):
        feature_names = list(dataset.feature_names)
    else:
        feature_names = ["F{0}".format(i)
                         for i in range(dataset.data.shape[1])]

    pca = PCA(n_components=2)
    X = pca.fit_transform(X)

    x_data = X[:,0]
    y_data = X[:,1]

    pc0 = pca.components_[0]
    pc1 = pca.components_[1]

    plt.figure()
    plt.scatter(x_data, y_data,
                c=dataset.target/len(set(dataset.target)),
                marker=".")

    for i in range(pc0.shape[0]):
        plt.arrow(0, 0, 
                  pc0[i]*arrow_mul, pc1[i]*arrow_mul,
                  color='r')
        plt.text(pc0[i]*arrow_mul*text_mul,
                 pc1[i]*arrow_mul*text_mul,
                 feature_names[i],
                 color='r')
    plt.show()

def main():
    iris = load_iris()
    wine = load_wine()

    biplot(iris, arrow_mul=2.5, scale=True)
    biplot(wine, arrow_mul=6, scale=True)

if __name__ == "__main__":
    main()

　今回はsklearnのデータセットを渡す形で関数にまとめました。ま、もしこのコードを流用したい人がいたら、必要なロジックだけ上手く切り出してください。

　結果は、こんな画像が出ます。

irisのバイプロット

wineのバイプロット

　上手く行ってる感じです。

　なお、上のコードでは変数をスケーリングしています（相関行列でPCAするのと等価）。スケーリングしなくてもできますが、やった方が矢印の長さが揃いやすいです（逆に変数のスケールを重視してPCAしたいときは、スケーリングしてはいけない。ケースバイケース）。

まとめ

　これくらい自作しなくても済めば良いのにと思いました。

*1:本当に手に取るようにわかるかはデータと見る人に依存しますが・・・

【python】numpyで主成分分析を実装してみた

2018-03-28T22:21:01+09:00

　numpyでPCA（principal component analysis：主成分分析）を実装してみました。自分の理解を深めるためです。

　sklearnに実装されているものと同じ結果を出すことを目標にしました。最終的には上手く行きました。

概要
実装
結果
まとめ

概要

　主成分分析のアルゴリズムの解説は他に譲ります。これは実装してみた記事です。

　実装のやり方は色々あるようですが、一番基本的な（だと思う）共分散行列の固有値と固有ベクトルを求める方法で行きます。

　やるべきこととしては、

データをセンタリングする（列ごとに平均を引く）
共分散行列を計算する
固有値と固有ベクトルを計算
データを固有ベクトルを使って写像する

　これらを実装すれば行けるはずです。というか、これで行くことにしました。

実装

　書いたソースコードを以下に示します。

# coding: UTF-8

import numpy as np

from sklearn.datasets import load_iris
from sklearn.decomposition import PCA

import matplotlib.pyplot as plt

class MyPCA:
    def __init__(self, n_components=2):
        self.n_components = n_components

    def fit_transform(self, X):
        """横着してfit_transformしか実装してない
        """

        # 平均を0にする
        X = X - X.mean(axis=0)

        # 共分散行列を作る
        self.cov_ = np.cov(X, rowvar=False)
        
        # 固有値と固有ベクトルを求めて固有値の大きい順にソート
        l, v = np.linalg.eig(self.cov_)
        l_index = np.argsort(l)[::-1]
        self.l_ = l[l_index]
        self.v_ = v[:,l_index] # 列ベクトルなのに注意

        # components_（固有ベクトル行列を途中まで取り出す）を作る
        self.components_ = self.v_[:,:self.n_components].T

        # データとcomponents_をかける
        # 上と下で二回転置してるのアホ・・・
        T = (np.mat(X)*(np.mat(self.components_.T))).A

        # 出力
        return T

def main():
    iris = load_iris()

    pca = PCA(n_components=2)
    sklearn_X = pca.fit_transform(iris.data)

    my_pca = MyPCA()
    my_X = my_pca.fit_transform(iris.data)

    print(pca.explained_variance_)
    print(my_pca.l_)

    print(pca.components_)
    print(my_pca.components_)

    plt.figure()
    plt.scatter(sklearn_X[:,0], sklearn_X[:,1], c=iris.target/3)
    plt.savefig("sklearn_resut.png")

    plt.figure()
    plt.scatter(my_X[:,0], my_X[:,1]*-1, c=iris.target/3)
    plt.savefig("my_result.png")

if __name__ == "__main__":
    main()

　numpyを使ったので簡単に書けました。アルゴリズム部分はコメントで解説を入れたので、それを読めばどんな感じかは理解して頂けると思います。

結果

　mainのテキスト出力を見ると、次のようになっていました。

# 固有値
[4.22484077 0.24224357]
[4.22484077 0.24224357 0.07852391 0.02368303]

# components_
[[ 0.36158968 -0.08226889  0.85657211  0.35884393]
 [ 0.65653988  0.72971237 -0.1757674  -0.07470647]]
[[ 0.36158968 -0.08226889  0.85657211  0.35884393]
 [-0.65653988 -0.72971237  0.1757674   0.07470647]]

　固有値が余計に出ちゃってますが、これは別に構いません。また、componentsの2次元目が符号反転していますが、これも特に問題ないこと（のはず）なので無視します。

　自作の方は第二主成分を反転させてプロットしてみました。

sklearnのPCAでirisを可視化

自作PCAでirisを可視化

　同じ図を2つ載せるなって怒られそうですが・・・とにかく上手く行ったようです。

まとめ

　numpyで実装してみたら思ったより簡単だったので、これで当分は「わかった気」になれそうです。

　ただ、今回は特異値分解やらなかったので、それはまた宿題ということで・・・。

【python】カーネル主成分分析を試してみる

2018-03-28T00:19:12+09:00

　カーネル主成分分析（Kernel PCA）はカーネル法と主成分分析を組み合わせて用い、データを非線形次元圧縮する方法です（こんな説明で良いのか・・・）。

　カーネル法のことは勉強中・・・というか正直勉強しようとしてもよくわからないで跳ね返されるのをこれまで4回くらい繰り返してきたのですが、とりあえず使ってみました。

試してみた

　非線形データが手元にあると良いのですが、あいにくありません。輪っか状のデータなどを生成してやってみるのは簡単にできますが、面白くなさそうです。だいたいsklearnの公式サンプルにすらあります。
Kernel PCA — scikit-learn 0.21.2 documentation

　そこで、分類問題での適用を考えます。これならいつものようにPCA+CLFとKPCA+CLFで比較するだけなので、簡単そうです。更に、カーネルのgammaはグリッドサーチして最適値を探すだけ・・・。

　ただし、irisやdigitsで散々色々試してみましたが、ぶっちゃけ普通にやるとなかなかPCAを上回る性能が得られませんでした。最終的に、「digitsを3次元に次元削減し、LDAで分類する」という問題でどうにかそれなりに性能が上回ることがわかりましたが、実用的な意味はあまりありません。

　たぶん、sklearnのtoy datasetは低次元で線形分離できるタチの良いデータばっかりなのだと思います。それはそれで良いことですが、ちょっとタチの悪いデータも混ぜておいてもらえると嬉しいところです（かといって20newsgroupsのBoWだとタチが悪すぎるし・・・2000データ400次元くらいのちょうど良いデータはどこかにないものか）。

　コードを以下に示します。

# coding: UTF-8

import numpy as np

from sklearn.datasets import load_digits
from sklearn.decomposition import PCA, KernelPCA as KPCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn.pipeline import Pipeline
from sklearn.model_selection import GridSearchCV, StratifiedKFold as SKF
from sklearn.model_selection import cross_val_score

def main():
    dataset = load_digits()
    print(dataset.data.shape)

    pca = PCA(n_components=3)
    kpca = KPCA(kernel="rbf", n_components=3)
    lda = LDA()
    pl_pca = Pipeline([("pca", pca), ("lda", lda)])
    pl_kpca = Pipeline([("kpca", kpca), ("lda", lda)])

    parameters = {"kpca__gamma" : np.arange(0.00001, 0.003, 0.0001)}

    clf = GridSearchCV(pl_kpca, parameters, verbose=0, n_jobs=-1)

    print(cross_val_score(pl_pca, dataset.data, dataset.target, 
                          cv=SKF(shuffle=True, random_state=0),
                          scoring="f1_macro").mean())                      
    print(cross_val_score(clf, dataset.data, dataset.target,
                          cv=SKF(shuffle=True, random_state=0),
                          scoring="f1_macro").mean())    

if __name__ == "__main__":
    main()

　PCAでは0.68らい、KPCAでは0.71くらいのF1値が得られました。

　だから？　って言われると、返す言葉は思いつきませんが・・・。

まとめ

　やってみた記事ですが、何かの参考になればと思います。意外と上手く使うのは難しいと感じました。というか分類の次元削減としてはたぶんそんなに適当ではないです。

　どんな問題に応用されてるんだろうか。やっぱり可視化？

追記

　文字列の編集距離の可視化に使ってみました。

www.haya-programming.com

　文字列カーネルというのもあるらしいのですが、sklearnで対応していないし、未確認。編集距離を使う分には無難に使えます。

【python】sklearnのPCAで相関行列を使う

2018-03-27T02:41:44+09:00

　主成分分析には共分散行列を用いる方法、相関行列を使う方法がある。

　sklearnのPCAを見ると、これに対応するオプションは存在しない。

sklearn.decomposition.PCA — scikit-learn 0.20.1 documentation

　ずっと不思議に思っていたが、ググってたらこんなものを見つけた。

Enhance: PCA options for using Correlation or covariance matrix · Issue #2689 · scikit-learn/scikit-learn · GitHub

　要約：特徴量をスケーリングしてPCAすれば相関行列でやったのと同じことになるよ。PipelineでStandardScalerと組み合わせてね。おわり。

本当か確認する

　確認してみる。

>>> import numpy as np
>>> from sklearn.datasets import load_iris
>>> from sklearn.preprocessing import StandardScaler
>>> from sklearn.decomposition import PCA
>>> from sklearn.pipeline import Pipeline
>>> iris = load_iris()
>>> pca = PCA(n_components=2)
>>> pca.fit(iris.data)
PCA(copy=True, iterated_power='auto', n_components=2, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)
>>> pca.get_covariance()
array([[ 0.67919741, -0.03258618,  1.27066452,  0.5321852 ],
       [-0.03258618,  0.18113034, -0.31863564, -0.13363564],
       [ 1.27066452, -0.31863564,  3.11934547,  1.28541527],
       [ 0.5321852 , -0.13363564,  1.28541527,  0.58961806]])
>>> ss = StandardScaler()
>>> p = Pipeline([("scaler", ss), ("pca", pca)])
>>> p.fit(iris.data)
Pipeline(memory=None,
     steps=[('scaler', StandardScaler(copy=True, with_mean=True, with_std=True)), ('pca', PCA(copy=True, iterated_power='auto', n_components=2, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False))])
>>> p.steps[1][1].get_covariance()
array([[ 0.9779242 , -0.10104477,  0.87069468,  0.86134879],
       [-0.10104477,  1.00395722, -0.41916911, -0.37286994],
       [ 0.87069468, -0.41916911,  1.04639367,  0.93676197],
       [ 0.86134879, -0.37286994,  0.93676197,  0.99857055]])
>>> np.corrcoef(iris.data, rowvar=False)
array([[ 1.        , -0.10936925,  0.87175416,  0.81795363],
       [-0.10936925,  1.        , -0.4205161 , -0.35654409],
       [ 0.87175416, -0.4205161 ,  1.        ,  0.9627571 ],
       [ 0.81795363, -0.35654409,  0.9627571 ,  1.        ]])

　違うじゃん。妥当そうなのはnumpyの結果だが（対角成分が1になってる）、とりあえずしょうがないのでスケーリングしたデータの共分散をnumpyで計算してみる。

>>> np.cov(ss.fit_transform(iris.data), rowvar=0, bias=1)
array([[ 1.00671141, -0.11010327,  0.87760486,  0.82344326],
       [-0.11010327,  1.00671141, -0.42333835, -0.358937  ],
       [ 0.87760486, -0.42333835,  1.00671141,  0.96921855],
       [ 0.82344326, -0.358937  ,  0.96921855,  1.00671141]])
>>> np.cov(ss.fit_transform(iris.data), rowvar=0, bias=1)
array([[ 1.        , -0.10936925,  0.87175416,  0.81795363],
       [-0.10936925,  1.        , -0.4205161 , -0.35654409],
       [ 0.87175416, -0.4205161 ,  1.        ,  0.9627571 ],
       [ 0.81795363, -0.35654409,  0.9627571 ,  1.        ]])

　標本分散はnp.corrcoefと等価だ。

　ここまでやったところでもう一回ドキュメントを読み、PCA.get_covariance()の結果が「Estimated covariance of data.」であり、厳密ではないことに気づいたので、問題は解決した。

　理論的にこうなる理由は、説明しようと思えばできるのだと思いますが、今回は大変なので触れません。

irisでやってみる

　irisの可視化にそれぞれを使ってみる。コードを以下に示す。

# coding: UTF-8

from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.pipeline import Pipeline

import matplotlib.pyplot as plt

def main():
    iris = load_iris()

    ss = StandardScaler()
    pca = PCA(n_components=2)
    p = Pipeline([("scaler", ss), ("pca", pca)])
    
    X = pca.fit_transform(iris.data)

    plt.figure()
    plt.scatter(X[:,0], X[:,1], c=iris.target/3)
    plt.savefig("iris_cov_pca.png")

    X = p.fit_transform(iris.data)

    plt.figure()
    plt.scatter(X[:,0], X[:,1], c=iris.target/3)
    plt.savefig("iris_corr_pca.png")

if __name__ == "__main__":
    main()

　結果は、

共分散行列で主成分分析したiris

相関行列で主成分分析したiris

　こうして見ると相関行列はあまりメリットがないように見えますが、実際には相関行列の方が良いタスクは色々あるようです。相関行列を使うことでbiplotが上手く行っているという例を出しているページを載せておきます。
PCA on correlation or covariance? - Cross Validated

まとめ

　とりあえずできることはわかったので良しとする。

　でも、「pipelineで出来るから要らねーよ」ってつもりらしいけど、ぶっちゃけオプション一つでできた方が親切だと思った（小並感）。

【python】sklearnのfetch_20newsgroupsで文書分類を試す(4)

2018-03-26T21:21:12+09:00

　前回は性能を追い求めると次元がでかくなりすぎて・・・というところで終わっていた。今回はもうちょっと頑張って次元を減らしてみる。

ストップワードの除去
PCA（主成分分析）とLDA（線形判別分析）
分類
ソースコード
結果とまとめ
次回
過去の回

ストップワードの除去

　とりあえずstop_wordsを指定していなかったので、指定してみる。

　stop_words="english"とすると、ストップワードを除去してくれる。

　結果だけ言うと、min_df=0.005のとき、

stop_words指定なし：3949次元
stop_words指定あり：3705次元

　だった。焼石に水。

PCA（主成分分析）とLDA（線形判別分析）

　PCAとLDAをかけ、次元削減をする。leakage怖いのでPipelineを使う（厳密なことを言い出すと、単語文書行列を作る段からPipelineに入れるべきなのだろうか？　きついのでパスさせて頂くが）。

　PCAは主にLDAの計算負荷削減と、変数の相関を除去することを意図してかける。1000次元まで落としてみたが、これでも累積寄与率は90%弱になる。まあ、正規化も何もしてないから、重要な情報を落としている可能性は否定できないのだが。

　LDAは次元削減に使う。有効性についてはこの前試してみたので、この記事を読んで欲しい。
【python】LDA（線形判別分析）で次元削減 - 静かなる名辞
　20newsgroupsは20クラスのデータなので、19次元に落とすことになる。相当早くなるだろうが、どこまで性能を維持できるかはデータの線形性にかかっている。

分類

　ランダムフォレストを使った。n_estimators=1000とし、他のパラメタはデフォルト。

ソースコード

　実験に使ったソースコードを以下に示す。

# coding: UTF-8

import numpy as np

from sklearn.datasets import fetch_20newsgroups
from sklearn.ensemble import RandomForestClassifier as RFC
from sklearn.feature_extraction.text import CountVectorizer as CV
from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import precision_recall_fscore_support as prf
from sklearn.pipeline import Pipeline

def main():
    news20 = fetch_20newsgroups()    
    
    cv = CV(min_df=0.005, max_df=0.5, stop_words="english")
    matrix = cv.fit_transform(news20.data).toarray()

    pca = PCA(n_components=1000, svd_solver="randomized")
    lda = LDA()
    rfc = RFC(n_estimators=1000, n_jobs=-1)

    clf = Pipeline([("pca", pca), ("lda", lda), ("rfc", rfc)])

    trues = []
    preds = []
    for train_index, test_index in StratifiedKFold().split(matrix, news20.target):
        clf.fit(matrix[train_index], news20.target[train_index])
        trues.append(news20.target[test_index])
        preds.append(clf.predict(matrix[test_index]))
    scores = prf(np.hstack(trues), np.hstack(preds), average="macro")[:3]
    print("p:{0:.6f} r:{1:.6f} f1:{2:.6f}".format(scores[0],
                                                  scores[1],
                                                  scores[2]))

if __name__ == "__main__":
    main()

結果とまとめ

p:0.764012 r:0.760731 f1:0.761510

　前回の0.8を超えるスコアには届かなかったが、とりあえずそれなりに軽くはなった。もうちょっと真面目に追い込めばという話はあるが、追求しない。次回はもうちょっと違うことをやってみたい。

次回

　このシリーズずっと放置していましたが、気が向いたので書きました。
www.haya-programming.com

過去の回

【python】sklearnのfetch_20newsgroupsで文書分類を試す(1) - 静かなる名辞
 【python】sklearnのfetch_20newsgroupsで文書分類を試す(2) - 静かなる名辞
 【python】sklearnのfetch_20newsgroupsで文書分類を試す(3) - 静かなる名辞

【python】matplotlibで3次元データを描画し、回転アニメーションにする

2018-03-24T22:16:53+09:00

　3次元くらいのデータを描画したいときがある。簡単に散布図にできると便利。

データの用意

　sklearnのload_irisなどで取得できるデータセットを入力にする前提の次のような関数を作った。

from sklearn.decomposition import PCA

def gen_3d_data(dataset):
    pca = PCA(n_components=3)
    return pca.fit_transform(dataset.data), dataset.target

　あとはirisなり何なりを入れてやる。

3次元プロット

　とりあえずプロットしたいだけならこれだけ。

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

def matplotlib_plt(X, y, filename):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection="3d")

    ax.scatter(X[:,0], X[:,1], X[:,2], c=y/len(set(y)))
    plt.savefig(filename)
    plt.show()

　点の色の指定がちょっとセコいが、まあ良いこととする。axes3dはパっと見使ってないように見えるが、importしないとエラーになるので必要と思われる。

　呼び出し元のmainも書く。

from sklearn.datasets import load_iris

def main():
    iris_X, iris_y = gen_3d_data(load_iris())
    matplotlib_plt(iris_X, iris_y, "iris.png")

　実行すると次のような画像が出力される。

irisの3次元描画

　あと、ぐりぐり回せるグラフのようなものが別ウィンドウで開く（plt.sow()に対応）。

回転させたアニメーションを表示

　このような例が公式で示されている。

for angle in range(0, 360):
    ax.view_init(30, angle)
    plt.draw()
    plt.pause(.001)

mplot3d example code: rotate_axes3d_demo.py — Matplotlib 2.0.2 documentation

　やると確かにぐるぐる回るアニメーションが表示される。こりゃあええわ、ということで、次はこれをgifアニメにすることを考える。

　matplotlibにもanimationというモジュールがあり、色々できるようだが使い方を覚えるのが大変そうだった。なので、「上のforループ内で一枚ずつ画像出力してffmpegで繋げば良いだろ」という手抜きの方針で行くことにする。

def matplotlib_rotate(X, y, dataname):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection="3d")

    ax.scatter(X[:,0], X[:,1], X[:,2], c=y/len(set(y)))

    for angle in range(0, 360):
        ax.view_init(30, angle)
        plt.savefig("figs/{0}_{1:03d}.jpg".format(dataname, angle))

　呼び方は、

matplotlib_rotate(iris_X, iris_y, "iris")

　こうするとfigs/以下に画像が360枚吐かれるので、ffmpegでつなぐ。

$  ffmpeg -r 10 -i figs/iris_%03d.jpg -pix_fmt rgb24 -f gif out.gif

　とりあえずこれで行けた。画質が悪い割に容量が重いので、どこか上手くない指定になってるのかもしれないけど。

回るiris

　上出来ではないだろうか。

【python】LDA（線形判別分析）で次元削減

2018-03-20T16:43:52+09:00

　一般によく使われる次元削減手法としてはPCA（主成分分析）がありますが、他にLDA（Linear Discriminant Analysis：線形判別分析）を使う方法もあります。

　これは本来は分類に使われる判別分析という古典的なアルゴリズムで、データが一番分離しやすくなる軸を求めていくものです。つまり教師ラベルを使います。教師ラベルを使うので、PCAのような教師なしの手法と比べて有利な可能性があります。

　線形判別分析の詳しい原理の説明などが欲しい方は、ググって出てくるwikipediaやqiitaなどを参考にしてください（投げやり）。この記事では、分類問題でこれを使ったとき、どのようなご利益があるのかを検証します。

実験

　sklearnのdigitsデータセットを使い、次元削減→分類というタスクを行って交差検証でスコアを出します。

　分類器は最初はSVMでやろうかと思ったけど、パラメタチューニングで幾らでも恣意的な結果になることに気づいたのでガウシアン・ナイーブベイズでやることにしました。

　実験に使ったコードは以下に示します。

# coding: UTF-8

import warnings
warnings.filterwarnings('ignore')

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn.naive_bayes import GaussianNB as GNB
from sklearn.pipeline import Pipeline
from sklearn.model_selection import StratifiedKFold as SKF
from sklearn.metrics import precision_recall_fscore_support  as prf

def main():
    digits = load_digits()

    gnb = GNB()

    df = pd.DataFrame([], columns=[
        "n_components",
        "pca-gnn precision", "pca-gnn recall", "pca-gnn f1",
        "lda-gnn precision", "lda-gnn recall", "lda-gnn f1"])
    for n_components in [5, 10, 15, 20, 25, 30, 40]:
        pca = PCA(n_components=n_components)
        lda = LDA(n_components=n_components)

        steps1 = list(zip(["pca", "gnb"], [pca, gnb]))
        steps2 = list(zip(["lda", "gnb"], [lda, gnb]))

        p1 = Pipeline(steps1)
        p2 = Pipeline(steps2)

        score_lst = []
        for decomp_name, clf in zip(["pca", "lda"], [p1, p2]):
            trues = []
            preds = []
            for train_index, test_index in SKF(
                    shuffle=True, random_state=0).split(
                    digits.data, digits.target):
                clf.fit(digits.data[train_index], 
                        digits.target[train_index])
                trues.append(digits.target[test_index])
                preds.append(clf.predict(digits.data[test_index]))
            scores = prf(np.hstack(trues), np.hstack(preds), average="macro")
            score_lst.extend(scores[:-1])
        df = df.append(pd.Series([n_components, *score_lst],
                                 index=df.columns),
                       ignore_index=True)
    print(df)
    plt.figure()
    df.plot(x="n_components", y=["pca-gnn f1", "lda-gnn f1"])
    plt.savefig("result.png")

if __name__ == "__main__":
    main()

結果

　次のようになりました。

　テキスト出力

   n_components  pca-gnn precision  pca-gnn recall  pca-gnn f1  \
0           5.0           0.847918        0.841684    0.841109   
1          10.0           0.915834        0.911346    0.912563   
2          15.0           0.926992        0.923032    0.924061   
3          20.0           0.934522        0.930192    0.931194   
4          25.0           0.941886        0.938611    0.939205   
5          30.0           0.946139        0.944251    0.944669   
6          40.0           0.945330        0.943644    0.943960   

   lda-gnn precision  lda-gnn recall  lda-gnn f1  
0           0.917464        0.917144    0.917031  
1           0.953751        0.952588    0.952950  
2           0.953751        0.952588    0.952950  
3           0.953751        0.952588    0.952950  
4           0.953751        0.952588    0.952950  
5           0.953751        0.952588    0.952950  
6           0.953751        0.952588    0.952950

結果（n_components対F1値）

　LDAを使った方が低い次元で、より高い分類性能が得られているようです。

まとめ

　LDAは良い。

おまけ

　ソースコードをちゃんと読んだ方は、最初に書かれた以下の記述に気づいたかと思います。

import warnings
warnings.filterwarnings('ignore')

　これを付けないとLDAはけっこうな警告（主に以下の2つ）を吐いてくれます。

UserWarning: Variables are collinear
UserWarning: The priors do not sum to 1. Renormalizing

　上の警告はPCAで説明変数の多重共線性を除去してやると消えます（本末転倒っぽいけど）。下の警告は、正直調べてもよくわかりませんでした。

　とりあえず、警告が出てもちゃんと動いてるみたいなので別に良いか・・・。

追記

　LDAのn_componentsには上限があり、クラス数-1以上のn_componentsは指定しても無意味です。

　実際にやってみても、クラス数-1以上にはなりません。

>>> from sklearn.datasets import load_digits
>>> from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
>>> lda = LDA(n_components=15)
>>> lda.fit(digits.data, digits.target)
>>> lda.explained_variance_ratio_
array([0.28912041, 0.18262788, 0.16962345, 0.1167055 , 0.08301253,
       0.06565685, 0.04310127, 0.0293257 , 0.0208264 ])

　決定境界をクラス数-1個引くので（SVMで言うところのone-versus-the-rest）、n_componentsも必然的にそれだけ必要になります（逆にそれ以上は必要になりません）。

　上のグラフはそのつもりで眺めてください。また、LDAはけっきょくのところ線形変換なので、クラス数-1次元の線形空間にうまく張り直せないような入力に対しては無力なことも覚えておく必要があるでしょう（PCAも非線形構造はダメだが・・・カーネルでも持ってくる必要がある）。

【python】sklearnのPCAでsvd_solverによる速度差を比較

2018-03-19T17:23:15+09:00

　sklearnのPCA（主成分分析）がやたら遅くて腹が立ちました。計算コストを下げるために次元削減してるのに、次元削減で計算コスト食ったら意味がありません。

　とにかくこのPCAを高速化したかったので、svd_solverを変えてどうなるか試しました。なお、腹が立つくらい遅かった理由は最終的にちゃんとわかったので、この記事の最後に載せます。

svd_solverとは
実験
結果
まとめ
おまけ：腹が立った理由

svd_solverとは

　PCAは内部で特異値分解（SVD）を使っています。この特異値分解がコンピュータにやらせるにはそれなりに計算コストの高い処理で、とりあえずアルゴリズムが何種類かあるようです。

　sklearnのPCAで使える（指定できる）アルゴリズムは次の4つです。

auto

　デフォルト値。500*500以下の入力データならfullを、それ以上ならrandomizedを使うそうです*1

full

　standard LAPACK solverを使うそうです。とりあえずぜんぶ丸ごと特異値分解してから、n_componentsで指定した次元数だけ取ってくるそうな

arpack

　Truncate SVDという手法を使う。一次元ずつ寄与率の大きい主成分から計算していくらしい。n_componentsが小さければ速いことが期待されるんだと思う

randomized

　randomized SVDという手法で計算する。乱数使って速くした。乱数なので厳密解ではない

　なお、以上の情報はすべて公式ドキュメントから得ました。
sklearn.decomposition.PCA — scikit-learn 0.20.1 documentation

　とりあえずautoはどうでも良いので、残りの3つを比較することにします。

実験

　PCAをかけたくなるような高次元データといえばBag of Words、ということでこのブログですでに何回も取り上げたことのある、sklearnのfetch_20newsgroupsとCountVectorizerの組み合わせを使います。前者はテキストのデータセット、後者はBoWを生成するクラスです。

　次のような実験用コードを書きました。

# coding: UTF-8

import time
from itertools import product

from sklearn.datasets import fetch_20newsgroups
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.decomposition import PCA

def main():
    news20 = fetch_20newsgroups()

    for min_df in [0.02, 0.01, 0.008, 0.005]:
        cv = CountVectorizer(min_df=min_df, max_df=0.5,
                             stop_words="english")
        X = cv.fit_transform(news20.data).toarray()

        print("min_df:{0} X.shape:{1}".format(min_df, X.shape))
        for n_components, svd_solver in product(
                [100, 500],
                ["full", "arpack", "randomized"]):
            pca = PCA(n_components=n_components, svd_solver=svd_solver)
            t1 = time.time()
            pca.fit_transform(X)
            t2 = time.time()
            print("n_components:{0}  solver:{1:>10}  "\
                  "time:{2:>6.2f}  CP:{3:.4f}".format(
                      n_components, svd_solver, t2-t1, 
                      pca.explained_variance_ratio_.sum()))
        print("")

if __name__ == "__main__":
    main()

　BoWの次元数をmin_dfで変えていき、n_componentsを100と500、svd_solverを上記3つで変化させてPCAをかけたときの速度と累積寄与率（CP：Cumulative Proportion）をそれぞれ測ります。

結果

　次のようになりました。

min_df:0.02 X.shape:(11314, 866)
n_components:100  solver:      full  time:  3.60  CP:0.7455
n_components:100  solver:    arpack  time:  3.90  CP:0.7455
n_components:100  solver:randomized  time:  1.72  CP:0.7443
n_components:500  solver:      full  time:  3.89  CP:0.9528
n_components:500  solver:    arpack  time: 19.42  CP:0.9528
n_components:500  solver:randomized  time:  8.91  CP:0.9516

min_df:0.01 X.shape:(11314, 1916)
n_components:100  solver:      full  time: 22.38  CP:0.8029
n_components:100  solver:    arpack  time:  8.41  CP:0.8029
n_components:100  solver:randomized  time:  4.86  CP:0.8028
n_components:500  solver:      full  time: 22.06  CP:0.9304
n_components:500  solver:    arpack  time: 53.73  CP:0.9304
n_components:500  solver:randomized  time: 13.47  CP:0.9293

min_df:0.008 X.shape:(11314, 2391)
n_components:100  solver:      full  time: 34.24  CP:0.7899
n_components:100  solver:    arpack  time: 10.42  CP:0.7899
n_components:100  solver:randomized  time:  5.75  CP:0.7897
n_components:500  solver:      full  time: 34.88  CP:0.9193
n_components:500  solver:    arpack  time: 63.37  CP:0.9193
n_components:500  solver:randomized  time: 15.18  CP:0.9182

min_df:0.005 X.shape:(11314, 3705)
n_components:100  solver:      full  time:100.52  CP:0.7701
n_components:100  solver:    arpack  time: 16.46  CP:0.7701
n_components:100  solver:randomized  time:  8.70  CP:0.7699
n_components:500  solver:      full  time:100.73  CP:0.9000
n_components:500  solver:    arpack  time: 94.33  CP:0.9000
n_components:500  solver:randomized  time: 20.04  CP:0.8988

　要約すると、

fullは基本的に遅い。入力の次元数が増えるとびっくりするくらい遅くなる
arpackは100次元に落とすときは威力を発揮している。500次元に落とすケースではかえって遅くなる。ヘタするとfullより遅い
randomizedは速い。ただし厳密解ではないことがCPからわかる（full、arpackとは微妙に違う数字になっている）

　こういう状況です。わかりやすいですね。

　それぞれの使い分けは、

入力次元数の小さい入力ではfullで良い。というかヘタにそれ以外を指定するとかえって遅いケースもある
入力次元数が大きく、入力次元数>>出力次元数で厳密解がほしければならarpackの使用を検討する
厳密解じゃなくても良いのでとにかく速いのを！　ってときはrandomized

　ってことになるかと思う・・・。

まとめ

　けっこう変わる。頑張って使い分けよう。

おまけ：腹が立った理由

　sklearnのPCAではn_componentsに小数を指定できます。そうすると累積寄与率がその数字になるように勝手に次元数を決めてくれるので、こりゃ便利だわいと思って私はよく使っていました。

　しかし、実はarpack、randomizedではこの小数での指定は使えません。そのことはドキュメントにもちゃんと書いてあります。無理矢理に指定すると次のようなエラーを吐かれます。

ValueError: n_components=0.95 must be between 1 and n_features=866 with svd_solver='arpack'

　ということは何が起こるか？　勝手にfullにされます。遅い訳です。なんてこった。

　わかってしまえば下らない話で、要するに私が使いこなせていなかっただけなのですが、このことは「ちゃんとドキュメントをよく読んで使おうね」という教訓を私に残したのでした。

*1:300*800だったりしたらどうなるんだろう？　それとも共分散行列のサイズなのだろうか？

sklearnのclassification_reportで多クラス分類の結果を簡単に見る

2018-03-12T21:35:24+09:00

はじめに

　多クラス分類をしていると、「どのクラスが上手く分類できてて、どのクラスが上手く行ってないんだろう」と気になることがままあります。

　そういった情報を簡単に要約して出力してくれるのがsklearnのclassification_reportで、簡単に使える割に便利なので実験中や開発中に威力を発揮します。

　※この記事はsklearn 0.19の時代に書きましたが、その後sklearn 0.20で使い方が変更されたので、2019/03/18に全面的に改稿しました。

使い方

　ドキュメントを見るととても簡単そうです。
sklearn.metrics.classification_report — scikit-learn 0.20.3 documentation

sklearn.metrics.classification_report(
    y_true, y_pred, labels=None, target_names=None,
    sample_weight=None, digits=2, output_dict=False)

　要するに真のラベルと予測ラベル、あとラベルに対応する名前を入れてあげればとりあえず使えます。文字列の返り値が出力になります。sample_weight, digitsはそれぞれサンプルの重みと結果に出力される桁数を表しますが、とりあえず入れなくても大した問題は普通はありません。output_dictはsklearn 0.20から追加された引数で、pandasデータフレームに変換可能な辞書を返します。

　さっそく使ってみましょう。

# coding: UTF-8

import numpy as np
from sklearn.datasets import load_iris
from sklearn.svm import SVC
from sklearn.metrics import classification_report
from sklearn.model_selection import StratifiedKFold as SKF

def main():
    # irisでやる
    iris = load_iris()

    # svmで分類してみる
    svm = SVC(C=3, gamma=0.1)

    # 普通の交差検証
    trues = []
    preds = []
    for train_index, test_index in SKF().split(iris.data, iris.target):
        svm.fit(iris.data[train_index], iris.target[train_index])
        trues.append(iris.target[test_index])
        preds.append(svm.predict(iris.data[test_index]))
        
    # 今回の記事の話題はここ
    print("iris")
    print(classification_report(np.hstack(trues), np.hstack(preds), 
                                target_names=iris.target_names))

if __name__ == "__main__":
    main()

　すると、次のような出力が得られます。

iris
              precision    recall  f1-score   support

      setosa       1.00      1.00      1.00        50
  versicolor       0.96      0.98      0.97        50
   virginica       0.98      0.96      0.97        50

   micro avg       0.98      0.98      0.98       150
   macro avg       0.98      0.98      0.98       150
weighted avg       0.98      0.98      0.98       150

　precision, recall, f1-scoreという代表的な評価指標と、support（＝y_trueに含まれるデータ数）が、クラスごとと全体の各種平均（後述）で出る、というのが基本的な仕組みです。

　まずクラスごとの結果を見ると、setasoは100%分類できていますが、versicolorとvirginicaはどうも混ざっているようです。以前の記事でirisを二次元にした画像を作ったので、再掲します。

irisをPCAで二次元にしたもの

　RGBの順でsetaso, versicolor, virginicaに対応しているはずです。ということはsetasoが綺麗に分離できてversicolorとvirginicaが混ざるというのは極めて妥当な結果ということになりそうです。

　また、下にあるmicro avg, macro avg, weighted avgは、それぞれマイクロ平均、マクロ平均、サンプル数で重み付けられた平均です。

　出る評価指標などの詳細については別途記事を書いたので、そちらを御覧ください。

www.haya-programming.com

output_dictを使って便利に集計する

　sklearn 0.20ではoutput_dictという引数がこの関数に追加されました。これを使うとデフォルトの文字列ではなく辞書形式で結果を得ることができ、結果をプログラム上で取り扱うことが容易になります。

　上のコードの出力部分を2行書き換えます。

    from pprint import pprint
    pprint(classification_report(np.hstack(trues), np.hstack(preds), 
                                 target_names=iris.target_names,
                                 output_dict=True))

　結果はこのようになります。

{'macro avg': {'f1-score': 0.97999799979998,
               'precision': 0.9801253834867282,
               'recall': 0.98,
               'support': 150},
 'micro avg': {'f1-score': 0.98,
               'precision': 0.98,
               'recall': 0.98,
               'support': 150},
 'setosa': {'f1-score': 1.0, 'precision': 1.0, 'recall': 1.0, 'support': 50},
 'versicolor': {'f1-score': 0.9702970297029702,
                'precision': 0.9607843137254902,
                'recall': 0.98,
                'support': 50},
 'virginica': {'f1-score': 0.9696969696969697,
               'precision': 0.9795918367346939,
               'recall': 0.96,
               'support': 50},
 'weighted avg': {'f1-score': 0.9799979997999799,
                  'precision': 0.980125383486728,
                  'recall': 0.98,
                  'support': 150}}

　この辞書の形式はpandasデータフレームに変換することも可能です。

    import pandas as pd
    d = classification_report(np.hstack(trues), np.hstack(preds), 
                              target_names=iris.target_names,
                              output_dict=True)
    df = pd.DataFrame(d)
    print(df)

　とすると、

            macro avg  micro avg  setosa  versicolor  virginica  weighted avg
f1-score     0.979998       0.98     1.0    0.970297   0.969697      0.979998
precision    0.980125       0.98     1.0    0.960784   0.979592      0.980125
recall       0.980000       0.98     1.0    0.980000   0.960000      0.980000
support    150.000000     150.00    50.0   50.000000  50.000000    150.000000

　のようにデータフレームとして見ることができます。ここからCSV, TeX, HTML, グラフなど任意のフォーマットに変換できるので、なにかと捗ると思います。

sklearn 0.20での変更点のまとめ

　別途記事を書いたので、そちらを御覧ください。

www.haya-programming.com

classification_reportを使わないとしたら

　このように大変便利なのですが、参考のためにこれを使わない方法も紹介しておきます。sklearn.metrics.precision_recall_fscore_supportを使います。

sklearn.metrics.precision_recall_fscore_support — scikit-learn 0.20.3 documentation

　使い方はこんな感じです。　

from sklearn.metrics import precision_recall_fscore_support
precision_recall_fscore_support(y_true, y_pred, average=None)

　結果はこんな感じになります（上のプログラムを対象に計算し、返り値をpprintしました）。

(array([1.        , 0.96078431, 0.97959184]),
 array([1.  , 0.98, 0.96]),
 array([1.        , 0.97029703, 0.96969697]),
 array([50, 50, 50]))

　numpy配列を格納したタプルが返ってますね。それぞれのnumpy配列がprecision, recall, fscore, supportに対応します。

まとめ

　簡単に使えるので、分類結果を見てみたいときはとりあえずこれに放り込むと良いかと思います。また、sklearn 0.20からはかなり便利になったので、汎用的な分類結果集計方法としても使えるようになりました。

【python】RandomForestの木の本数を増やすとどうなるか？

2018-03-08T20:14:00+09:00

はじめに

　RandomForest（ランダムフォレスト）には木の本数という重要なパラメータがある。slearnのデフォルトは10だが、実際に使うときは1000以上にしてやらないと良い性能が得られないということをよく経験する。

　これを大きくすることで、一体どんな効果が得られるのだろうか？

予想1：より複雑な形状の分離超平面を学習できるようになる
予想2：汎化性能が向上する

　予想1の効果は恐らく木の本数が相対的に少ないとき（100本以下）に顕著に現れると考えられる。その後、木の本数が増えていくに従ってモデルのバリアンスが下がり、予想2の通り汎化性能は向上する方向に向かうと考えられる。

　ここで思い浮かぶ疑問は、「とにかく木の本数を増やしさえすれば、SVMみたいに高い汎化性能が得られるのか？」という点である。RandomForestは決定木なので、基本的にデータの次元軸に直交する決定境界しか引けないという弱点がある。それでも、とにかく木を増やしていけば、丸くてぬるぬるした決定境界になったりするのだろうか？

　実際にどうなるのか、やってみよう。

実験：2次元データで木の本数を変えながら予測確率を評価する

　コードは流し読みしてください。結果の画像だけ見ればわかります。

# coding: UTF-8

import numpy as np

from sklearn.ensemble import RandomForestClassifier as RFC
from sklearn.svm import SVC

from matplotlib import pyplot as plt

def make_circle(a=1, b=1, xy=(0,0), phi=0, n=100, random_s=0.1):
    theta = np.arange(0, 2*np.pi, 2*np.pi/n)
    X = a*np.cos(theta)
    Y = b*np.sin(theta)
    data_mat = np.matrix(np.vstack([X, Y]))
    phi_d = np.deg2rad(phi)
    rot = np.matrix([[np.cos(phi_d), -np.sin(phi_d)],
                     [np.sin(phi_d), np.cos(phi_d)]])
    rot_data = rot*data_mat
    X = rot_data[0].A
    Y = rot_data[1].A

    rand1 = np.random.normal(scale=random_s, size=theta.shape)
    rand2 = np.random.normal(scale=random_s, size=theta.shape)

    return X+rand1+xy[0], Y+rand2+xy[1]

def gen_data():
    n = 150
    X1, Y1 = make_circle(a=6.5, b=4.5, n=n, random_s=0.4)
    X2, Y2 = make_circle(a=5, b=3, n=n, random_s=0.4)
    X3, Y3 = make_circle(a=2.5, b=1.5, n=n, random_s=0.4)

    X = np.hstack([X1, X2, X3])
    Y = np.hstack([Y1, Y2, Y3])
    data = np.vstack([X, Y]).T
    target = np.array([0]*n + [1]*n + [0]*n)

    return data, target

def plot_proba(data, target, clf, filename="fig.png", text=""):
    plt.figure()
    ax = plt.subplot()
    ax.set_xlim((-9, 9))
    ax.set_ylim((-6, 6))

    x = np.arange(-9, 9, 18/250)
    y = np.arange(-6, 6, 12/250)

    X, Y = np.meshgrid(x, y)
    points = np.c_[X.ravel(), Y.ravel()]

    probs = clf.predict_proba(points)
    probs = probs[:,1].reshape(X.shape)

    plt.pcolormesh(X, Y, probs, cmap=plt.cm.RdBu)
    plt.scatter(data[:,0], data[:,1], c=["rb"[t] for t in target], edgecolors='k')
    plt.title(text)

    plt.savefig(filename)

def main():
    data, target = gen_data()

    for i, n_estimators in enumerate([5, 10, 30, 50, 100, 300, 500, 1000]):
        rfc = RFC(n_estimators=n_estimators, oob_score=True, n_jobs=-1)
        rfc.fit(data, target)
        print("{0}:{1}".format(n_estimators, rfc.oob_score_))
        plot_proba(data, target, rfc,
                   "rfc_{0}.png".format(n_estimators),
                   "RFC n_estimators={0}".format(n_estimators))

    svm = SVC(C=20, gamma=0.05, probability=True)
    svm.fit(data, target)
    plot_proba(data, target, svm, "svm.png", "SVM C=20, gamma=0.05")

if __name__ == "__main__":
    main()

　昨日の記事の楕円形データ生成を使っています。
hayataka2049.hatenablog.jp

結果

　そこそこ良くなるが、100以上では改善度合いは微妙かもしれない。OOB errorは、

5:0.8288888888888889
10:0.8822222222222222
30:0.9044444444444445
50:0.9133333333333333
80:0.9155555555555556
100:0.9222222222222223
300:0.92
500:0.9177777777777778
1000:0.9222222222222223

　やはり100以上の改善は微妙という、画像を見て思う感覚を裏付けるものになっている。

　では、SVMだとどんな画像が得られるだろうか？

　これは勝てない。RandomForestだとどうしてもカクカクが残るのに。

考察

　この結果の妥当性は率直に言って判断しづらい。

　そもそも、2次元データを入力している以上、ランダムフォレストはデフォルトで（）の特徴量を使って木を作ってくれている訳で、つまり1次元だけで判断してくれている。ちょっとあんまりなので、max_features=2としたのが次の画像。データが（ランダム絡みで）変わってるのでそこだけ注意。

　まあ、これを見てもSVMみたいに滑らかな決定境界が引けてるとは言い難いものがあるけど・・・（考えようによっては1次元でやった上の画像の方が汎化性能は高い、ような気もしてこなくはない。全体的にもやっとしてて、相対的に滑らかに見える）。

　でも、SVMもSVMで、パラメタ次第ではどんな複雑怪奇な決定境界だって引けるといえば引ける。

　こういう問題（けっきょく汎化性能が得られるかどうかはパラメタ次第）があるので、SVMの方が良いと一概に言えるかはけっこう微妙。

　もっと言えば、「それぞれの軸に意味がある」「ノイズもけっこう混ざってる」「スパース」「高次元」という性質のデータを対象とする場合、SVMの汎化性能（滑らかな決定境界を引ける）はかえって邪魔になるのではないだろうか？　そのようなデータでは、とにかく重要な軸を見つけてきて、そこで判断するRandomForestの方が良い性能が得られることが多いと経験的に感じる。ちなみに自然言語処理で使うBoW(Bag of Words)はその典型例である。

　逆に、軸に意味がなくて相対的に低次元でデンスな空間を相手にする場合、SVMの方が良い結果を産むということもよく経験することである。PCAで低次元に落としてしまったデータとか、word2vecで生成される単語の空間とかが割とそんな感じである。

　なんだか話が脱線してきたのでこれくらいにするけど、けっきょく「滑らかな決定境界を引く能力はどうやってもSVMの方が高い（あたりまえ）」「滑らかだから良いというものでもない」「使い分けが重要」という当たり障りのない結論に落ち着いてしまった。

　あと、木の本数は無尽蔵に増やすわけにはいかない。ランダムフォレストは計算量は軽いけど意外とリソース消費の激しいアルゴリズムで、増やしすぎると効率が悪化する。

ランダムフォレストはサンプル数が多いとメモリ消費量が大きい - 静かなる名辞

　汎化性能はできるだけ他の方法で確保したいところ。

まとめ

　SVMの方が滑らかでした（小並感）。

【python】混合ガウスモデル (GMM)でハード・ソフトクラスタリング

2018-03-06T04:39:50+09:00

はじめに

　先日はFuzzy c-meansによるソフトクラスタリングを行いました。

【python】skfuzzyのFuzzy c-meansでソフトクラスタリング - 静かなる名辞

　ソフトクラスタリングの有名な手法としてはFuzzy c-meansの他に、混合ガウスモデル（混合正規分布モデル）を使った手法があります。この手法はデータが「複数の正規分布から構成されている」と仮定し、その正規分布のパラメタ*1をEMアルゴリズム（expectation–maximization algorithm）という手法を使って最尤推定します。

　ごちゃごちゃと書きましたが、要するに「3つのクラスタにクラスタリングしたければ、（各クラスタのデータの分布が正規分布に従うと仮定して）3つの正規分布が重なりあってると思ってGMMを使って解く」という乱暴なお話です。正規分布が重なりあっているとみなすということは、どの分布に属するかも確率でわかる訳で、これがソフトクラスタリングに使える理由です。ハードクラスタリングに使いたいときは、確率最大のクラスタラベルに振ることになるかと思います。

　このGMM、pythonではsklearnに入っているので簡単に使えます。

sklearn.mixture.GaussianMixture — scikit-learn 0.20.1 documentation

　ということで、他のクラスタリング手法と比較してみることにしました。

実験の説明

　先日の記事でやったのと同様、irisをPCAで二次元に落としたデータに対してクラスタリングを行います。クラスタリング結果（所属するクラスタの確率）はirisが3クラスのデータなのを利用し、色（RGB）で表現します。

　比較するクラスタリング手法はk-means（ハード）、Fuzzy c-means（ソフト）、GMM（ハード・ソフト）です。

　前回はFuzzy c-meansのパラメタmを動かして結果を見たりしましたが、今回これは2で決め打ちにします。

　実験用ソースコードは次のものです。走らせるにはいつもの定番ライブラリ以外にscikit-fuzzyというライブラリを入れる必要があります（あるいはFuzzy c-means関連の部分をコメントアウトするか。でもskfuzzyはpipで一発で入るし、入れておいても別に損はない）。

# coding: UTF-8

import numpy as np

from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans as KM
from sklearn.mixture import GaussianMixture as GMM
from matplotlib import pyplot as plt

from skfuzzy.cluster import cmeans

def target_to_color(target):
    if type(target) == np.ndarray:
        return (target[0], target[1], target[2])
    else:
        return "rgb"[target]

def plot_data(data, target, filename="fig.png"):
    plt.figure()
    plt.scatter(data[:,0], data[:,1], c=[target_to_color(t) for t in target])
    plt.savefig(filename)

def gen_data():
    iris = load_iris()
    pca = PCA(n_components=2)
    return pca.fit_transform(iris.data), iris.target

def main():
    data, target = gen_data()
    plot_data(data, target, filename="origin.png")

    km = KM(n_clusters=3)
    km_target = km.fit_predict(data)
    plot_data(data, km_target, filename="kmeans.png")

    cm_result = cmeans(data.T, 3, 2, 0.003, 10000)
    plot_data(data, cm_result[1].T, filename="cmeans_2.png")

    gmm = GMM(n_components=3, max_iter=1000)
    gmm.fit(data)
    gmm_target = gmm.predict(data)
    gmm_target_proba = gmm.predict_proba(data)
    plot_data(data, gmm_target, filename="gmm.png")
    plot_data(data, gmm_target_proba, filename="gmm_proba.png")

if __name__ == "__main__":
    main()

結果

オリジナルデータ

元データ

　これが元のデータです。できるだけこれに近いようなクラスタリング結果を得ることを目標とします。

k-means

　図の左側のクラスタは分離できていますが、右側は割と悲惨です。クラスタ同士が隣接していて細長い形だったりすると上手く行かないことが多いのがk-meansの特徴です。

c-means

Fuzzy c-means

　こうして見るとc-meansは「ファジー理論を入れて境界を曖昧にしたk-means」という気がしてきます。実際アルゴリズムもそんな感じなんですけど。

GMM

GMM-based clustering (hard)

GMM-based clustering (soft)

　一見して「おお」って感じですね。k-means、c-meansと比較して、元データのラベルに近いクラスタリング結果が得られています（図の右側の2つのクラスタの境界が右肩上がりになっている）。まあ、ちょっと元データのラベルとはずれているんですが（右下の方はかなり怪しい）、普通はこちらの方がk-meansやc-meansより「良い」クラスタリング結果だ、と判断されることが多いでしょう。

　どうしてこうなるのかというと、「irisのデータが正規分布していた」ということに尽きます。ま、アヤメの花びらの大きさとかのデータですから、正規分布しているんでしょう、きっと。

　こうして見るとGMMの方が良さそうな気もしますが、「ちゃんと正規分布してるか」が怪しいとちょっと適用するのを躊躇うのと、あと計算コスト自体はk-meansより高いはずなので*2、いまいちk-meansと比べて使われていない、というのが実情に近いかもしれません。

まとめ

　GMMを使ってみたらけっこう良かったです。

*1:一次元なら平均と分散、多次元なら共分散みたいな話になってくるのだろうか？

*2:Fuzzy c-meansとどっちが良いかは未調査

【python】skfuzzyのFuzzy c-meansでソフトクラスタリング

2018-03-03T20:25:58+09:00

はじめに

　Fuzzy c-meansはソフトクラスタリングの手法です。

　pythonではscikit-fuzzyというライブラリで利用できるようです。ということで、実際に使ってみました。

　基本的な理論はこちらのサイトなどを参考にしてください。

ファジィc-means法 - 機械学習の「朱鷺の杜Wiki」

はじめに
sklearnではなかった
基本的な使い方
- パラメタ
- 返り値
実験
結果
まとめ
- ※追記
付録

sklearnではなかった

　ドキュメントはここです。
Module: cluster — skfuzzy v0.2 docs

　親切なことに使用例のページもあります。
Fuzzy c-means clustering — skfuzzy v0.2 docs

　一読してわかるとおり、sklearnライクなインターフェースがある訳ではないようです。それも一つの方針として悪くはないと思いますが、ちょっとカルチャーショックを受けました。

基本的な使い方

　sklearnライクではないので、使い方を理解するまで手間取りました。簡単に（コメントで）説明しておきます。

# cmeans『関数』をimportする
# 間違ってもクラスなどではないことを理解すること
from skfuzzy.cluster import cmeans

# dataは(データ数,次元数)のいつもの特徴量
# が、転置して渡してやる必要がある
# 5つのパラメタは省略できない。詳しくは後述
cm_result = cmeans(data.T, 3, 2, 0.003, 10000)
# 返り値も実際は7つくらい返っている
# 必要なものだけ使う

　率直な感想は「作りが古いなぁ・・・」です。「sklearn？　なにそれ食えるの。オブジェクト指向？　知らんよ」みたいな。別に使う側としてはどんなインターフェースでも、解析アルゴリズムがちゃんと走れば良いので、文句言うつもりはないんですが。それでも正直なところ、ちょっと引きました。

　とりあえずパラメタと返り値について説明します。上記ドキュメントを適当に（省略しつつ）和訳しただけなので、ちゃんと使いたい人はドキュメントを読んで下さい。

パラメタ

data : 2d array, size (S, N)

　データ。Sはデータの次元数、Nはデータ数です。

c : int

　クラスタ数

m : float

　c-meansのパラメタ。詳しくは上の理論のページを読んで下さい。変えたときの変化は後で見せます。

error : float

　fuzzy c-meansは処理を繰り返して収束させていく系のアルゴリズムなので、こういうものが必要になります。エラー率が一定以下になったらループを抜ける訳です。理論をちゃんと勉強していないので、エラー率の計算方法とか私は知りませんけれど。

maxiter : int

　最大繰り返し回数。

init : 2d array, size (S, N)

　これは省略できます。計算に使う初期値です。省略するとランダム配列にされます。普通「よい初期値」なんて持ってないので、省略して使うことになると思います。

seed : int

　initが省略されたときのランダム初期化で使うseedです。

返り値

cntr : 2d array, size (S, c)

　各クラスタの中心。

u : 2d array, (S, N)

　各データがどのクラスタに割り当てられたかを比率で示した行列です。要するにこれがc-meansの結果です。あと、ドキュメントには(S, N)って書いてあるけど、これはどう考えても(c, N)の間違いです。

u0 : 2d array, (S, N)

　uの初期値みたいなことが書いてある。

d : 2d array, (S, N)

　「Final Euclidian distance matrix.」。クラスタ中心に対する距離ってことか？

jm : 1d array, length P

　「Objective function history.」だって。さっぱりわからない。

p : int

　ループが回った回数

fpc : float

　「Final fuzzy partition coefficient.」

　まあ、要するに適当な引数で呼び出してあげて、実質的に使う返り値はuだけということです・・・。

　とにかく使い方はわかったので、実験してみましょう。

実験

　irisでやってみます。irisは4次元データですが、4次元だと見づらいのでPCAで2次元に落とし、2次元空間上で3クラスタ（irisが3種類のアヤメの花のデータなので）にクラスタリングすることにします。

　結果として、

元のデータを色分けした絵
k-meansでクラスタリングした結果（比較用）
mを1.5, 2, 3, 4, 5で変化させたときの結果それぞれ

　を得ることにし、考察します。

　実験に使ったソースコードは記事の最後に付録として載せます。次の章で先に結果を見せます。

結果

　1つずつ見せます

元のデータ

　似たような絵を何回も見たことがあります。いつも感じることですが、irisは綺麗なデータで良い子だと思います。

k-means

　上とはクラスタの色が違いますが、気にしないでください。以下のデータでも（出力するたびにクラスタ番号変わるので）色はバラバラです。

　右側の境界をうまく捉えられていませんが（クラスタリングだから当たり前）、それを除けば真っ当な結果に見えます。

Fuzzy c-means

　上から順にm=1.5, 2, 3, 4, 5です。各クラスタへの帰属の度合いをそのままRGBにしました。

m=1.5

m=2

m=3

m=4

m=5

　少しずつ色が暗くなっていく傾向にありますが、これは「どのクラスタにも同程度に帰属する」みたいな結果になっていることを表します。m=1.5のときは[0.7, 0.2, 0.1]みたいな結果が、m=5のときは[0.4, 0.3, 0.3]みたいな結果が得られていると解釈してください。

　どの程度のmを設定するとちょうど良いのか・・・？　私にはわかりません。要するにmを小さくするとシャープに切れ、大きくするとゆるーく移り変わるってことなんですが。適当に決めるしかないってことなのかな。

　それはそうとして、よく見ると各クラスタの中心付近が（相対的に）鮮明な色になっていて、クラスタ中心から外れると（正確には他のクラスタ中心に接近していくと）色が混ざり合うのがお分かりいただけたでしょか？　これがソフトクラスタリングの成果です。

まとめ

　とりあえず、使える（意図したとおり動かせる）のはわかりました。

　何に使えるかは・・・ごめんなさい、思いつきそうで思いつかなかったです。

　一応フォローしておくと、ソフトクラスタリング自体は「高次元の特徴量を相対的に低次元に持ってくる」とか「元の空間の情報を保持したまままったく別の空間に移す」みたいな用途のために使われることがあるようです（各クラスタへの帰属度合いを新たな特徴量とする）。たとえばこんな応用があるらしいです。

qiita.com

　ただ、上のページでは混合正規分布モデルという別の手法を使って実装していますけど・・・。どっちが良いんだろう？　ぶっちゃけ謎（たぶん既存の知見があるんだろうけど）。そのへんも含めて今後検討していきたいところです。

※追記

　混合正規分布モデルでもやりました。こちらの記事をご参照ください。
hayataka2049.hatenablog.jp

付録

　実験に使ったソースコードです。python3系で必要なパッケージを入れればそのまま動くと思います。

▶クリックで展開

# coding: UTF-8

import numpy as np

from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans as KM

from matplotlib import pyplot as plt

from skfuzzy.cluster import cmeans

def target_to_color(target):
    if type(target) == np.ndarray:
        return (target[0], target[1], target[2])
    else:
        return "rgb"[target]

def plot_data(data, target, filename="fig.png"):
    plt.figure()
    plt.scatter(data[:,0], data[:,1], c=[target_to_color(t) for t in target])
    plt.savefig(filename)

def gen_data():
    iris = load_iris()
    pca = PCA(n_components=2)
    return pca.fit_transform(iris.data), iris.target

def main():
    data, target = gen_data()
    plot_data(data, target, filename="origin.png")

    km = KM(n_clusters=3)
    km_target = km.fit_predict(data)
    plot_data(data, km_target, filename="kmeans.png")

    for m in [1.5, 2, 3, 4, 5]:
        cm_result = cmeans(data.T, 3, m, 0.003, 10000)
        plot_data(data, cm_result[1].T, filename="cmeans_{0}.png".format(m))

if __name__ == "__main__":
    main()

【python】sklearnのPipelineを使うとできること

2018-02-22T23:40:11+09:00

　機械学習では、何段もの前処理をしてから最終的な分類や回帰のアルゴリズムに入力するということがよくあります。

　前処理にはけっこう泥臭い処理も多く、leakageの問題なども絡んできます。はっきり言って自分で書こうとすると面倒くさいです。

　こういう問題を（ある程度）解決できるのがsklearnのPipelineです。これについては、以前から「何かあるらしいな」というのは知っていましたが、実際に使ったことはありませんでした。でも、このたび使ってみたら「すげえ」となったので、こうして記事を書いている訳です。

　この記事ではPipelineのコンセプトと使い方を簡単に説明します。雰囲気は伝わるかと思いますが、細かい使い方はライブラリの公式ドキュメントを参照してください。

sklearn公式

sklearn.pipeline.Pipeline — scikit-learn 0.21.3 documentation

問題設定
実装
結果
参考サイト

問題設定

　今回は例として、sklearnのdigits（load_digits）を対象データにして説明します。

sklearn.datasets.load_digits — scikit-learn 0.21.3 documentation

　これは0～9の数字を分類する問題で、特徴量は8*8の画像データをflattenして64次元にしたものです。このデータの分類は割と簡単な方で、直接SVMにでもかけてパラメタを追い込めばF1値にして0.95以上のスコアが得られたりするのですが、もうちょっとタチの悪いデータのつもりで扱います。

　具体的には、以下の処理をしてやることにします。

RandomForestの特徴重要度を使って特徴選択
PCAで累積寄与率pに次元削減
SVMで分類
グリッドサーチでパラメタチューニング
交差検証して性能評価

　大変そうです。でも、Pipelineを使えばすぐできます。

実装

　まず、上の1～3について、それぞれの部品を作ります。それからPipelineのインスタンスで一つにまとめます。グリッドサーチと交差検証は自分で書くことにします。

特徴選択

　これにはSelectFromModelを使うと良さそうです。分類器のfeature_importances_に基づき、重要度の高い特徴だけ残してくれます。

sklearn.feature_selection.SelectFromModel — scikit-learn 0.21.3 documentation

　何次元残すかの指定ができると使いやすかったのですが、実際はmean, medianとそれらのfloat倍、そして重要度の下限を指定できるようです。とりあえずmeanとmedianのどちらかにしてみましょう。

　分類器には今回はRandomForestを使います。つまり、次のようなコードになります。

from sklearn.ensemble import RandomForestClassifier as RFC
from sklearn.feature_selection import SelectFromModel

rfc = RFC(n_estimators=100, n_jobs=-1)
fs = SelectFromModel(rfc)

　パラメタはあとでGridsearchSVを使ってチューニングするので、今の段階で指定する必要はありません。

次元削減

　次元削減にはPCAを使います。普通にやるだけなので説明は割愛します。

from sklearn.decomposition import PCA

pca = PCA()

分類

　分類にはSVMを使います。これもインスタンスを作っておきます。

from sklearn.svm import SVC

svm = SVC()

パイプライン化

　今回の記事のキモです。パイプラインを使って上記「特徴選択」「次元削減」「分類」をすべてまとめてしまいます。

　書き方はこんな感じです。

from sklearn.pipeline import Pipeline

estimators = zip(["feature_selection", "pca", "svm"], 
                 [fs, pca, svm])
pl = Pipeline(estimators)

　とてもあっさりしていますが、これで特徴選択をし、次元削減して、分類するという一連の流れをまとめて行うインスタンスができました。

パラメタチューニング

　パラメタチューニングはGridsearchCVを使うと簡単？にできます。

sklearn.model_selection.GridSearchCV — scikit-learn 0.21.3 documentation

from sklearn.model_selection import GridSearchCV

parameters = {"feature_selection__threshold" : ["mean", "median"],
              "pca__n_components" :[0.8, 0.5],
              "svm__gamma" : [0.001, 0.01, 0.05],
              "svm__C": [1, 10]}

clf = GridSearchCV(pl, parameters)

　ちょっとパラメタ指定周りが面倒くさいですが、とにかくこうすれば後は全部自動でやってくれます。パラメタは「モデルにつけた名前__（アンダーバー２つ）パラメータ名」という形で書いてください。

交差検証

　これはStratifiedKFoldを使い、後は自分で書きます。クロスバリデーションをやってくれる関数もsklearnにはありますが、今回はちょっと複雑な制御（1回目のFoldでパラメタチューニングして2回目以降はそのパラメタを使いまわしたい）をするので使いません。

結果

　最終的なソースコードはこうなりました。

# coding: UTF-8

import numpy as np

from sklearn.datasets import load_digits
from sklearn.pipeline import Pipeline
from sklearn.model_selection import StratifiedKFold as SKF, GridSearchCV
from sklearn.ensemble import RandomForestClassifier as RFC
from sklearn.svm import SVC
from sklearn.feature_selection import SelectFromModel
from sklearn.decomposition import PCA
from sklearn.metrics import precision_recall_fscore_support as prf

def main():
    rfc = RFC(n_estimators=100, n_jobs=-1)
    fs = SelectFromModel(rfc)
    pca = PCA()
    svm = SVC()
    estimators = zip(["feature_selection", "pca", "svm"], 
                     [fs, pca, svm])
    pl = Pipeline(estimators)
    
    parameters = {"feature_selection__threshold" : ["mean", "median"],
                  "pca__n_components" :[0.8, 0.5],
                  "svm__gamma" : [0.001, 0.01, 0.05],
                  "svm__C": [1, 10]}

    gclf = GridSearchCV(pl, parameters, n_jobs=-1, verbose=2)

    digits = load_digits()
    X = digits.data
    y = digits.target
    first_fold = True
    trues = []
    preds = []
    for train_index, test_index in SKF().split(X, y):
        if first_fold:
            gclf.fit(X[train_index], y[train_index])
            clf = gclf.best_estimator_
            first_fold = False
        clf.fit(X[train_index,], y[train_index])
        trues.append(y[test_index])
        preds.append(clf.predict(X[test_index]))

    true_labels = np.hstack(trues)
    pred_labels = np.hstack(preds)
    print("p:{0:.6f} r:{1:.6f} f1:{2:.6f}".format(
        *prf(true_labels, pred_labels, average="macro")))

if __name__ == "__main__":
    main()

　実行するとGridsearchCVのverboseがたくさん出力された後、スコアが出てきます。スコアは今回は

p:0.948840 r:0.948205 f1:0.948379

　でした。

　かなり面倒くさい処理なのに、Pipelineのおかげでシンプルに書けているのがお分かりいただけたでしょうか。

　便利なのでこれから色々使っていこうと思います。

参考サイト

Scikit-learnのpipeleine.Pipelineが便利 - KAZ log TechMemo
Scikit learnよりグリッドサーチによるパラメータ最適化 - Qiita

【python】正準相関分析（Canonical Correlation Analysis）を試してみる

2018-02-16T02:13:14+09:00

　正準相関分析を使うと、2つの多次元データ同士の関連性を分析できるらしい。

　面白そうなので試してみた。ちなみに正準相関はsklearn.cross_decomposition.CCAで使える。正準相関自体の解説はほとんどしないので、文中のリンクを参考にして欲しい*1。

一応概要だけ
ノイズに埋もれた波形を取り出す
もうちょっとデータ分析っぽいことをしてみる

一応概要だけ

　代表的な多変量解析の手法（といって良いのかどうか少し悩むけど）として、主成分分析や重回帰分析が存在する。

主成分分析：一つの多変量データを直交するより少ない変数に縮約する
重回帰分析：一つの多変量データを一つの単変量データに変換する

　主成分分析にしろ重回帰分析にしろ、変換の係数だったり行列だったりを求めてそれで変換するのが実際にやることである。

　さて、正準相関は上の流れで説明すると、

正準相関分析：二つの多変量データをそれぞれ直交するより少ない変数に縮約して、かつ二つの変換されたデータの間で相関を最大化する

　という目的の分析である。主成分分析と重回帰を混ぜた感じ。

　気づいた人もいると思うけど、多変量vs多変量のデータでどちらかを単変量に分解して個別に重回帰で解くことも可能である。それに対するメリットとしては、

個別に重回帰するより全体の構造みたいなものを捉えられる可能性がある
個別に重回帰すると係数の数が全体でとても多くなるので解釈が面倒くさいが、一度次元を下げて直交した空間に持っていくことでそこが楽になる

　というあたりがあり、要するに解釈性がいいということ。

　この説明でもよくわからん、という人は、ニューラルネットのオートエンコーダーとか思い浮かべていただくとかえってわかりやすいかもしれない。

ノイズに埋もれた波形を取り出す

　参考URLの通りにやることにする。

　単一の信号源に複数のプローブを当てていて、それぞれに独立のノイズが乗って信号が埋もれてしまった・・・みたいな状況から元の信号を取り出そうとしているらしい。脳波計測とかで使えるのだろうか？

　参考URL：https://www.jstage.jst.go.jp/article/jnns/20/2/20_62/_pdf

　とりあえずこのようなコードを書き、

# coding: UTF-8
import numpy as np
import matplotlib.pyplot as plt
from sklearn.cross_decomposition import CCA

def plot_wave(data, filename):
    fig, [ax1,ax2] = plt.subplots(2,1,figsize=(16,9))

    ax1.plot(data[:,0], color="b")
    ax2.plot(data[:,1], color="r")

    plt.savefig(filename)

def gen_pulse_data():
    common_pulse = np.array([-1]*50 + [0]*50 + [1]*50 + [0]*50 + 
                            [-1]*50 + [0]*50 + [1]*50 + [0]*50, dtype=np.float64)
    common_pulse += (np.random.random(common_pulse.shape) - 0.5)*0.1

    noise1 = (np.random.random(common_pulse.shape) - 0.5)*50
    noise2 = (np.random.random(common_pulse.shape) - 0.5)*50

    data1 = np.vstack([common_pulse + noise1,  common_pulse - noise1]).T
    data2 = np.vstack([common_pulse + noise2,  common_pulse - noise2]).T

    return data1, data2

def main():
    X1, X2 = gen_pulse_data()
    plot_wave(X1, "X1.png")
    plot_wave(X2, "X2.png")

    cca = CCA(n_components=2)
    cca.fit(X1, X2)

    Y1 = cca.transform(X1)
    Y2 = cca.transform(X2)
    plot_wave(Y1, "Y1.png")
    plot_wave(Y2, "Y2.png")

if __name__ == "__main__":
    main()

　実行する。

　まずは元の信号。

　わかる訳ねえな、という感じ。

　CCAでX1とX2の相関が最大になるような変換を計算し、その変換に基いてX1を変換したものをY1とすると、

　こんな感じになった。Y2も同じようなものなので省略。ここでは2次元出しているが、元のパルス信号が1次元なので2次元目（下）はノイズだけ出ている。

　まあ、上手く動いているのではないだろうか。

もうちょっとデータ分析っぽいことをしてみる

　如何せん、「ノイズに埋もれた信号を取り出せる！」というだけでは、データ分析っぽくなくて（個人的には）面白くない。正準相関自体はもっと色々なことができるはず。

　ここで足を引っ張るのは「正準相関向きのサンプルデータが見つからない」ということである。

正準相関向きのデータを探すのは困難

　2つの多次元データが対応しているようなデータで、適当にわかりやすいものがあれば良いのだが・・・なかなか良いデータがない。上に挙げた解説論文でも、「知名度は低い」とか書かれちゃってるし、正準相関自体、ニッチな感じがする。そこが素敵なのだが。

　一応、ネット上にある解説例だと、

統計学入門−第19章

　医学の分野で、肝機能の検査値（複数）と腎機能の検査値（複数）の対応を見るとか、

http://ogasun.la.coocan.jp/hanbetsubunseki.pdf

　中学生の体格（身長、体重、座高とか）と運動能力（50m走、走り幅跳びとか）の対応を見るとか、

　そういう感じのことをやっているのだが、この手のデータを探してくるのがまず面倒くさいし、見つけてもプログラムに流し込めるようにするまでがまた苦行だろうな、ということは容易に想像できるのである。

　この点で悩んで、この記事も一週間くらい出すか出さないか迷ってたんだけど、やることにした。ただ、結局良いデータは見つからなかったので、それっぽくでっちあげることにした。

作成したデータ

　ある架空の中学校で集計したという設定の、20人の生徒のデータである。「学外での勉強や取り組み」と「学校の成績」が対応付いている。

　「学外での勉強や取り組み」には、

一ヶ月に何冊読書するか
一年に何回博物館に行くか
毎週何日塾に通っているか
毎日何時間自習しているか

　の4つの変数がある。一方、「学校の成績」は、

国語
数学
社会
理科
英語

　の5つの科目があり、5段階評価で成績が付く。

　本来であれば適当に線形モデルでも作ってあげて数字を作るべきところだが、面倒くさいので私の想像で適当に埋めた（ツッコミポイント）。

　一応、次のような方針を考え、それに沿ったデータになるようにでっちあげた・・・つもり。

読書量と国語の成績は比例する
博物館に行った回数と社会、理科の成績は比例する
塾に通う頻度、自習時間は成績全体に影響を及ぼす

　よって、こういう結果が出てくるか、という勝負になる。

実験と結果

　こういうプログラムを書いた。

# coding: UTF-8
import numpy as np
from scipy.stats import pearsonr
from sklearn.cross_decomposition import CCA

def gen_data():
    # X1:
    # 毎月何冊の本を読むか,
    # 一年に何回博物館に行くか,
    # 塾に週何日通うか,
    # 毎日何時間自習するか
    X1 = np.array([[1,0,2,1],
                   [3,2,4,2],
                   [0,0,2,0],
                   [9,4,2,1],
                   [1,1,3,1],
                   [8,1,6,3],
                   [0,9,7,8],
                   [2,2,4,1],
                   [5,0,0,1],
                   [2,0,4,0],
                   [0,0,7,8],
                   [4,4,2,2],
                   [5,1,2,1],
                   [1,1,5,2],
                   [8,6,2,1],
                   [0,0,0,1],
                   [6,1,3,1],
                   [2,0,3,1],
                   [4,8,5,3],
                   [5,0,1,1]])

    # X2:
    # 国語,数学,社会,理科,英語の成績
    X2 = np.array([[3,3,3,3,3],
                   [4,3,4,4,5],
                   [2,2,3,3,2],
                   [5,4,3,3,3],
                   [3,3,4,4,4],
                   [5,5,5,4,5],
                   [3,5,5,4,5],
                   [4,4,4,5,3],
                   [5,3,3,3,3],
                   [3,4,3,4,3],
                   [5,5,4,5,5],
                   [4,4,5,5,3],
                   [4,3,3,3,3],
                   [4,4,5,4,5],
                   [5,3,5,5,3],
                   [2,2,2,1,2],
                   [5,3,4,4,4],
                   [3,4,3,4,3],
                   [5,5,5,5,5],
                   [5,3,3,3,3]])
    return X1, X2

def main():
    X1, X2 = gen_data()

    cca = CCA(n_components=4)
    cca.fit(X1, X2)

    print("Correlation Coefficient")
    for i in range(4):
        print("{0}:{1:.4f}".format(i, pearsonr(cca.x_scores_[:,i], cca.y_scores_[:,i])[0]))

    print("")
    np.set_printoptions(formatter={'float': '{: 0.4f}'.format})
    print("X1 loadings")
    print(cca.x_loadings_.T)
    print("")
    print("X2 loadings")
    print(cca.y_loadings_.T)

if __name__ == "__main__":
    main()

　「学外での勉強や取り組み」=X1と「学校の成績」=X2を4次元の空間上に写像して相関を最大化する、という問題を解かせる。軸同士は直交していて無相関なので、写像したデータの各軸の値同士の相関だけ見てやれば良い。写像したデータは、cca.x(or y)_scores_かcca.transform(X1(or X2))で取得できる*2。

　あとはX1とX2の各成分が、写像先の各軸にどれだけ寄与しているかがわかれば良い。そのためにはcca.x(or y)_loadings_を見る。転置した方が見やすいのでそうしている。

　こうして見ると、PCAに似ている。実際、CCAとPCAは親戚らしい。ま、あまり理論的な話に深入りしてもボロが出るので、これくらいにしておく。

　さて、結果はこのようになった。

Correlation Coefficient
0:0.9558
1:0.8978
2:0.5980
3:0.2927

X1 loadings
[[-0.4224  0.4244  1.0047  0.7353]
 [ 0.9326  0.4315  0.2511  0.3450]
 [ 0.2350  0.9269 -0.1714 -0.2514]
 [-0.4357  0.4062 -0.0634  0.8007]]

X2 loadings
[[-0.0558  0.7802  0.6574  0.6206  0.8044]
 [ 0.9392  0.5254  0.4702  0.3366  0.4271]
 [-0.2184 -0.0705  0.9353  0.5442 -0.3528]
 [-0.4310  0.0676 -0.2136 -0.8751  0.0008]]

　まず見るべきはCorrelation Coefficientで、写像先の空間の軸にどれだけ相関（＝やった意味）があるかを示している。0,1,2次元目はまあまあ強い相関だが、3次元目は相関係数0.3じゃ大した意味はなさそうだな、という風に解釈しておく。

　次にX1 loadingsとX2 loadingsを見る。X1 loadingsは4*4、X2 loadingsは4*5で、つまり行が写像先の軸、列が元の空間の軸に対応するように表示している。

　X1 loadingsの各行を見ていくと、

1行目

　塾と自習に熱心

2行目

　読書

3行目

　博物館

4行目

　自習と博物館だけ？

　なんとか解釈できる。数字がでかいところだけ重視するのがこつ。X2 loadingsも同様にやると、

1行目

　国語以外のすべて。国語にはほぼ中立。特に強いのは英語

2行目

　国語。他もそれなりに

3行目

　社会と理科

4行目

　理科にとてもネガティブ。全体的にネガティブな感じ

　ここまで出揃えば後はなんとかなる。このデータを作った方針を再掲する。

読書量と国語の成績は比例する
博物館に行った回数と社会、理科の成績は比例する
塾に通う頻度、自習時間は成績全体に影響を及ぼす

　0次元目は「塾に通う頻度、自習時間は成績全体に影響を及ぼす」とに、1次元目は「読書量と国語の成績は比例する」に、2次元目は「博物館に行った回数と社会、理科の成績は比例する」に対応していることがわかり、まあ妥当な結果なんじゃないの、という気はする。相関係数の低い3次元目はそこまで重視する必要はない。

　今回は先に方針を決めてデータをでっち上げたのであまり感動がないような気もするが、実際はデータにどんな構造があるのかは分析してみないとわからない。その構造を理解する上で正準相関が役に立つことは、上の例でなんとなく理解できた。

*1:正準相関でググって1ページ目に出てくるようなページばかり・・・

*2:今回はどちらも同じ値が返るが、transformだと学習時とは違うデータも入れられる

主成分分析 - 静かなる名辞

SVMのsupport vectorを可視化してみた

はじめに

見方

コード

まとめ

sklearnとmatplotlibでiris（3クラス）の予測確率を可視化した話

はじめに

まずやる

他の分類器も試す

まとめ

【python】高次元の分離境界をなんとか2次元で見る

はじめに

方法

実験

まとめ

【python】PCAと非負値行列因子分解のバイプロットを見比べる

はじめに

比較実験

メリット

まとめ

本当は怖いSVMと交差検証

概要

コード

結果

怖くない線形SVM

怖いかどうか悩むランダムフォレスト

怖くない気がする多層パーセプトロン

SVMも怖くない！

現実的な話

まとめ

【python】sklearnのOneClassSVMを使って外れ値検知してみる

はじめに

実験

結果

まとめ

【python】sklearnのFeatureAgglomerationを使ってみる

はじめに

使い方

実験

考察

まとめ

【python】sklearnのSparsePCAを使ってみる

はじめに

SparsePCAとは？

sklearnの実装

実験

まとめ

GridSearchCV『の』パラメータ・チューニング 高速化中心に

はじめに

下準備とベースライン

cvを指定する（効果：大）

return_train_score=Falseする（効果：それなり）

まとめ

それでも時間がかかりすぎるときは

【python】MeanShiftのbandwidthを変えるとどうなるか実験してみた

プログラム

結果

結論

【python】sklearnのMeanShiftクラスタリングを試してみる

はじめに

使い方

実験

プログラム

結果

結論

【python】複数の特徴をまとめるFeatureUnion

複数の目的変数で回帰を行う方法

はじめに

目的変数の数だけ回帰モデルを作る方法

複数の目的変数に対応したモデルを使う

正準相関分析

ランダムフォレスト回帰

多層パーセプトロン（ニューラルネットワーク回帰）

まとめ

【python】SOMのライブラリSomocluはかなりおすすめ

SOMの概要

ライブラリがない

それでも頑張ってググった

使ってみた

GridSearchCV『の』パラメータ・チューニング高速化中心に