静かなる名辞

pythonとプログラミングのこと


日本語モダリティ解析器 Zundaを試す

自然言語処理 構文解析

　日本語のモダリティを解析できるらしい。「文中のイベント（動詞や形容詞など）に対して、その真偽判断（イベントが起こったかどうか）、仮想性（仮定の話かどうか）などを解析します」とのこと。

　公式ページはたぶんここ。

jmizuno.github.io

環境

　ubuntu14.04

インストール

　CaboCha (>=0.60)と、Boost Library (>=1.41)を予め入れおく必要がある。CaboChaは入ってたけどBoost Libraryはなかったので、apt-getで入れた。

$ sudo apt-get install libboost-all-dev

　後はtarballを落としてきてmakeで入れろと公式に書いてある。どんなエラーが出てくるかとびくびくしながらやったけど、まったく問題なく入った。

$ ./configure
$ make
$ sudo make install

試してみる

　こうやって使うらしい。

$ echo -e "次郎は大阪に行った。\n太郎は東京には行かず地元に残ろうとした" | zunda
#EVENT0	4	wr:筆者	非未来	0	叙述	成立	0	0
* 0 2D 0/1 -2.249829
次郎	名詞,固有名詞,人名,名,*,*,次郎,ジロウ,ジロー
は	助詞,係助詞,*,*,*,*,は,ハ,ワ
* 1 2D 0/1 -2.249829
大阪	名詞,固有名詞,地域,一般,*,*,大阪,オオサカ,オーサカ
に	助詞,格助詞,一般,*,*,*,に,ニ,ニ
* 2 -1D 0/1 0.000000
行っ	動詞,自立,*,*,五段・カ行促音便,連用タ接続,行く,イッ,イッ
た	助動詞,*,*,*,特殊・タ,基本形,た,タ,タ
。	記号,句点,*,*,*,*,。,。,。
EOS

#EVENT0	5	wr:筆者	未来	0	叙述	不成立	0	0
#EVENT1	9	wr:筆者	未来	0	意志	高確率	ポジティブ	0
#EVENT2	12	wr:筆者	非未来	0	叙述	成立	0	0
* 0 4D 0/1 -1.650377
太郎	名詞,固有名詞,人名,名,*,*,太郎,タロウ,タロー
は	助詞,係助詞,*,*,*,*,は,ハ,ワ
* 1 2D 0/2 0.320510
東京	名詞,固有名詞,地域,一般,*,*,東京,トウキョウ,トーキョー
に	助詞,格助詞,一般,*,*,*,に,ニ,ニ
は	助詞,係助詞,*,*,*,*,は,ハ,ワ
* 2 4D 0/1 -1.650377
行か	動詞,自立,*,*,五段・カ行促音便,未然形,行く,イカ,イカ
ず	助動詞,*,*,*,特殊・ヌ,連用ニ接続,ぬ,ズ,ズ
* 3 4D 0/1 -1.650377
地元	名詞,一般,*,*,*,*,地元,ジモト,ジモト
に	助詞,格助詞,一般,*,*,*,に,ニ,ニ
* 4 -1D 3/4 0.000000
残ろ	動詞,自立,*,*,五段・ラ行,未然ウ接続,残る,ノコロ,ノコロ
う	助動詞,*,*,*,不変化型,基本形,う,ウ,ウ
と	助詞,格助詞,引用,*,*,*,と,ト,ト
し	動詞,自立,*,*,サ変・スル,連用形,する,シ,シ
た	助動詞,*,*,*,特殊・タ,基本形,た,タ,タ
EOS

　公式ページには改行すれば別の文として扱われると書いてあるが、echoに-eオプションを渡さないと改行文字は改行として出力されないので注意。

　それはそうと、この結果の見方だが、

#EVENT0	4	wr:筆者	非未来	0	叙述	成立	0	0」

　という結果は4番目の形態素がどんなモダリティなのかを表す。つまり「行っ 動詞,自立,*,*,五段・カ行促音便,連用タ接続,行く,イッ,イッ」に対応する。単純だけどわかりやすいかどうかは微妙かもしれない。

色々なことを試す

　とりあえず、もうちょっと色々な文を入れてみる。

$ echo "遊びに行きたいな" | zunda
#EVENT0	2	wr:筆者	未来	0	欲求	0	ポジティブ	0
* 0 1D 0/1 0.000000
遊び	名詞,一般,*,*,*,*,遊び,アソビ,アソビ
に	助詞,格助詞,一般,*,*,*,に,ニ,ニ
* 1 -1D 0/2 0.000000
行き	動詞,自立,*,*,五段・カ行促音便,連用形,行く,イキ,イキ
たい	助動詞,*,*,*,特殊・タイ,基本形,たい,タイ,タイ
な	助詞,終助詞,*,*,*,*,な,ナ,ナ
EOS

　なるほど。

$ echo -e "吾輩は猫である。\n名前はまだない。" | zunda
* 0 1D 0/1 0.000000
吾輩	名詞,代名詞,一般,*,*,*,吾輩,ワガハイ,ワガハイ
は	助詞,係助詞,*,*,*,*,は,ハ,ワ
* 1 -1D 0/2 0.000000
猫	名詞,一般,*,*,*,*,猫,ネコ,ネコ
で	助動詞,*,*,*,特殊・ダ,連用形,だ,デ,デ
ある	助動詞,*,*,*,五段・ラ行アル,基本形,ある,アル,アル
。	記号,句点,*,*,*,*,。,。,。
EOS

#EVENT0	3	wr:筆者	非未来	0	叙述	成立	0	0
* 0 2D 0/1 -2.377508
名前	名詞,一般,*,*,*,*,名前,ナマエ,ナマエ
は	助詞,係助詞,*,*,*,*,は,ハ,ワ
* 1 2D 0/0 -2.377508
まだ	副詞,助詞類接続,*,*,*,*,まだ,マダ,マダ
* 2 -1D 0/0 0.000000
ない	形容詞,自立,*,*,形容詞・アウオ段,基本形,ない,ナイ,ナイ
。	記号,句点,*,*,*,*,。,。,。
EOS

　ふーん。

　まあ、勝手はなんとなくわかった。それはそうと、zundaにはpythonバインディングなんて親切なものはないらしい。zundaオリジナルなのは#から始まる行だけで、その下はcabochaの-f 1フォーマットがそのまま出てるだけなので、プログラミングで使う側としてはテキスト処理でゴリ押して使うことになるだろう。

　その際は、コマンドラインオプションで何百kbも受け渡しするのは流石にアレなので、まとめてファイルを処理させるか、立ち上げておいて標準入出力でやりとりするかのどちらかになる。

【python】nltkで英語の形態素解析

python nltk 自然言語処理

概要

　形態素解析、いわゆるPOS taggingと呼ばれるようなタスクをnltkを使うと簡単に行なえます。日本語の解析ではmecabやjumanを使うと思いますが、英語だとnltkに入っているものが使えるので（インストールとかが）楽です。

　目次

スポンサーリンク



使い方

　凝ったやり方は幾らでもあるのですが、簡単のために次の二つを使います。

  • nltk.word_tokenize
  • nltk.pos_tag

　word_tokenizeは文を入力に受け取り、単語のリストを作ります。pos_tagはその単語のリストを受け取って、POSタグを推測してくれるようです。

　では、実際に試してみます。なお、最初に使うときは「解析に必要なリソースがないよ」という感じのエラーが出ると思います。親切なことにリソースの入れ方がエラーメッセージで出るので（nltk内にインストーラーが実装されててメソッド呼ぶだけで入る）、その通り対処してください。インターネット接続だけは必要です。

　リソースをダウンロードしたら、対話的インタプリタで以下のように打ち込んでみます。

>>> import nltk
>>> string = "I AM A CAT. As yet I have no name."
>>> words = nltk.word_tokenize(string)
>>> words
['I', 'AM', 'A', 'CAT', '.', 'As', 'yet', 'I', 'have', 'no', 'name', '.']
>>> nltk.pos_tag(words)  # 可読性のため結果を改行しています
[('I', 'PRP'), ('AM', 'VBP'), ('A', 'DT'), ('CAT', 'NNP'), ('.', '.'), 
 ('As', 'IN'), ('yet', 'RB'), ('I', 'PRP'), ('have', 'VBP'),
 ('no', 'DT'), ('name', 'NN'), ('.', '.')]

　このように解析が行えました。

各タグの意味

　各タグの意味は次のコマンド（というかメソッド呼び出し）で出ます。

>>> nltk.help.upenn_tagset()

　ものすごくたくさん出てくるし、英語なので、必要になったら皆さんで使ってください（NLPガチ勢の人的には「こんなの当たり前にわかるよ！」でそもそも必要ないのかもしれないが・・・）。

　なお、上の例で出てきたタグの説明だけ抜粋すると、

  • PRP

　pronoun, personal

  • VBP

　verb, present tense, not 3rd person singular

  • DT

　determiner

  • NNP

　noun, proper, singular

  • .

　sentence terminator

  • IN

　preposition or conjunction, subordinating

  • RB

　adverb

  • NN

　noun, common, singular or mass

　抜粋したのは説明の表題だけで、実際は3行くらいの説明が一緒に付いてます。なのでちゃんと読めばどのタグが何を意図してるかわからなくて困ることは、ないはずです・・・。

まとめ

　nltkさえ入れれば解析できる英語って言語は良いな・・・と思いました。

関連記事

www.haya-programming.com
　単語に分割する前にまず文に区切りたい、という場合はこちらを使います。

www.haya-programming.com
　ステミングとかです。

【python】文字列を一文字ずつのリストにする

python Tips

　文字列はそもそもiterableなので、これが必要なことは滅多にないんだけど・・・。

>>> for c in "hogehoge": # そのままfor文に渡してオッケー。
...     print(c)
... 
h
o
g
e
h
o
g
e

　この前必要になったので（ライブラリの引数でlist型を要求された）考えてみる。

　目次

スポンサーリンク



考えられる方法

空リストを用意してfor文で一文字ずつappend

　さすがに面倒くさいので論外。

リスト内包表記

　pythonは何でもリスト内包表記で書ける。よってこの処理もリスト内包表記で書ける。

>>> [c  for c in "hogehoge"]
['h', 'o', 'g', 'e', 'h', 'o', 'g', 'e']

　メリットは、

  • 他の方法と比べると可読性（不可誤読性）と簡潔さのバランスが良い

　デメリットは、

  • 冗長
  • ぶっちゃけforループなので遅い

　ぶっちゃけ長いので、積極的に使いたくはない。

list(str)する

　list(str)で一文字ずつのリストになる。

>>> list("hogehoge")
['h', 'o', 'g', 'e', 'h', 'o', 'g', 'e']

　メリットは、

  • 簡潔
  • 組み込みなので速いと思われる

　デメリットは、

  • 誤読されるかも？

　誤読されるというのは難癖に近いが、

>>> list("hogehoge") # -> 他人に["hogehoge"]だと解釈される可能性がある！

　ということ。ないとは言い切れない。

starred expression

　こんな書き方もできる。

>>> [*"hogehoge"]
['h', 'o', 'g', 'e', 'h', 'o', 'g', 'e']

　メリットは、

  • 恐らく一番簡潔

　たった3文字の追加でリストにできる！　list()だと6文字なので半分の記述量で済む

  • 見た目がカッコいい

　デメリットは、

  • 古いバージョンだと確か使えない
  • 可読性は良くないかも。というか始めてみた人は面食らう気がする
  • 速度性能が未知数

　しょうじき見た目は一番好きなので積極的に使っていきたいのだが、遅いと困る・・・。

時間を計測してみた

　計測しました。

# coding: UTF-8

import time

def time_measure(string, splitter, n=1000):
    time_list = []
    for i in range(n):
        t1 = time.time()
        splitter(string)
        t2 = time.time()
        time_list.append(t2-t1)
    return sum(time_list)/n

def main():
    string = "hoge"*100000
    
    lm = lambda s:None
    f1 = lambda s:[c for c in s]
    f2 = lambda s:list(s)
    f3 = lambda s:[*s]

    # 念の為空のラムダも回す
    print(time_measure(string, lm))
    print(time_measure(string, f1))
    print(time_measure(string, f2))
    print(time_measure(string, f3))

if __name__ == "__main__":
    main()

　結果は、

1.652240753173828e-07
0.008317208766937256
0.00259556245803833
0.0026724901199340822

　意外なことに、list(str)と[*str]は互角だった（プログラム走らせる度に誤差で優劣が変わる程度）。

　個人的には古いバージョンのpythonは使わないし、可読性もどっこいどっこいだと思うので、今後は[*str]で行くことにしました。

　追記：上のように思っていた時期もありました。でもやっぱりトリッキーな感じがするので、たいたいlist(str)ばっかり使っています。この方が白python的な気もするし……。[*str]は黒python寄り。

【python】scipy.statsのzscoreで警告が出るときの対策

python scipy 統計

概要

　z得点を計算しようとしたとき、このような警告を見かけることがあります。

RuntimeWarning: invalid value encountered in true_divide

　これが出た場合、結果にはnanが含まれています。なので後段の分析で落ちたりします。

>>> import numpy as np
>>> from scipy.stats import zscore
>>> a = np.array([[1,2,3,4], [1,2,3,4],[1,3,4,5]])
>>> zscore(a, axis=0)
stats.py:2248: RuntimeWarning: invalid value encountered in true_divide
  return (a - mns) / sstd
array([[        nan, -0.70710678, -0.70710678, -0.70710678],
       [        nan, -0.70710678, -0.70710678, -0.70710678],
       [        nan,  1.41421356,  1.41421356,  1.41421356]])

　どうしてエラーになるかというと、z得点は標準偏差で割るので、標準偏差が0だと0除算エラーが発生するからです。標準偏差0の列が含まれるようなゴミだらけの汚い疎行列をそのまま入れると、これが出来ます。

対策

　どうせ標準偏差0の軸とか要らないので、予め消し飛ばしておく。

>>> a[:,np.std(a, axis=0) != 0]
array([[2, 3, 4],
       [2, 3, 4],
       [3, 4, 5]])
>>> zscore(a[:,np.std(a, axis=0) != 0], axis=0)
array([[-0.70710678, -0.70710678, -0.70710678],
       [-0.70710678, -0.70710678, -0.70710678],
       [ 1.41421356,  1.41421356,  1.41421356]])

　めでたしめでたし。

【python】numpyでデータをランダムサンプリング

python numpy Tips 機械学習

　機械学習に使うデータをランダムサンプリングしたいときがある。簡単そうなのにやり方が見つからないから自分で書く。

　目次

スポンサーリンク



実装方針

重複ありランダムサンプリング

　これはnp.random.randintでn個のindexを作り、Fancy Indexingするだけ。

重複なしランダムサンプリング

　ちょっと面倒くさいというかやり方が幾つか思いつくが、とりあえず今回方針としては、

  1. オリジナルデータのindexの配列を作る
  2. shuffleする
  3. 先頭からn個取り出してsample_indexとする
  4. そのsample_indexを使ってオリジナルデータに対してFancy Indexingし、取り出す

　これで行くことにする。

実装と結果

　以上を踏まえて、次のようなソースを書いた。

# coding: UTF-8

import numpy as np
from sklearn.datasets import load_iris

def sampling_Xy(X, y, sample_n=50, dup=False):
    n_data = X.shape[0]

    if dup:
        samples_index = np.random.randint(0, n_data, sample_n)
    else:
        original_index = np.arange(n_data)
        np.random.shuffle(original_index)
        samples_index = original_index[:sample_n]
    return X[samples_index], y[samples_index], samples_index

def main():
    iris = load_iris()
    print(iris.data.shape, iris.target.shape)
    sample_X, sample_y, samples_index= sampling_Xy(
        iris.data, iris.target, sample_n=80, dup=False)
    print(sample_X.shape, sample_X.shape, len(set(samples_index)))
    sample_X, sample_y, samples_index = sampling_Xy(
        iris.data, iris.target, sample_n=80, dup=True)
    print(sample_X.shape, sample_X.shape, len(set(samples_index)))

if __name__ == "__main__":
    main()

　結果は、

(150, 4) (150,)
(80, 4) (80, 4) 80
(80, 4) (80, 4) 59

　となり正しくランダムサンプリングできていることが確認できた。

そもそもなにに使いたかったの？

　手軽に実験を回してフィードバックを得て改良していくにはでかすぎるデータセットとか重すぎるアルゴリズムが世の中にはたくさんあるので、そういうのに対して使いたかったんです・・・。

裏技

　この方法だと教師ラベルの比率にまでは気を遣っていないので、特に少ないデータ数では問題（特定クラスに偏る）が発生する可能性がある。

　手軽に回避する方法として、StratifiedKFold.splitの返り値の先頭の先頭だけ使うことが考えられる。本来の使い方ではないけど、コード量は圧縮できそう。

　ただし、元データの教師ラベルの比率が均等な保証はないので、厳密に元データと同じ比率でサンプリングしたい・・・とか考え出すと無駄だし、そこまでやる意義も少ない。けっきょく、ランダムサンプリングで良いじゃん、ということになる。

※追記（参照することをおすすめします）

　numpyでできるもっと良い方法を見つけました。

　こっちを使うと、自分で書く量を減らして簡単にできそうです。

www.haya-programming.com

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

python sklearn 次元削減 主成分分析 線形判別分析 Pipeline 機械学習


　一般によく使われる次元削減手法としては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値）
結果（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による速度差を比較

python sklearn 主成分分析 次元削減 機械学習 20newsgroups CountVectorizer 自然言語処理 Tips 速度計測シリーズ

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

　とにかくこのPCAを高速化したかったので、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とは微妙に違う数字になっている）

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

　それぞれの使い分けは、

  1. 入力次元数の小さい入力ではfullで良い。というかヘタにそれ以外を指定するとかえって遅いケースもある
  2. 入力次元数が大きく、入力次元数>>出力次元数で厳密解がほしければならarpackの使用を検討する
  3. 厳密解じゃなくても良いのでとにかく速いのを！　ってときは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だったりしたらどうなるんだろう？　それとも共分散行列のサイズなのだろうか？

【python】tfidfは分類精度を向上させるのか？→向上しなかった

python 自然言語処理 sklearn 20newsgroups CountVectorizer TfidfVectorizer tf-idf 特徴抽出 機械学習 ランダムフォレスト

　目次

はじめに――長年の疑問

　自然言語処理でテキスト分類などに、よくtf-idfが使われます（最近はそうでもないのかもしれないが）。一般には、tf-idfを使うことで分類精度の向上効果があると認識されているようです。

　このことを長年疑問に思っていました。tf-idfのうち、tfは文書中の単語の出現回数（あるいは相対頻度）ですから、単なるBag of Wordsと変わりません。また、idfは文書全体でのその単語の出現する文書数の対数みたいなものですから、文書集合全体で各単語に1つのidfが定まります。

　けっきょく、tfi-dfはBoWにidfの列ベクトルをかけたものとみなせそうです。ということは、とても単純な線形変換ですから、こんなもので本当に分類精度が上がるんかいな？　という疑問をずっと抱いてました。分類器のアルゴリズムによってはある程度効果は期待できるかもしれないが（特に単純なものなら：k近傍法とか）、たとえば確率分布として取り扱うナイーブベイズや、線形変換をかけまくって分類できる軸を探すLDA(Linear Discriminant Analysis：線形判別分析）、あるいは決定木で分類に有効な特徴を探し出すRandomForestのような手法ではまったく効かないんじゃないの、という仮説をずっと考えていました。

　せっかくなので検証してみます。

スポンサーリンク



検証

　検証用のデータは、sklearnのdatasetsから使える20newsgroupsにしました。fetch_20newsgroupsで使えます。

　ただし、このデータは量が多くて（1.1万件ほど）処理を回すのが大変なので、全体のだいたい40%をランダムサンプリングすることにしました。

　また、min_df=0.03, max_df=0.5, stop_words="english"を指定し、予め次元数を501次元にしています。ここがスタートラインです。

　このデータから4種類の方法で特徴量を作りました。

  1. 単語の出現回数（CountVectorizerで直接作成）
  2. 1をl2ノルムで割って正規化したもの
  3. tf-idf（TifdfVectorizerで生成。norm=Noneを指定して正規化なしの条件でやる）
  4. 正規化tf-idf（norm="l2"を指定）

　これらに対し、以下の分類器で交差検証を回して分類スコアを計算しました。

  • ナイーブベイズ（Gaussian Naive Bayes）
  • k近傍法（K Nearest Neighbors）
  • 線形判別分析（Linear Discriminant Analysis）
  • SVM（Support Vector Machine）
  • ランダムフォレスト（RandomForest Classifier）

　以下に検証に使ったソースコードを載せておきます。

結果

　次のような結果になりました。

f:id:hayataka2049:20180319124128p:plain

　分類器別に端的にまとめると、

  • ナイーブベイズ：正規化効果あり。tf-idf効果なし
  • k近傍：正規化、tf-idfともに効果あり
  • 線形判別分析：正規化効果あり。tf-idf効果なし
  • SVM：正規化効果あり。正規化なしtf-idfまったくダメ。正規化+tf-idfは効いてる可能性あり
  • ランダムフォレスト：正規化、tf-idfともに恐らく効果なし。効果があるとしても1%以下

　という結果になりました。要するに、

  1. tf-idfは基本的に無力
  2. tf-idfをする暇があったらl2ノルムで割る。こっちの方が効く
  3. ただし一番精度を出せてるRandomForestでは、l2ノルムで割る正規化すら大して効いていないので、どこまで意味があるかは正直微妙

　という結論です。

tf-idfは死んだのか？

　少なくとも文書分類における特徴抽出手法としては死んだ、と言って構わないでしょう。

　このことは仮説の通りだったので、驚きはあまりないです。tf-idfは『分類精度を上げる目的』ではほとんど使えないというのが結論です。

　ではなぜtf-idfがこれまでもてはやされてきたのか？　相対頻度への変換や正規化によって、生BoW（単語の出現回数数えただけ）より良い結果が得られてきたためではないでしょうか。肝心のidfによる重み付けは「ちっとも意味がない」と言わざるを得ないと思います。

　では、tf-idfは使えない子なのか？　分類に使う特徴量としては上記の通り無意味ですが、tf-idfは特徴語抽出に使えます。tf-idf（の文書集合内の平均）が高すぎず低すぎない単語を抜き出すことで、文書の特徴をよく表す単語を抽出するという使い方です。これは教師なしで計算の軽い特徴選択手法として利用できますから、そっちでは役に立つでしょう、たぶん。

まとめ

　思った通り向上しませんでした。

ランダムフォレストとSVMの使い分け

機械学習 ランダムフォレスト SVM Tips

はじめに

　ランダムフォレスト（RandomForest）とSVM(Support Vector Machine)はよく比較される分類器です。でも、様々なシチュエーションで、けっきょくどちらを使うべきなのか、という指針はあまり見かけません。

　私は研究などで*1両者を使ってきて、それなりに両者のクセもわかるようになってきたので、「こんな感じで使い分ければ良いんじゃない？」という情報を簡単にまとめておきます。

　結論を先に言うと、

  • 次元数とデータ数が多いか？

　Yesならばランダムフォレスト
　NoならばSVM

  • データの軸に明確な意味があるか？

　Yesならば（どちらかといえば）ランダムフォレスト
　Noならば（どちらかといえば）SVM

  • 最後は両方試して選ぼう
  • ランダムフォレストは多芸だけどSVMはそうでもないよ

　こんな感じです。一つずつ説明していきます。

　目次

スポンサーリンク



次元数とデータ数

　低次元でデータが少なければSVM、高次元でデータが多ければランダムフォレストです。これは精度云々より計算コストの問題です。

　ランダムフォレストは実用的な精度を得ようとすると、ある程度木の本数を増やす必要があります（500～1000くらいが相場、データによってはもっと）。それだけ多量の計算を回す必要がある上、決定木の処理というのはそれなりに重くてGPGPUなどによる並列化の恩恵も受けづらい（なにせ本質的にIF文の森）ので、どうしてもコストが高く付きます。なので、簡単な問題であればあるほどSVMの方が安上がりです。

　ただし、SVMは計算量（オーダー）が重く、高次元で大規模データであればあるほど重くなります。厳密な定義は難しいようですが（カーネルのアルゴリズムによっても変わるだろうし）、たとえば下記ページを見ると、 O(n\_samples^2 * n\_features) なんて計算量が出ています（あくまでも一例として示されています。もっと速い場合もあるとか書いてあるし）。

https://www.quora.com/What-is-the-computational-complexity-of-an-SVM

　その点ランダムフォレストはもうちょっと色々マシなので*2、データや次元数が増えていくとどこかでSVMと逆転します。

軸の意味

　これについてはっきり断言してくれた文献をこれまでほとんど（まったくかも）見たことがないんですが、ぶっちゃけ一番重要なんじゃないですかね。計算コストのことは置いておいて性能の優劣を論じる場合、だいたい軸に意味があるかどうかで決まってくる印象です。

　軸の意味といっても釈然としない方もいるかもしれませんが、たとえば二次元空間上に適当な正規分布を二個置いて決定境界を描く、みたいな例がよくありますが、ああいうのが「軸に意味がない」ケースです。とにかく空間があって、その上でクラスタを形成しているのを選り分けるようなケース。

　一方、軸にもうちょっと明快な意味がある場合もあります。典型的な例はテキスト分類で使われるBoW(Bag of Words)で、これは単語の出現回数をそのままベクトルにしたものです。100000単語のボキャブラリがある文書集合は100000次元のベクトルになります。こういう例では、軸の意味は非常にはっきりします（BoWならその軸に対応する単語の出現回数）。つまり、いわゆる頻度データと呼ばれるものでは、ランダムフォレストの有効性が高くて、SVMはあまり効かない印象があります。

　ランダムフォレストは決定木ですから、つまるところ「この軸の値がこの範囲なら……」という処理を繰り返して分類を行います。だから軸に直行する分離境界しか引けないとか言われたりします。下のページを見ると、そのことはよくわかります。

Classifier comparison — scikit-learn 0.21.3 documentation

　だけど、それが即欠点という訳ではなく、かえってよく働くケースは多々あると思います。上述したBoWは実はランダムフォレストを使うと上手く分類できるデータの好例で、私があるタスクで使ったときは、BoW*3のデータを分類させたランダムフォレストはSVMに対して10%くらいの優位を叩き出していました。そしてSVMをいくらチューニングしてもランダムフォレストには追いつきませんでした。

　考えてみるとわかりますが、SVMは基本的に滑らかな分離超平面を引くようにできています。そう言うとなんとなく良さそうですが、千次元とか二千次元の滑らかな超平面を学習データから完全に学習するのは、はっきり言って無理難題です。データが少なかったりノイジーだったりするとなおさら無理で、どうやっても破綻します。それより、さっさと決定木を作ってしまった方が無理がない訳で。

　ただし、だからといってランダムフォレスト最強でSVMが駄目だ、ということを言いたい訳ではないのに注意してください。滑らかな分離超平面が必要なときはSVMの方が適しています。たとえば、テキスト分類でも、BoWという巨大な疎行列を扱うのではなく、word2vecやらdoc2vecやらを使って低次元で密な分散表現を得るというアプローチもあります。そして、そっちで作った特徴量だとSVMの方がランダムフォレストより有利になったりもします*4。

　要するに、欲しいものは「滑らかな（軸に対して斜めだったり曲がってたりする）分離超平面」なのか、「軸に直交する分離超平面」なのかをちゃんと考えてあげる必要があるということです。もちろんある程度は融通が効くので、この通りやらなくてもある程度の性能は出ますが、最終的に性能を追い込もうとするとこの辺りのことがとても響いてきます。

最後は両方試して選ぼう

　まあ、正直言って「どちらかを選ぶ」というシチュエーションなら両方試して（交差検証でF値出して）軽くて性能が高い方選べば良いと思います。pythonならsklearnのインターフェースは統一されていますから、そんなに難しいことじゃありません。というか、そういう目的で便利に使う（色々試して最適な方策を決める）ためにsklearnなんかのライブラリが整備されているのだし、データサイエンス・機械学習は最終的にはトライアル＆エラーになるのが普通なんじゃないでしょうか。もっと言えば、選択肢に載せるアルゴリズムをランダムフォレストとSVMに限定する必要もなく、行けそうなものは全部試せば良いと思います。

　ぶっちゃけ、アルゴリズムを選ぶ理論的な根拠なんてカッコつけ以上の意味はないんです。分類問題は評価指標が出るんだから、「やってみたらこれが一番良かったです」で根拠としては十分でしょう。

　ただし、アルゴリズムの特性がわかっていると、

  • 明らかに無駄な試行は減らせる。どう見てもでかすぎるデータにSVM使って一晩回し続けるとか
  • ある程度は挙動（結果）が予想できるようになる。やる前からだいたいこうなるだろう、と見当がつくので、その通りにならなかったらなにか失敗してるかもしれないとか、未知の現象が介在してるかもしれないとか、そういことがわかる*5
  • 色々な難癖を付けたがる人は多いので、理論武装は色々しておくに越したことはない

　こういうご利益はあるので、この記事も皆様の何らかのお役に立てば良いな、と思います。

ランダムフォレストは多芸

　ランダムフォレストには次のような能力があります。

  • OOB誤り率が計算できる。交差検証の代わりに使える。パラメータチューニングには十分
  • 特徴重要度が計算できる。これを特徴選択に使える
  • predict probabilityが自然に計算できる

　この辺が面白いからランダムフォレストが好きだ、という人はけっこう多いんじゃないでしょうか。実用的にも便利なので、これらを期待して使うシチュエーションもあると思います。

まとめ

　どっちが良いかはデータ次第、というのが結論です。機械学習やデータサイエンスは、データの性質（分布）さえわかっていれば、後は知識と経験とカンでやっても場合によってはなんとかなります*6。

　まあ、そういう「知識と経験とカン」に依存する要素を減らそうということで今はニューラルネットワークが流行っている訳ですが、ランダムフォレストやSVMのようなアルゴリズムも実際にはまだまだ使われていると思うので、当分の間は「使い分け」を意識してあげることも大切だと思います。

*1:ったって卒研ですが

*2:一本の木を構築するのに選ぶデータ数や次元数をどうするかにもよるが、通常は O(\sqrt{n})くらいで済む。

*3:本当は幾らでもあるBoW風の特徴量のうちの幾つか。生BoWではない

*4:これも実際に経験しました。ちなみにBoW+RandomForestより分散表現+SVMの方が性能面では当然有利です。分散表現はここ数年流行ってるので、ランダムフォレストの重要性は相対的に低下したと思います

*5:ただし、これに頼りすぎると「よしよし思い通り上手く行ってるな」→実は上手く行ってないけどなかなか気づかない、というパターンにたまにハマる。なので気をつける必要がある。コンピュータは間違えないけど使う人間は普通に間違えるので、眉にたくさん唾を付けておくに越したことはない

*6:慣れるとパラメータもダロカンで決めて、一発で最適パラメタを探り当てちゃうとか（笑）

【python】クラスタリング結果をエントロピーで評価する

python sklearn クラスタリング 評価指標 機械学習

はじめに

　クラスタリング結果の良し悪しを評価したいことがあります。

　正解ラベルがないデータに対してクラスタリングを行った場合（つまり本当に教師なし学習でやる場合）、基本的にクラスタ内距離二乗和やクラスタ中心間の距離などを使ってやる以外の方法はありません*1。

　しかし、正解ラベルがあるデータに対してクラスタリングをかける、というシチュエーションもあります。たとえばクラスタリング手法の性能評価だとか、相手にするデータセットが教師なしでやってどれくらいまとまる性質があるか調べるとか、そんな感じのときです。

　こういうとき、評価指標がないと困ります。そこで、評価指標を考えます。何通りかあるようですが、簡単に計算できるのはエントロピーです。

説明

　エントロピー（厳密にはシャノン情報量）は次式で定義されます。

　 H(P) = - \sum_{A\in\Omega} P(A) \log_2 P(A)
　 \Omegaは対象とする事象を表します。

　厳密な理論や定義などはwikipediaでも読んでいただくとして、エントロピーは乱雑さの指標になります。scipyで計算できますからやってみましょう。

>>> from scipy.stats import entropy
>>> entropy([0.5, 0.5],base=2)
1.0
>>> entropy([1, 0],base=2)
0.0
>>> entropy([0.25, 0.25, 0.25, 0.25],base=2)
2.0
>>> entropy([1, 0, 0, 0],base=2)
0.0

　渡しているリストの乱雑さを測っていると見ることができます。中身が一様（=乱雑）なら取り得る値の \log_2となり、乱雑でない（=どれかが1で他が0）なら0となります。

　ということは、クラスタごとにエントロピーを計算してやると、クラスタの中身が特定の正解ラベルに偏っているならば値は0、まんべんなくすべての正解ラベルを含んでいれば \log_2 正解ラベル数となり、0に近いほど良いという性質の評価指標になります。エントロピー自体はクラスタごとに計算できますが、クラスタに割り振られたデータ数で重み付き平均を取ることで全体の結果としてまとめます。

実装

　ライブラリがありそうだけど見つけられなかったので、結局自前実装することに。sklearn.metricsにあると思ったら、0～1のレンジに収まるような改良されまくった指標ばっかりでそのものズバリのエントロピー加重平均はありませんでした（と、思う。見落としてたら恥ずかしい）。

　2.3. Clustering — scikit-learn 0.20.1 documentation

def entropy_score(true_y, cl_y):
    df = pd.crosstab(cl_y, true_y)
    arr = df.as_matrix()

    weight = arr.sum(axis=1)
    weight = weight/(weight.sum())

    cl_entropy = np.array([entropy(cl, base=2) for cl in arr])
    return (cl_entropy*weight).sum()

　もうちょっと綺麗な書き方は色々あると思いますが、こんなので動きました。scipyのentropyは引数の和が1になるように勝手に揃えてくれるので、それに頼り切ってコードを書いています。

実験

　irisでKMeansとGMM(混合ガウスモデル)によるクラスタリングを比較してみましょう。GMMの方が良さそうなのは以前の記事で確認しているので、その定量的評価を今回はするということです。

# coding: UTF-8

import numpy as np
import pandas as pd

from scipy.stats import entropy

from sklearn.datasets import load_iris
from sklearn.cluster import KMeans
from sklearn.mixture import GaussianMixture as GMM

def entropy_score(true_y, cl_y):
    df = pd.crosstab(cl_y, true_y)
    arr = df.as_matrix()

    weight = arr.sum(axis=1)
    weight = weight/weight.sum()

    cl_entropy = np.array([entropy(cl, base=2) for cl in arr])
    return (cl_entropy*weight).sum()

def main():
    iris = load_iris()
    km = KMeans(n_clusters=3, n_init=30, max_iter=1000)
    cluster_labels = km.fit_predict(iris.data)
    print("km")
    print(entropy_score(iris.target, cluster_labels))
 
    gmm = GMM(n_components=3)
    cluster_labels = gmm.fit(iris.data).predict(iris.data)

    print("gmm")
    print(entropy_score(iris.target, cluster_labels))

if __name__ == "__main__":
    main()

　結果は、

km
0.39388631839664884
gmm
0.1611488952045549

　やっぱりGMMが良かったです。

蛇足

　「それにしてもsklearnにないのは変だなぁ」と思い、上記資料を読んだ感じだとhomogeneity_scoreが今回やったことに近いようですが、

  • 条件付きエントロピーできっちり定義
  • 0~1になるように元データのデータ比率を log_2したもので割り、1から引く

　という感じでした。そのうちちゃんと調べます。あと、実務的に使う場合はsklearnのhomogeneity_score、completeness_score（homogeneityのクロス集計を逆転させて計算しようという発想に近い）、v_measure_score（homogeneityとcompletenessの調和平均。要するにF1値みたいなもん）でやれば良いので、今回書いた方法の出番はないと思います。まあ、何かの参考にしてください。

*1:こちらの記事が参考になります。クラスタリング結果の評価の尺度基準 - froglog

【python】クラスタリング結果を積み上げ棒グラフで可視化する

python pandas matplotlib クラスタリング Tips 機械学習 可視化

はじめに

　ラベル付きデータをクラスタリングすることがよくあります（そんな頻繁にあるか？　まあ、クラスタリングの使い方次第でたまにはあるからこうして記事にしている訳ですが）。

　各クラスタの中身がどんなラベルで構成されているのか、知りたくなります。積み上げ棒グラフで出てくれると嬉しいですね（嬉しさがわからない方も読み進めて頂ければわかるので大丈夫）。

　pythonでの積み上げ棒グラフの描きをググると、matplotlibを駆使した怖い（大変そうな）描き方がいくらでも出てくるのですが、そんなことで苦労したくないので簡単なやり方でやります。

やりかた

　こちらを参考にしました。
python - pandasのデータフレームから積み上げ棒グラフを作りたい - スタック・オーバーフロー
　ぜんぶpandasの機能でできるらしいです。素敵。要するにクロス集計してplotしてやれば良い、ということのようです。

やってみた

　irisとdigitsをクラスタリングし、上の方法を参考にグラフ化してみます。

　ソースコード

# coding: UTF-8

import numpy as np
import pandas as pd

from sklearn.datasets import load_digits, load_iris
from sklearn.cluster import KMeans

import matplotlib.pyplot as plt

def visualize(y_true, y_cl, y_names, filename="result.png"):
    y_true = np.array([y_names[y] for y in y_true])

    df = pd.crosstab(y_cl, y_true, 
                     rownames=["cluster number"], colnames=["true label"])
    print(df)

    plt.figure()
    df.plot(kind='bar',stacked=True, legend=False)
    plt.legend(bbox_to_anchor=(1.13, 1), loc="upper right")
    plt.savefig(filename)

def main():
    iris = load_iris()
    km = KMeans(n_clusters=3, n_init=30, max_iter=1000)
    cluster_labels = km.fit_predict(iris.data)
    print("iris")
    visualize(iris.target, cluster_labels, 
              iris.target_names, filename="iris.png")

    digits = load_digits()
    km = KMeans(n_clusters=10, n_init=30, max_iter=1000)
    cluster_labels = km.fit_predict(digits.data)
    print("digits")
    visualize(digits.target, cluster_labels, 
              digits.target_names, filename="digits.png")

if __name__ == "__main__":
    main()

　とりあえず、クロス集計した結果のDFをテキストで吐いてみたので、そっちから見せます。

iris
true label      setosa  versicolor  virginica
cluster number                               
0                   50           0          0
1                    0          48         14
2                    0           2         36
digits
true label        0    1    2    3    4    5    6    7    8    9
cluster number                                                  
0               177    0    1    0    0    0    1    0    0    0
1                 0   24  148    0    0    0    0    0    3    0
2                 0    0    3    7   10    0    0  177    5    8
3                 1    0    0    0  166    2    0    0    0    0
4                 0    1   13  155    0    2    0    0    2    6
5                 0    1    0    2    0  136    0    0    4    5
6                 0    2    0    0    0    1  177    0    2    0
7                 0    0    2   12    0   41    0    0   51  140
8                 0  100    8    7    2    0    3    2  101    1
9                 0   54    2    0    3    0    0    0    6   20

　なるほど。まあ、なんとなくどんな状態になっているのかはこの表からもわかります。

　そしてお待ちかねのグラフはこちら。

irisの結果
iris
digitsの結果
digits
　わかりやすいですね。この記事で言いたいことは、この絵を簡単に得られるpandasは便利、ということだけです。

　一応グラフの説明をすると、このグラフは各クラスタに割り振られたデータのラベルの件数を表しています。そして、たとえばirisのグラフからは、setosaは完全に一つのクラスタを形成していますが、versicolorとviriginicaは綺麗にクラスタには分かれず混ざる、ということがわかります。versicolorとviriginicaが似ている、というかベクトル空間上で近くにいる、という知見が得られる訳です。

まとめ

　pandasは独自の世界観があって正直苦手なんですが、たまに便利に使えることがあるなぁと思いました（小並感）。

【python】分類タスクの評価指標の解説とsklearnでの計算方法

python 評価指標 機械学習 sklearn

はじめに

　分類結果の評価指標として、混同行列（confusion matrix）、適合率（precision）、再現率（recall）、F1値（F1-measure）*1などがあります。

　分類の評価をやるときはとりあえずこれらを出せば良い、ということで日常的に用いられるものですが、意外とまとまった解説をネット上で見かけません。私もこれまでなんとなく使っていましたが、それじゃいかんなぁ、とずっと思っていました。

　この記事はこれらの評価指標について解説します。ついでにsklearnでの計算方法も書いておきます。

　目次

スポンサーリンク



理論の解説編

基本編（二値分類）

混同行列の話

　混同行列は分類結果を表、というか行列の形で表したものです。混同行列はすべての評価指標の基本となります。

　一つの数値として指標が得られる訳ではない、という意味では厳密には評価指標には含まれないかもしれませんが、これを理解しておくと、適合率・再現率・F1値もとても理解しやすくなります（というか混同行列がわからない状態で、この3つを理解するのは困難）。なので、最初は混同行列から理解するのが良いでしょう。

　混同行列とは、こんな奴のことを指します。
f:id:hayataka2049:20180314112047p:plain
（Excelで5分で描いた画像なので見栄えが寂しいのは許してください）

　重要なのは中央のTP, TN, FP, FNという4つの記号です。これはそれぞれTrue Positive, True Negative, False Positive, False Negativeの略です。

　TP, TNは「正しくPositive or Negativeに分類されました（本当にPositive or Negative）」という意味です。同様に、FP, FNは「間違ってPositive or Negativeになっちゃいました（本当はNegative or Positive）」という意味だと覚えておけば良いと思います*2。

　さて、実際にはTP, TN, FP, FNには、それぞれ対応する数字（条件に当てはまるデータの件数）が入ります。400件のデータを分類して混同行列を得たら、TP, TN, FP, FNをぜんぶ足し合わせるとちゃんと400になる、という具合です。

　評価指標を用いなくても、この混同行列を見て分類結果の良し悪しを判断することも（主観的になりそうですが）可能です。

　たとえば、400件の分類結果に対して、TP=150, TN=100, FP=50, FN=100という結果を得たとしましょう。まずTP+FN、TN+FPで真値がPositive, Negativeなデータの件数を計算すると、250、150という数字が出ます。400件の中にはPositiveが多め、ということがこれでわかります。次に正解率（accuracyと言います）を出してみましょう。これはTP+TNを全データ数400で割れば良く、つまり250/400=0.625という数字になります。あまり良くなさそうです。ではどんな風に良くないのかというと、FPは50、FNは100ということで、間違ってPositiveになってしまうことは相対的に少ないけど、間違ってNegativeを出している率が相対的に高いことがこれでわかります。

　このように、混同行列を見ることでだいたいどんな分類結果なのかを理解できます。また、他の評価指標には含まれない情報も含むので（データ、予測結果に実際にどんな比率でP/Nが含まれているのか、多クラス分類の場合はどのクラス同士で混同が起きやすいかなど）、けっこう重宝します。

評価指標の概念の説明

　混同行列から分類結果の良し悪しを判断できると言っても、上のようなことを一々考えるのはけっこう大変そうですね。そこで、数字一つで済む評価指標の出番になる訳です。

　先に言っておくと、適合率・再現率・F1値が表す情報はすべて混同行列に含まれています。実際、これらの評価指標は混同行列から計算できます。よくある説明は、混同行列の記号を用いた数式を描いて、それで終わり、というものです。

　しかし概念を理解することも重要なので、この記事では数式を出す前に言葉だけで説明してみます。数式は後ほど出します。

  • 適合率

　適合率（precision、精度と訳されることもあります）は、いわば「どれだけ嘘の分類結果が少ないか」の指標です。端的にいうと、Positiveと判定されたものの中にどれだけPositiveが含まれているかの割合です（ちなみにNegativeでも理屈の上では同様に計算できる。これについてはずっと後ろで後述する）。

　インフルエンザ検査で例えると、検査が陽性になったけど実はインフルエンザにかかっていませんでした、という人がどれだけいるかの指標になります。

  • 再現率

　再現率(rcall)は、「どれだけ取りこぼさないで分類できたか」の指標です。真値がPositiveの人の何割をPositiveと判定できたかを表します。

　インフルエンザ検査では、インフルエンザにかかってるのに検査しても陽性にならなかった、という人がどれだけいるかの指標になります。

  • F1値

　適合率や再現率はわかりやすい指標ですが、ではこれらを高めれば良いのか、というと一概にそうは言えない面があります。

　たとえば、インフルエンザ患者が50人、ただの風邪の人が50人いたとしましょう。

　100人中一番インフルエンザっぽい人に『だけ』陽性を出す検査というものがあったとして、その適合率はほぼ100%くになるに違いないでしょう。しかし、残りのインフルエンザ患者49人は陰性になってしまうのだから、ちょっと困ったものです。というか使い物になりません。

　同様に、100人全員インフルエンザ患者、という結果を出す検査があったとして、再現率は100%です。しかし、これは当然何の役にも立ちません。

　この適合率と再現率には、実はトレードオフがあります。上の例で、それぞれ適合率がほぼ100%のときの再現率、再現率が100%のときの適合率を考えてみると理解できるかと思います。そこで、中間くらいを取ってやれ、という指標を考えることができます。平均で良いじゃん、と思う訳ですが、実際には調和平均を使います。それがF1値（F1 measure）です。

　調和平均を使う理論的な根拠の説明はけっこう難しいようですが、こちらの記事が参考になります。

F値に調和平均を使う理由（再） - あらびき日記

　直感的な理解としては、調和平均は低い方に引っ張られる（それも高低差が大きければ大きいほど引っ張られる）平均なので、適合率、再現率どちらかに特化しすぎた結果はポイされることになり、とても都合が良い、と思っておけば良いかと思います。なお、調和平均については下のページが参考になります。

http://www004.upp.so-net.ne.jp/s_honma/mean/harmony2.htm

　F1値は分類結果の直接の良し悪しの指標になります。ですから、論文などではこれの改善*3を以て提案手法の有効性をアピールする、といったことがよく行われます。

評価指標を数式で書く

　お待ちかね（？）の数式です。といっても、基本的には上で書いたことを数式で表現するだけです。ここで混同行列を先に説明して、記号を導入しておいたことが生きてきます。

  • 適合率

Precision = \frac{TP}{TP+FP}

  • 再現率

Recall = \frac{TP}{TP+FN}
　※FNは「間違ってNになっちゃってるけど本当はP」なので TP+FNで「本当にP」のデータの総数になる。

  • F1値

F1-measure = \frac{2\cdot Precision \cdot Recall}{Precision+Recall}
　これは調和平均の定義通り式を書いていますが、PrecisionとRecallを展開して式を整理することができ、最終的には次のようになります。
F1-measure = \frac{2TP}{2TP+FP+FN}
　なんとなく適合率と再現率の式を混ぜた感じに見える・・・というか分子分母をそれぞれ足した値ですね。どうしてこうなるかというと、これは調和平均を計算したので、
\frac{1}{F1-measure} = \frac{1}{2}(\frac{TP+FP}{TP} + \frac{TP+FN}{TP})
　という式になる訳で、あとは逆数にしてやると上の式になることが理解できるかと思います。適合率、再現率でともに分子がTPなので、式が綺麗になるのがミソ。

どの指標を使えば良い？

　適合率・再現率・F1値を紹介しましたが、どれを重視するかはケースバイケースです。単純に考えるとF1値が良さそうですが、実際にはタスクによってどれを重視するかが変わってきます。

  • 例1：致死率が数十%の危険な感染症が海外で流行している。入国者を空港で検査し、検査が陽性なら隔離したい

　→この場合、FPを出す（偽陽性と言います）より、FNを出してすり抜けられてしまうことの方が遥かに大問題です。ですからとにかく再現率の高い検査を持ってこい、ということになると思います。適合率が10%とか、場合によっては1%であっても許容されるかもしれません（検査に引っかかった人全員を潜伏期間が終わると予想される時期まで隔離するとか、再現率99.9%、適合率1%の簡易検査で引っ掛けて精密検査でFPを弾くといった戦略もあり得る）。

  • 例2：単語で検索すると関連するWebページを出してほしい（Web検索）

　→この場合、その単語に関連するWebページがすべて表示されなくても、とりあえず5ページとか10ページ表示されれば実用上そんなに問題はないでしょう。しかし、無関係なページが結果に紛れ込んでしまうとユーザーエクスペリエンスを大きく損ないそうです。この場合、どちらかといえば再現率より適合率を重視すべき（FPを減らすべき）です。

  • 例3：ミカンとオレンジが混ざってしまったので分類したい

　→この場合、例1や例2とは違って、適合率と再現率のどちらが重要ということは特にありません。むしろ変にどちらかに偏った分類結果になる方が厄介なので、F1値で評価しておけば良いシチュエーションということになります。

　要するにFPとFNのどちらを重視するか、それともどちらも同程度に重視するかで着目する指標が変わってきます。もちろんすべて高いに越したことはないのですが、適合率と再現率にはトレードオフがある以上、最終的には「タスクに相応しい指標を重視する」ことも重要になってくると思います。

多クラス分類編

　これまではPositive or Negativeだけの二値分類について考えてきました。しかし、実際の分類問題では多クラス分類という、三値以上のクラスに分類するケースもけっこうあります。

　この場合、評価指標の計算方法が二値分類とは変わってきます。しかもややこしいことに、結果全体に対して評価指標を計算するか、複数のクラスそれぞれで評価指標を計算するか、それぞれで計算した場合はどう平均して全体の結果を出すかといった問題があります。

　このあたりの問題をちゃんと説明したWebページや書籍は、これまでほとんど見たことがありません。ある意味では適当にやってもなんとかなる*4世界だからかもしれませんが、はっきり言って問題だと思います。この記事ではできるだけ丁寧に説明してみるつもりです。

多クラス分類の混同行列

　とりあえず一例として、A,B,Cの3つのラベルを考えることにします。別にミカン、オレンジ、デコポンでも、巨峰、マスカット、デラウェアでも、なんでも構わないと言えば構わないのですが。

　多クラス分類でも二値分類の場合と同じように混同行列を描くことができます。混同行列は真のクラスラベルと予測先がわかっていれば描けます。行列の概形としてはこうなります。

f:id:hayataka2049:20180314112221p:plain

　ここでFA(B)といった記号を導入していますが、一般的に使われる記号ではないので注意してください。私が勝手に記号を付けただけです。これは「間違ってAになっちゃった（本来はB)」という意味です。以下ではこの記号を使って説明していきます。

　なお、上述の通り多クラス分類の場合では、評価指標の出し方に何通りかの方法があります。それだけ「本当のところどうなってるのよ」という疑問も出てきやすいので*5、二値分類の場合と比較して多クラス分類では混同行列による評価がより重視される傾向にあるかと思います。論文などを読み書きするときは、その辺りへの配慮も重要になってくるかもしれません*6。

クラスごとの適合率・再現率・F1値

　基本的に、多クラス分類の評価指標はクラスごとに計算するのが一番簡単です。クラスAに対する評価指標、Bに対する評価指標、Cに対する評価指標という形で計算することができます。

　つまり、以下の式のようになります。

P_A = \frac{TA}{TA+FA(B)+FA(C)}
R_A = \frac{TA}{TA+FB(A)+FC(A)}
F_A = H(P_A, R_A)
(P:Precision, R:Recall, F:F1-measure, H:harmonic mean)

　基本的な考え方は二値分類のときと変わりませんが、上の混同行列ではABCの3クラスがありますから、当然3つのクラスに対してそれぞれ適合率・再現率・F1値が求まります。

マクロ平均

　上の方法で各クラスごとに評価指標を計算した場合、各クラスごとの分類の状況はわかるものの、分類結果全体の良し悪しの判断は難しくなります。そこで、なんとか全体で一つの評価指標にまとめたい、というニーズが出てきます。

　もっとも単純なのがマクロ平均（macro average）です。大げさな名前が付いていますが、単なる全クラスの結果の平均に他なりません。適合率、再現率、F1値をそれぞれのクラスごとに求め、それを平均するだけです。

　あと、これは私の長年の疑問で、蛇足なので読み飛ばして頂いても構わないんですが、「各クラスごとにF1値を求めてからF1値のマクロ平均を計算する」方法だと、F1値のマクロ平均は「適合率のマクロ平均と再現率のマクロ平均の調和平均」と一致しない気がするんですが、それで良いんですかね・・・。で、恐らくこの事態を避けるために「適合率のマクロ平均と再現率のマクロ平均の調和平均」をF1値として評価指標に用いている論文も読んだことがあるんですが、どっちが妥当なんでしょうか。ご存知の方は教えて頂けると助かります。ちなみに後述するsklearnの関数はaverage="macro"を指定すると「各クラスごとにF1値を求めてからF1値のマクロ平均を計算する」でやってくれるみたいですが、そっちの方が良いのかな・・・。

重み付き平均

　マクロ平均には各クラスごとのデータ件数の偏りが反映されないという欠点、というか特性があります。1000件のデータがあるクラスも10件のデータがあるクラスも一緒くたに同じ重みで扱われる、ということで、1000件の方の評価指標は0.95なのに10件の方が0.6になってくれたせいで0.775くらいにされちゃったとか、逆に1000件の方が0.3でも10件の方が0.9なら0.6扱いになるなど、理不尽というか「それで良いの？」という結果が出ることがあります。

　これを改善したのが重み付き平均で、これも基本的にはデータ件数の比率で重み付き平均を計算してやるというわかりやすい処理をするだけですから、数式で説明するほどのことはありません。

　ただ、真のクラスの比率に基づいた重みを使うか、予測ラベルに基づいた重みを使うかという微妙な問題が、あるといえばあるような気はします。まあ、普通は真のクラスの方でやるんだと思いますが。

　ちなみにこの重み付き平均が論文などで使われているのは、ほぼ見たことがありません。上述のマクロ平均か、後述のマイクロ平均が使われていることが多いという印象です。

マイクロ平均

　マイクロ平均（micro average*7）は、「クラスごとに計算してどうにか平均する」という方針を取りません。

　ではどうやるかというと、まず混同行列全体でTP, FP, FNを集計し、それに基いて二値分類と同様に評価指標を算出します。

TP = TA+TB+TC
（TNはどうせ使わないので省略）
FP = FP_A+FP_B+FP_C = FA(B)+FA(C)+FB(A)+FB(C)+FC(A)+FC(B)
FN = FN_A+FN_B+FN_C = FB(A)+FC(A)+FA(B)+FC(B)+FA(C)+FB(C)

　まあ、簡単な話ですね。

　これには当然データ件数が効いてくるので、重み付き平均と似たような性質があります。じゃあ重み付き平均と一致するのかというと、上でちろっと書いた「真のクラスの比率に基づいた重みを使うか、予測ラベルに基づいた重みを使うか」の問題が首をもたげてきます。要するに「適合率に対しては予測ラベルの比率」、「再現率に対しては真のラベルの比率」で重み付き平均を計算すればマイクロ平均と同じ結果になるんだと思いますが*8、そんな面倒くさいことするくらいなら素直にマイクロ平均って言えば済む話ですね。なんとなくこちらの方が理論的に妥当そうな気もする、ということがあり、重み付き平均よりマイクロ平均が積極的に使われる結果を招いています。

どれを使うか

　クラスごと、マクロ平均、マイクロ平均の3択で考えることにします。重み付き平均は、選択肢に載せないことにしましょう。

  • クラスごとの評価指標

　ある意味これを見るのが一番間違いがない。ただし、クラス数と同じ個数出てくるので（適合率、再現率、F1値をそれぞれ計算するとクラス数*3個ですね！）、クラス数が多いと評価が困難になる。

  • マクロ平均

　クラスごとのデータ件数の比率を考慮しないのは駄目では？　と思われがちだが、クラスごとにデータ件数にばらつきがないケースではこれを使っても問題は起きない。また、データ件数が少ないクラスでひどい結果になっていたりした場合、マイクロ平均では埋もれてしまうがマクロ平均では敏感に（1/クラス数の比率で）評価指標に反映されるので、一概にそれが悪いとは言えない。むしろ「どのクラスも同じ価値がある」とみなすマクロ平均は、分類結果の良し悪しをストレートに反映するという考え方もできる。ただし、データ件数10件のクラスとかあったらそのクラスの評価指標はざっと0.1刻みになってしまう訳で、1000件のクラスと10件のクラスでマクロ平均するとほぼ10件のクラスに引っ張られて使い物にならん・・・という考え方もできる。要するにケースバイケース。

  • マイクロ平均

　クラスごとのデータ件数の比率が、実際に使ったりするときとだいたい同じであれば、全体的な良し悪しの指標として妥当そう。ただし、小さいクラスで結果が悪い場合は普通に埋もれるので、そういった観点で細かくケアしてあげないといけない気がする。

　正直、間違いがないのはクラスごと、マクロ平均とマイクロ平均は率直に言って甲乙つけがたいというか、どっちもどっちという感想です。紙幅に制限がなければ、マクロ平均、マイクロ平均、クラスごとの評価指標、混同行列をすべて示すのが無難というか誠実な気がしますが、そうも行かない場合が多いので、実際は「えいや」でマクロ平均かマイクロ平均のどちらかを選ぶ、といった感じになることが多いと思います。「結果が酷くない」「クラスごとのデータ件数に酷いばらつきがない」という条件なら、マクロ平均とマイクロ平均はどちらを選んでも大差はないですし・・・。

二値分類の取り扱いについて

　これまで多クラス分類について説明してきましたが、実は二値分類を「2クラスの多クラス分類」と解釈することも可能です。そしてPositiveとNegativeそれぞれで評価指標を計算したり、マクロ平均やマイクロ平均を計算することも当然可能です。

　二値分類のところで示した評価指標の計算式はけっきょくPositiveについてだけ計算していると解釈することができますが、Negativeについても同様に計算すると当然異なった数字が出てきます。そしてPositiveとNegativeでマクロ平均やマイクロ平均を取ってやると、当然Positiveだけで出した数値とは異なる結果になります。

　何が言いたいかというと、ミカンとオレンジを分類するときなどは本来「2クラスの多クラス分類」とみなしてマクロ平均やマイクロ平均で評価するのが妥当なはず、ということです。インフルエンザ検査は「ポジネガの結果になる二値分類」とみなしても構わないかもしれませんが・・・。

　この辺りのことをちゃんと説明してくれた資料にはこれまで出会ったことがありませんし、実用上これをどこまで厳密に使い分ける必要性があるのかもなんとも言えませんが、理屈の上では上で書いた通りになると思っています。ですから、二値分類の評価指標の算出は意外と気楽にはできないと思っています。「本当にその評価指標は妥当なの？」という質問は、実験結果を根底から否定されかねないとても厄介なものです。皆さんにはしっかり対策して挑んでほしいところです。

スポンサーリンク


実践編（sklearnでの計算方法）

　おまたせしました、pythonで計算する方法です。手っ取り早くやるためにsklearnを使います。

　評価指標の計算にはsklearn.metricsというモジュールを使います。便利な評価指標計算関数がたくさん入っているので、これを活用して計算していくことになります。

API Reference — scikit-learn 0.21.3 documentation

　基本的な使い方はどの関数も大体同じで、真のラベル、予測ラベルの順に引数を渡してあげると評価指標が返ります（返り値の形式は関数によって色々ある）。

混同行列の計算

　次のように計算することができます。

>>> true = [0,0,0,0,0,1,1,1,1,1,2,2,2,2,2,3,3,3,3,3]
>>> pred = [0,0,1,1,1,0,1,1,1,1,2,1,1,2,2,3,2,1,3,3]
>>> from sklearn.metrics import confusion_matrix
>>> confusion_matrix(true, pred)
array([[2, 3, 0, 0],
       [1, 4, 0, 0],
       [0, 2, 3, 0],
       [0, 1, 1, 3]])

　結果はnumpy配列で返りますので、画像形式で出力したい場合は別途matplotlibなどで描画する必要があります。

　seabornを使うのが簡単です。これについては過去に記事を書いたので、ご参考にどうぞ。
hayataka2049.hatenablog.jp

評価指標の計算

　適合率、再現率、F1値を計算する関数はそれぞれ存在しています。

  • sklearn.metrics.precision_score
  • sklearn.metrics.recall_score
  • sklearn.metrics.f1_score

　しかし、これらを利用するより、一括で結果を計算してくれる「precision_recall_fscore_support」を利用した方が簡単です。

>>> true = [0,0,0,0,0,1,1,1,1,1,2,2,2,2,2,3,3,3,3,3]
>>> pred = [0,0,1,1,1,0,1,1,1,1,2,1,1,2,2,3,2,1,3,3]
>>> from sklearn.metrics import precision_recall_fscore_support
>>> precision_recall_fscore_support(true, pred)
(array([0.66666667, 0.4       , 0.75      , 1.        ]),
 array([0.4, 0.8, 0.6, 0.6]),
 array([0.5       , 0.53333333, 0.66666667, 0.75      ]),
 array([5, 5, 5, 5])) # 結果は見やすいよう整形しています

　ここでは返り値はtupleになっています。このtupleの中身は、先頭から順に適合率（precision）、再現率（recall）、F1値（fscore）、クラスごとのデータの件数（support）を各分類クラスごとに計算して格納したnumpy配列となっています。つまり「クラスごとの評価指標」がこれで得られたことになります。

　この関数の重要な引数としてaverageがあります。これを指定することで、理論の解説編で説明したマクロ平均、重み付き平均、マイクロ平均を計算できます*9。

>>> precision_recall_fscore_support(true, pred, average="macro")
(0.7041666666666666, 0.6000000000000001, 0.6124999999999999, None)
>>> precision_recall_fscore_support(true, pred, average="weighted")
(0.7041666666666666, 0.6, 0.6124999999999999, None)
>>> precision_recall_fscore_support(true, pred, average="micro")
(0.6, 0.6, 0.6, None)

　つまり何も悩まなくても結果は出してくれるということで、非常に心強いです。ただし、内部でどんなことをしているのかは使う人間がしっかり把握する必要があります。何もソースを読め等と言うつもりはありませんが、解説編で書いた程度の内容は理解して使うべきものです。

　また、多クラス分類の結果を簡単に見たいときは、classification_reportも有用です。

>>> from sklearn.metrics import classification_report
>>> print(classification_report(true, pred))
              precision    recall  f1-score   support

           0       0.67      0.40      0.50         5
           1       0.40      0.80      0.53         5
           2       0.75      0.60      0.67         5
           3       1.00      0.60      0.75         5

    accuracy                           0.60        20
   macro avg       0.70      0.60      0.61        20
weighted avg       0.70      0.60      0.61        20

　このように見やすく結果をまとめてくれます。また、クラスラベルを表示用に当てたり、結果を辞書で返したりといった機能もサポートされているので、割と役に立ちます。ざっくり全体の傾向を把握したいときにおすすめです。

　参考：
sklearnのclassification_reportで多クラス分類の結果を簡単に見る - 静かなる名辞

まとめ

　軽い気持ちで書き始めた記事ですが、ちゃんと説明しようとすると意外と内容が膨れ上がってしまい、長い記事になってしまいました。評価指標はそれだけ濃いテーマだ、ということを実感しました。

　結果の評価というのは非常に大切な部分で、この数字（評価指標）を良くすることを目標に頑張るシチュエーションというのはとても多いです。なので、その基本的な仕組みはしっかりと理解しておくべきでしょう。

　なお、間違いの指摘や質問などは大歓迎です。遠慮なくコメントして頂けると幸いです。

*1:F値と表記されることも多いですが、この記事では敢えてF1値で通します

*2:ところで、PositiveとNegativeというのは、日本語に訳すと陽性、陰性になります。インフルエンザ検査のときのあの陽性、陰性です。

*3:改善前と改善後で統計的有意差があることを検定で示せば良い。ということは、当然データを入れ替えたりして複数の試行を行う必要がある。余談だが、この改善がせいぜい数%の差しかなかったりすると検定のために何十回も試行してデータを取る必要があり、処理が重かったりするとけっこう大変。私が文書分類のテーマで卒論を書いたときは二日二晩パソコンをぶん回したが、それでも軽い方という業界もあるくらいだと思う（ディープラーニングで交差検証を10回やるのに二週間かかるとか・・・そういう業界はどうやってるんだろう？）

*4:最悪比較対象にする数字で評価方法が一貫していればなんとかなる。更にクラス数に極端な偏りがなければどの方法で計算しても問題は起きづらい

*5:有利な指標の計算方法を恣意的に選んでないか？　とか問題になる可能性がある

*6:書く場合は、評価指標だけ示すのではなく混同行列も載せるようにする。読むときは混同行列もじっくり見て、著者の導いている結論が妥当かどうかを確認する

*7:日本語ではマイクロ平均と読み、ミクロ平均と呼ばれることはありません。ドイツ語読みだとミクロで学術用語ではこちらを使うことが多い（日本の学問の多くは明治期以降にドイツから入ってきた）気がするんですが、統計とかデータサイエンスは英語圏から入ってきたので、英語読みが定着してるのか・・・

*8:数式で確認すれば良いんだろうけどやってない。もし間違ってたらごめんなさい

*9:なお、precision_scoreなどの関数も同様にaverageを指定して使うことができます。まあ、普通はprecision_recall_fscore_supportを使えば事足りるでしょう。

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

python sklearn 機械学習 評価指標 Tips SVM 主成分分析

はじめに

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

　そういった情報を簡単に要約して出力してくれるのが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で二次元にしたもの
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の木の本数を増やすとどうなるか？

python sklearn 機械学習 ランダムフォレスト SVM 主成分分析

はじめに

　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

結果

f:id:hayataka2049:20180308194151p:plain
f:id:hayataka2049:20180308194158p:plain
f:id:hayataka2049:20180308194207p:plain
f:id:hayataka2049:20180308194223p:plain
f:id:hayataka2049:20180308194232p:plain
f:id:hayataka2049:20180308194241p:plain
f:id:hayataka2049:20180308194248p:plain

　そこそこ良くなるが、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だとどんな画像が得られるだろうか？
f:id:hayataka2049:20180308194626p:plain
　これは勝てない。RandomForestだとどうしてもカクカクが残るのに。

考察

　この結果の妥当性は率直に言って判断しづらい。

　そもそも、2次元データを入力している以上、ランダムフォレストはデフォルトで（ \sqrt{2}）の特徴量を使って木を作ってくれている訳で、つまり1次元だけで判断してくれている。ちょっとあんまりなので、max_features=2としたのが次の画像。データが（ランダム絡みで）変わってるのでそこだけ注意。
f:id:hayataka2049:20180308200130p:plain
　まあ、これを見てもSVMみたいに滑らかな決定境界が引けてるとは言い難いものがあるけど・・・（考えようによっては1次元でやった上の画像の方が汎化性能は高い、ような気もしてこなくはない。全体的にもやっとしてて、相対的に滑らかに見える）。

　でも、SVMもSVMで、パラメタ次第ではどんな複雑怪奇な決定境界だって引けるといえば引ける。
f:id:hayataka2049:20180308200634p:plain

　こういう問題（けっきょく汎化性能が得られるかどうかはパラメタ次第）があるので、SVMの方が良いと一概に言えるかはけっこう微妙。

　もっと言えば、「それぞれの軸に意味がある」「ノイズもけっこう混ざってる」「スパース」「高次元」という性質のデータを対象とする場合、SVMの汎化性能（滑らかな決定境界を引ける）はかえって邪魔になるのではないだろうか？　そのようなデータでは、とにかく重要な軸を見つけてきて、そこで判断するRandomForestの方が良い性能が得られることが多いと経験的に感じる。ちなみに自然言語処理で使うBoW(Bag of Words)はその典型例である。

　逆に、軸に意味がなくて相対的に低次元でデンスな空間を相手にする場合、SVMの方が良い結果を産むということもよく経験することである。PCAで低次元に落としてしまったデータとか、word2vecで生成される単語の空間とかが割とそんな感じである。

　なんだか話が脱線してきたのでこれくらいにするけど、けっきょく「滑らかな決定境界を引く能力はどうやってもSVMの方が高い（あたりまえ）」「滑らかだから良いというものでもない」「使い分けが重要」という当たり障りのない結論に落ち着いてしまった。

　あと、木の本数は無尽蔵に増やすわけにはいかない。ランダムフォレストは計算量は軽いけど意外とリソース消費の激しいアルゴリズムで、増やしすぎると効率が悪化する。

ランダムフォレストはサンプル数が多いとメモリ消費量が大きい - 静かなる名辞

　汎化性能はできるだけ他の方法で確保したいところ。

まとめ

　SVMの方が滑らかでした（小並感）。

【python】numpyで楕円形のデータを生成する

python numpy

　楕円形のデータを生成したいと思った。

理論

　楕円の式は、基本的に次の定義でいく。

　 x=a\cos (\theta)\\y=b\sin (\theta)
　 0 < \theta < 2\piで、a,\ bが楕円のパラメタである（長半径と短半径）。

　よって、適当なスパンで0～2\piのデータを生成して定義通りxy座標を求めれば良い。簡単そう。

　データを生成する観点からは色々な楕円がほしいので、ちょっと余計な真似をする。

  • 回転角を指定して回転できるようにする

　線形代数の基本。回転行列をかけるだけ。参考：回転行列 - Wikipedia

  • 中心を指定できるようにする

　xyそれぞれの座標に足すだけ。回転させた後に足さないと変なところに飛ばされるので、そこだけ注意する

  • ノイズを乗せる

　理論的には二次元正規分布にするべき・・・とか考えると面倒くさいので、とりあえず各点のXY座標にそれぞれ独立なノイズを乗せて誤魔化すことにする。手抜き。二次元正規分布にするのも然るべき分散共分散行列を与えればそう難しくないはずだけど、やらない（必要な人は自分でやってみて）

実装

# coding: UTF-8

import numpy as np
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 main():
    plt.figure()
    X, Y = make_circle(a=3, b=2, xy=(0,0), phi=0, n=100, random_s=0.2)
    plt.scatter(X, Y, c="r")
    X, Y = make_circle(a=10, b=4, xy=(2,2), phi=45, n=100, random_s=0.5)
    plt.scatter(X, Y, c="g")
    X, Y = make_circle(a=20, b=2, xy=(-3,-5), phi=-70, n=100, random_s=0.5)
    plt.scatter(X, Y, c="b")
    plt.savefig("circle.png")

if __name__ == "__main__":
    main()

　結果は、
f:id:hayataka2049:20180307220454p:plain
　こんな感じで、とりあえず不満のないものが得られました。

　ただ、媒介変数（角度）の変化＝一定とし、楕円の円弧の長さの変化＝一定としなかったので、割と楕円の長半径近くに点がまとまる。これを改善するのは大変そうだけど、やってみるのも面白そうだ。