静かなる名辞

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

2019/03/22:TechAcademyがteratailの質問・回答を盗用していた件
2019/03/26:TechAcademy盗用事件 公式発表と深まる疑念



【python】numpyでバイナリサーチをするsearchsorted

 numpy.searchsortedを使うとnumpyでソート済み配列に対するバイナリサーチ(二分探索)を行えます。

numpy.searchsorted — NumPy v1.16 Manual

 ただし、1次元配列しか受け付けません。まあ、いいか。

 次のように使えます。

>>> import numpy as np
>>> a = np.linspace(0, 100, 1000)
>>> a.shape
(1000,)
>>> np.searchsorted(a, 64)
640
>>> a[640]
64.06406406406407

 つまり、ソート済みの一次元配列を第一引数に、探したい値を第二引数に渡すと、順序を狂わせないで第二引数を挿入できる場所のインデックスを返します。そういう意味では、普通の二分探索とは少し異なります。もしぴったり一致するものがほしいのであれば、返り値のインデックスを得てから自分でチェックしたほうが良いでしょう。

 optionalな引数としてはsideとsorterがあります。

 sideは、"right"と"left"という文字列を指定できます。デフォルトは"left"です。お察しの方もおられると思いますが、配列内のある要素とぴったり一致する値を指定したとき、左右どちらが等号なしの不等号になるかを示します。

side returned index i satisfies
left a[i-1] < v <= a[i]
right a[i-1] <= v < a[i]

>>> np.searchsorted(a, 64.06406406406407, side="left")
640
>>> np.searchsorted(a, 64.06406406406407, side="right")
641

 もう一つの重要な引数はsorterです。名前からはソートを行うcallableでも渡すのかな? という気がしますが、実際に渡さないといけないのはargsortした結果です。

>>> a = np.random.rand(1000)
>>> idx = a.argsort()
>>> np.searchsorted(a, 0.5, sorter=idx)
539
>>> a[idx][539]
0.5011315014064937

 しょうじき使いづらいと思った。

 他に注意しておくべきことは、渡した配列の範囲内に値が入れられなかった場合の返り値についてです。

>>> a = np.random.rand(1000) # aは上でrandで生成したので[0, 1)の区間に収まる
>>> np.searchsorted(a, -1, sorter=idx)
0  # 0になる
>>> np.searchsorted(a, 2, sorter=idx)
1000  # 配列の要素数n+1になる

 範囲内に収まる保証がないときは、何らかのチェックが必要になるでしょう。

 最後に、一応素直に線形探索する場合と比べた速度を見ておきます。

>>> a = np.arange(10**4)
>>> np.nonzero(a == 5000)[0][0]
5000
>>> np.searchsorted(a, 5000)
5000
>>> import timeit
>>> timeit.timeit(lambda : np.nonzero(a == 5000)[0][0], number=10**5)
1.3569234880014847
>>> timeit.timeit(lambda : np.searchsorted(a, 5000), number=10**5)
0.30346718700093334

 少なくとも数倍の高速化が見込めるでしょう。もちろん計算量が変わってくるので、実際のところどれだけ速くなるかはケースバイケースです。

【python】np.randomの関数の配列サイズの渡し方をいつも忘れるのでメモった

 np.random以下には色んな乱数生成関数があるのですが、毎回「生成される配列のサイズの指定方法がわからない、なんだっけ?」と思っているので、この際備忘録として残しておきます。

Random sampling (numpy.random) — NumPy v1.16 Manual

 を見るとわかりますが、

  • randおよびrandnだけが
rand(d0, d1, …, dn)	
randn(d0, d1, …, dn)

 というフォーマットを取る。

  • それ以外は主要なものはすべてsize引数で渡す(array.shapeと同じフォーマットのtuple)

 ということになっており、大変覚えやすいのでした。

 これまでは、

np.random.normal(2,3,4)

 とか、

np.random.rand(size=(3,4,5))

 とか、あとはなぜか、

np.random.randn(shape=(4,5,6))

 みたいに書いたりして怒られるということがしょっちゅうありましたが、この備忘録を書いたことで今後は間違えないで済みそうです。

【python】Pool.mapをProcessPoolExecutor.mapに置き換えてみる

はじめに

 concurrent.futures.ProcessPoolExecutorは便利そうなので、Poolの代わりに使ってみようと思います。

17.2. multiprocessing — プロセスベースの並列処理 — Python 3.6.5 ドキュメント
17.4. concurrent.futures – 並列タスク実行 — Python 3.6.5 ドキュメント

相違点

 非同期で使えることを除くと、以下のような違いがあります。

  • ProcessPoolExecutor.mapは複数の引数をサポートする。つまり普通のmapみたいに使える
  • ただし返り値がitertools.chainになる。リストでほしければ自分で変換する

 利便性の観点からすると、まあまあよさげです。

実験1:使えるかどうか確認する

 とりあえず同等に使えることを確認するため、次のようなコードを書いてみました。

import time
import numpy as np
from multiprocessing import Pool
from concurrent.futures import ProcessPoolExecutor

def f1(x):
    a, b = x
    return a * b

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

def main():
    p1 = Pool(4)
    p2 = ProcessPoolExecutor(4)
    a, b = np.random.rand(2, 10**2, 10**2)
    alst = [a+x for x in range(1000)]
    blst = [b+x for x in range(1000)]

    t1 = time.time()
    result1 = p1.map(f1, zip(alst, blst))
    t2 = time.time()
    print("Pool:{:.6f}".format(t2 - t1))

    t1 = time.time()
    result2 = list(p2.map(f2, alst, blst))
    t2 = time.time()
    print("ProcessPoolExecutor:{:.6f}".format(t2 - t1))
    print(np.array_equal(result1, result2))

if __name__ == "__main__":
    main()

# 結果
""" =>
Pool:0.586397
ProcessPoolExecutor:1.137316
True
"""

 動いてるけど、おっそ。ドキュメントには「長いデータ渡すときはchunksizeを指定すると速くなるよ」と書いてあったので、そうしてみます。

# 中略

    result2 = list(p2.map(f2, alst, blst, chunksize=128))

# 中略
""" =>
Pool:0.637173
ProcessPoolExecutor:0.655285
True
"""

 何回か回しているとなんとなく微妙に遅い気もするのですが(たぶん5%~10%くらい?)、とりあえず一応ほぼ互角。

実験2:メソッドを試す

 インスタンスメソッドでも行けるかどうか見ます。

import time
import numpy as np
from multiprocessing import Pool
from concurrent.futures import ProcessPoolExecutor

class F:
    def f1(self, x):
        a, b = x
        return a * b

    def f2(self, a, b):
        return a * b

def main():
    p1 = Pool(4)
    p2 = ProcessPoolExecutor(4)
    a, b = np.random.rand(2, 10**2, 10**2)
    alst = [a+x for x in range(1000)]
    blst = [b+x for x in range(1000)]
    f = F()

    t1 = time.time()
    result1 = p1.map(f.f1, zip(alst, blst))
    t2 = time.time()
    print("Pool:{:.6f}".format(t2 - t1))

    t1 = time.time()
    result2 = list(p2.map(f.f2, alst, blst, chunksize=128))
    t2 = time.time()
    print("ProcessPoolExecutor:{:.6f}".format(t2 - t1))
    print(np.array_equal(result1, result2))

if __name__ == "__main__":
    main()

 普通に走ったので問題なさそう。古めのPythonだとメソッドは渡せなかったはずなのですが、最近は渡せるみたいですね。

実験3:if __name__ == "__main__":を外す

 multiprocessingを使っているとたまにこれがないと落ちることがあるみたいなので、外してみます。

import time
import numpy as np
from multiprocessing import Pool
from concurrent.futures import ProcessPoolExecutor

def f1(x):
    a, b = x
    return a * b

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

p1 = Pool(4)
p2 = ProcessPoolExecutor(4)
a, b = np.random.rand(2, 10**2, 10**2)
alst = [a+x for x in range(1000)]
blst = [b+x for x in range(1000)]

t1 = time.time()
result1 = p1.map(f1, zip(alst, blst))
t2 = time.time()
print("Pool:{:.6f}".format(t2 - t1))

t1 = time.time()
result2 = list(p2.map(f2, alst, blst, chunksize=128))
t2 = time.time()
print("ProcessPoolExecutor:{:.6f}".format(t2 - t1))
print(np.array_equal(result1, result2))

 これも問題なかったです。

結論

 概ね同等に使えます。なにしろ引数を複数渡せるのが便利なところで、wrapper関数をわざわざ書かなくて済みます。あとはlistに変換するまでは非同期で動いている訳で、それを期待して使うのも状況次第ではありでしょう。

 ただしchunksizeは明示的に大きめに指定することです。

【python】str.findとstr.indexの違い

はじめに

 pythonで文字列に対して部分文字列の検索を行うメソッドとして、str.findとstr.indexがあります。ほとんど同じものなのですが、どのような違いがあるのでしょうか?

 str.findとstr.indexはどちらも文字列のメソッドで、引数に渡した文字列の位置を返します。

>>> "hoge".find("og")
1
>>> "hoge".index("og")
1

 「一体なにが違うんだっけ」とふと思って調べてしまったので、メモします。

違い

 大きな違いはありません。ヒットする場合はほとんど同じです。しかし、ヒットしなかった場合の挙動だけ異なります。

 端的に言えばstr.findは見つからなかった場合-1を返しますが、str.indexは例外を送出します。

>>> "hoge".find("fuga")
-1
>>> "hoge".index("fuga")
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: substring not found

str.index(sub[, start[, end]])(原文)
find() と同様ですが、部分文字列が見つからなかったとき ValueError を送出します。

組み込み型 — Python 3.7.3 ドキュメント

 それだけの違いのようです。単純ですね。

速度差

 まさかこんなのに速度差なんてないだろう、そう思っていた時期が私にもありました。

>>> s = "hoge"*1000 + "fuga"
>>> timeit.timeit(lambda : s.find("fuga"))
3.8523642780000955
>>> timeit.timeit(lambda : s.index("fuga"))
3.803682490000938

 ヒットする場合は互角です。

>>> def f1():
...     ret = s.find("piyo")
...     if ret == -1:
...         pass
...     else:
...         pass
... 
>>> def f2():
...     try:
...         ret = s.index("piyo")
...     except:
...         pass
... 
>>> timeit.timeit(f1)
3.907751840999481
>>> timeit.timeit(f2)
4.395755831001225

 例外処理は単純な比較より重いため、少し時間を食うようです。なんてこった。

 もちろん、処理目的にもよるため、「速いからfindの方が良い」とは一概には言えません。意図に応じて選ぶべきです。

 ただ、例外処理を書くのは少し面倒だし、行数も増えるので、findの方が便利に使えるケースが多いかな?

まとめ

 見つからなかったとき、str.findは-1を、str.indexはValueErrorを出すと覚えましょう。

【python】scipyのpdistとsquareformの使い方と仕組み

はじめに

 scipyで距離行列を扱うときはscipy.spatial.distanceのpdist, squareformなどを主に使いますが、長年よくわからないままに使っていたので、整理してまとめておきます。

 なお、以下のドキュメントを参考にします。
scipy.spatial.distance.pdist — SciPy v1.2.1 Reference Guide
scipy.spatial.distance.squareform — SciPy v1.2.1 Reference Guide

pdist

 pdistを使うと様々な距離の距離行列を計算できます。デフォルトはユークリッド距離ですが、他にもいろいろ指定できる距離の種類があります。詳細はドキュメントを見てください。

 そしてpdistは「condensed distance matrix」という表現方式の距離行列を返します。この方式にはあまり馴染みがありませんが、このように使えます。

>>> import numpy as np
>>> a = np.random.random(size=(5,3))
>>> from scipy.spatial.distance import pdist
>>> pdist(a)
array([0.58791701, 0.71055164, 0.46750397, 0.3389748 , 0.39126501,
       0.28675164, 0.79462635, 0.39543688, 0.93094202, 0.59666923])

 結果は常に1次元配列で帰ります。5つのベクトルの距離を計算してみた結果、要素数は10になっています。5*4/2=10なので、距離行列の上下どちらか半分の三角行列がflattenされて返っているのだとわかります。上下の三角行列は対称ですし、対角要素(自分同士の距離)は常に0なので、最低限これだけあれば表現できる訳です。

 結果として省メモリな表現になりますが、人間にとってはわかりづらくなります。また、関数やライブラリによってはこの形式を受け付けないこともままあります。そこで、次のsquareformが登場します。

squareform

 squareformはpdistで生成された配列と、馴染みのある距離行列(「square-form distance matrix」と称しています)の相互変換を担います。

 上の実行例から続けてこう打つと、距離行列が得られます。

>>> from scipy.spatial.distance import squareform
>>> squareform(pdist(a))
array([[0.        , 0.58791701, 0.71055164, 0.46750397, 0.3389748 ],
       [0.58791701, 0.        , 0.39126501, 0.28675164, 0.79462635],
       [0.71055164, 0.39126501, 0.        , 0.39543688, 0.93094202],
       [0.46750397, 0.28675164, 0.39543688, 0.        , 0.59666923],
       [0.3389748 , 0.79462635, 0.93094202, 0.59666923, 0.        ]])

 なお、squareformは1次元配列か2次元配列を引数として受け取ります。1次元配列を受け取ると「condensed distance matrix」→「square-form distance matrix」の変換、2次元配列を受け取ると「square-form distance matrix」→「condensed distance matrix」の変換として動くという仕組みです。要するに何が言いたいかというと、二回通すと元に戻ります。

>>> squareform(squareform(pdist(a)))
array([0.58791701, 0.71055164, 0.46750397, 0.3389748 , 0.39126501,
       0.28675164, 0.79462635, 0.39543688, 0.93094202, 0.59666923])

 便利かどうかはちょっと即答しかねますが、けっこう面白い気はします。

インデックスの変換

 さて、ここまで読んできた方は「squareformのインデックスとpdistのインデックスはどう対応付けられるの?」という興味を持ったかと思います。

 pdistのドキュメントには意味深な一文があります。

Returns a condensed distance matrix Y. For each and (where ),where m is the number of original observations. The metric dist(u=X[i], v=X[j]) is computed and stored in entry ij.

 結論から先にいうと、この文に深い意味はないので、無視して構いません。真に受けるとかえって混乱します。

 日本人で英語が苦手な私には上の一文の意図するところはよくわかりませんが、英語圏の方もわからないようで、stackoverflowでこんなスレッドを見つけました。

python - How does condensed distance matrix work? (pdist) - Stack Overflow

 基本的にはここの回答に書いてあることが参考になります。ただし、python2時代のスレッドなので、整数の割り算が/で記載されていることには注意する必要があります(python3では/を使うとfloatの結果になり、整数の結果を得たければ//を使う必要があります)。

 とはいえ、stackoverflow見てねというのも不親切なので、この記事でも説明しておきます。

 とりあえず4つの要素の距離の「square-form distance matrix」の側の三角行列のインデックスについて考えることにすると、以下のようになります。

上三角行列のインデックス
上三角行列のインデックス

 pdistで作る「condensed distance matrix」の場合は以下のようなものです。

pdistで作った配列のインデックス
pdistで作った配列のインデックス

 対応をまとめると、

0,1→0
0,2→1
0,3→2
1,2→3
1,3→4
2,3→5

 となります。

 以下は上記stackoverflowの回答に基づいた記述ですが、pythonの数式で変換するとすると、

n*j - j*(j+1)//2 + i - 1 - j

 こんなですが、覚える必要はないです。

 また、itertools.combinationsの結果と同じ順で返ったりします。

>>> from itertools import combinations
>>> list(combinations(range(4), 2))  # 4つから2つ取り出す
[(0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3)]

 概念的にはこちらの方がわかりやすいですね。

まとめ

 scipyで距離行列の絡む処理をしようとすると避けて通れないpdistとsquareformですが、実はやっていることは単純です。慣れればそれほど困惑することもないので、気楽に使えるようになると良いと思います。

【python】scipyで階層型クラスタリングするときの知見まとめ

はじめに

 scipyの階層型クラスタリングを使う機会がありましたが、使い方がわかりづらいと思ったのでまとめておきます。

 目次

関数がいっぱいある

 いっぱいあるんですよ。

Hierarchical clustering (scipy.cluster.hierarchy) — SciPy v1.2.0 Reference Guide

 私の数え間違えがなければ31個。多いですね。

 とはいえ、本質的なもの(実際によく使うもの)は以下くらいです。

  • linkage

 実際に階層型クラスタリングを行う。これがないと始まらない。

  • fcluster

 任意の深さで木を切り、クラスタに分割する。

  • cophenet

 yを渡すとコーフェン相関係数なる評価指標を出してくれるらしい。

  • dendrogram

 デンドログラムを描画する。

 他に、こんな関数があります。

  • fclusterdata, leaders

 fclusterの処理と絡む便利関数。

  • single, complet , average,...など

 クラスタリングアルゴリズムに対応。すべてlinkageの引数で文字列を使って指定できるので実際にこれらの関数を使うことはない。

 他にも各種の数値を出せたり、MATLABのフォーマットと相互変換してくれたり、ポインタのリンクで表現された構造化データに変換するクラスだったり、いろいろなものが押し込まれています。必要に応じてリファレンスから探してくれば良いので、それぞれ述べることはしません。

使い方

 とりあえず上で挙げた実際によく使う関数を中心に説明していきますよ。

linkage

 データをlinkageに通すことで階層型クラスタリングが行えます。返り値として木の情報を表す配列が返ります。それに対して用意されている関数であれこれ処理していくというのが基本的な流れです。

scipy.cluster.hierarchy.linkage(y, method='single', metric='euclidean', optimal_ordering=False

scipy.cluster.hierarchy.linkage — SciPy v1.2.0 Reference Guide

 yはデータですが、基本的にはscipy.spatial.distance.pdistで作れる距離行列のフォーマット(一般的な正方距離行列とは異なるので注意)で渡してくれ、ということになっています。

scipy.spatial.distance.pdist — SciPy v1.2.0 Reference Guide

 でもshape=(サンプル数, ベクトル次元数)のnumpy配列も受け取ってくれるので、それほど気を配る必要はありません。

 気を配らないといけないのはむしろ距離行列を入れてクラスタリングしたい場合で、その場合は正方行列として入れようとすると正しく処理されません。pdistのフォーマットに変換する必要があります。これはscipy.spatial.distance.squareformで変換できるので、そうしてください。

scipy.spatial.distance.squareform — SciPy v1.2.0 Reference Guide

 あとは距離(距離行列で入力した場合は無視される)とクラスタリング方法を指定できます。選べるものの一覧はリファレンスを見てください。選択肢はいろいろ実装されているので、ここで「欲しいものがなくて困る」というシチュエーションはそう滅多にないと思います。

 linkageの返り値はnumpy配列です。細かいフォーマットについてはリファレンスを(ry。自分でこれをいじってどうこうしようという機会はあまりないと思うので、知らなくてもなんとなく使えます。リファレンスの説明ではだいたいZという変数に代入しているので、それに倣うと良いと思います。

fcluster

 クラスタリングしてくれます。

scipy.cluster.hierarchy.fcluster(Z, t, criterion='inconsistent', depth=2, R=None, monocrit=None)

scipy.cluster.hierarchy.fcluster — SciPy v1.2.0 Reference Guide

 いちいち細かく述べることはしませんが、Zにlinkageの返り値を入れ、tはスレッショルドでクラスタの分割/結合の基準として渡します。あとはcriterionに好きなアルゴリズムを選ぶ……という仕組みですね。

cophenet

 「cophenetic distances」なるものを計算します。

scipy.cluster.hierarchy.cophenet(Z, Y=None)

scipy.cluster.hierarchy.cophenet — SciPy v1.2.0 Reference Guide

 
 Y(クラスタリングの元になった距離行列。当然pdistフォーマット)を渡すとクラスタリング全体の評価指標を出してくれるので、その目的で使うのが良いと思います。返り値は2つあり得るんですが、Yを渡さなければ1つめ(全体の評価指標)は省略されて2つめだけ返ります。

 MATLABから輸入された関数らしく、MATLABのドキュメントを読んだ方がわかりやすいと思います。

コーフェン相関係数 - MATLAB cophenet - MathWorks 日本

 あとで書く実践編でちょっと触れます。

dendrogram

 デンドログラムを描きます。デンドログラム出しとけばとりあえず納得感があるので、これだけ出してみてからあれこれすれば良いと思います。

scipy.cluster.hierarchy.dendrogram(Z, p=30, truncate_mode=None, color_threshold=None, get_leaves=True, orientation='top', labels=None, count_sort=False, distance_sort=False, show_leaf_counts=True, no_plot=False, no_labels=False, leaf_font_size=None, leaf_rotation=None, leaf_label_func=None, show_contracted=False, link_color_func=None, ax=None, above_threshold_color='b')

scipy.cluster.hierarchy.dendrogram — SciPy v1.2.0 Reference Guide

 引数が泡吹いて倒れるほどいっぱいある。到底説明なんてできないので、代表的なものだけピックアップします。

  • Z

 説明不要。linkageの結果を渡す。

  • p

 木を省略して表示するときのパラメータ。次のtruncate_modeと絡みます。

  • truncate_mode

 どのように木を省略するか。これはデータ数がある程度多いときに威力を発揮します。今回は使いません。

  • color_threshold

 色分けに絡むしきい値。木の中の最大ノード間距離との比率で色分けを決めます。

  • labels

 葉のラベル。

  • ax

 matplotlibのAxesを渡せます。渡すとそこに描画されます。

 あとはまあ、いろいろ。基本的にはどれも表示フォーマットに絡む引数なので、見た目を変えたくなったら合う引数を探すという感じです。

実践編

 ではこれから実践してみます。

データを作る

 階層型クラスタリングは50>データ数くらいが適しています。多くても見づらいし計算量が嵩むからです。これくらいでのサイズの使いやすいデータはすぐに思い浮かびませんでしたが、sklearnのdigitsでラベルごとに平均を取ると良さげなことに気づいたので、そうします。

 ついでに可視化もする。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_digits

def gen_data():
    digits = load_digits()
    label_uniq = np.unique(digits.target)
    result = []
    for label in label_uniq:
        result.append(digits.data[digits.target == label].mean(axis=0))
    return result, label_uniq

def visualize():
    X, y = gen_data()
    fig, axes = plt.subplots(nrows=2, ncols=5)
    for ax, x, label in zip(axes.ravel(), X, y):
        ax.set_title(label)
        ax.imshow(x.reshape(8, 8))
    plt.savefig("data.png")

if __name__ == "__main__":
    visualize()

data.png
data.png

 よさげなので、この10件のデータでやってみます。

手法を選ぶ

 とりあえずmetric="euclidean"に固定してmethodを変化させ、評価指標を出してみましょう。

from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import linkage, cophenet

def clustering_score():
    X, y = gen_data()
    methods = ["single", "complete", "average", "weighted",
               "centroid", "median", "ward"]
    for method in methods:
        S = pdist(X)
        Z = linkage(S, method=method)
        c, d = cophenet(Z, S)
        print("{0} {1:.3f}".format(method, c))

if __name__ == "__main__":
    clustering_score()
single 0.722
complete 0.752
average 0.769
weighted 0.766
centroid 0.681
median 0.730
ward 0.720

 average(UPGMA、いわゆる群平均法)がベストっぽいので、これを使うことにします。

 どれを選んでもだいたい0.75くらいのコーフェン相関係数なので、25%くらいはデータの性質が狂っていると思って結果を解釈する必要がある、ということだと思います。あまり細かいところを見ても無意味です。

クラスタに分ける

 fclusterを使うと適当な数のクラスタに分割できます。4つのクラスタに分割してみましょう。

from collections import defaultdict
from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import linkage, fcluster

def clustering_fcluster():
    X, y = gen_data()
    S = pdist(X)
    Z = linkage(S, method="average")
    result = fcluster(Z, t=4, criterion="maxclust")
    d = defaultdict(list)
    for i, r in enumerate(result):
        d[r].append(i)
    for k, v in d.items():
        print(k, v)

if __name__ == "__main__":
    clustering_fcluster()
1 [1, 2, 3, 5, 8, 9] なんとなく形が似てるかな・・・?
2 [7]  # 独特の位置
3 [4, 6] # 似ているかと言うと微妙なものがある
4 [0]  # 独特の位置

 いまいち納得感の少ない結果になりました。

デンドログラムを描く

 デンドログラムを描いてみましょう。

from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import linkage, dendrogram

def clustering_dendrogram():
    X, y = gen_data()
    S = pdist(X)
    Z = linkage(S, method="average")
    dendrogram(Z)
    plt.savefig("dendro1.png")

if __name__ == "__main__":
    clustering_dendrogram()

dendro1.png
dendro1.png

 上の方でまとまっているのはあまりあてにならないので、実質的に意味があるのはとりあえず1,8と3,9あたりです。

 1,8は縦の真ん中あたりに集中するクラスタ。3,9は9の上の丸の左側を切り取ればほぼ3みたいな形になります。

 上のfclusterでやったのと同じ4つのクラスタに分かれるようにする方法がないのか? 基本的に違う考え方に基づいて着色されるので厳しいものがありますが、縦軸を見ながら4つに分かれるあたりを狙ってcolor_thresholdを決め打ちすると一応それに近いことはできます。

from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import linkage, dendrogram

def clustering_dendrogram2():
    X, y = gen_data()
    S = pdist(X)
    Z = linkage(S, method="average")
    dendrogram(Z, color_threshold=31)
    plt.savefig("dendro2.png")

if __name__ == "__main__":
    clustering_dendrogram2()

dendro2.png
dendro2.png

 このcolor_thresholdは距離の近さがスレッショルド未満のものを1つにまとめる、という挙動です。スレッショルドより大きい距離はabove_threshold_color(デフォルトは"b"で青)になります。今回はスレッショルドより上で1サンプルで別れてそのまま1クラスタを形成するという厄介な子がいるので微妙な結果になってしまいますが、もう少し性質の良いデータだとうまく合わせることはできると思います。

遊ぶ

 せっかくmethodが7種類あるので、それぞれでやってみます。

import matplotlib.pyplot as plt
from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import linkage, cophenet, dendrogram

def clustering_dendro_each_method():
    X, y = gen_data()
    methods = ["single", "complete", "average", "weighted",
               "centroid", "median", "ward"]

    fig, axes = plt.subplots(nrows=2, ncols=4, figsize=(10, 7))
    for ax, method in zip(axes.ravel(), methods):
        S = pdist(X)
        Z = linkage(S, method=method)
        c, d = cophenet(Z, S)
        dendrogram(Z, color_threshold=31, ax=ax)
        ax.set_title("{0} {1:.3f}".format(method, c))
    plt.savefig("dendro_each.png")

if __name__ == "__main__":
    clustering_dendro_each_method()

結果
結果

 けっこう結果が変わるけど、細部はそれなりに揃っている感じ。コメントはとりあえず差し控えたいと思います。

まとめ

 このように使えます。便利関数や引数が多いのがややこしいだけで、基本的な使い方そのものは難しくないです。

 階層型クラスタリングは結果の良し悪しを客観的に評価するのが難しいので、割り切って使うと良いと思います。デンドログラムの細かい枝分かれの順序どうこうより、各クラスタの成因を判別分析や検定などでちゃんと調べてあげる方が本当は大切です。

【python】dictの集合演算を辞書ビューオブジェクトで行う

はじめに

 pythonのdictは便利なデータ型ですが、複数のdictに対して重複を除去する、逆に共通部分のみを抜き出すといった集合のような演算を行いたいときがあります。

 dictそのものは集合演算をサポートしていませんが、辞書ビューオブジェクトというものを使うと集合演算が容易に行なえます。

辞書ビューオブジェクトとは?

 普段よく使うdict.keys(), dict.values(), dict.items()はすべて辞書ビューオブジェクトのことを指します。

4. 組み込み型 — Python 3.6.5 ドキュメント

 たいていはforでイテレータを取り出すためだけに使われるので、あまり詳細には知られてはいませんが、このように振る舞います。

>>> d = {x:str(x) for x in range(3)}
>>> d
{0: '0', 1: '1', 2: '2'}
>>> items = d.items()
>>> items
dict_items([(0, '0'), (1, '1'), (2, '2')])
>>> d["hoge"] = "fuga"
>>> items
dict_items([(0, '0'), (1, '1'), (2, '2'), ('hoge', 'fuga')])

 要するに、元の辞書をずっと参照していて、元の辞書を変更すれば変更が波及するという性質があるので、ビューという名前になっています。

辞書は集合演算できないが、辞書ビューオブジェクトはできる

 さて、辞書に対して集合演算をかけようとするとエラーになります。

>>> d1 = {x:str(x) for x in range(0, 3)}
>>> d2 = {x:str(x) for x in range(1, 4)}
>>> d1 & d2
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: unsupported operand type(s) for &: 'dict' and 'dict'

 これはもうどうにもなりませんね。サポートされていないんですから。

 しかし、辞書ビューオブジェクトは集合演算をサポートします。意外な仕様ですが、ドキュメントにもそう書いてあります。 

キーのビューは、項目が一意的でハッシュ可能であるという点で、集合に似ています。すべての値がハッシュ可能なら、 (key, value) 対も一意的でハッシュ可能なので、要素のビューも集合に似ています。(値のビューは、要素が一般に一意的でないことから、集合に似ているとは考えられません。) 集合に似ているビューに対して、抽象基底クラス collections.abc.Set で定義されている全ての演算 (例えば、 ==、<、^) が利用できます。
4. 組み込み型 — Python 3.6.5 ドキュメント

 早速やってみましょう。

>>> d1.keys() & d2.keys()
{1, 2}

 できました。ちなみに返る型はsetです。

keysはできるがitemsには制約がある。valuesはできない

 上のドキュメントの引用にも書いてある通り、集合演算が使えるためには

  • 項目が一意的でハッシュ可能である

 という条件があります。考えればわかることですが、それぞれ整理しておくと、

  • keysの場合、上の条件を常に満たすので集合演算可能
  • itemsの場合、値がすべてハッシュ可能なら集合演算可能(listなどのimmutableが値に入っていると駄目)
  • valuesは駄目({"hoge":"a", "fuga":"a"}みたいな辞書があり得るため)

 ということになります。

ちょっとした応用

 itemsで辞書の共通部分を抜き出した新しい辞書を作りたい場合、集合型で返りますから辞書に再変換する必要があります。

>>> dict(d1.items() & d2.items())
{1: '1', 2: '2'}

 値にimmutableが含まれるがキーに基づいて処理できれば良い、という場合は、keysで処理してから辞書内包表記で辞書を組み立てることが出来ます。ただし演算の種類によっては、ややこしくなります。

>>> d1 = {x:[x] for x in range(0, 3)}
>>> d2 = {x:[x] for x in range(1, 4)}
>>> {k:d1[k] if k in d1 else d2[k] for k in d1.keys() | d2.keys()}
{0: [0], 1: [1], 2: [2], 3: [3]}

 これは他の方法でやるべきですかね。

まとめ

 このように、辞書ビューオブジェクトを使うと複雑な辞書同士の操作がときに簡単に行なえます。もし必要になったら使ってみてください。