pandasで指数平滑移動平均

昨日の記事が移動平均だったので、今日は指数平滑移動平均を扱います。
初めて知った日は衝撃だったのですが、pandasには指数平滑移動平均を計算する専用の関数が用意されています。
(pythonを使い始める前はExcel VBAでいちいち実装していたので非常にありがたいです。)

馴染みがない人もいると思いますので軽く紹介しておきます。
元のデータを${x_t}$とし、期間$n$に対して指数平滑移動平均${EWMA_t}$は次のように算出されます。
$$
\begin{align}\alpha &= \frac{2}{1+n}\\
EWMA_0 &= x_0\\
EWMA_t &= (1-\alpha)*EWMA_{t-1} + \alpha * x_t
\end{align}
$$

3番目の式を自分自身に逐次的に代入するとわかるのですが、
$EWMA_t$は、$x_t$から次のように算出されます。
$$
EWMA_t = \alpha\sum_{k=0}^{\infty}(1-\alpha)^k x_{t-k}
$$
$(1-\alpha)$の絶対値は1より小さいので、この無限級数の後ろの方の項は無視できるほど小さくなります。
結果的に、過程${x_t}$の最近の値に重みを置いた加重平均と見做せます。

さて、早速ですが計算してみましょう。

pandasのDataFrameおよび、Seriesに定義されているewm関数を使います。
pandas.DataFrame.ewm


import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# データ作成
data = pd.Series(np.random.normal(0, 100, 200).cumsum() + 20000)
# 指数平滑移動平均の計算
data_ewm = data.ewm(span=10).mean()
# 可視化
plt.rcParams["font.size"] = 14
fig = plt.figure(figsize=(12, 7))
ax = fig.add_subplot(1, 1, 1)
ax.plot(data, label="元データ")
ax.plot(data_ewm, label="指数平滑移動平均")
plt.legend()
plt.show()

出力がこちら。

ここで一つ注意する点があります。
data_ewm = data.ewm(span=10).mean()
という風に、spanという変数名で期間$10$を渡しています。
ドキュメントを読んでいただくとわかるのですが、span=をつけないと、
comという別の変数に値が渡され、$\alpha$の計算が、
$\alpha=1/(1+com)$となり、結果が変わります。

また、spanやcomを使う以外にも、alpha=で$\alpha$のあたいを直接指定することも可能です。

pandasで移動平均や高値線、安値線を計算する

前回がローソク足だったので今回も市場データでよく使われるテクニックから移動平均を取り上げてみたいと思います。
ついでにHLバンド(ドンチャンチャンネル/高値線,安値線)も同様にもとまるので紹介します。

技術としては、window関数と呼ばれる種類の関数を使って算出します。

ドキュメントはこの辺り。
Window
pandas.DataFrame.rolling
pandas.Series.rolling

DataFrameとSeries両方に実装されていて、同じように使うことができます。
rolling() で 指定期間ごとに区切ったデータを作り、それに対して、 meanやmax,minなどの関数を適用して
平均や最大値、最小値を算出して配列として返します。

実際に見た方が早いと思うので、乱数でランダムウォークデータを生成し、
計算して可視化してみましょう。


import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# データ作成
data = pd.DataFrame(np.random.normal(0, 100, 200).cumsum() + 20000)
# 移動平均、期間の高値/安値線の計算
ma = data.rolling(10).mean()
high_band = data.rolling(20).max()
low_band = data.rolling(20).min()
# 可視化
fig = plt.figure(figsize=(12, 7))
ax = fig.add_subplot(1, 1, 1)
ax.plot(data, label="元データ")
ax.plot(ma, label="移動平均")
ax.plot(high_band, label="高値線")
ax.plot(low_band, label="安値線")
plt.legend()
plt.show()

結果がこちら。

期間の最初の方はデータ不足により線が途切れています。
この辺りの制御は min_periods などの引数で細かく調整できるので、
データ量やその時の目的に応じて調整して使いましょう。

pythonでローソク足を描く

以前(pythonを勉強し始めた頃)は、matplotlibでローソク足をかけたはずなのですが、最近は方法が変わってしまったようなのでそのメモです。
なお、ここでサンプルに使うデータはすでにcsvフィアルか何かで保存されているものとします。

以前は matplotlib.finance というのをimport でき、これを使ってかけたのですが、version 2.0 からなくなってしまったようです。
matplotlib.finance
This module is deprecated in 2.0 and has been moved to a module called mpl_finance.

そしてさらに良くないことに、移動先の mpl_finance ですが、あまりしっかり保守されてない様子。

githubのリポジトリに下記の文言があります。
The code is provided as is and is basically un-maintained.

ただ、一応動くようなので動かしてみましょう。
anacondaには含まれていないようなので、インストールから必要です。
pip install mpl_finance
これで、
mpl-finance==0.10.0
が入りました。

さて、使い方ですがun-maintainedの宣言通り、 mpl-finance の公式ドキュメントらしきものは見当たらず、
上の、matplotlib.finance時代のドキュメントを読んで使わないといけないようです。

ローソク足を書く関数は次の4つあり、それぞれデータの渡し方が違います。
.candlestick2_ochl(ax, opens, closes, highs, lows, width=4, colorup=’k’, colordown=’r’, alpha=0.75)
.candlestick2_ohlc(ax, opens, highs, lows, closes, width=4, colorup=’k’, colordown=’r’, alpha=0.75)
.candlestick_ochl(ax, quotes, width=0.2, colorup=’k’, colordown=’r’, alpha=1.0)
.candlestick_ohlc(ax, quotes, width=0.2, colorup=’k’, colordown=’r’, alpha=1.0)

今回は手元のデータと相性が良いので .candlestick_ohlc を使います。
quotes に 日付、始値、高値、安値、終値、の5列のデータがデータ件数行だけ並んだ配列を渡してあげる必要があります。
ここで面倒なのは日付の渡し方で、float型で渡す必要があります。
ドキュメントに time must be in float days format – see date2numとある通り、専用の関数があるのでそれを使います。

matplotlib.dates.date2num(d)

また、この関数は引数がdatetime型なので、元々が2019-05-07 のような文字列になっているならば、
datetime型に変換しておく必要があります。
それにはpandasの to_datetimeを使います。
pandas.to_datetime
(いつもならそれぞれ1記事使ってるようなテクニックですね。to_datetimeの方は便利なのでそのうち専用記事書くかも。)

前置きが長くなりましたが、ここまでの情報でできるので日経平均のcsvデータからローソク足を書いてみましょう。


import pandas as pd
import mpl_finance
import matplotlib.pyplot as plt
from matplotlib.dates import date2num

# データの読み込み
df = pd.read_csv("./日経平均データ.csv")
print(df.head())

'''
        date      open      high       low     close
0  2019-3-11  21062.75  21145.94  20938.00  21125.09
1  2019-3-12  21361.61  21568.48  21348.81  21503.69
2  2019-3-13  21425.77  21474.17  21198.99  21290.24
3  2019-3-14  21474.58  21522.75  21287.02  21287.02
4  2019-3-15  21376.73  21521.68  21374.85  21450.85
'''

# dateの型変換
# まずdatetime型にする
print(df["date"].dtypes)  # object
df["date"] = pd.to_datetime(df["date"])
print(df["date"].dtypes)  # datetime64[ns]
# 続いて float型へ
df["date"] = matplotlib.dates.date2num(df["date"])
print(df["date"].dtypes)  # float64

print(df.head())
'''
       date      open      high       low     close
0  737129.0  21062.75  21145.94  20938.00  21125.09
1  737130.0  21361.61  21568.48  21348.81  21503.69
2  737131.0  21425.77  21474.17  21198.99  21290.24
3  737132.0  21474.58  21522.75  21287.02  21287.02
4  737133.0  21376.73  21521.68  21374.85  21450.85
'''

# 可視化
fig = plt.figure(figsize=(13, 7))
ax = fig.add_subplot(1, 1, 1)
mpl_finance.candlestick_ohlc(ax, df.values)
plt.show()

こうして出来上がるチャートが次です。

正直これ単体では手間の割に可視化するメリットがないなーというのが正直なところです。
ただ、matplotlibの仕組みに乗っかっているので、
自分のオリジナルの指標などを追加していくことができます。

matplotlibのpcolorとpcolormesh

先日のトピックモデルの記事中で、試しにヒートマップでの可視化を試みた時、
matplotlibのpcolorって関数を使用しました。

参考:pythonでトピックモデル(LDA)

matplotlibでヒートマップを描こうと思うと少々無理やりな実装になるものも含めて、
imshowや、contourf、pcolorなど複数の方法が考えられ、結構迷いますが一番自然に書ける気がして採用しました。
しかしどうやらこのpcolor、あまり評判がよろしくないようです。

公式ドキュメントを見ても次ように記載があります。

matplotlib.axes.Axes.pcolor

Hint

pcolor() can be very slow for large arrays. In most cases you should use the similar but much faster pcolormesh instead. See there for a discussion of the differences.

要するに pcolormesh を使う方が良いようです。

ドキュメントはこちら。
matplotlib.axes.Axes.pcolormesh

Differences between pcolor() and pcolormesh()
の部分を読んでも、あまり pcolorにメリットを感じないので、
いっそのこと pcolor 自体を pcolormesh のエイリアスか何かに変えてしまっても良さそうなのですが、
戻り値の型が違うこともありますし、何か同じように使えない事情もあるようですね。

とりあえず今後はこれまでpcolorを使っていた場面ではpcolormeshを使うようにしようと思います。
ちなみに、x軸y軸が等間隔の場合はimshowの方がさらに速いそうです。

imshow
If X and Y are each equidistant, imshow can be a faster alternative.

標準化レーベンシュタイン距離

以前の記事で、レーベンシュタイン距離を計算できるライブラリを紹介しました。
参考:pythonで編集距離(レーベンシュタイン距離)を求める

どちらかというと、ライブラリよりもアルゴリズム側の説明の続きなのですが、
標準化されたレーベンシュタイン距離(normalized Levenshtein distance)というものも提案されています。
これは、二つの文字列のレーベンシュタイン距離を、文字数が多い方の文字数で割った値として定義されます。
固有名詞の名寄せなどでレーベンシュタイン距離を使う場合、
こちらを使った方がうまく行くことが多いようです(個人的な経験から。)

以前紹介した、 python-Levenshtein
実装されているんじゃないかと期待してしばらく調べていたのですが、どうやらこれは実装されてないようです。
特にLevenshtein.ratio という関数に期待したのですがこれは全然違いました。

ということで自分で実装しましょう。


import Levenshtein


def normalized_distance(text0, text1):
    dist = Levenshtein.distance(text0, text1)
    max_len = max(len(text0), len(text1))
    return dist / max_len

関数名は normalized_levenshtein_distance にしたかったが流石に長すぎるので少し短縮。
ただ、これでも長いですね。(pep8的につらい)

これによって、例えば通常のレーベンシュタイン距離では、
“バニラ” と “アイス” の組み合わせも “チョコレート” と “チョコレートアイス” もどちらも距離は3でしたが、
標準化した距離を使うことで、
前者は距離1、後者は 1/3(=0.333…)と算出されるようになり、前者の方が離れてると見なせます。

使ってみた結果がこちら。


print(Levenshtein.distance("アイス", "バニラ"))  # 3
print(Levenshtein.distance("チョコレートアイス", "チョコレート"))  # 3
print(normalized_distance("アイス", "バニラ"))  # 1.0
print(normalized_distance("チョコレートアイス", "チョコレート"))  # 0.3333333333333333

pythonでトピックモデル(LDA)

この前の記事で、scikit-learnのニュース記事のジャンルをロジスティック回帰で予測するというモデルを作ってみました。

参考:scikit-learnのニュースデータをロジスティック回帰で分類

今回はアプローチを変えて、トピックモデルを試してみようと思います。
どちらかというと、20newsのデータセットでもう少し何かやりたいというのが主目的で、
トピックモデルの理論的な説明については今回は省略します。
興味のあるかたへは、講談社から出ている岩田具治先生の、 トピックモデル (機械学習プロフェッショナルシリーズ)
が非常にわかりやすかったのでおすすめです。ページ数も少なめでありがたい。
(数式が多くて書くのが大変なのですがゆくゆくは時系列分析みたいにこのブログでも説明したい。)

さて、pythonでトピックモデルを実装するには gensim を使うのが一般的のようです。
gensim topic modelling for humans
ただ、今回はいつも使っているscikit-learnでやってみました。
(gensimはword2vec等で使ってるのですがscikit-learnに比べると少し苦手。)

scikit-learnでトピックモデルを実装するために読むドキュメントはこちら。
sklearn.decomposition.LatentDirichletAllocation
Topic extraction with Non-negative Matrix Factorization and Latent Dirichlet Allocation

サンプルコードと同じことをしてもしょうがないので、少し工夫をしています。
・サンプルデータのカテゴリーを前回の記事同様に5個に絞る(その代わりそのカテゴリの全データを使用)
・カテゴリーごとに各文章のトピッックを可視化

前置きが長くなりましたが、やってみましょう。
必要ライブラリーのインポートとデータの読み込み


from sklearn.feature_extraction.text import CountVectorizer
from sklearn.decomposition import LatentDirichletAllocation
from sklearn.datasets import fetch_20newsgroups
import matplotlib.pyplot as plt
import numpy as np

remove = ('headers', 'footers', 'quotes')
categorys = [
        "rec.sport.hockey",
        "soc.religion.christian",
        "sci.med",
        "comp.windows.x",
        "talk.politics.mideast",
    ]
twenty_news = fetch_20newsgroups(
                                subset='all',
                                remove=remove,
                                categories=categorys
                            )
X = twenty_news.data

続いて、単語の出現頻度を数え、LDAのモデルを構築して学習します。
トピック数は カテゴリー数と同じ5でも試したのですが、
どうやら6か7にして、あまり重要でない単語を引き受けるトピックを作った方が納得性の高いものになりました。
サンプルコードは20カテゴリーを10トピックでうまく処理できているのに何故だろう?


# 単語の出現頻度データを作成
tf_vectorizer = CountVectorizer(max_df=0.90, min_df=5, stop_words='english')
tf = tf_vectorizer.fit_transform(X)
len(tf_vectorizer.get_feature_names())

# LDAのモデル作成と学習
lda = LatentDirichletAllocation(
                            n_components=7,
                            learning_method='online',
                            max_iter=20
                        )
lda.fit(tf)

それでは、学習した7個のトピックについて、それぞれの頻出語をみてみます。


features = tf_vectorizer.get_feature_names()

for tn in range(7):
    print("topic #"+str(tn))
    row = lda.components_[tn]
    words = ', '.join([features[i] for i in row.argsort()[:-20-1:-1]])
    print(words, "\n")

出力は下記の通りです。(乱数の影響で、モデルの学習をやり直すと結果は変わります。)

topic #0
god, people, think, don, know, just, like, does, say, believe, jesus, church, time, way, did, christ, things, good, christian, question

topic #1
25, 10, 11, 12, 14, 16, 15, 17, 20, 13, 18, 19, 55, 30, la, period, 24, 21, pit, 92

topic #2
armenian, armenians, turkish, people, turkey, armenia, turks, greek, genocide, russian, azerbaijan, government, history, muslim, university, soviet, war, 000, ottoman, killed

topic #3
game, don, said, team, just, didn, hockey, like, know, went, year, time, games, think, got, people, going, did, ll, came

topic #4
israel, jews, jewish, israeli, arab, state, people, world, right, public, arabs, rights, human, war, anti, peace, adl, states, country, palestinian

topic #5
medical, health, disease, cancer, patients, use, new, hiv, doctor, season, good, treatment, years, aids, high, drug, number, time, information, vitamin

topic #6
edu, use, file, window, com, server, program, dos, windows, available, motif, using, version, widget, sun, set, display, mit, x11, information

#1があまり意味のない数値を引き受けてくれていますが、
それ以外は、トピックごとに、宗教や国際的な話題、スポーツに医療に、コンピューターなどの単語が分類されています。

最後に、元の各テキストが、カテゴリーごとに妥当なトピック(話題)を持つと判定さているのか可視化してみてみましょう。
どんな可視化方法が一番わかりやすいか色々試したのですが、カラーマップが比較的良さそうでしたので紹介します。
(このほか箱ひげ図などもそこそこ綺麗に特徴が出ましたが。)


topic_data = lda.transform(tf)
fig = plt.figure(figsize=(6, 25))
for i in range(5):
    ax = fig.add_subplot(6, 1, 1+i)
    im = ax.pcolor(topic_data[twenty_news.target == i], vmax=1, vmin=0)
    fig.colorbar(im)
    # 軸の設定
    ax.set_xticks(np.arange(7) + 0.5, minor=False)
    ax.set_xticklabels(np.arange(7))
    ax.set_title(twenty_news.target_names[i])
plt.show()

出力がこちら。

概ね、カテゴリーごとに別のトピックに分類されているのがみて取れます。

プライバシーポリシーページを更新しました

このブログのフッターからアクセスできるプライバシーポリシーページを実態に合わせて更新しました。
もともとwordpressがデフォルトで作成した英語のページがあったのですが、
今回日本語で作っています。

手順

  1. wordpressの管理画面の左ペインの設定内にあるプライバシーを選択
  2. プライバシーポリシーページを変更する の 新規ページを作成
  3. 作成されたページを実態に合わせて編集

初期設定のページは英語だったのですが、サイト設定を日本語にしているせいか、新規作成したページは日本語でした。
ただ、内容がかなり充実しすぎていて関係ない記載も多々あったのでそれらは削除し、
コメントに関する記載だけ残しています。

そこに、アクセス解析のため入れているGAのことを記載して完了です。

Google アナリティクス利用規約
7. プライバシー に明記されていルルールにちゃんと対応できてないのを気にしていたのですが、これで一安心です。

お客様はプライバシー ポリシーを公開し、そのプライバシー ポリシーで、お客様がデータ収集のために Cookie を使用していることを必ず通知するものとします。

scikit-learnのニュースデータをロジスティック回帰で分類

以前書いた、ニュース記事のテキストのサンプルデータを読み込む記事ですが、読み込んだ後何かした記事をまだ書いてなかったのでちょっとやってみようと思います。
参考:20ニュースグループのテキストデータを読み込んでみる
といっても、凝ったモデルを作るのではなく、対した前処理もせずに単純なBoWと ロジスティック回帰だけでどの程度の性能が出るものなか見てみます。

まずはライブラリのインポートと、データの読み込みです。
20種類全部使うとデータが多く時間がかかるのと、miscのカテゴリーも入りどうやってもある程度以上の性能は出ないので、
適当に5カテゴリー選択しました。(categoriesという引数で指定しています。)

また、テキストにヘッダーやフッターが含まれると、ほとんどそこだけの情報で分類できてしまうので、
removeを使って除去しています。
この辺りの仕様はドキュメント参照。
sklearn.datasets.fetch_20newsgroups


# 必要ライブラリのインポート
from sklearn.datasets import fetch_20newsgroups
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score
from sklearn.metrics import classification_report
from sklearn.metrics import confusion_matrix
from hyperopt import hp,  fmin, tpe, Trials, space_eval

# データの読み込み

remove = ("headers", "footers", "quotes")
categorys = [
        "rec.sport.hockey",
        "soc.religion.christian",
        "sci.med",
        "comp.windows.x",
        "talk.politics.mideast",
    ]

twenty_train = fetch_20newsgroups(
                                subset="train",
                                remove=remove,
                                categories=categorys
                            )
twenty_test = fetch_20newsgroups(
                                subset="test",
                                remove=remove,
                                categories=categorys
                            )

X_train = twenty_train.data
y_train = twenty_train.target
X_test = twenty_test.data
y_test = twenty_test.target

# パラーメータチューニングのため訓練データを2つに分ける
X_train_A, X_train_B, y_train_A, y_train_B = train_test_split(
                                                            X_train,
                                                            y_train,
                                                            test_size=0.2,
                                                            stratify=y_train
                                                        )

見ての通り、X_train, X_test をさらに二つのグループに分けています。
ハイパーパラメーターをチューニングする際に、ここで作ったAグループで学習して、Bグループで評価するようにし、
デフォルトで用意されているテストデータは最後の評価時まで触らずに取っておきます。
本当はクロスバリデーションなどを真面目にやった方がいいのですが、今回はこの方針で。

続いて、ハイパーパラメーターの決定に進みます。
グリッドサーチでもいいのですが、せっかくなのでhyperoptを使ってみます。

spaceは最初はもっと広い範囲で探索していましたが、何度かトライして絞り込みました。
(lr__penalty の l1なども試していたのですが、 Cが大きい時に非常に時間がかかったのと、
今回のデータではl2の方が性能が良かったので探索範囲からも除外。)


# ハイパーパラメーターの探索の準備
space = {
        'cv__min_df': 1 + hp.randint('min_df', 20),
        'cv__max_df': hp.uniform('max_df', 0.5, 0.9),
        'lr__penalty': hp.choice('penalty', ['l2']),
        'lr__C': hp.loguniform('C', -5, 5),
    }


def create_model(args):
    clf = Pipeline(
        [
            ("cv", CountVectorizer(min_df=args["cv__min_df"], max_df=args["cv__max_df"])),
            ("lr", LogisticRegression(C=args["lr__C"], penalty=args["lr__penalty"]))
        ]
    )
    return clf


def objective(args):
    clf = create_model(args)
    clf.fit(X_train_A, y_train_A)
    return - accuracy_score(y_train_B, clf.predict(X_train_B))


trials = Trials()
best = fmin(
            fn=objective,
            space=space,
            algo=tpe.suggest,
            max_evals=100,
            trials=trials,
        )

# 結果の表示
print(best)
print(space_eval(space, best))
# 以下出力
'''
100%|██████████| 100/100 [01:42<00:00,  1.02it/s, best loss: -0.9016949152542373]
{'C': 0.11378962521452059, 'max_df': 0.5208423432316657, 'min_df': 0, 'penalty': 0}
{'cv__max_df': 0.5208423432316657, 'cv__min_df': 1, 'lr__C': 0.11378962521452059, 'lr__penalty': 'l2'}
'''

これでパラーメーターが決まりましたので、モデルを作って評価します。
ここでの評価には最初にとっておいたテストデータを使います。


# 最良のパラメーターでモデルを構築し学習
clf = create_model(space_eval(space, best))
# 訓練データはすべて使う
clf.fit(X_train, y_train)

# 評価
print("正解率:", accuracy_score(y_test, clf.predict(X_test)))
print(classification_report(y_test, clf.predict(X_test), target_names=twenty_test.target_names))

# 以下出力
'''
正解率: 0.8615071283095723
                        precision    recall  f1-score   support

        comp.windows.x       0.90      0.90      0.90       395
      rec.sport.hockey       0.79      0.95      0.86       399
               sci.med       0.87      0.80      0.84       396
soc.religion.christian       0.88      0.82      0.85       398
 talk.politics.mideast       0.87      0.83      0.85       376

           avg / total       0.86      0.86      0.86      1964
'''

正解率約86%。
ほとんど何も工夫していないロジスティック回帰にしてはそこそこの結果だと思います。

NLTKを使う準備をする

普段、文章の形態素解析にはMeCabを使用しているのですが、
とあるサンプルコードを動かそうとした時に、その中でNLTKが使われており、思ったように動かなかったのでそのメモです。

ちなみに、 Anaconda で環境を作ったので、 nltk自体はインストールされていました。


~$ python --version
Python 3.6.8 :: Anaconda, Inc.
~$ pip freeze | grep nltk
nltk==3.3

サンプルコードを動かした時にデータエラーがこちら


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

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

ドキュメントを読むと、何かダウンロードをやらないといけないようです。

Installing NLTK Data

こういう、基本的な処理であってもエラーになります。


>>> import nltk
>>> nltk.word_tokenize('hello nltk!')
Traceback (most recent call last):
    (略)
LookupError:
**********************************************************************
  Resource punkt not found.
  Please use the NLTK Downloader to obtain the resource:

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

  Searched in:
    (略)
**********************************************************************

エラーメッセージを見る限りでは、 ‘punkt’ってのだけで良さそうですが、一気に入れてしまっておきましょう。


~$ python
>>> import nltk
>>> nltk.download()
showing info https://raw.githubusercontent.com/nltk/nltk_data/gh-pages/index.xml

CLI内で完結すると思ったら windowが立ち上がったので少しびっくりしました。

all を選んで Downloadします。
結構時間がかかります。

これで動くようになりました。


>>> import nltk
>>> nltk.word_tokenize("hello nltk!")
['hello', 'nltk', '!']

scipyで距離行列を計算する

前の記事でちらっと pdist関数が登場したので、scipyで距離行列を求める方法を紹介しておこうと思います。

距離行列の説明はwikipediaにあります。
距離行列 – Wikipedia

要するに、N個のデータに対して、(i, j)成分がi番目の要素とj番目の要素の距離になっているN*N正方行列のことです。

これは次の二つの関数を使って計算できます。
scipy.spatial.distance.pdist
scipy.spatial.distance.squareform

numpyで適当に5点取ってやってみましょう。


import numpy as np
from scipy.spatial.distance import pdist
from scipy.spatial.distance import squareform

# 出力する桁数を抑える
np.set_printoptions(precision=3)
# 乱数生成
X = np.random.randint(-5, 6, size=(5, 2))
print(X)
'''
[[-2 -4]
 [-3 -4]
 [ 2 -1]
 [ 4 -2]
 [-1 -2]]
'''

y = pdist(X)
print(y)
'''
[1.    5.    6.325 2.236 5.831 7.28  2.828 2.236 3.162 5.   ]
'''

M = squareform(y)
print(M)
'''
[[0.    1.    5.    6.325 2.236]
 [1.    0.    5.831 7.28  2.828]
 [5.    5.831 0.    2.236 3.162]
 [6.325 7.28  2.236 0.    5.   ]
 [2.236 2.828 3.162 5.    0.   ]]
'''

変数Mに距離行列が入っているのが確認できますね。
(1,2)成分が [-2 -4] と [-3 -4] の距離の1,
(1,3)成分が [-2 -4] と [ 2 -1] の距離の5など、きちんと距離が入っています。

なお、 pdistの戻り値自体は、正方行列の形をしておらず、squareformで正方行列に変形する必要があるので注意です。
pdist の戻り値(コード中では変数y)には、上三角行列部分の値が、1行めから順番に入っています。
配列の長さは$(N*(N-1)/2) = 5*4/2 = 10$ です。

pdist 関数は metric という引数に様々な距離関数の値をとることができ、いろんな種類の距離行列を作ることができます。
(デフォルトはユークリッド距離)

ここまで紹介した、pdistは1つのデータセットに対して距離行列を生成しましたが、
cdistという関数を使うと、2個のデータセットからそれぞれ1個ずつ要素を抽出して距離行列を作ることもできます。

こちらはsquareformを使わなくても初めから行列の形で結果が得られます。


from scipy.spatial.distance import cdist
XA = X[:3]
XB = X[3:]
print(cdist(XA, XB))

# 出力
'''
[[6.325 2.236]
 [7.28  2.828]
 [2.236 3.162]]
'''

最初の例においても、
pdistとsquareform を使うよりも、 cdist(X, X) した方が便利かもしれません。
結果は同じになります。