Pythonで等高線プロット(matplotlib.pyplot.contour)

こんにちは.
Pythonでmatplotlib.pyplot.contourを使った等高線プロットです.
ここではf(x, y)=0を描くことで円を表現してみました.
hこれは分類問題の境界を描く時とかに使えると思います.
単にcontourのlevelsパラメータを変更すれば良いみたいです.

# ~*- coding: utf-8 -*-
"""
等高線プロットのテスト
"""

import numpy as np
import matplotlib.pyplot as plt

n = 100
x = np.linspace(-3, 3, n)
y = np.linspace(-3, 3, n)
X, Y = np.meshgrid(x, y)

# 原点中心,半径1の円
Z = X**2 + Y**2 - 1
plt.contour(X, Y, Z, levels=[0])
plt.show()

# 中心(2, 1),半径1の円
Z = (X-2)**2 + (Y-1)**2 - 1
plt.contour(X, Y, Z, levels=[0])
plt.show()

# 関数で表現
def circle(x, y):
    return (x-2)**2 + (y-1)**2 - 1
Z = circle(X, Y)
plt.contour(X, Y, Z, levels=[0])
plt.show()

# パラメータR,cを与えることのできる関数から作成
# デフォルトは中心(0,0)半径1の円
def circle(x, y, c=[0, 0], R=1):
    return (x-c[0])**2 + (y-c[1])**2 - R**2
Z = circle(X, Y, c=[2, 1], R=2)
plt.contour(X, Y, Z, levels=[0])
plt.show()

広告

微分せずに微分する

数学を勉強していたら、微分は記号処理のように感じられて仕方ない。
違うんだ、ということを意識するためにコードを書いてみた。
非常に素朴な実装だ。
2階微分までは行けそうだ。つまりニュートン法にも使えるということか。

Source

import matplotlib.pyplot as plt
import numpy as np

def f(x):
    return x**3

def gradient(f, h=1./10**2):
    def fx(x):
        return (f(x + h) - f(x)) / h
    return fx

if __name__ == '__main__':
    x = np.linspace(-1, 1, 100)

    fig = plt.figure(1)
    #1st
    plt.subplot(2, 2, 1)
    plt.plot(x, f(x), "b-")

    #2nd
    fx = gradient(f)
    plt.subplot(2, 2, 2)
    plt.plot(x, fx(x), "b-")

    #3rd
    fxx = gradient(fx)
    plt.subplot(2, 2, 3)
    plt.plot(x, fxx(x), "b-")

    #4th
    fxxx = gradient(fxx)
    plt.subplot(2, 2, 4)
    plt.plot(x, fxxx(x), "b-")
    fig.savefig("floattest.png")


Result

file:///Users/hiroyukitanaka/Dropbox/py/thinkstatdata/floattest.png

並列計算(python multiprocessing)

並列計算を勉強しようと思いまして、pythonのmultiprocessingを触っております。 ただ、あまり並列計算用にコードを書き換えるということをしたくないので、何か簡単な方法はないかと探しておりました。 そんな時に目についたのがこの記事。

Poolのmapメソッドを使えば、簡単に並列計算ができるらしいが、複数の引数を取る関数(例えば以下)に適用しようとするとエラーが出る。 そこで、これを解決しようという内容だ。

def myfunc(a, b):
    return a*b

これはリスト内包表記とどちらが早いのか、やってみた。 コードは以下

def argwrapper(args):
    return args[0](*args[1:])

def myfunc(a, b):
    return a*b

if __name__ == '__main__':
    from time import time

    i = 100
    j = 100

    #usual way
    start = time()
    results = [myfunc(a, b) for a in xrange(1, i) for b in xrange(1, j)]
    end = time()
    print "comprehension:", end-start

    #multi processing

    from multiprocessing import Pool
    p = Pool(8)
    func_args = [(myfunc, a, b) for a in xrange(1, i) for b in xrange(1, j)]
    start = time()
    results = p.map(argwrapper, func_args)
    end = time()
    print "multiprocessing:", end-start

結果が以下

comprehension:   0.00245380401611
multiprocessing: 0.00831413269043

おいおい、普通に内包表記の方が早いではないか…
しかも引数リストを作る時間まで合わせるととても使える代物ではない…
どうしたことか。処理に時間のかかるメソッドなら早いのかな?
というわけで以下のメソッドで再計測

def myfunc(a, b):
    for i in xrange(100000):
        pass

結果は以下

comprehension: 18.3397388458
multiprocessing: 5.15106010437

なんと!3分の一(の純情な感情)やないか!これは驚き。
処理の思いメソッドなら、並列計算の恩恵を受けることができるんですね。
いやぁ久々に震えました笑。どんどん並列化しよう。
今回はここまで。

tfidfでお手軽キーワード抽出(2)

前回:tfidfでキーワード抽出
これと同じことをnltk.TextCollectionを使ってやってみました。 非常に簡単に書けました。前回の記事で時間と空間のトレードオフとか言ってましたが、実は全然メモり削減できてませんでしたね…
環境:
OSX 10.8.4
Python 2.7.2

Implementation

# -*- coding: utf-8 -*-
from nltk.corpus import reuters
import nltk
import math

if __name__ == '__main__':

    target_class = "acq"
    N = 200
    docs = [reuters.words(fileid) for fileid in reuters.fileids(categories=target_class)]
    docs = docs[:N]
    doc = docs[0]
    collection = nltk.TextCollection(docs)
    candidates = set(doc)
    for candidate in sorted(candidates, key=lambda x:collection.tf_idf(x, doc), reverse=True)[:20]:
        print "%s: %f" % (candidate, collection.tf_idf(candidate, doc))




Result

Sumitomo: 0.091463
Komatsu: 0.058798
bank: 0.032187
Sogo: 0.026132
Heiwa: 0.026132
business: 0.023753
securities: 0.019674
We: 0.018469
": 0.018430
,": 0.017610
will: 0.017266
But: 0.016395
analysts: 0.014775
Gottardo: 0.013066
cases: 0.013066
network: 0.013066
Smithson: 0.013066
deregulation: 0.013066
quickly: 0.013066
Kleinwort: 0.013066

Summary

nltk.TextCollectionを使えば簡単にtfidfを求めることができることがわかった。
しかし、扱う文書を増やすとToo many open filesというエラーが出るのは解決していない。
というわけでこれを解決したいと思う。もう少しロイターコーパスで遊ぶ。

tfidfでお手軽なキーワード抽出

そんな大それたものではありません。 単語をtfidfで重み付けするという手法はかなりtraditionalな手法だと思いますが、これを利用してキーワード抽出をやってみました。
環境: python 2.7.2
OSX 10.8.4

tfidf

tfidfというのは以下で定義される量です。自分なりに書くのでwikiを観て頂いた方が良いかと。

\displaystyle  tfidf(w, d, D) = tf(w, d) \cdot idf(w, D) \\ tf(w, d) = \frac{t(w, d)}{N(d)} \\ idf(w, D) = \frac{1}{df(w, D)}

ただし、
t(w, d)は文書dにおける単語wの出現回数、
N(d)は文書dの総単語数(異なり語数ではない)、
df(w, D)は文書集合D=\{d_{1}, \ldots, d_{n}\}の内、単語wが出現する文書の数、
df(w, D)が0の時はidf(w, D)=0 (org2blogでの条件式の書き方がわかりませんでした。
とする。
なぜキーワード抽出にtfidfを使うのか。一つの理由にストップワードを明示的に指定しなくても、ストップワードがキーワードになることを防ぐことかできるという理由がある。(自分なり
言い切ってしまった。 例えば”the”だとか”a”だとかは当然tfが高くなるが、ほとんどの文章で出現するのでidfが0になる。

Implementation

目的は
文書dが与えられたとき、そのキーワードを抽出すること
です。 今回はnltkのロイターコーパス(nltk.corpus.reuters)を使います。 中でも”acq”クラスに属する文書の内200件だけ使います。それ以上使うとnltk too many open filesというエラーが発生したので(僕のMacだけでしょうか? コードを以下に示します。 これは非常に効率の悪いコードです。 「時間と空間のトレードオフ」という言葉がありますが、かなり空間を優先しています。つまり時間を犠牲にしています。 何回も文書を走査するので、当然遅いです。でかいコーパスを使うと、こういう問題にも遭遇するのだな、という良い勉強になりました。

# -*- coding: utf-8 -*-
from nltk.corpus import reuters
import math

def tf(word, words, lower=False):
    if lower:
        words = [word.lower() for word in words]
        word = word.lower()
    if words:
        return float(words.count(word)) / len(words)
    return 0.0

def df(word, docs, lower=False):
    if lower:
        docs = [[word.lower() for word in doc] for doc in docs]
        word = word.lower()
    return float(sum([word in doc for doc in docs])) / len(docs)

def idf(word, docs, lower=False):
    if lower:
        docs = [[word.lower() for word in doc] for doc in docs]
        word = word.lower()
    df_val = df(word, docs)
    if df_val == 0.0:
        return 0.0
    return math.log(1.0 / df_val)

def tfidf(word, doc, docs, lower=False):
    return tf(word, doc, lower=lower) * idf(word, docs, lower=lower)

target_class = "acq"
N = 200
docs = [reuters.words(fileid) for fileid in reuters.fileids(categories=target_class)]

docs = docs[:N]
doc = list(docs[0])
doc_set = set(doc)
doc_set = sorted(doc_set, key=lambda x:tfidf(x, doc, docs), reverse=True)

for word in doc_set[:20]:
    print "%s: %f" % (word, tfidf(word, doc, docs, lower=False))

Result

結果は以下のようになりました。(上位20件) 印象はまぁ悪くなさそうです。”とか,”が入っているのが何とも言えませんが笑。

Sumitomo: 0.091463
Komatsu: 0.058798
bank: 0.032187
Sogo: 0.026132
Heiwa: 0.026132
business: 0.023753
securities: 0.019674
We: 0.018469
": 0.018430
,": 0.017610
will: 0.017266
But: 0.016395
analysts: 0.014775
Gottardo: 0.013066
cases: 0.013066
network: 0.013066
Smithson: 0.013066
deregulation: 0.013066
quickly: 0.013066
Kleinwort: 0.013066    

Summary

nltkのロイターコーパスからtfidfを使って、お手軽なキーワード抽出を実装してみた。
最近ロイターコーパス使って遊んでいるので、その一環としてのポスト。 そういえばtfidfを使ったことがあんまりないなぁと思ったので、その復習もかねました。 今回思ったこと
1.でかいコーパスを扱うと、メモリの問題なども出てくる
2.dfの計算に結構時間がかかる
次回はいろいろ思いつくのですが、とりあえず並列計算をやってみたい。 あと、コーパスのあるカテゴリによく出てくる単語を抽出してみたい。 これにはまたtfidfを使うことになるかも。