圧縮行格納方式(CRS)の疎行列のデータ構造

今回も疎行列のお話です。
前回の記事で登場したcrsとcscについて、具体的にどのようなデータ構造なのかを紹介します。
ちなみにcrsとcscはそれぞれ、
圧縮行格納方式 (Compressed Sparse Row) と、
圧縮列格納方式 (Compressed Sparse Column) の略です。
ほぼ同じ処理を行方向に行うか列方向の違いしかなく、転置を取るとそれぞれ入れ替わるので、 CSRの方を紹介します。

ちなみに、 wikipediaの説明で理解したので、それをみながら記事を書いています。
例として取り上げる行列はこれ。(wikipediaの例と同じ。)
$$
\left(
\begin{matrix}
1 & 2 & 3 & 0 \\
0 & 0 & 0 & 1 \\
2 & 0 & 0 & 2 \\
0 & 0 & 0 & 1 \\
\end{matrix}
\right)
$$

まず、csr形式のデータで作りましょう。
今回はarrayで作って変換するのではなく、お作法にしたがい、lil形式で生成してから変換します。


from scipy import sparse
# 成分が全て0の 4 * 4 の lil形式の疎行列を作成。
M_lil = sparse.lil_matrix((4, 4))
# 各成分を代入。
M_lil[0, 0] = 1
M_lil[0, 1] = 2
M_lil[0, 2] = 3
M_lil[1, 3] = 1
M_lil[2, 0] = 2
M_lil[2, 3] = 2
M_lil[3, 3] = 1
# M_csr形式に変換
M_csr = sparse.csr_matrix(M_lil)
# 確認
print(M_csr)
# 出力
'''
  (0, 0)    1.0
  (0, 1)    2.0
  (0, 2)    3.0
  (1, 3)    1.0
  (2, 0)    2.0
  (2, 3)    2.0
  (3, 3)    1.0
'''

これで、csr形式の変数、M_csrに例の行列が格納されました。
printすると整形されて表示されるのですが、実際のデータ構造はこうはなっていません。
wikipediaの説明と、ドキュメントをみながら確認します。
まず、 実際のデータは、次の3つの属性に格納されています。

data ・・・ CSR format data array of the matrix
indices ・・・ CSR format index array of the matrix
indptr ・・・ CSR format index pointer array of the matrix

具体例を見てから説明します。


print(M_csr.data)
# [1. 2. 3. 1. 2. 2. 1.]
print(M_csr.indices)
# [0 1 2 3 0 3 3]
print(M_csr.indptr)
# [0 3 4 6 7]

まず、data が 疎行列の0では無い要素の値を、左上から行方向(右側へ)に順番に並べたものです。
(csrのrが対応。 cscの場合はここで列方向(下向き)に並べたものになります。)
そして、indices が、それぞれの要素が、何列目の要素なのかを示す配列です。

明らかにわかるように、data と indices の要素の数はその疎行列の0では無い成分の個数です。
あとは、dataの各要素が何行目なのかがわかれば、行列を復元できますが、
それを担っているのが、indptr です。
これだけ、wikipediaの説明と異なっていて非常にわかりにくいですが、次のように解釈できます。


# 行列の最初の行のデータは、indptrの最初の2個のデータで作ったスライスの値
print(M_csr.data[0: 3])
# [1. 2. 3.]
# 次の行のデータは、indptrの一つずらした2個のデータで作ったスライスの値
print(M_csr.data[3: 4])
# [1.]
# 以下繰り返し
print(M_csr.data[4: 6])
# [2. 2.]
print(M_csr.data[6: 7])
# [1.]

明らかにわかる通り、 indptr の要素の個数は行の数より1つ大きくなります。

これで、csr_matrixの中のデータの構造がわかりました。
また、data属性の中に行単位でデータが固まって存在してて、
行単位の取り出しや演算が得意なことにも納得できると思います。

csc_matrixとcsr_matrix

前回の記事に続いて疎行列の話です。
参考:scipyにおける疎行列

今回は疎行列を実装する型(クラス)のうち、csc_matrixcsr_matrixの紹介です。
前回紹介した、lil_matrixには効率的に疎行列を構成できるメリットがありますが、
計算が遅いというデメリットがあります。
ドキュメントに次のように書いてある通り、演算するときはcsrかcsc形式に変換した方が良いです。
arithmetic operations LIL + LIL are slow (consider CSR or CSC)

とりあえず本当に遅いのか時間を測っておきましょう。
今回は 10000*10000の大きめの疎行列を作り、実験します。
(arrayの場合の時間も測定したいので、最初にarrayで作成して変換していますが、本来は最初からlil型で生成してメモリを節約の恩恵を受けた方が良いです。)
また、csrと、cscはそれぞれ行方向、列方向の取り出しに有利という特徴があるので、実験用データによって有利不利無いように対称行列を作りました。


from scipy import sparse
import numpy as np
# 疎行列作成
M_array = np.random.choice(
    [0, 1, 2, 3, 4, 5],
    p=[0.99, 0.002, 0.002, 0.002, 0.002, 0.002],
    size=[10000, 10000]
)
# 自身の転置行列と足し合わせて対称行列にする
M_array = M_array + M_array.T

# それぞれの型に変換する
M_lil = sparse.lil_matrix(M_array)
M_csc = sparse.csc_matrix(M_lil)
M_csr = sparse.csr_matrix(M_lil)

# 和をとるのにかかる時間を計測
%timeit M_array + M_array
# 406 ms ± 2.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit M_lil + M_lil
# 618 ms ± 7.04 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit M_csc + M_csc
# 7.92 ms ± 17.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit M_csr + M_csr
# 7.94 ms ± 51.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

ご覧の通り、 lil 同士の演算は通常のarray形式よりも遅く、一方でcscやcsr同士の演算は圧倒的に速くなりました。

あとは、cscとcsrの違いですが、それぞれ、列方向のスライスと行方向のスライスが効率的に行えます。
(あまり指摘されないようですが、スライスをとるだけであればarrayが一番早い。)

これもそれぞれみておきましょう。


# 列方向のスライス
%timeit M_array[:, 400]
# 225 ns ± 3.06 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
%timeit M_lil[:, 400]
# 11.4 ms ± 201 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit M_csc[:, 400]
# 151 µs ± 2.73 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit M_csr[:, 400]
# 2.86 ms ± 127 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

列方向のスライスについては、array が最速ではありますが、 csc が csr に比べてかなり早いことが確認できます。

次は行方向。


# 行方向のスライス
%timeit M_array[400, :]
# 235 ns ± 10.1 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
%timeit M_lil[400, :]
# 36.6 µs ± 1.37 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit M_csc[400, :]
# 2.9 ms ± 112 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit M_csr[400, :]
# 66 µs ± 3.39 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

今度は、csr が csc より速くなりました。
とはいえ、lil が csrよりも速く処理できています。(array最速は変わらず。)

これは、前回の記事で説明した通り、lilがrows属性に、行ごとの値を保有しているからのようです。

csrの方が演算が速く、わざわざ疎行列を作ったうえでなんの演算もせずにスライスで値を取り出すだけという場面もあまり無いので、
csrかcscを選んで使うことが多いと思いますが、行方向のスライスだけならlilが速いことは一応覚えておこうと思います。

scipyにおける疎行列

自然言語処理やレコメンドシステムの開発などをやっていると、サイズは非常に大きいが成分のほとんどは0という行列が頻繁に登場します。
そのような行列を疎行列と言います。
疎行列 – Wikipedia
これをそのまま扱うと、メモリは無駄に消費するし無駄な計算は増えるしで非常に効率の悪いことになります。
そのため、scipyでは疎行列を専門に扱う scipy.sparse というモジュールが用意されています。

ドキュメント: Sparse matrices (scipy.sparse)

リンク先を見たらわかる通り、疎行列を表す型は結構色々な種類があります。(Sparse matrix classes 参照)
それぞれ、他の型との変換が高速とか、行方向/列方向の取り出しが早いとか個別にメリットを持っていますが、
共通しているのはデータ量を大幅に節約できる点です。

今後の記事でいくつか紹介する予定ですが、とりあえずデータ量を削減できるって特徴だけでも確認しておきましょう。

疎行列のメリットを感じるにはかなり小さいのですが、乱数で10*10の疎行列(array型)を作って、
それを lil_matrix に変換して、中のデータを見てみましょう。

まずデータ作成。
(今回はサンプルとして、array型でデータを作って変換していますが、
省メモリの恩恵をきちんと受けるには最初から疎行列でデータを生成するべきである点には注意してください。)


from scipy import sparse
import numpy as np
M_array = np.random.choice(
    [0, 1, 2, 3, 4, 5],
    p=[0.9, 0.02, 0.02, 0.02, 0.02, 0.02],
    size=[10, 10]
)
M_array

# 出力例
array([[0, 3, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 1, 0, 4, 0, 0, 0, 0],
       [0, 0, 0, 0, 5, 0, 3, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 2, 0, 0],
       [0, 5, 0, 0, 0, 3, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 2],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])

これを型変換します。方法は、strとintの変換と似たような感じで、lil_matrix(D)のようにするとできます。


M_lil = sparse.lil_matrix(M_array)
M_lil
# 出力
<10x10 sparse matrix of type ''
	with 9 stored elements in LInked List format>

これで変換できました。

printすると、0ではない成分の座標と値を保有していることが確認できます。
(ただし、実際のデータ構造は結構違います。print用に整形しているようです。)


print(M_lil)
# 以下出力
  (0, 1)	3
  (1, 3)	1
  (1, 5)	4
  (2, 4)	5
  (2, 6)	3
  (3, 7)	2
  (4, 1)	5
  (4, 5)	3
  (8, 9)	2

0行1列目が3, 1行3列目が1,という風に読んでいくと、確かに0以外の成分が記録されているのがわかります。

なお、実際のデータ構造は型ごとに違うのですが、lil_matrixの場合は次ように、
dataとrowsとうい二つの属性を使ってデータを保有しています。


print(M_lil.data)
'''
[list([3]) list([1, 4]) list([5, 3]) list([2]) list([5, 3]) list([])
 list([]) list([]) list([2]) list([])]
'''

print(M_lil.rows)
'''
[list([1]) list([3, 5]) list([4, 6]) list([7]) list([1, 5]) list([])
 list([]) list([]) list([9]) list([])]
'''

data方は、各行の0ではない要素の値、rowsにそれぞれの要素が何番目に入っているのかの情報を保有しています。
合計18個の数字で、10*10=100個の要素を持つ行列を表現できています。

あと、lil_matrixの重要なメリットして、スライスで値を取り出せるし、代入もできるという点があります。


print(M_lil[0, 1])
# 3

(できて当然に思えますが、他の疎行列のデータ型によってはこれができないものがある。)
そのため、 lil_matrix で疎行列を作って、それを(その他のメリットを持つ)他の型に変換して使うというやり方をよくやります。
(Intended Usageにそう書いてあるので従っています。)

最後に、arrayに戻す方法は toarray()です。


M_lil.toarray()
# 出力
array([[0, 3, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 1, 0, 4, 0, 0, 0, 0],
       [0, 0, 0, 0, 5, 0, 3, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 2, 0, 0],
       [0, 5, 0, 0, 0, 3, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 2],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=int64)

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) した方が便利かもしれません。
結果は同じになります。

scipyの階層的クラスタリングで使える距離関数について

再び階層的クラスタリングの記事の続きです。
参考:scipyで階層的クラスタリング

この記事中では、元のデータの点と点の距離はユークリッド距離を採用しました。
metric='euclidean'と指定しているところがそれです。


# ユークリッド距離とウォード法を使用してクラスタリング
z = linkage(X, metric='euclidean', method='ward')

しかし実際は、これ以外にも多くの距離関数が利用できます。
ドキュメントを読むと、 pdistのページを読めと書いてあるのでそちらをみてみましょう。
scipy.cluster.hierarchy.linkage

metric : str, optional
The distance metric to use. See the distance.pdist function for a list of valid distance metrics.

ちなみに method のほうで指定できる値とその意味は、ちゃんとlinkageのページに載っています。
’single’/’complete’/’average’/’weighted’/’centroid’/’median’/’ward’

scipy.spatial.distance.pdist

metric : str or function, optional
The distance metric to use. The distance function can be ‘braycurtis’, ‘canberra’, ‘chebyshev’, ‘cityblock’, ‘correlation’, ‘cosine’, ‘dice’, ‘euclidean’, ‘hamming’, ‘jaccard’, ‘jensenshannon’, ‘kulsinski’, ‘mahalanobis’, ‘matching’, ‘minkowski’, ‘rogerstanimoto’, ‘russellrao’, ‘seuclidean’, ‘sokalmichener’, ‘sokalsneath’, ‘sqeuclidean’, ‘yule’.

この他、自分で定義した関数も使えるようです。

Y = pdist(X, f)
Computes the distance between all pairs of vectors in X using the user supplied 2-arity function f. For example, Euclidean distance between the vectors could be computed as follows:

dm = pdist(X, lambda u, v: np.sqrt(((u-v)**2).sum()))

それぞれの距離関数はこちらのページから辿っていくと確認しやすいです。
Distance computations (scipy.spatial.distance)

scipyのlinkage関数の結果について

前回の記事の続きです。
参考:scipyで階層的クラスタリング

前回の記事で階層的クラスタリングを実行し可視化するところまで紹介しましたが、
今回は一歩戻ってlinkage関数の戻り値の中身を見てみます。
とりあえず、 linkage matrix をprintして結果を見てみましょう。


from sklearn.datasets import load_iris
from scipy.cluster.hierarchy import linkage
X = load_iris().data[::10, 2:4]
print(X.shape)  # (15, 2)
# ユークリッド距離とウォード法を使用してクラスタリング
z = linkage(X, metric='euclidean', method='ward')
print(z.shape)  # (14, 4)
print(z)
# 以下出力
[[ 2.          3.          0.1         2.        ]
 [ 0.          1.          0.1         2.        ]
 [12.         14.          0.14142136  2.        ]
 [ 4.         16.          0.2081666   3.        ]
 [ 6.          8.          0.31622777  2.        ]
 [ 5.          9.          0.36055513  2.        ]
 [ 7.         11.          0.36055513  2.        ]
 [15.         18.          0.39072582  5.        ]
 [10.         17.          0.43969687  3.        ]
 [13.         23.          0.73598007  4.        ]
 [20.         21.          1.0198039   4.        ]
 [19.         25.          2.00831604  6.        ]
 [24.         26.          3.72312593 10.        ]
 [22.         27.          9.80221064 15.        ]]

前回の記事で可視化したのと同じデータなので、以降の説明は前回の記事中の図と見比べながら読むとわかりやすいと思います。
結果のlinkage matrixは、z.shape の値から分かる通り、14行4列の行列の形をしています。
で、この14という値は、元のデータの個数15個から1減らした値です。
階層的クラスタリングのプロセスの中で、1個ずつグルーピングして集約し、もともと15個あったグループを1つにまとめるのでこうなってます。

そして、列ですが、pythonのインデックスいうと 0列目と1列目はあたらに同じグループに含まれるデータのインデックス、
2列目はそれらの要素orクラスタ間の距離、3列めはそこで新たに作られたクラスタに含まれれるデータの個数を意味します。

具体的には、次のデータは、X[2]とX[3]の距離が0.1で、この二つをまとめて要素が2個のクラスタを作ったことを意味します。
[ 2. 3. 0.1 2. ]
そして、明示はされていませんが、その新しいクラスタには、インデックス15が振られます。(元のデータが0~14の15個なので。)

同様に、次のデータで0と1がまとめられてインデックス16のクラスタが作られます。
[ 0. 1. 0.1 2. ]
で、このインデックス16のクラスタは次のデータで4番目の要素とグルーピングされて、要素数3個のクラスタになります。
[ 4. 16. 0.2081666 3. ]
前回の記事のデンドログラムで確かに0と1でできたクラスタに4が合流しているのが描かれていますね。

このようにして、 linkage matrix の中身を直接読むことができます。

scipyで階層的クラスタリング

今回紹介するのは階層型クラスタリングをscipyで実施する方法です。

階層型クラスタリングの各種アルゴリズム自体は、まだエンジニアを本職にしてたころに知り、その時はこれは面白い手法だと感心していたのですが、
いざデータサイエンティストに転職してからはあまり使ってきませんでした。
クラスタリングする時はk-meansなど、クラスタ数を指定する手法を使うことが多いし、
どれとどれが近いとか言った分析は距離行列眺めたり、次元削減してプロットしたりすることが多かったので。
ただ、他の職種のメンバーに説明するときの樹形図(デンドログラム)のわかりやすさを活用したくなり最近使い始めています

さて、本題に戻ります。

階層型クラスタリングを雑に説明すると、一旦個々のデータを全部別々のクラスタに分類し、
その中から近いものを順番に一つのクラスタにまとめるという操作を繰り返し、最終的に全データを1個のクラスタにまとめる方法です。
この操作を途中で打ち切ることで、任意の個数のクラスタにここのデータを分類することができます。
この際、個々の要素の距離をどう定義するのか、またクラスタと要素、クラスタとクラスタの距離をどのように定義するかによって、
手法が複数存在し、結果も変わります。

この階層型クラスタリングを行う関数が、scipyに用意されています。
Hierarchical clustering (scipy.cluster.hierarchy)

非常に多くの関数がありますが使うのは次の3つです。

scipy.cluster.hierarchy.linkage ・・・ 主役。これが階層型クラスタリングの実装。
scipy.cluster.hierarchy.fcluster ・・・ 各要素について、クラスタリングの結果どのクラスタに属するのかを取得する。
scipy.cluster.hierarchy.dendrogram ・・・ 樹形図(デンドログラム)を描く。

ここまで読んであまりイメージがつかめないと思うので、とりあえずやってみましょう。
データは何でもいいのですが、いつものirisでやります。
(そんなに多くの件数も必要ないので、1/10の件数のデータに絞って使い、4つある特徴量のうち2個だけ使います)


# データ取得。 (15件だけ取得し、特徴量も petal lengthとpetal width に絞る
from sklearn.datasets import load_iris
import matplotlib.pyplot as plt
X = load_iris().data[::10, 2:4]

# データを可視化。
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(1, 1, 1, title="iris (sample)")
plt.scatter(X[:, 0], X[:, 1])
for i, element in enumerate(X):
    plt.text(element[0]+0.02, element[1]+0.02, i)
plt.show()

出力がこちら。これをクラスタリングしていきます。
散布図には番号を振っておきましたが、この番号が結果に出てくる樹形図(デンドログラム)内の番号に対応します。

つぎに階層型クラスタリングを実行して可視化します。


# 階層型クラスタリングに使用する関数インポート
from scipy.cluster.hierarchy import linkage
from scipy.cluster.hierarchy import dendrogram

# ユークリッド距離とウォード法を使用してクラスタリング
z = linkage(X, metric='euclidean', method='ward')

# 結果を可視化
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(1, 1, 1, title="樹形図")
dendrogram(z)
plt.show()

結果がこちら。

0,1,2,3,4 が早い段階で一つのクラスタにまとまっていたり、 12と14、6と8が早々にまとまっていたりと納得性のある形になっています。

あとは、もっと扱いやすい形で、何番のデータが何番目のクラスタに所属するのかのリストを作りましょう。
ここで fcluster 関数を使います。
詳しくはドキュメントにありますが、 criterion に ‘maxclust’を指定して、最大クラスタ数で決めたり、
criterion に’distance’ を指定して、距離で閾値を指定したりできます。

ドキュメントを読むよりやってみたほうがわかりやすいと思うのでそれぞれやってみます。
(この辺の話は専用に別記事を用意して取り上げるかも。)


from scipy.cluster.hierarchy import fcluster

# クラスタ数を指定してクラスタリング
clusters = fcluster(z, t=3, criterion='maxclust')
for i, c in enumerate(clusters):
    print(i, c)

# 以下出力
0 1
1 1
2 1
3 1
4 1
5 3
6 3
7 3
8 3
9 3
10 2
11 3
12 2
13 2
14 2

0,1,2,3,4 と 5,6,7,8,9,11 と 10,12,13,14 の3グループにきちんと別れました。

この他、樹形樹に横線を引いてその位置で分けることもできます。
距離=3くらいで分ければ同じ3グループに分かれるのですが、せっかくなので別のところで切りましょう。
距離=1.7を閾値にして、4グループに分かれるのを確認します。


# 距離の閾値を決めてクラスタリング
clusters1 = fcluster(z, 1.7, criterion='distance')
for i, c in enumerate(clusters1):
    print(i, c)

# 以下出力
0 1
1 1
2 1
3 1
4 1
5 4
6 3
7 4
8 3
9 4
10 2
11 4
12 2
13 2
14 2

クラスタ1や、クラスタ2は先ほどと同じですが、6,8 と 5,7,9,11 が別のクラスタに別れました。
これらが分かれる理由は上の樹形図を見ていただければ理解できると思います。

scipyで手軽にカイ二乗検定

業務でカイ二乗検定を行う場面は、割と多くあります。
自分の場合は、ABテストの評価で使うことが多いです。

Excelでも自分でコードを書いてでもさっとできるのですが、
scipyに専用の関数が用意されているのでその紹介です。

scipy.stats.chi2_contingency

これを使うと、引数にデータのクロス表を放り込むだけで、
カイ二乗値、p値、自由度、期待度数表を一発で算出してくれます。
3つの値と、期待度数がarray形式ではいったタプルが戻るので、引数を4つ用意して受け取ると便利です。


import pandas as pd
from scipy.stats import chi2_contingency

# 検定対象のデータ(サンプルなので値はダミー)
df = pd.DataFrame([[300, 100], [800, 200]])

chi2, p, dof, expected = chi2_contingency(df, correction=False)

print("カイ二乗値:", chi2)
print("p値:", p)
print("自由度:", dof)
print("期待度数:", expected)

# 以下出力
カイ二乗値: 4.242424242424244
p値: 0.03942583262944296
自由度: 1
期待度数: [[314.28571429  85.71428571]
 [785.71428571 214.28571429]]

このようにめんどな手間がなく一瞬で検定ができました。
ちなみに、
correction=False
はイエーツの修正を”行わない”設定です。
デフォルトがTrueなので注意してください。

データから確率分布のパラメーターを推定する

データから、そのデータを生成した背景にある確率分布を推定したいことはよくあります。
正規分布やポアソン分布を仮定するのであれば、簡単ですが、多くの分布では結構面倒です。
そこで、scipyのstatsにある、fitとという便利な関数を使って最尤推定します。

今回はベータ分布を例に取り上げます。
公式ドキュメントはここです。
scipy.stats.rv_continuous.fit
ここ、ベータ関数を使ったサンプルも乗ってるんですよね。
初めて読んだ時はもっと早く読めばよかったと思いました。

それでは、真の分布を設定して、そこからデータを生成し、パラメーターを推定してみます。


# モジュールのインポート
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta

# これから推定したい真の分布
frozen_beta_true = beta.freeze(a=3, b=7, loc=-2, scale=4)
# 真の分布に従うデータを生成
data = frozen_beta_true.rvs(500)

# データから最尤推定 (全パラメーター)
fit_parameter = beta.fit(data)
print(fit_parameter)
# 出力
# (2.548987294857196, 4.380552639785355, -1.946453704152459, 3.1301112690818194)

bの値と、scaleがちょっと乖離が大きいかなと感じられるのですが、
結構妥当な値が推定できました。

経験上、ベータ分布を使いたい時は、取りうる値の範囲が決まっていることが多いです。
そのため、locやscaleは固定して推定を行いたいのですが、
その時は、パラメーターにfをつけて、fitに渡すと、
それらのパラメーターは固定した上で残りを推定してくれます。


# データから最尤推定 (loc と scaleは指定する)
fit_parameter = beta.fit(data, floc=-2, fscale=4)
print(fit_parameter)
# 出力
# (3.1998198349509672, 7.4425425953673505, -2, 4)

かなり真の値に近い結果が出ました。
最後に推定した確率分布の確率密度関数を可視化してみましょう。


# 推定したパラメーターで確率分布を生成
frozen_beta = beta.freeze(*fit_parameter)

# 可視化
plt.rcParams["font.size"] = 14
x = np.linspace(-2, 2, 51)
fig = plt.figure(figsize=(8, 5))
ax = fig.add_subplot(1, 1, 1, xlim=(-2, 2), title="scipyによる最尤推定")
ax.plot(x, frozen_beta_true.pdf(x), label="真の分布")
ax.plot(x, frozen_beta.pdf(x), label="推定した分布")
ax.hist(data, bins=30, alpha=0.7, density=True, label="サンプルデータの分布")
ax.legend()
plt.show()

出力されたのがこちらの図です。
うまく推定されているように見えますね。

pythonを触り始めたばかりの頃は、scipyをうまく使えず、
確率分布はnumpyでスクラッチで書いて、この種の推定もゴリゴリ自分で実装していました。
(かなり効率の悪いアルゴリズムで)
fitを知ってからも、しばらくは4つの戻り値のどれがaでどれがlocなのかよくわからなかったり、
locやscaleを固定する方法を知らず長いこと敬遠していたのですが、
ちゃんとドキュメントを読めば全部書いてあるものです。