【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で作った配列のインデックス

　対応をまとめると、

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

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

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

　この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]}

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

まとめ

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

本当は怖いSVMと交差検証

概要

　SVMと交差検証を組み合わせて使うと、たとえ交差検証で高いスコアが出て汎化性能確保できた！　と思っても想像とかけ離れた分離超平面になっていることがままある。

　なのでこの組み合わせは少し怖いということを説明する。

コード

　irisを分類します。二次元で決定境界を可視化するために、irisを主成分分析を使って2次元に落としておきます。

　GridSearchCVを使って交差検証し、ベストパラメータを探ります。その後、ベストパラメータの分類器で、平面上で散布図と決定境界を可視化してみます。

　ちょっと長いけど、ざっくり読んでみてください。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.svm import SVC
from sklearn.model_selection import GridSearchCV
from sklearn.decomposition import PCA
from sklearn.datasets import load_iris

def main():
    # データの準備
    iris = load_iris()
    X, y = PCA(2).fit_transform(iris.data), iris.target

    # GridSearchCVの準備
    svm = SVC()
    params = {"C":10**np.arange(-3, 3, dtype=float),
              "gamma":10**np.arange(-3, 3, dtype=float)}
    gscv = GridSearchCV(svm, params, cv=8, iid=False, 
                        return_train_score=False,
                        verbose=1, n_jobs=-1)
    gscv.fit(X, y)

    # GridSearchCVの結果を表示する
    result_df = pd.DataFrame(gscv.cv_results_)
    print(result_df[["param_C", "param_gamma",
                     "rank_test_score", "mean_test_score"]
                ][result_df["rank_test_score"]==1])

    # 最良推定器をclfに代入
    clf = gscv.best_estimator_
    
    # 可視化の準備
    xmin, xmax, ymin, ymax = (X[:,0].min()-1, X[:,0].max()+1,
                              X[:,1].min()-1, X[:,1].max()+1)    
    x_ = np.arange(xmin, xmax, 0.01)
    y_ = np.arange(ymin, ymax, 0.01)
    xx, yy = np.meshgrid(x_, y_)

    # 予測
    zz = clf.predict(np.stack([xx.ravel(), yy.ravel()], axis=1)
                 ).reshape(xx.shape)

    # 可視化
    plt.pcolormesh(xx, yy, zz, cmap="winter", alpha=0.1, shading="gouraud")
    plt.scatter(X[:, 0], X[:, 1], c=y, edgecolors='k', cmap="winter")
    plt.savefig("result.png")

if __name__ == "__main__":
    main()

結果

　printされた結果は良さげなものでした。

   param_C param_gamma  rank_test_score  mean_test_score
22       1          10                1         0.972222

　0.972222なら悪くないaccuracyです。

　しかし、出力された画像は思っていたものとは違いました。

SVMの結果

　こわっ。どう考えてもこういうデータではないと思いますが、こういうことが現実に起こります。(-1.5, 0)あたりにデータが来たら、青かねずみ色のどちらかの色のグループに分類されてほしいところですが、実際は緑になってしまう訳です。

　ちょっと信頼ならないですね、SVM。

怖くない線形SVM

　ついでにいろいろな分類器を見てみます。

　LinearSVCをimportし、「# GridSearchCVの準備」以下「# 最良推定器をclfに代入」の上までを次のように書き換える。

    # GridSearchCVの準備
    svm = LinearSVC()
    params = {"C":10**np.arange(-3, 3, dtype=float)}
    gscv = GridSearchCV(svm, params, cv=8, iid=False, 
                        return_train_score=False,
                        verbose=1, n_jobs=-1)
    gscv.fit(X, y)

    # GridSearchCVの結果を表示する
    result_df = pd.DataFrame(gscv.cv_results_)
    print(result_df[["param_C",
                     "rank_test_score", "mean_test_score"]
                ][result_df["rank_test_score"]==1])

　結果。

  param_C  rank_test_score  mean_test_score
4      10                1         0.960317

　accuracyはわずかに下がるだけ。

線形SVMの結果

　はるかに納得感の高い結果になっています。こういうことがあるので、ほぼ線形分離可能なことがわかっているデータは、まずは線形なモデルで試すことがおすすめです。非線形なモデルで1%とかaccuracyをあげられるとしても、それで未知のデータに対して良い推定ができるかどうかは交差検証ではわからないのです*1。

怖いかどうか悩むランダムフォレスト

　LinearSVCと同様にRandomForestClassifierをimportし、同じ箇所を書き換えます。

    # GridSearchCVの準備
    clf = RandomForestClassifier(n_jobs=-1)
    params = {"n_estimators":10**np.arange(3),
              "min_samples_leaf":[1,2]}
    gscv = GridSearchCV(clf, params, cv=8, iid=False, 
                        return_train_score=False,
                        verbose=1, n_jobs=-1)
    gscv.fit(X, y)

    # GridSearchCVの結果を表示する
    result_df = pd.DataFrame(gscv.cv_results_)
    print(result_df[["param_n_estimators", "param_min_samples_leaf",
                     "rank_test_score", "mean_test_score"]
                ][result_df["rank_test_score"]==1])

　性能。SVMと比べると少し低下するかな。理由はよくわからないけど、木の本数10本で葉の最小サンプル数2のときが最高性能（本当になんで？　基本的に木が多い方が性能が高いはずなのだが・・・）。

  param_n_estimators param_min_samples_leaf  rank_test_score  mean_test_score
4                 10                      2                1         0.953373

　可視化。

ランダムフォレストの結果

　この状態ではrefitしているので、accuracyは1になるはずです。さて、ランダムフォレストの特徴は、決定木なので軸と90度の直線の組み合わせで決定境界が表現されることです。また、少ない分割でエントロピーが下がるような決定境界を追求していくので、データの全体の傾向はあまり見てくれません。

　それでもSVMよりはマシな感じでしょうか。

怖くない気がする多層パーセプトロン

　疲れてきたので最後。やり方は上と同じ。

    # GridSearchCVの準備
    clf = MLPClassifier(max_iter=3000)
    params = {"hidden_layer_sizes":[(5*x,) for x in range(1, 5)]}
    gscv = GridSearchCV(clf, params, cv=8, iid=False, 
                        return_train_score=False,
                        verbose=1, n_jobs=-1)
    gscv.fit(X, y)

    # GridSearchCVの結果を表示する
    result_df = pd.DataFrame(gscv.cv_results_)
    print(result_df[["param_hidden_layer_sizes",
                     "rank_test_score", "mean_test_score"]
                ][result_df["rank_test_score"]==1])

　同率一位が3つもありました。ま、どれでも良いか。

  param_hidden_layer_sizes  rank_test_score  mean_test_score
1                    (10,)                1          0.96627
2                    (15,)                1          0.96627
3                    (20,)                1          0.96627

多層パーセプトロンの結果

　まともな感じ。右側の決定境界はちょっと甘いかなー、という気はします。

　ただしパラメータ数の少ない隠れ層1層の多層パーセプトロンだからこれくらいの素直な結果に落ち着くのであって、深層学習は油断するとすぐ過学習するので注意が必要です。

SVMも怖くない！

　汎化重視のパラメータにすれば変なことにはならないので安心してください。基本的にCもgammaも低ければ低いほど過学習しづらくなります。

　今までと同じ箇所をこう書き換える。

    # SVMをfit
    clf = SVC(C=1, gamma=0.2)
    clf.fit(X, y)

汎化重視のパラメータのSVMの結果

　別に問題ないですよね。SVMそのものはこのようにパラメータで自由に分離超平面の複雑さを調整できる、優れた分類器です。

　ただし、交差検証で機械的にパラメータを決めてしまうとあまり良くない結果を招く可能性がある訳です。

現実的な話

　二次元で見ているから違和感があるのであって、高次元空間はサクサクメロンパン問題があるのでまた異なった挙動になります。

　↓サクサクメロンパン問題

windfall.hatenablog.com

　直感的な理解は難しいと思いますが、こういうことが問題になるケースは少ないということはなんとなくわかります。

　とはいえ、「交差検証して最高性能だったSVMを投入したらトンチンカンな結果を出してくる」みたいなことが実際に起こらないとも限らないので、ある程度は注意した方が良いでしょう。

まとめ

　SVMはパラメータ設定を過学習するようにしたとき、他の手法と比べても群を抜いて変な結果になるのですが、何が駄目なんでしょうね。やっぱりカーネル使ってるからか。

　ちょっと注意が要るよなー、というお話でした。

*1:手持ちのデータの分布に依存する話ですが。未知のデータなんてないぜ！　といえるくらい豊富なデータがあればそれはそれで構いません

tkinterで遅い処理を別スレッドに投げ画面が固まらないようにする

はじめに

　pythonで書いた何かの処理にGUIのwrapperを付けたいので、お手軽にtkinterで作ろう！　みたいなことをする人が果たしてどれだけいるのかは知りませんが、tkinterのコールバックで遅い処理をすることがあります。まあ、あるいは、練習でタイマーとかを作るときにwaitさせたくなるというケースの方が多いのかもしれませんが。

　この状況を再現した下のようなコードについて考えてみます。

import time
import tkinter as tk

def f(event):
    print("start callback")
    s.set("start callback")
    time.sleep(2)
    print("end callback")
    s.set("end callback")
    time.sleep(2)
    print("(空白に戻す)")
    s.set("")

def main():
    global s
    root = tk.Tk()

    s = tk.StringVar()
    s.set("")
    label = tk.Label(root, textvariable=s)
    label.pack()

    button = tk.Button(root, text="Button")
    button.bind("<1>", f)
    button.pack()
    
    root.mainloop()

if __name__ == "__main__":
    main()

　ボタンが押されたとき、コールバック関数の中で動的にラベルを書き換えようという訳ですね。

　でも実行してみるとわかりますが、stdoutに対するprintは出てきても、画面はうまく書き換わりません。

空白の画面のまま・・・

原因

　コールバック関数がreturnするまでの間、呼び出し元の処理は止まっています。よって、画面の更新処理などもすべて止まってしまい、StringVarへの書き換えも反映されません。

　ちなみに、この状態のときに操作しようとすると「応答なし」になると思います。

対策

　コールバック関数で別スレッドにやりたい処理の実行を投げ、コールバック関数そのものはさっさとreturnします。

import time
import tkinter as tk
import threading

def f(event):
    print("start callback")
    s.set("start callback")
    time.sleep(2)
    print("end callback")
    s.set("end callback")
    time.sleep(2)
    print("(空白に戻す)")
    s.set("")

def callback(event):
    th = threading.Thread(target=f, args=(event,))
    th.start()

def main():
    global s
    root = tk.Tk()

    s = tk.StringVar()
    s.set("")
    label = tk.Label(root, textvariable=s)
    label.pack()

    button = tk.Button(root, text="Button")
    button.bind("<1>", callback)
    button.pack()
    
    root.mainloop()

if __name__ == "__main__":
    main()

　今度は問題なく意図通りの結果が得られます。スレッドをjoinしていませんが、コールバック関数で待つわけには行きませんからこうするしかありません。このスレッドは放っておけば（実行が終われば）勝手に消えるので、たぶん大丈夫です。

このように画面が変化する

多重実行対策

　上の方法だとスレッドが幾つでも同時に走り得ます。end callbackになったあたりでもう一度ボタンを押すとカオスな切り替わりが観察できます。

　一つの考え方としては、実行中はボタンをDISABLEDにしてしまえば良いというものがあります。ただし、これはcommandでコールバックを指定した場合しか効きません。bindを使った場合は別の手立てが必要になります。

import time
import tkinter as tk
import threading

def f():
    print("start callback")
    s.set("start callback")
    time.sleep(2)
    print("end callback")
    s.set("end callback")
    time.sleep(2)
    print("(空白に戻す)")
    s.set("")
    button.configure(state=tk.NORMAL)

def callback():
    button.configure(state=tk.DISABLED)
    th = threading.Thread(target=f)
    th.start()

def main():
    global s, button
    root = tk.Tk()

    s = tk.StringVar()
    s.set("")
    label = tk.Label(root, textvariable=s)
    label.pack()

    button = tk.Button(root, command=callback, text="Button")
    button.pack()
    
    root.mainloop()

if __name__ == "__main__":
    main()

このような表示になりクリックが効かなくなる

　bindを使うので別の方法で多重実行を防ぎたいというような場合、ロックを使うのが正攻法だと思います。

import time
import tkinter as tk
import threading

def f(event):
    print("start callback")
    s.set("start callback")
    time.sleep(2)
    print("end callback")
    s.set("end callback")
    time.sleep(2)
    print("(空白に戻す)")
    s.set("")
    button.configure(state=tk.NORMAL)
    lock.release()

def callback(event):
    if lock.acquire(blocking=False):
        button.configure(state=tk.DISABLED)
        th = threading.Thread(target=f, args=(event,))
        th.start()
    else:
        print("(すでに実行中)")

def main():
    global s, button, lock

    lock = threading.Lock()
    root = tk.Tk()

    s = tk.StringVar()
    s.set("")
    label = tk.Label(root, textvariable=s)
    label.pack()

    button = tk.Button(root, text="Button")
    button.bind("<1>", callback)
    button.pack()
    
    root.mainloop()

if __name__ == "__main__":
    main()

　threading.Lock.acquireメソッドで、ロックが獲得されていない場合はロックを獲得してTrueを返します。獲得されていた場合の挙動ですが、第一引数blockingをTrueにしていればロックが解放されるまで待ち（デフォルトの挙動）、FalseにすればすぐにFalseを返してくれます。今回はTrueにしたら台無しなのでFalseにし、ifで動作を切り分けることにしています。

17.1. threading — スレッドベースの並列処理 — Python 3.6.5 ドキュメント

　このようにロックを使うことでも多重実行を防止できます。

まとめ

　何もしなければシングルスレッドで動く、というのはわかっていれば当たり前のことなのですが、それに気づかないと困惑すると思います。

　固まるのを妥協してしまうのも一つの考え方ですが、遅い処理を走らせているときでも画面を動かすのはそんなに難しい処理という訳でもないので、よかったら皆さんやってみてください。

【python】クラスでデコレータ！

はじめに

　デコレータといえば関数で作るものだと思っている人も大勢いると思いますが、クラスでも__call__メソッドを実装すればクラスインスタンスはcallableになり、呼び出しできるので、デコレータたりえます。

import time

class deco:
    def __init__(self, f):
        self.f = f
    def __call__(self, *args, **kwargs):
        t1 = time.time()
        ret = self.f(*args, **kwargs)
        t2 = time.time()
        print("deco:", t2-t1)
        return ret

@deco
def f(message):
    print("f:", message)

f("hoge")
""" =>
f: hoge
deco: 6.508827209472656e-05
"""

　このことには先に書いた記事で気づいたのですが、もう少し追求してみます。

　目次

デコレータの種類ごとに実装してみる

　こちらの記事によると、デコレータは4種類に分類できます。
qiita.com

　以下の2項目の組み合わせによる

• 引数の有無
• ラッパー関数を返すか否か

　ただし、引数なしでラッパー関数を返さないデコレータは簡単すぎて意味がない、としているので3種類考えることにします。

　これらすべてを実装できたらすごいですね。ということで、実装していきます。なお、以下で上記記事からの引用コードを一部示します。

引数なしでラッパー関数を返すデコレータ

def args_logger(f):
    def wrapper(*args, **kwargs):
        f(*args, **kwargs)
        print('args: {}, kwargs: {}'.format(args, kwargs))
    return wrapper

（引用）

　これはクラスのインスタンス自身をラッパー関数として取り扱えば簡単です。冒頭のコードもそれです。もうやったので除外。

引数ありでラッパー関数を返さない場合

funcs = []
def appender(*args, **kwargs):
    def decorator(f):
        # args や kwargsの内容によって処理内容を変えたり変えなかったり
        print(args)
        if kwargs.get('option1'):
            print('option1 is True')

        # 元の関数をfuncsに追加
        funcs.append(f)
    return decorator

（引用）

　簡単そうです。

funcs = []
class deco:
    def __init__(self, cond):
        self.cond = cond

    def __call__(self, f):
        print(self.cond)  # 条件によって処理を変えたり
        funcs.append(f)

@deco("aaa")
def hoge(message):
    print("hoge:", message)

@deco("bbb")
def fuga(message):
    print("fuga:", message)

for f in funcs:
    f("called")

""" =>
aaa
bbb
hoge: called
fuga: called
"""

　できた。

引数ありでラッパー関数を返す場合

def args_joiner(*dargs, **dkwargs):
    def decorator(f):
        def wrapper(*args, **kwargs):
            newargs = dargs + args  # リストの結合
            newkwargs = {**kwargs, **dkwargs}  # 辞書の結合 (python3.5以上で動く)
            f(*newargs, **newkwargs)
        return wrapper
    return decorator

（引用）

　ちょっとむずかしそうですかね。デコレータ自身の引数は__init__で受け取ることになります。次に__call__が関数を受け取るしかないのですが、更にwrapperがないと関数そのものの引数が取り扱えません。

　まあ、wrapperもメソッドで作れるので、実は難しいと見せかけて楽勝なのですが。

class deco:
    def __init__(self, decomessage):
        self.decome = decomessage

    def __call__(self, f):
        self.f = f
        return self.wrapper

    def wrapper(self, *args, **kwargs):
        print("deco:", self.decome)
        return self.f(*args, **kwargs)

@deco("hoge")
def func(x):
    print("func:", x)
func("foo")
    
@deco("fuga")
def func(x):
    print("func:", x)
func("bar")

@deco("piyo")
def func(x):
    print("func:", x)
func("baz")

""" =>
deco: hoge
func: foo
deco: fuga
func: bar
deco: piyo
func: baz
"""

　ということで、ぜんぶできました。よって、クラスもデコレータになれますね。

考察

　思いついたことをつらつらと考察します。

どうしてこの書き方が一般的ではないの？

　単純に__call__を使うのがキモいからだと思います。。。

利点は？

　幾つか考えられます。

• 関数定義のネストが深くならない。
• 通常の関数によるデコレータでは状態をコールスタック上で管理するが（クロージャ）、クラスで書くと明示的にselfの状態として取り扱える。

　以上2点によって、人によっては可読性が高いと感じる可能性があります。

• 各種拡張が行いやすい

　通常のデコレータはいかんせんクロージャなので、あまり大規模なものの実装は避けると思います。この方法だと何しろクラスなので、その気になればなんでも作れるでしょう。まあ、クラスを他に定義しておいて普通のデコレータでそれのインスタンスを返しても良い訳ですが。

　まったく利点がない訳ではないということです。

欠点は？

• とりあえずキモい。
• 更に一般的でもない。
• そもそも濫用しすぎ。
• 決定的な利点がある訳ではないので、上記欠点がむしろ目立つ。

　うーん、苦しいかな。

総評

　面白いけど、これを使ったコードを書いて他人に渡すのは嫌がらせに近いかもしれない。
（いやでも、普通のデコレータもなかなか読みづらいけどなぁ）

まとめ

　とにかくできることはわかりました。「ここが駄目！」とか「こんなふうに活用できるね！」などの意見がありましたら、遠慮なくコメントに書いてください。

追記

　pandasのコードを読んでいたら、実際に使っている例を見かけました。

pandas/_decorators.py at master · pandas-dev/pandas · GitHub

　やるときはけっこうやるんですね。

【python】使いやすい関数の呼び出し回数カウンタを考える

はじめに

　関数の呼び出し回数を数えたい、というシチュエーションはたまにあります。

　その都度、場当たり的にカウンタ変数を増やしたりして対処するのも、まあ、ありといえばありですが、使いやすいものを作るとしたらどうなるかな？　というのを興味本位で書いてみました。

　あくまでも興味本位なので実用的な頑強性までは保証しません。

　なお、ある程度は以前書いた記事の

www.haya-programming.com

　を踏まえた内容になっているので、よろしければこちらも御覧ください。

要件

　「使いやすい」がコンセプトなので、最低限以下のような要件を抑えようと思います。

• 関数内を書き換えない
• 追加が容易
• 関数でもクラスのメソッドなどでも同様に扱える
• カウント回数を取得する、リセットするなどの機能があると嬉しい

　これは独自クラスのインスタンスを返すデコレータですかね。

作ってみる

　とりあえずカウントするためのクラスを作ります。

class Count:
    def __init__(self, func):
        self.count = 0
        self.func = func

    def __call__(self, *args, **kwargs):
        self.count += 1
        return self.func(*args, **kwargs)

    def reset(self):
        self.count = 0

　最低限これだけあれば良いや。カウントは直接count属性を見て取得します。「ゲッターを作れ」「いやそこはプロパティで」といった指摘が飛んできそうですが、python的ではない気がするのと面倒くさいので（こちらが本音）、やりません。また、resetメソッドでカウンタを0にできます。

　クラスなので、必要なら幾らでも処理を追加できます（オーバーヘッドの許す限り）。

　さて、次にデコレータにするために、普通はこんなコードを書いてみるのではないでしょうか。

def count(f):
    return Count(f)

　……要らないんじゃね？　と思ってクラスをそのままデコレータにできないか試したら、いけました。なので、クラス名をcountにしてそのままデコレータとして使うことにしました。

動かしてみる

　最終的にこんなコードになり、望み通りの結果が得られました。

class count:
    """
    カウントのためのクラス
    """
    def __init__(self, func):
        self.count = 0
        self.func = func

    def __call__(self, *args, **kwargs):
        self.count += 1
        return self.func(*args, **kwargs)

    def reset(self):
        self.count = 0

@count
def f(x, y=None):
    """
    実験用関数
    """
    print("f:", x, y)

# 以下で確認
for i in range(2):
    for j in range(3+i):
        f(i, y=j)
        print("count:", f.count)
    f.reset()

""" 結果 =>
f: 0 0
count: 1
f: 0 1
count: 2
f: 0 2
count: 3
f: 1 0
count: 1
f: 1 1
count: 2
f: 1 2
count: 3
f: 1 3
count: 4
"""

　問題なさげ。できました。

まとめ

　以前の記事と見比べると、カウンタそのものを実装するのと比べて大した手間はかかりません。それでいて汎用的です。高階関数が扱えてデコレータのある言語は良いですね（ってそれpythonだけじゃん）。

　また、クラスが直接デコレータに使えるのは新たな発見でした。これを使っているものはあるのだろうか？（クラスだけだと引数を取れないデコレータになるので、いろいろと厳しいかもしれません）　これについてはもう少し追求してみたいです。