静かなる名辞

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


【python】collections.ChainMapの使い方を理解する

はじめに

 pythonで複数の辞書をマージするにはどうしたらいいのでしょうか。forループ? 辞書内包表記を使う? updateメソッド? 実は、ChainMapというものもあります。

 その使い方について説明します。

先に結論

 記事本文でだらだらと説明していますが、要約すると、

  • 複数の辞書を一つにまとめる方法
  • 元辞書への参照を張ってごにょごにょするだけなので、新しい辞書を作るより(ケースバイケースだが)速い。ただしマージが速くても参照が遅いという微妙な性質がある
  • 当然元辞書を変更すると反映されてしまう

 これだけです。

基本的な使い方

 リストを渡すのかと思ったら、違うようです。

>>> from collections import ChainMap
>>> d1 = {"hoge":"ほげ"}
>>> d2 = {"fuga":"ふが"}
>>> cm = ChainMap(d1, d2)
>>> cm["hoge"]
'ほげ'
>>> cm["fuga"]
'ふが'

 このように使えます。ただ、辞書のリストを動的に生成して渡したい場合もあるでしょう。その場合、starred expressionを活用して、

>>> cm = ChainMap(*[d1, d2])
>>> cm
ChainMap({'hoge': 'ほげ'}, {'fuga': 'ふが'})

 このようにやることになるかと思います。

参照先を変更してみる

>>> d1["HOGE"] = "ホゲ"
>>> cm
ChainMap({'hoge': 'ほげ', 'HOGE': 'ホゲ'}, {'fuga': 'ふが'})

 このようにけっこう恐ろしい性質があります。これが嫌なときは他の方法で辞書をマージするか、deepcopyを使うことになるはずです。

 参考:複数の辞書のマージ方法いろいろ - Qiita

速いらしいという噂があるので測る

 測ってみました。上記参考サイトの辞書内包表記版と比較します。

 予想では、マージは速いものの探索は恐らくhashのリストを線形探索していく上、重複を取り除く処理などもあるので、遅い要素がありそうです。

 items()と一つのキーへのアクセスで近似的に探索を表現します。

# coding: UTF-8

import time
from collections import ChainMap

def time_measure(f, lst):
    time_lst = []
    for i in range(10000):
        t1 = time.time()
        result = f(lst)
        t2 = time.time()
        time_lst.append(t2-t1)
    return result, sum(time_lst)/10000

def main():
    f1 = lambda dicts:{k: v for dic in dicts for k, v in dic.items()}
    f2 = lambda dicts:ChainMap(*dicts)

    d1 = dict(zip(range(100), range(100)))
    d2 = dict(zip("ho", "ge"))
    d3 = dict(zip(range(100), "p"*100))

    result, time = time_measure(f1, [d1, d2, d3])
    print("辞書内包マージ", time)
    f = lambda _:[(k,v) for k,v in result.items()]
    _, time = time_measure(f, None)
    print("辞書内包items()", time)
    f = lambda _:result["h"]
    _, time = time_measure(f, None)
    print("辞書内包キーアクセス", time)
    result, time = time_measure(f2, [d1, d2, d3])
    print("ChainMapマージ", time)
    f = lambda _:[(k,v) for k,v in result.items()]
    _, time = time_measure(f, None)
    print("ChainMap items()", time)
    f = lambda _:result["h"]
    _, time = time_measure(f, None)
    print("ChainMapキーアクセス", time)

if __name__ == "__main__":
    main()

 結果は

辞書内包マージ 1.2324166297912597e-05
辞書内包items() 7.125544548034668e-06
辞書内包キーアクセス 1.9309520721435546e-07
ChainMapマージ 1.3034582138061524e-06
ChainMap items() 6.01630687713623e-05
ChainMapキーアクセス 1.183009147644043e-06

 マージは若干遅いが探索は速い辞書、マージは速いものの肝心の探索が遅いChainMapという結論。

 探索速度(実タスクでの探索速度ではないが)に1桁の差が出てる訳で、安直に「ChainMap使おう」とはならないと思う・・・。普通はアクセスの速さを享受したいから辞書にすると思うのだが。

まとめ

 使い所は難しいと思いました。

【python】外部プロセスと標準入出力で通信する

 pythonで外部プロセス(subprocess)と、標準入出力を介したやりとりをしたいときがある。

 目次

やってみること

 今回は問題例として、形態素解析器MeCabをpythonから呼んでみる。MeCabはpythonバインディングがあるので、実用的には外部プロセスとして起動する必要はないのだが、わかりやすい例として敢えて挙げる。

 なお、MeCabを知らない方のために説明すると、自然言語処理で使う形態素解析器というもので、次のような使い方をするものである。

$ mecab
吾輩は猫である。 # 入力行。この下の行から解析結果
吾輩	名詞,代名詞,一般,*,*,*,吾輩,ワガハイ,ワガハイ
は	助詞,係助詞,*,*,*,*,は,ハ,ワ
猫	名詞,一般,*,*,*,*,猫,ネコ,ネコ
で	助動詞,*,*,*,特殊・ダ,連用形,だ,デ,デ
ある	助動詞,*,*,*,五段・ラ行アル,基本形,ある,アル,アル
。	記号,句点,*,*,*,*,。,。,。
EOS

 
 pythonスクリプトから外部プロセスとしてMeCabを起動し、入力をpythonスクリプトから送り込み、解析結果を得ることを目標にする。

簡単な方法

とりあえずやる

 標準ライブラリのsubprocessを用いる。次のように書く。

# coding: UTF-8

import subprocess

def main():
    input_string = "吾輩は猫である。"

    mecab_p = subprocess.Popen(["mecab"], stdin=subprocess.PIPE,
                               stdout=subprocess.PIPE)
    output_string = mecab_p.communicate(input_string.encode())[0].decode()
    print(output_string)

if __name__ == "__main__":
    main()

 結果

吾輩	名詞,代名詞,一般,*,*,*,吾輩,ワガハイ,ワガハイ
は	助詞,係助詞,*,*,*,*,は,ハ,ワ
猫	名詞,一般,*,*,*,*,猫,ネコ,ネコ
で	助動詞,*,*,*,特殊・ダ,連用形,だ,デ,デ
ある	助動詞,*,*,*,五段・ラ行アル,基本形,ある,アル,アル
。	記号,句点,*,*,*,*,。,。,。
EOS

 できた。

解説

 subprocess.Popenでプロセスを立ち上げパイプをつなぐ。後はmecab_p.communicateの引数で入力を、返り値で出力を得る。なお、返り値は

 注意すべき点としては、端末から打ち込むのと基本的に同じなので文字コードを意識する必要がある。python内部で使われている文字列型は、(python3では)unicode文字列(str型)という「とりあえず文字コードを意識しなくても使える型」になっているのだが、外部プロセスに投げるときはencodeしてbytes型という要するにただのバイト列にしてやる必要があり、また外部プロセスから受け取った文字列もまたbytes型なのでdecodeしてstr型に戻してやらないといけない。

 今回はlinux上の環境で実験しているので、文字コードはpython、システム、mecabすべてutf-8で統一されており、その場合は上のコードのように簡単に書ける。

 windows環境だったりすると(けっきょくはpythonとmecabの設定に依存するが)、encodingが統一されていないこともあり得るので、そういうときはちゃんと自分で環境を設定した通りにencodingを取り扱ってやらないと、たぶん動かない。

応用:複数回入出力を送りつけたい

 subprocess.Popen.communicateを2回呼ぶと、次のようなエラーが出る。

ValueError: Cannot send input after starting communication

 つまり1回しかやり取りできない。送りつけて出力を受け取った瞬間、パイプは閉じられてしまう。安全といえば安全だが、これは不便である。

 この回避のためには、

  • 複数回の入力を一つにまとめて送りつける
  • stdin, stdoutを直接叩く

 といった方法がある、らしい。

 参考:python - Communicate multiple times with a process without breaking the pipe? - Stack Overflow


 前者の方法はそれで済めば良いだろうけど、たまには対応しきれないときもある。後者の方法は、そんな見るからにおっかない(制御ミスるとすぐ死ねそうな)方法は使いたくない。というか、実際やってみても上手く動かせなかった。

 なので、Pexpectというライブラリに頼る。これは擬似端末を立ち上げてプロセスとおしゃべりする系のもの。コマンドの自動化などに使うものらしい。

 結論から言うと、これを使うのはけっこう辛かった。挙動のクセが強すぎる。正規表現に慣れててexpectの世界観も理解してる人には良いんだろうけど、僕みたいな素人にはただただ望む出力を得るまで試行錯誤を重ねる苦行だった。

def main2():
    inputs = ["吾輩は猫である。", "名前はまだない。"]

    p = pexpect.spawn("/bin/bash")
    #p.logfile_read = sys.stdout.buffer

    p.sendline("mecab")
    p.readline()
    p.sendline("")
    p.expect("EOS")

    results = []
    for s in inputs:
        p.sendline(s)
        p.expect("EOS")
        results.append(p.before.decode())

    p.terminate()
    p.expect(pexpect.EOF)

    for r in results:
        print(r)

 うまい正規表現が思いつかなかったのでMeCabのEOSを取ることに。

 なんだか微妙な結果になってしまったけど、Pexpect自体はたぶん良いライブラリなので、然るべき場面で活用してあげてください。標準入出力を自動化したいときはたぶん使えます。ただ、今回の用途にはいまいち合ってなかった気がする。

まとめ

 パイプで投げつけるだけだろと思ってたら、けっこう色々なことに気を遣う羽目になりました。

【python】nltkで英語のStemmingとLemmatization

 Stemming(ステミング)は単語の語幹を取り出したいとき、Lemmatization(レンマ化、敢えてカタカナ表記するとレンマタイゼーション)はカテゴリごとにグルーピングしたりしたいときに使う。

 公式ドキュメントはここ。
nltk.stem package — NLTK 3.4 documentation

 目次

スポンサーリンク



Stemming

概要

 nltkでStemmingに使えるクラスはたくさんある(ように見える)。nltk.stemに実装されているものだけでも、

  • ARLSTem Arabic Stemmer*1

 アラビア語用

  • ISRI Arabic Stemmer*2

 アラビア語用

  • Lancaster Stemmer*3

 英語用。古い(1990年)

  • Porter Stemmer*4

 英語用。古い(1980年) 

  • Regexp Stemmer

 正規表現で実装されている。どんな正規表現にするかは使うときに指定する

  • RSLP Stemmer

 ポルトガル語用

  • Snowball Stemmers

 なんか色々な言語に対応している

 これだけある。注釈はぜんぶドキュメントからコピペした出典です。

 とりあえず英語で何も考えずに使えるのは、Lancaster、Porter、Snowball Stemmersの3択。LancasterとPorterは基本的にルールベース、Snowball StemmersはPorter StemmerのPorterさんの作ったフレームワークで、多言語対応していることが特徴なのだけど肝心の英語版のアルゴリズムはPorterらしいので(ドキュメントにそう書いてある)、選択肢から外すと、けっきょくLancaster、Porterの2択になる。

 なお、Stemmingの研究がどこまで進んでいるかはこのドキュメントを見て勉強した。
http://kenbenoit.net/assets/courses/tcd2014qta/readings/Jivani_ijcta2011020632.pdf

 統計的言語処理に走って、けっきょくこれは難しいタスクというか、一単語ずつ見てもどうしようもないよね、という壁にぶち当たった・・・って感じか。今はどんな感じになってるんだろうか(というかstemming自体、今はどこまで重要視されているのだろうか)。一応、nltkに実装されているアルゴリズムも、ある程度の完成度には達しているのだろうけど。

 nltk.stemで使えるアルゴリズムはそういった難しい問題が立ちはだかる以前のものなので、とりあえず気楽に使えると言えば使えそう。n-gramとかHMMとかを応用したStemmerも提案はされてるっぽいのだが、nltkで使えなさそうなので没。

Porterを使う

 とりあえずPorterから先に使ってみる。

>>> from nltk.stem.porter import PorterStemmer as PS
>>> ps = PS()
>>> ps.stem("dogs")
'dog'
>>> ps.stem("cats")
'cat'
>>> ps.stem("ceiling")
'ceil'

 これだけ。最後のceilingは「eg. ceil- is not the stem of ceiling」という例なので、解析がいまいち上手く行っていないということになる(この例はnltk.stemの公式ドキュメントの最後に出ている)。

Lancasterを使ってみる

>>> from nltk.stem.lancaster import LancasterStemmer as LS
>>> ls = LS()
>>> ls.stem("dogs")
'dog'
>>> ls.stem("cats")
'cat'
>>> ls.stem("ceiling")
'ceil'

 まったく同様に使えるようだ。

Lemmatizing

 nltk.stemのページの一番最後にしれっと載ってる。wordNetを使って単語をカテゴリに落としてくれる。一体どんな粒度でやってくれるのかよくわからないのだが、公式の例を見るとstemmingと同様に機能している気がする。

 なお、初回呼び出し時は親切なことに、

LookupError: 
**********************************************************************
  Resource wordnet not found.
  Please use the NLTK Downloader to obtain the resource:

  >>> import nltk
  >>> nltk.download('wordnet')

 こんなエラーを吐いてくれる。downloadすると使えるようになる。

>>> from nltk.stem.wordnet import WordNetLemmatizer as WNL
>>> wnl = WNL()
>>> wnl.lemmatize("dogs")
'dog'
>>> wnl.lemmatize("cats")
'cat'
>>> wnl.lemmatize("ceiling")
'ceiling'

 Lemmatizer特有のメリットとしては、

>>> wnl.lemmatize("good", pos="a")
'good'
>>> wnl.lemmatize("better", pos="a")
'good'

 こういうことができる。この例はpos="a"(adjective)を指定しないと名詞として取り扱われて上手く行かないみたいだけど。このposキーワードについては、そもそもPOSタグ判定が成功しないと上手く使えない上、POS taggerの吐くPOSタグと指定できるPOSタグに互換性がなさそうなので、使い方はかなり難しそう。無理を重ねている感じ。

 なお、公式曰く、「Returns the input word unchanged if it cannot be found in WordNet.」とのこと。つまり「元々語幹の形で入力されたのでchangeする必要がなかった」のか「見つからないからとりあえずそのまま返した」だけなのか区別がつかない。ふざけんなと言いたい。せめて検索して見つかったかどうかくらいはT/Fでも1/0でも良いから返してくれ。

 ま、そういうことになっているので、実際どう使うべきなのかは率直に言って悩む。

結論

 研究者・開発者の皆さんには申し訳ないけど、どれもどんぐりの背比べ感が否めない・・・。

*1:K. Abainia, S. Ouamour and H. Sayoud, A Novel Robust Arabic Light Stemmer , Journal of Experimental & Theoretical Artificial Intelligence (JETAI‘17), Vol. 29, No. 3, 2017, pp. 557-573.

*2:Taghva, K., Elkoury, R., and Coombs, J. 2005. Arabic Stemming without a root dictionary. Information Science Research Institute. University of Nevada, Las Vegas, USA.

*3:Paice, Chris D. “Another Stemmer.” ACM SIGIR Forum 24.3 (1990): 56-61.

*4:Porter, M. “An algorithm for suffix stripping.” Program 14.3 (1980): 130-137.

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

 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次元描画
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
回るiris

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

【python】flymakeのエラー

 次のようなエラーを吐かれた。

Flymake: Configuration error has occurred while running(***/flymakes ***.py).Flymake will be switched OFF

原因

 色々あるらしいけど、今回は文字コード指定をタイポした瞬間エラーが出た。

# coding: UTF-9

 たぶん、これを書いた瞬間にflymakeが走って「そんな文字コードないよ!解釈できない!」→落ちる、という流れだったんだと思う。

対策

 これも色んなサイトに書いてあることだが、次の設定をemacsの設定ファイルに付け足す。

(defadvice flymake-post-syntax-check (before flymake-force-check-was-interrupted)
  (setq flymake-check-was-interrupted t))
(ad-activate 'flymake-post-syntax-check)

 参考にしたページ:
flymakeのsyntax-checkが異常終了しても無視するようにする - すぎゃーんメモ

 とりあえずflymakeは落ちなくなったが、UTF-9とか書いちゃってもエラー表示してくれる訳ではないので、微妙っちゃ微妙。でもそれに対応させる気概もないのでこれで行く。

【python】関数オブジェクトは辞書のキーに使える!

 「まさかできねーだろw」と思ってやったらできたのでびっくりしたよ!

>>> def hoge():
...     return None
... 
>>> {hoge:1}
{<function hoge at 0x7f025f18cf28>: 1}
||<

  えぇ…。

 lambdaでもできる。

>|python|
>>> fuga = lambda:None
>>> {fuga:2}
{<function <lambda> at 0x7f0243880510>: 2}

 hogeとfugaは関数として同値である(敢えてそう言えば)。ということは・・・!?

>>> hoge in {hoge:1}
True
>>> fuga in {hoge:1}
False

 よかった、ちょっとホッとした。

 ん、待てよ・・・

>>> def piyo():
...     return None
... 
>>> piyo in {hoge:1}
False

 こっちもちゃんとFalseだった。よかった。

どうしてこうなるのか。

 こちらの記事を読んだ。
Python における hashable - Qiita
 __hash__メソッドがありさえすれば辞書のキーに使えるらしい。試してみよう。

>>> hoge.__hash__()
-9223363308844643086
>>> fuga.__hash__()
8727981228113

 あるんだ・・・。

>>> hoge.hogehoge = "hogee^~"
>>> hoge.__hash__()
-9223363308844643086

 定義時に決まるimmutableだと思う。まあ、それは予想できてた。hashが更新できたら大騒ぎだ(関数の処理内容に応じて決まるのが本来のhashではないか? という気がしなくはないが。でも関数オブジェクトの処理内容は後から変更できる訳ではないので(高階関数で渡す等は除く)別に構わないのだろう)。

悪用方法

 国際難読化pythonコードコンテストがあれば使えそう(ないけど)。たとえば、15秒くらいで思いついたちょっと格好良いコードはこんな感じ。

>>> d = None
>>> d = {lambda :d:d}

 lambdaのコロンとdictのコロンが混ざって非常に見た目が素敵。こんなのでもちゃんと構文解析されるので、pythonは偉い。

 それ以外の使い道・・・keyを呼び出せるというのがすべてですが、それが活きそうな場面は思いつきそうで思いつかない。微妙な感じです。

 30秒ほど悩んだところで「ヘタに使い道があると使うバカが出てきて危険」ということに思い至り、それ以上考えることをやめました。

【python】cabochaのpythonバインディングの変な挙動

環境

 ubuntu 14.04
 cabocha 0.69
 cabocha-python 0.69

問題の概要

 変な挙動だった。というか率直に言ってバグなのでは?

>>> import CaboCha
>>> cparser = CaboCha.Parser()
>>> tree1 = cparser.parse("吾輩は猫である。")
>>> print(tree1.toString(CaboCha.FORMAT_TREE))
    吾輩は-D
  猫である。
EOS

>>> tree2 = cparser.parse("吾輩は猫ではない。")
>>> print(tree1.toString(CaboCha.FORMAT_TREE))
      吾輩は-D
  猫ではない。
EOS

 これはおかしい。CaboChaはこわれている。

 いや、「これで仕様通り動いてる。おまえの使い方が間違ってるんだ」て言われたら反論できないけど。詳細なドキュメントを見かけたことがないので、もしかしたらアホなこと(Parserのインスタンスの使い回し)をやっているのかもしれない。

回避するために試したこと

  • CaboCha.Parser("")する(コンストラクタに空文字列を渡す)

 効果なし。

  • 一度空文字列に対してparseを呼ぶ

 ここを参考に「もしかしたら効くかも」と思ってやってみた。
MeCabのparseToNodeのひどいバグ - 北野坂備忘録
 効果なし。

  • 文字列を変数に入れる、encodeする

 ここを参考に(以下略)。
MeCabをPythonから使う注意点とか
 効果なし。encodeに至ってはやったら落ちた。

  • 仕方がないのでCaboCha.Parser()を毎回作る
>>> tree = CaboCha.Parser().parse("吾輩は猫である。")
>>> tree.toString(CaboCha.FORMAT_TREE)
Segmentation fault (コアダンプ)

 たぶん本体のメモリ管理とpythonの接合が上手く行っていないのだろうけど、さて困った。

>>> parser1 = CaboCha.Parser()
>>> parser2 = CaboCha.Parser()
>>> tree1 = parser1.parse("吾輩は猫である。")
>>> tree2 = parser2.parse("吾輩は猫ではない。")
>>> tree1.toString(CaboCha.FORMAT_TREE)
'    吾輩は-D\n  猫である。\nEOS\n'

 一応回避できることはわかった。これで書くと極めて非python的なプログラミングを強いられるという問題はあるが、たぶんなんとかなる。

 ちなみに、ParserのインスタンスがGCに回収されると treeだけ残っててもtoStringできないようです(Segmentation faultを吐いてくれる)。どんな作りになってるのかなんとなくわかってきたけど、率直に言って○○。

def parse(sentences):
    """
    sentencesは一文ずつに区切られた文のリストとして扱う
    """
    trees = []
    plist = []
    for s in sentences:
        parser = CaboCha.Parser()
        trees.append(parser.parse(s))
        plist.append(parser)

 このようなものを書いてあげれば、後からtreeを使うことができることがわかった。率直に言ってまったく嬉しくない。

問題原因の切り分け

 は、できてないです。
 

  • うちの環境固有の問題
  • cabocha-pythonの特定のバージョンの問題
  • cabocha-python固有の問題
  • cabocha固有の問題

 とりあえず逃げれることはわかったので、僕はやらない(明言)。

対策

 たぶん解析結果のtreeオブジェクトを使いまわそうという発想が間違っていて、cabochaのtreeオブジェクトを使わないで初手でXMLか何かに変換して取り扱うのが楽だと思います。そんなことを強いるバインディングって何よ? って気がしますが。

 もう面倒くさいから、JUMAN/KNPに鞍替えしようかなと思う今日この頃。

【python】区切り文字を含めてsplitする

 
 正規表現によるsplitで区切り文字(あるいは区切り文字列)を含めたいときがある。デフォルトでは区切り文字は消える。

 たとえば、文を句点で分割する場合。

>>> import re
>>> string = "吾輩は猫である。名前はまだない。"
>>> re.split("。", string)
['吾輩は猫である', '名前はまだない', '']

 実は「区切り文字(列)をリストの要素にする」のか「区切り文字(列)を前の要素に(あるいは後の要素に)くっつける」のかという問題もあるので、一筋縄ではいかない。検索すると前の方ばっかり出てくるので注意。

注意

 この記事の内容は一部古いです。おそらく大半の読者が求めているであろうものは追記4にあります。そちらをご覧ください。

区切り文字をリストの要素に入れる

 これは検索してたくさん出てくる方。正規表現のグループ化というものを使えば簡単にできる。

>>> re.split("(。)", string) # 半角丸括弧でくくる
['吾輩は猫である', '。', '名前はまだない', '。', '']

 今回の例では、それが欲しいんじゃないんだけどなぁ・・・と相成る。

# こうなってほしい(空文字列はご愛嬌)
['吾輩は猫である。', '名前はまだない。', '']

区切り文字を前の要素に付ける。

 できるだけ簡単なやり方を考える。

後処理で連結する

 真っ先に思いつくのはこの方法。forループで処理する。

>>> lst = []
>>> for elem in re.split("(。)", string):
...     if elem != "。":
...         lst.append(elem)
...     else:
...         lst[-1] = lst[-1] + elem
... 
>>> lst
['吾輩は猫である。', '名前はまだない。', '']

 書き方は他にも幾らでもある。この書き方だと、区切り文字(列)が連続したときは連続が途切れるまで連結される。

 ただ、これは遅いかもしれない。

正規表現で書く

 rubyでやってるページは見つけた。
[Ruby]後読みを使ったsplitで句点(。)を残したまま文に分割 - Qiita

 ここの通りにやろうと思ったら、できませんでした・・・。

 rubyの正規表現のsplitとは仕様が違うのかなぁ。たぶん『長さゼロの文字列』へのマッチを何らかの方法で書けば同様にできると思われる。どうやるんだろうか、というかできるんだろうか。詳しい人に教えてほしい。

まとめ

 正規表現とかよくわからないけど、それに学習コスト投入するのは面倒くさいので騙し騙し適当に処理してプログラム書くことにします・・・。

追記

 findallすれば等価のことができた。

>>> re.findall(".*?。", s)
['吾輩は猫である。', '名前はまだない。']

 ただし厳密に言えばこれはsplitではない。なのでこうなる。

>>> s2 = "吾輩は猫である。名前はまだない。どこで生まれたかとんと見当がつかぬ"
>>> re.findall(".*?。", s2)
['吾輩は猫である。', '名前はまだない。']

 邪道だと思うが、

>>> s2 = "吾輩は猫である。名前はまだない。どこで生まれたかとんと見当がつかぬ"
>>> re.findall(".*?。|.*$", s2)
['吾輩は猫である。', '名前はまだない。', 'どこで生まれたかとんと見当がつかぬ', '']

 こんな方向性で工夫すればなんとかなる可能性はある。正規表現が得意な人なら上手く書く方法を思いつくんだろうなぁ。

追記2

 コメントでこういう場合の書き方を教えていただきました。

>>> s2 = "吾輩は猫である。名前はまだない。どこで生まれたかとんと見当がつかぬ"
>>> import re
>>> re.findall("[^。]+。?",s2)
['吾輩は猫である。', '名前はまだない。', 'どこで生まれたかとんと見当がつかぬ']

 やっぱりsplitとは違う気がするけど、実用的には良さそうです。

追記3

 ※この追記3は以前のバージョンの記述に基づいた内容です。バージョン3.7以降で事情が変わっているので、追記4をご参照ください。

 ドキュメントをつらつらと眺めていたら、関連する記述を見つけたのでメモ。

現在、split() は空のパターンマッチでは文字列を分割しません。例えば、次のようになります:

>>> re.split('x*', 'axbc')
['a', 'bc']

'x*' は 'a' の前、 'b' と 'c' との間、 'c' の後の 0 個の 'x' にもマッチしますが、現在これらのマッチは無視されます。正しい動作 (空のマッチでも文字列を分割し、['', 'a', 'b', 'c', ''] を返す) は、Python の将来のバージョンで実装されます。これは、後方互換生のない変更であるため、移行期間中は FutureWarning が送出されます。

空の文字列のみとマッチするパターンは、現在文字列を全く分割しません。これは望ましい動作ではないため、Python 3.5 から ValueError が送出されます:

>>> re.split("^$", "foo\n\nbar\n", flags=re.M)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  ...
ValueError: split() requires a non-empty pattern match.

re --- 正規表現操作 — Python 3.7.4rc1 ドキュメント

 将来的にはできるようになるけど、今のpythonは空のマッチはできないよ、と書いてある。

 たとえば、下記の記事のようにすると、

qiita.com

>>> import re
>>> re.split("(?<=。)", 'おはよう。こんにちは。こんばんは。')
ValueError: split() requires a non-empty pattern match.

 こうなるからできない、と。

 ただしその下に記述があるように、findall, finditerは空のマッチを含む。finditerなら位置の情報が得られる。

 なので、こんな感じで表題の通りの動作は実現できる。

>>> s = 'おはよう。こんにちは。こんばんは。'
>>> result = []
>>> before_start = 0
>>> for mobj in re.finditer("(?<=。)", 'おはよう。こんにちは。こんばんは。'):
...     result.append(s[before_start:mobj.start()])
...     before_start = mobj.start()
... 
>>> result
['おはよう。', 'こんにちは。', 'こんばんは。']

 これがさえたやり方かどうかはわからないけど、最初にやりたかったことには一番近いのかな・・・? という気がする。

追記4

 python3.7から空のマッチに対応し、上のような面倒なことが不要になりました。

バージョン 3.7 で変更: 空文字列にマッチしうるパターンでの分割をサポートするようになりました。

re --- 正規表現操作 — Python 3.7.4rc1 ドキュメント

>>> import re
>>> re.split("(?<=。)", 'おはよう。こんにちは。こんばんは。')
['おはよう。', 'こんにちは。', 'こんばんは。', '']

 ここにたどり着くまでは長い道のりでした。ようやく理想的な結果が得られました。

日本語モダリティ解析器 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で英語の形態素解析

概要

 形態素解析、いわゆる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】文字列を一文字ずつのリストにする

 文字列はそもそも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で警告が出るときの対策

概要

 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でデータをランダムサンプリング

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

 目次

スポンサーリンク



実装方針

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

 これは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(線形判別分析)で次元削減


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

 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だったりしたらどうなるんだろう? それとも共分散行列のサイズなのだろうか?