クラスタ分析

今回は、教師なし学習のクラスタリング分析について見ていく
クラスタ分析では、正しい答えが事前にわかっていないデータから隠れた構造を見つけ出すことができる
クラスタリングの目標は、データを自然なグループにまとめる方法を見つける事

k-means法を使った類似度によるオブジェクトのグループ化

k-means法は、産業界でも学術界でも広く使用されている
クラスタリングを利用すれば、類似したオブジェクトをグループにまとめる事ができる

これを利用すれば音楽、映画、商品のレコメンドなどができる様になる

k-means法のメリット
他のクラスタリングアルゴリズムと比べて計算効率が良いところ
実装するのが簡単な所

k-means法のデメリット
クラスタの個数kを指定しなければならない事
kの値が不適切である場合には、クラスタリングの性能が悪くなる


from sklearn.datasets import make_blobs

#クラスタリング用に正規分布でサンプルを作る関数
X, y = make_blobs(n_samples=150, 
                  n_features=2, 
                  centers=3, 
                  cluster_std=0.5, 
                  shuffle=True, 
                  random_state=0)

#作成したサンプルを散布図にする
import matplotlib.pyplot as plt

plt.scatter(X[:, 0], X[:, 1], c='white', marker='o', s=50)
plt.grid()
plt.tight_layout()
#plt.savefig('./figures/spheres.png', dpi=300)
plt.show()

このデータを分類する

k-means法の手続き
1,クラスタの中心の初期値として、サンプル点からk個のセントロイドをランダムに選ぶ
2,各サンプルの最も近いセントロイドに割り当てる
3,セントロイドに割り当てられたサンプルの中心にセントロイドを移動する
4,サンプル点へのクラスタの割り当てが変化しなくなるが、ユーザー定義の許容値またはイテレーションの最大値まで2,3を繰り返す

次はどれだけ類似しているかを測定する
類似度は距離が離れていないほど高い

連続値の特徴量を持つサンプルにおいてはユークリッド距離の二乗が用いられる

j次元の特徴量を持つ点Xと点yの類似度は


d(x,y)^2 = \sum_m^{j=1}(x_j - y_j)^2 = ||x-y||^2_2

実際に最適化していくにあたって、k-means法では、クラスタ内誤差平方和(SSE)を用いる


SSE = \sum_{i=1}^n\sum_{j=1}^kw^{(i,j)}||x^i - μ^j||^2_2

μ:セントロイド
サンプル点がクラスタ内に存在する場合はw = 1,存在しない場合はw = 0

ステップ2,3ではこのSSEを小さくしていく

k-means++法というk-means法を改良したものもある
k-means法では初期のセントロイドをランダムで指定したが
クラスタリングがうまくいかなかったり、収束に時間がかかる事がある

そこで、k-means++では、初期のセントロイドを互いに離れた位置に配置する
これにより、従来のk-meansよりも一貫性のある結果が得られる

k-means法の初期化のステップ
k個のセントロイドを格納するために、空のデータセットMを初期化する
入力サンプルから初期のセントロイドμをランダムに選択し、Mに割り当てる
Mに含まれていないサンプル毎にMのセントロイドに対して最小となるセントロイドを求める
次のセントロイドを選択するためには、各サンプルの重みを等しくした以下の確率分布を用いる


\frac{d(μ^p , M)^2}{\sum_id(x^i , M)^2}

k個のセントロイドが選択されるまで3,4を繰り返す

from sklearn.cluster import KMeans

km = KMeans(n_clusters=3, 
            init='random', 
            n_init=10, 
            max_iter=300,
            tol=1e-04,
            random_state=0)
y_km = km.fit_predict(X)

plt.scatter(X[y_km == 0, 0],
            X[y_km == 0, 1],
            s=50,
            c='lightgreen',
            marker='s',
            label='cluster 1')
plt.scatter(X[y_km == 1, 0],
            X[y_km == 1, 1],
            s=50,
            c='orange',
            marker='o',
            label='cluster 2')
plt.scatter(X[y_km == 2, 0],
            X[y_km == 2, 1],
            s=50,
            c='lightblue',
            marker='v',
            label='cluster 3')
plt.scatter(km.cluster_centers_[:, 0],
            km.cluster_centers_[:, 1],
            s=250,
            marker='*',
            c='red',
            label='centroids')
plt.legend()
plt.grid()
plt.tight_layout()
#plt.savefig('./figures/centroids.png', dpi=300)
plt.show()

このデータセットではうまくグループ化が行われたように思える

しかし、k-means法には注意しなければならない課題がいくつかある
それはクラスタの個数kを指定しなければいけない事である
現実のデータセットではkがわからない事が多い

エルボー法を使ってクラスタの最適な個数を求める

教師なし学習の問題は、正解がわからないことである
そのため、クラスタリングの性能を数値化するには、本章で説明したクラスタ内SSEのような指標を用いて様々なk-meansクラスタリングの性能を比較する必要がある

SSEであればKMeansモデルを適合したあとにはinertia_で確かめられる

print('Distortion: %.2f' % km.inertia_)

Distortion: 72.48

エルボー法では、最も歪みの増え初めるkを特定し、最適な個数とする
ひとまずプロットしてみる


distortions = []
for i in range(1, 11):
    km = KMeans(n_clusters=i, 
                init='k-means++', 
                n_init=10, 
                max_iter=300, 
                random_state=0)
    km.fit(X)
    distortions.append(km.inertia_)
plt.plot(range(1, 11), distortions, marker='o')
plt.xlabel('Number of clusters')
plt.ylabel('Distortion')
plt.tight_layout()
#plt.savefig('./figures/elbow.png', dpi=300)
plt.show()

このグラフからk=3が適当であるとわかる

シルエット図を使ってクラスタリングの性能を数値化する

クラスタリングを評価するもう一つの指標はシルエット分析である
シルエット分析では、クラスタ内のサンプルがどの程度凝集しているかの目安となるグラフをプロットできる

サンプルについてシルエット係数を計算するには以下のステップをふむ
1,サンプルx^iと同じクラスタ内のその他すべてのサンプルとの平均距離として
凝集度a^iを計算する
2,サンプルx^iと最も近くにあるクラスタ内の全てのサンプルとの平均距離として、その次に近いクラスタからの乖離度b^iを計算する
3,クラスタの凝集度と乖離度の差をそれらのうち大きい方の値でわり、シルエット係数sを計算する


s^i = \frac{b^i - a^i}{max(b^i , a^i)}

a=bの時がシルエット係数0となり最も小さい

このシルエット係数を点毎に計算し、その値をグラフにする

import numpy as np
from matplotlib import cm
from sklearn.metrics import silhouette_samples

km = KMeans(n_clusters=3, 
            init='k-means++', 
            n_init=10, 
            max_iter=300,
            tol=1e-04,
            random_state=0)
y_km = km.fit_predict(X)

cluster_labels = np.unique(y_km)
n_clusters = cluster_labels.shape[0]
silhouette_vals = silhouette_samples(X, y_km, metric='euclidean')
y_ax_lower, y_ax_upper = 0, 0
yticks = []
for i, c in enumerate(cluster_labels):
    c_silhouette_vals = silhouette_vals[y_km == c]
    c_silhouette_vals.sort()
    y_ax_upper += len(c_silhouette_vals)
    color = cm.jet(float(i) / n_clusters)
    plt.barh(range(y_ax_lower, y_ax_upper), c_silhouette_vals, height=1.0, 
             edgecolor='none', color=color)

    yticks.append((y_ax_lower + y_ax_upper) / 2.)
    y_ax_lower += len(c_silhouette_vals)

silhouette_avg = np.mean(silhouette_vals)
plt.axvline(silhouette_avg, color="red", linestyle="--") 

plt.yticks(yticks, cluster_labels + 1)
plt.ylabel('Cluster')
plt.xlabel('Silhouette coefficient')

plt.tight_layout()
# plt.savefig('./figures/silhouette.png', dpi=300)
plt.show()

このシルエット図を用いてさまざまなクラスタのサイズを細かく調べる事で、外れ値を含んでいるクラスタを特定できる

次に、クラスタリングがうまくいかなかった場合のシルエット図を見てみる
クラスタの個数kを2つでやってみる

import numpy as np
from matplotlib import cm
from sklearn.metrics import silhouette_samples

km = KMeans(n_clusters=3, 
            init='k-means++', 
            n_init=10, 
            max_iter=300,
            tol=1e-04,
            random_state=0)
y_km = km.fit_predict(X)

cluster_labels = np.unique(y_km)
n_clusters = cluster_labels.shape[0]
silhouette_vals = silhouette_samples(X, y_km, metric='euclidean')
y_ax_lower, y_ax_upper = 0, 0
yticks = []
for i, c in enumerate(cluster_labels):
    c_silhouette_vals = silhouette_vals[y_km == c]
    c_silhouette_vals.sort()
    y_ax_upper += len(c_silhouette_vals)
    color = cm.jet(float(i) / n_clusters)
    plt.barh(range(y_ax_lower, y_ax_upper), c_silhouette_vals, height=1.0, 
             edgecolor='none', color=color)

    yticks.append((y_ax_lower + y_ax_upper) / 2.)
    y_ax_lower += len(c_silhouette_vals)

silhouette_avg = np.mean(silhouette_vals)
plt.axvline(silhouette_avg, color="red", linestyle="--") 

plt.yticks(yticks, cluster_labels + 1)
plt.ylabel('Cluster')
plt.xlabel('Silhouette coefficient')

plt.tight_layout()
# plt.savefig('./figures/silhouette.png', dpi=300)
plt.show()

一つのセントロイドが真ん中に落ちてしまっている


cluster_labels = np.unique(y_km)
n_clusters = cluster_labels.shape[0]
silhouette_vals = silhouette_samples(X, y_km, metric='euclidean')
y_ax_lower, y_ax_upper = 0, 0
yticks = []
for i, c in enumerate(cluster_labels):
    c_silhouette_vals = silhouette_vals[y_km == c]
    c_silhouette_vals.sort()
    y_ax_upper += len(c_silhouette_vals)
    color = cm.jet(float(i) / n_clusters)
    plt.barh(range(y_ax_lower, y_ax_upper), c_silhouette_vals, height=1.0, 
             edgecolor='none', color=color)

    yticks.append((y_ax_lower + y_ax_upper) / 2.)
    y_ax_lower += len(c_silhouette_vals)

silhouette_avg = np.mean(silhouette_vals)
plt.axvline(silhouette_avg, color="red", linestyle="--") 

plt.yticks(yticks, cluster_labels + 1)
plt.ylabel('Cluster')
plt.xlabel('Silhouette coefficient')

plt.tight_layout()
# plt.savefig('./figures/silhouette_bad.png', dpi=300)
plt.show()

シルエット図長さと幅が全く異なるため、最適なクラスタリングとは言えない

クラスタを階層木として構成する

このアルゴリズムの利点は樹形図をプロットできる事である
樹形図をみる事で結果を解釈するのに役立つ

もう一つの利点は、クラスタ数を指定する必要がない事である

階層的クラスタリングには凝集型、分割型の二つがあるが今回は凝集型に着目する
凝集型ではまず個々のサンプルを一つのクラスタとして扱い、クラスタが一つ残った状態になるまで
最も近くにある二つのクラスタをマージしていく

凝集型クラスタリングには単連結方と完全連結方の二つがある
単連結法：最も類似度の高いメンバー同士の距離を計算、最も類似度の高いどうしの距離が最小になるよう二つのクラスタをマージする

完全連結法：最も類似度の低いものの距離を計算する

ここでは完全連結法に着目する

完全連結法のステップとしては
1,全てのサンプルの距離行列を計算する
2,各データ点を単一のクラスタとみなして表現する
3,最も類似度の低いものの距離に基づき、二つの最も近いクラスタをマージする
4,距離行列を更新する
5,クラスタが一つだけ残った状態になるまで、ステップ3,4を繰り返す

まずサンプルを作成する


import pandas as pd
import numpy as np

np.random.seed(123)

variables = ['X', 'Y', 'Z']
labels = ['ID_0', 'ID_1', 'ID_2', 'ID_3', 'ID_4']

X = np.random.random_sample([5, 3])*10
df = pd.DataFrame(X, columns=variables, index=labels)
df

Scipyのspatialdistanceサブモジュールのpdist関数を使用し、距離行列を計算する



from scipy.spatial.distance import pdist, squareform

row_dist = pd.DataFrame(squareform(pdist(df, metric='euclidean')),
                        columns=labels,
                        index=labels)
row_dist

次に完全連結法に基づく凝集型階層的クラスタリングを適用する
これにはScipyのcluster.hierarchyサブモジュールのlinkage関数を使用する


from scipy.cluster.hierarchy import linkage

row_clusters = linkage(df.values, method='complete', metric='euclidean')
pd.DataFrame(row_clusters,
             columns=['row label 1', 'row label 2',
                      'distance', 'no. of items in clust.'],
             index=['cluster %d' % (i + 1)
                    for i in range(row_clusters.shape[0])])

1列目：最も類似度の低い（高い？）メンバー
2列目：最も類似度が低い（高い？）メンバー
3列目：それらのメンバーの距離
4列目：各クラスタのメンバーの個数

実際の樹形図を見てみる

from scipy.cluster.hierarchy import dendrogram

# make dendrogram black (part 1/2)
# from scipy.cluster.hierarchy import set_link_color_palette
# set_link_color_palette(['black'])

row_dendr = dendrogram(row_clusters, 
                       labels=labels,
                       # make dendrogram black (part 2/2)
                       # color_threshold=np.inf
                       )
plt.tight_layout()
plt.ylabel('Euclidean distance')
#plt.savefig('./figures/dendrogram.png', dpi=300, 
#            bbox_inches='tight')
plt.show()

次にヒートマップがでてきますが何に使うんですかねこれ

scikitlearnでもやはり用意されてます
AgglomerativeClusteringというクラス
今回はn_clunsterに2を指定する事でk=2としてクラスタリングしてみる


rom sklearn.cluster import AgglomerativeClustering

ac = AgglomerativeClustering(n_clusters=2, 
                             affinity='euclidean', 
                             linkage='complete')
labels = ac.fit_predict(X)
print('Cluster labels: %s' % labels)

Cluster labels: [0 1 1 0 0]

となり、先ほどの結果と合致している

DBSCANを使って高密度の領域を特定する

Density-based- Spatial Clustering of Application with Noise
ここでいうdensityとは、指定された半径ε以内の点の個数のこと

指定された半径ε以内に少なくとも指定された個数の隣接点があればコア点
指定された半径ε以内に指定された個数の隣接点がなければボーダー点
コア点でもボーダー点でもない点はノイズ点と見なされる

DBSCANのメリット
・k-meansのようにクラスタが球状である事を前提としない
・ノイズ点を除去する能力がある

具体的な例として
半月状の構造を持つデータセットを分類する事で
k-means,階層的クラスタリング,DBSCANを比較してみる


from sklearn.datasets import make_moons

X, y = make_moons(n_samples=200, noise=0.05, random_state=0)
plt.scatter(X[:, 0], X[:, 1])
plt.tight_layout()
# plt.savefig('./figures/moons.png', dpi=300)
plt.show()


f, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 3))

km = KMeans(n_clusters=2, random_state=0)
y_km = km.fit_predict(X)
ax1.scatter(X[y_km == 0, 0], X[y_km == 0, 1],
            c='lightblue', marker='o', s=40, label='cluster 1')
ax1.scatter(X[y_km == 1, 0], X[y_km == 1, 1],
            c='red', marker='s', s=40, label='cluster 2')
ax1.set_title('K-means clustering')

ac = AgglomerativeClustering(n_clusters=2,
                             affinity='euclidean',
                             linkage='complete')
y_ac = ac.fit_predict(X)
ax2.scatter(X[y_ac == 0, 0], X[y_ac == 0, 1], c='lightblue',
            marker='o', s=40, label='cluster 1')
ax2.scatter(X[y_ac == 1, 0], X[y_ac == 1, 1], c='red',
            marker='s', s=40, label='cluster 2')
ax2.set_title('Agglomerative clustering')

plt.legend()
plt.tight_layout()
#plt.savefig('./figures/kmeans_and_ac.png', dpi=300)
plt.show()

k-means:二つのクラスタを識別できていない
階層的クラスタリング:複雑な形状になった

次にDBSCANでやってみる


from sklearn.cluster import DBSCAN

db = DBSCAN(eps=0.2, min_samples=5, metric='euclidean')
y_db = db.fit_predict(X)
plt.scatter(X[y_db == 0, 0], X[y_db == 0, 1],
            c='lightblue', marker='o', s=40,
            label='cluster 1')
plt.scatter(X[y_db == 1, 0], X[y_db == 1, 1],
            c='red', marker='s', s=40,
            label='cluster 2')
plt.legend()
plt.tight_layout()
#plt.savefig('./figures/moons_dbscan.png', dpi=300)
plt.show()

うまくできた

しかしデメリットもある
DBSCANを使用するには二つのハイパーパラメータ(MinPts, ε)を最適化する必要がある
データセットの密度の違いが大きい場合には、効果的なハイパーパラメータを見つけるのは難しいかもしれない

［読書会］Python機械学習プログラミング_11章