np.random以下には色んな乱数生成関数があるのですが、毎回「生成される配列のサイズの指定方法がわからない、なんだっけ？」と思っているので、この際備忘録として残しておきます。
Random sampling (numpy.random) — NumPy v1.16 Manual
を見るとわかりますが、
rand(d0, d1, …, dn) randn(d0, d1, …, dn)
というフォーマットを取る。
ということになっており、大変覚えやすいのでした。
これまでは、
np.random.normal(2,3,4)
とか、
np.random.rand(size=(3,4,5))
とか、あとはなぜか、
np.random.randn(shape=(4,5,6))
みたいに書いたりして怒られるということがしょっちゅうありましたが、この備忘録を書いたことで今後は間違えないで済みそうです。
concurrent.futures.ProcessPoolExecutorは便利そうなので、Poolの代わりに使ってみようと思います。
17.2. multiprocessing — プロセスベースの並列処理 — Python 3.6.5 ドキュメント
17.4. concurrent.futures – 並列タスク実行 — Python 3.6.5 ドキュメント
非同期で使えることを除くと、以下のような違いがあります。
利便性の観点からすると、まあまあよさげです。
とりあえず同等に使えることを確認するため、次のようなコードを書いてみました。
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%くらい？）、とりあえず一応ほぼ互角。
インスタンスメソッドでも行けるかどうか見ます。
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だとメソッドは渡せなかったはずなのですが、最近は渡せるみたいですね。
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は明示的に大きめに指定することです。
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 を送出します。
それだけの違いのようです。
まさかこんなのに速度差なんてないだろう、そう思っていた時期が私にもありました。
>>> 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
例外処理は単純な比較より重いため、少し時間を食います。
どちらが適しているかは処理の意図などにもよるので一概には言えません。
見つからなかったとき、str.findは-1を、str.indexはValueErrorを出すと覚えましょう。
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は「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は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」の場合は以下のようなものです。
対応をまとめると、
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ですが、実はやっていることは単純です。慣れればそれほど困惑することもないので、気楽に使えるようになると良いと思います。
scipyの階層型クラスタリングを使う機会がありましたが、使い方がわかりづらいと思ったのでまとめておきます。
目次
いっぱいあるんですよ。
Hierarchical clustering (scipy.cluster.hierarchy) — SciPy v1.2.0 Reference Guide
私の数え間違えがなければ31個。多いですね。
とはいえ、本質的なもの（実際によく使うもの）は以下くらいです。
実際に階層型クラスタリングを行う。これがないと始まらない。
任意の深さで木を切り、クラスタに分割する。
yを渡すとコーフェン相関係数なる評価指標を出してくれるらしい。
デンドログラムを描画する。
他に、こんな関数があります。
fclusterの処理と絡む便利関数。
クラスタリングアルゴリズムに対応。すべてlinkageの引数で文字列を使って指定できるので実際にこれらの関数を使うことはない。
他にも各種の数値を出せたり、MATLABのフォーマットと相互変換してくれたり、ポインタのリンクで表現された構造化データに変換するクラスだったり、いろいろなものが押し込まれています。必要に応じてリファレンスから探してくれば良いので、それぞれ述べることはしません。
とりあえず上で挙げた実際によく使う関数を中心に説明していきますよ。
データをlinkageに通すことで階層型クラスタリングが行えます。返り値として木の情報を表す配列が返ります。それに対して用意されている関数であれこれ処理していくというのが基本的な流れです。
scipy.cluster.hierarchy.linkage(y, method='single', metric='euclidean', optimal_ordering=Falsescipy.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という変数に代入しているので、それに倣うと良いと思います。
クラスタリングしてくれます。
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に好きなアルゴリズムを選ぶ……という仕組みですね。
「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 日本
あとで書く実践編でちょっと触れます。
デンドログラムを描きます。デンドログラム出しとけばとりあえず納得感があるので、これだけ出してみてからあれこれすれば良いと思います。
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
引数が泡吹いて倒れるほどいっぱいある。到底説明なんてできないので、代表的なものだけピックアップします。
説明不要。linkageの結果を渡す。
木を省略して表示するときのパラメータ。次のtruncate_modeと絡みます。
どのように木を省略するか。これはデータ数がある程度多いときに威力を発揮します。今回は使いません。
色分けに絡むしきい値。木の中の最大ノード間距離との比率で色分けを決めます。
葉のラベル。
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の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です。
上のドキュメントの引用にも書いてある通り、集合演算が使えるためには
という条件があります。考えればわかることですが、それぞれ整理しておくと、
ということになります。
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と交差検証を組み合わせて使うと、たとえ交差検証で高いスコアが出て汎化性能確保できた！ と思っても想像とかけ離れた分離超平面になっていることがままある。
なのでこの組み合わせは少し怖いということを説明する。
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です。
しかし、出力された画像は思っていたものとは違いました。
こわっ。どう考えてもこういうデータではないと思いますが、こういうことが現実に起こります。(-1.5, 0)あたりにデータが来たら、青かねずみ色のどちらかの色のグループに分類されてほしいところですが、実際は緑になってしまう訳です。
ちょっと信頼ならないですね、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はわずかに下がるだけ。
はるかに納得感の高い結果になっています。こういうことがあるので、ほぼ線形分離可能なことがわかっているデータは、まずは線形なモデルで試すことがおすすめです。非線形なモデルで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層の多層パーセプトロンだからこれくらいの素直な結果に落ち着くのであって、深層学習は油断するとすぐ過学習するので注意が必要です。
汎化重視のパラメータにすれば変なことにはならないので安心してください。基本的にCもgammaも低ければ低いほど過学習しづらくなります。
今までと同じ箇所をこう書き換える。
# SVMをfit clf = SVC(C=1, gamma=0.2) clf.fit(X, y)
別に問題ないですよね。SVMそのものはこのようにパラメータで自由に分離超平面の複雑さを調整できる、優れた分類器です。
ただし、交差検証で機械的にパラメータを決めてしまうとあまり良くない結果を招く可能性がある訳です。
二次元で見ているから違和感があるのであって、高次元空間はサクサクメロンパン問題があるのでまた異なった挙動になります。
↓サクサクメロンパン問題
直感的な理解は難しいと思いますが、こういうことが問題になるケースは少ないということはなんとなくわかります。
とはいえ、「交差検証して最高性能だったSVMを投入したらトンチンカンな結果を出してくる」みたいなことが実際に起こらないとも限らないので、ある程度は注意した方が良いでしょう。
SVMはパラメータ設定を過学習するようにしたとき、他の手法と比べても群を抜いて変な結果になるのですが、何が駄目なんでしょうね。やっぱりカーネル使ってるからか。
ちょっと注意が要るよなー、というお話でした。
*1:手持ちのデータの分布に依存する話ですが。未知のデータなんてないぜ！ といえるくらい豊富なデータがあればそれはそれで構いません