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

はじめに

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

はじめに
関数がいっぱいある
使い方
- linkage
- fcluster
- cophenet
- dendrogram
実践編
- データを作る
- 手法を選ぶ
クラスタに分ける
- デンドログラムを描く
- 遊ぶ
まとめ

関数がいっぱいある

　いっぱいあるんですよ。

Hierarchical clustering (scipy.cluster.hierarchy) — SciPy v1.3.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.3.0 Reference Guide

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

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

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

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

scipy.spatial.distance.squareform — SciPy v1.3.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.3.0 Reference Guide

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

cophenet

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

scipy.cluster.hierarchy.cophenet(Z, Y=None)
scipy.cluster.hierarchy.cophenet — SciPy v1.3.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.3.0 Reference Guide

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

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

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

truncate_mode

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

color_threshold

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

labels

　葉のラベル。

　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()

　よさげなので、この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()

　上の方でまとまっているのはあまりあてにならないので、実質的に意味があるのはとりあえず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()

　この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とプログラミングのこと