変分混合ガウスモデル(VBGMM:Variational Bayesian Gaussian Mixture)でのクラスタリングを試してみます。

VBGMMは、クラスタ数が不明な場合に有効な手法です。

複数のガウス分布を仮定して、各データがどのガウス分布に所属するのかを決定してクラスタ分析を行います。

MeanShiftとは違って、ベイズ推定に基づいて確率分布を計算しながらクラスタ数や分布の形状を決定します。

VBGMMでクラスタリング

データは前前前回使用したワインの分類データセットを使います。（前前前回記事を参照）

比較のため、最初にk-meansのクラスタリング結果を可視化し(1～5行目)、その後にGMMでクラスタリングを行ってから同じく可視化しています(7～15行目)。

VBGMMモデルはmixture.BayesianGaussianMixture関数で作成し、クラスタ数の上限としてn_componentsに10を指定しています。(8行目)

[Google Colaboratory]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

plt.figure(figsize=(10,3))
plt.scatter(x,y, c=z_km.labels_)
plt.scatter(z_km.cluster_centers_[:,0],z_km.cluster_centers_[:,1],s=250, marker="*",c="red")
plt.suptitle("k-means")
plt.show
 
from sklearn import mixture
vbgmm = mixture.BayesianGaussianMixture(n_components=10, random_state=0)
vbgmm=vbgmm.fit(X_norm)
labels=vbgmm.predict(X_norm)
 
plt.figure(figsize=(10,3))
plt.scatter(x,y, c=labels)
plt.suptitle("vbgmm")
plt.show

[実行結果]

VBGMMで最適なクラスタ数を探す

VBGMMでは、weightsを参照するとクラスタごとの各データ分布を確認することができます。

縦軸がweightsで、横軸がクラスタ番号の棒グラフを作成します。(3行目)

[Google Colaboratory]

1
2
3
4
5

x_tick =np.array([1,2,3,4,5,6,7,8,9,10])
plt.figure(figsize=(10,2))
plt.bar(x_tick, vbgmm.weights_, width=0.7, tick_label=x_tick)
plt.suptitle("vbgmm_weights")
plt.show

[実行結果]

クラスタ数の上限(n_components)を10で実行しましたが、1つも分類されていないクラスタがあります。

クラスタ数の変更

グラフで見ると、だいたい３つ(クラスタ番号１～２と４と７)に分類されているので、n_componentsに3を設定して、もう一度クラスタリングを行います。(1行目)

[Google Colaboratory]

1
2
3
4
5
6
7
8
9
10
11
12
13
14

vbgmm = mixture.BayesianGaussianMixture(n_components=3, random_state=0)
vbgmm=vbgmm.fit(X_norm)
labels=vbgmm.predict(X_norm)
 
plt.figure(figsize=(10,3))
plt.scatter(x,y, c=labels)
plt.suptitle("vbgmm")
plt.show
 
x_tick =np.array([1,2,3])
plt.figure(figsize=(10,2))
plt.bar(x_tick, vbgmm.weights_, width=0.7, tick_label=x_tick)
plt.suptitle("vbgmm_weights")
plt.show

[実行結果]

今回実施したように、VBGMMではweightsを確認して低い割合のクラスタがある場合は、n_componentsをさげて再度クラスタリングを行うという手順でクラスタ数を調整することが可能です。

混合ガウスモデル(GMM:Gaussian Mixture Model)でのクラスタリングを試してみます。

GMMは各データが、どのガウス分布に所属している確率が最も高いかを求めてラベリングを行います。

GMMでクラスタリング

データは前前回使用したワインの分類データセットを使います。（前前回記事を参照）

比較のため、最初にk-meansのクラスタリング結果を可視化し(1～5行目)、その後にGMMでクラスタリングを行ってから同じく可視化しています(7～15行目)。

GMMモデルはmixture.GaussianMixture関数で作成し、クラスタ数(n_components)には3を指定しています。(8行目)

[Google Colaboratory]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

plt.figure(figsize=(10,3))
plt.scatter(x, y, c=z_km.labels_)
plt.scatter(z_km.cluster_centers_[:,0], z_km.cluster_centers_[:,1], s=250, marker="*", c="red")
plt.suptitle("k-means")
plt.show

from sklearn import mixture
gmm = mixture.GaussianMixture(n_components=3, covariance_type="full")
z_gmm = gmm.fit(X_norm)
z_gmm = z_gmm.predict(X_norm)
 
plt.figure(figsize=(10,3))
plt.scatter(x, y, c=z_gmm)
plt.suptitle("gmm")
plt.show

[実行結果]

k-meansとGMMでのクラスタリング結果はほとんど同じように見えます。

ただk-meansでは同心円状での分類であり、GMMは傾いた楕円形の分類であるために多少の違いがあるようです。

またGMMではガウス分布を仮定するので、各データがどのクラスタに所属するのかという確率を求めることができます。

今回はX-meansというクラスタリングを試します。

前回試したMeanShiftと同様に、X-meansもクラスタ数を自動で設定してくれます。

X-meansは、k-meansのアルゴリズムに加えて、あるクラスタが正規分布２つで表されるのと１つで表されるのとでどちらが適切かを判定し、２つが適切な場合はクラスタ数を２つに分けるというアルゴリズムになります。

X-meansでクラスタリング

データは前回使用したワインの分類データセットを使います。（前回記事を参照）

まずx-meansライブラリをインストールします。（1行目）

クラスタ数の初期値として2を指定（5行目）していますが、適切なクラスタ数でクラスタリングを実施してくれます。

kmaxパラメータには最大クラスタ数を指定しています。（6行目）

最後にprocess関数をコールしてX-meansでクラスタリングを実行します。（7行目）

[Google Colaboratory]

1
2
3
4
5
6
7

!pip install pyclustering

from pyclustering.cluster.xmeans import xmeans
from pyclustering.cluster.center_initializer import kmeans_plusplus_initializer
xm_c = kmeans_plusplus_initializer(X_norm, 2).initialize()
xm_i = xmeans(data=X_norm, initial_centers=xm_c, kmax=20, ccore=True)
xm_i.process()

[実行結果]

クラスタリング結果の可視化

クラスタリング結果を可視化します。

[Google Colaboratory]

1
2
3
4
5
6
7
8

z_xm = np.ones(X_norm.shape[0])
for k in range(len(xm_i._xmeans__clusters)):
    z_xm[xm_i._xmeans__clusters[k]] = k+1
plt.figure(figsize=(10,3))
plt.scatter(x,y, c=z_xm)
centers = np.array(xm_i._xmeans__centers)
plt.scatter(centers[:,0],centers[:,1],s=250, marker="*",c="red")
plt.show

[実行結果]

クラスタ数を指定せずに、期待するクラスタ数3となっていることが確認できました。

MeanShiftを使ったクラスタリングを行います。

MeanShiftはクラスタ数が分からない場合に、データをクラスタを分類する方法です。

複数のガウス分布（正規分布）を仮定して、各データがどのガウス分布に所属するのかを決定し、クラスタ分析を行います。

ワインの分類データセット

まず、オープンデータであるワインの分類データセットを準備します。

[Google Colaboratory]

1
2
3
4

df_wine_all = pd.read_csv("https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data", header=None)
df_wine = df_wine_all[[0,10,13]]
df_wine.columns = [u"class", u"color", u"proline"]
pd.DataFrame(df_wine)

今回はワインの品種(0列目)、色(10列目)、プロリン量(13列目)を使用します。

[実行結果(一部略)]

データセットを可視化

抽出したワインのデータセットを可視化します。

[Google Colaboratory]

1
2
3
4
5
6
7
8
9

X = df_wine[["color","proline"]]
sc = preprocessing.StandardScaler()
X_norm = sc.fit_transform(X)
x = X_norm[:,0]
y = X_norm[:,1]
z = df_wine["class"]
plt.figure(figsize=(10,3))
plt.scatter(x,y, c=z)
plt.show

[実行結果]

３品種ごとに色分けされたデータが確認できます。

k-meansでクラスタリング

k-meansでクラスタリングを行います。

品種ごとに３つに分類したいのでクラスタ数は3に設定します。

[Google Colaboratory]

1
2
3
4
5
6
7

km = cluster.KMeans(n_clusters=3)
z_km = km.fit(X_norm)

plt.figure(figsize=(10,3))
plt.scatter(x,y, c=z_km.labels_)
plt.scatter(z_km.cluster_centers_[:,0],z_km.cluster_centers_[:,1],s=250, marker="*",c="red")
plt.show

[実行結果]

中心点から同心円状に広がって分類されていることが分かります。

MeanShiftでクラスタリング

MeanShift関数(1行目)を使って、MeanShiftでクラスタリングを行います。

[Google Colaboratory]

1
2
3
4
5
6
7
8
9
10
11

ms = cluster.MeanShift(seeds=X_norm)
ms.fit(X_norm)
labels = ms.labels_
cluster_centers = ms.cluster_centers_
print(cluster_centers)

plt.figure(figsize=(10,3))
plt.scatter(x,y, c=labels)
plt.plot(cluster_centers[0,0], cluster_centers[0,1], marker="*",c="red", markersize=14)
plt.plot(cluster_centers[1,0], cluster_centers[1,1], marker="*",c="red", markersize=14)
plt.show

[実行結果]

MeanShiftでは２つに分類されて実際の分類とはだいぶ違う結果となりました。

MeanShiftではクラスタ数を指定しなくてもクラスタリングが実施可能なので、パラメータはseeds（乱数シード）のみ設定しています。（1行目）

MeanShiftはk-meansをベースとして近いクラスタをまとめていき、既定の距離より近くなったクラスタをまとめて１つにします。

今回は、SpectralClusteringでクラスタリングを行います。

SpectralClusteringはk-meansとは異なり、データ密度でクラスタを作成するため、同心円状になっていないデータもうまくクラスタリングを行うことができます。

k-meansではうまくクラスタリングできなくても、SpectralClusteringであればきちんとクラスタリングできるケースを確認していきます。

ムーンデータの生成

scikit-learnのサンプルデータセットであるムーンデータと呼ばれる三日月状のデータセットを生成します。

また、データ加工用・可視化用のライブラリもインポートしておきます。

[Google Colaboratory]

1
2
3
4
5
6
7
8
9
10

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn import cluster, preprocessing
from sklearn import datasets

X,z = datasets.make_moons(n_samples=200, noise=0.05, random_state=0)
sc = preprocessing.StandardScaler()
X_norm = sc.fit_transform(X)
display(X_norm)

[実行結果（途中略）]

ムーンデータの可視化

生成したデータを可視化します。

[Google Colaboratory]

1
2
3
4
5

x = X_norm[:,0]
y = X_norm[:,1]
plt.figure(figsize=(10,3))
plt.scatter(x,y, c=z)
plt.show

[実行結果]

三日月の形をしたデータが上下に２つあることが分かります。

この２つにクラスタリングができれば良さそうです。

k-meansでクラスタリング

クラスタ数を２に設定してk-meansでクラスタリングを行ってみます。(1行目)

[Google Colaboratory]

1
2
3
4
5
6
7

km = cluster.KMeans(n_clusters=2)
z_km = km.fit(X_norm)

plt.figure(figsize=(10,3))
plt.scatter(x,y, c=z_km.labels_)
plt.scatter(z_km.cluster_centers_[:,0],z_km.cluster_centers_[:,1],s=250, marker="*",c="red")
plt.show

[実行結果]

k-meansでは中心点からの距離でクラスタが決まるため、同心円状に広がっていないムーンデータではうまく分類できていないようです。

SpectralClusteringでクラスタリング

今度は、クラスタ数を同じく２に設定してSpectralClusteringでクラスタリングを行ってみます。(1行目)

パラメータのaffinity(親和性)は、クラスタリングを実施する際に作成するグラフ行列の作成法を設定します。

グラフ行列とは、「データがそれぞれどのようにつながっているか」を表す行列です。

設定している“nearest_neighbors”は、最近傍法を意味し「あるデータに対し、もっとも近いデータが属するクラスタに分類する」ことを意味します。

[Google Colaboratory]

1
2
3
4
5
6

spc=cluster.SpectralClustering(n_clusters=2, affinity="nearest_neighbors")
z_spc=spc.fit(X_norm)

plt.figure(figsize=(10,3))
plt.scatter(x,y, c=z_spc.labels_)
plt.show

[実行結果]

k-means法とは違い、データ密度でクラスタリングを作成するため、上下の三日月ごとにうまくクラスタリングできていることが分かります。

今回は群平均法を使ってクラスタリングをします。

群平均法は、クラスタ同士のすべての点同士の距離の平均をクラスタの距離とする方法です。

データセットの準備

前回と同様にscikit-learnのアヤメのデータセットを読み込みます。

比較しやすくするために、サンプル数を１０分の１に減らします。(4行目)

最後に散布図に表示します。(10行目)

[Google Colaboratory]

1
2
3
4
5
6
7
8
9
10

from sklearn.datasets import load_iris
import matplotlib.pyplot as plt

X = load_iris().data[::10]
fig = plt.figure(figsize=(6, 3))
ax = fig.add_subplot(1, 1, 1, title="iris")
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()

[実行結果]

群平均法とウォード法の比較

群平均法とウォード法の両方でクラスタリングを行い、それぞれの樹形図を表示します。

群平均法ではパラメータmethodに“average”を指定(1行目)し、ウォード法では“ward”を指定(7行目)します。

[Google Colaboratory]

1
2
3
4
5
6
7
8
9
10
11

Z = linkage(X, method="average", metric="euclidean")
fig2, ax2 = plt.subplots(figsize=(6,3))
ax2 = dendrogram(Z)
fig2.suptitle("average")
fig2.show()

Z = linkage(X, method="ward", metric="euclidean")
fig2, ax2 = plt.subplots(figsize=(6,3))
ax2 = dendrogram(Z)
fig2.suptitle("ward")
fig2.show()

[実行結果]

距離の違いはありますが、同じクラスタリング結果となりました。

群平均法は、鎖効果や拡散現象を起こさないため、ウォード法と同様に用いられることが多い手法です。

今回は最長距離法を使ってクラスタリングをします。

最長距離法は完全リンク法（complete_linkage）と呼ばれることもあり、各クラスタにおいて一番遠い点同士の距離をクラスタの距離とする方法です。

データセットの準備

前回と同様にscikit-learnのアヤメのデータセットを読み込みます。

比較しやすくするために、サンプル数を１０分の１に減らします。(4行目)

最後に散布図に表示します。(10行目)

[Google Colaboratory]

1
2
3
4
5
6
7
8
9
10

from sklearn.datasets import load_iris
import matplotlib.pyplot as plt

X = load_iris().data[::10]
fig = plt.figure(figsize=(6, 3))
ax = fig.add_subplot(1, 1, 1, title="iris")
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()

[実行結果]

最長距離法とウォード法の比較

最長距離法とウォード法の両方でクラスタリングを行い、それぞれの樹形図を表示します。

最長距離法ではパラメータmethodに“complete”を指定(1行目)し、ウォード法では“ward”を指定(7行目)します。

[Google Colaboratory]

1
2
3
4
5
6
7
8
9
10
11

Z = linkage(X, method="complete", metric="euclidean")
fig2, ax2 = plt.subplots(figsize=(6,3))
ax2 = dendrogram(Z)
fig2.suptitle("complete")
fig2.show()

Z = linkage(X, method="ward", metric="euclidean")
fig2, ax2 = plt.subplots(figsize=(6,3))
ax2 = dendrogram(Z)
fig2.suptitle("ward")
fig2.show()

[実行結果]

クラスタ数３として比較すると、まとまっている位置が変わっていますが同じように分類されているようです。

最長距離法は計算量が少なく、分類感度も比較的高いですが、拡散現象が生じることがあります。

拡散現象とはクラスターが大きくなるにつれ、他のデータと最長距離を多く持つようになり、次のクラスターの形成の候補に選ばれにくくなる現象です。

前回までは、ウォード法でクラスタリングを行ってきましたが、今回は最短距離法を使ってクラスタリングをします。

データセットの準備

まずscikit-learnのアヤメのデータセットを読み込みます。

比較しやすくするために、サンプル数を１０分の１に減らします。(4行目)

最後に散布図に表示します。(10行目)

[Google Colaboratory]

1
2
3
4
5
6
7
8
9
10

from sklearn.datasets import load_iris
import matplotlib.pyplot as plt

X = load_iris().data[::10]
fig = plt.figure(figsize=(6, 3))
ax = fig.add_subplot(1, 1, 1, title="iris")
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()

[実行結果]

最短距離法とウォード法の比較

最短距離法とウォード法の両方でクラスタリングを行い、それぞれの樹形図を表示します。

最短距離法ではパラメータmethodに“single”を指定(1行目)し、ウォード法では“ward”を指定(7行目)します。

[Google Colaboratory]

1
2
3
4
5
6
7
8
9
10
11

Z = linkage(X, method="single", metric="euclidean")
fig2, ax2 = plt.subplots(figsize=(6,3))
ax2 = dendrogram(Z)
fig2.suptitle("single")
fig2.show()

Z = linkage(X, method="ward", metric="euclidean")
fig2, ax2 = plt.subplots(figsize=(6,3))
ax2 = dendrogram(Z)
fig2.suptitle("ward")
fig2.show()

[実行結果]

結果はほとんど変わっていないようですが、5,13,7,11の分類が少し異なっています。

最短距離法のメリットは計算量が少ないことです。

ただ鎖効果（ある１つのクラスタに対象が１つずつ順番に吸収されながらクラスタが形成されること）のために、クラスタが帯状になり分類感度が低いことがデメリットになります。

樹形図(デンドログラム)の表示

前回実行した階層クラスタリングの結果を、樹形図(デンドログラム)に表示します。

[Google Colaboratory]

1
2
3
4
5

from scipy.cluster.hierarchy import dendrogram
import matplotlib.pyplot as plt
fig2, ax2 = plt.subplots(figsize=(20,5))
ax2 = dendrogram(Z)
fig2.show()

[実行結果]

2,3,4,0,1が１つのクラスタにまとまっていたり、10,12,14がまとまっていたり、5,9がまとまっていることなどが分かります。

クラスタ数を指定

fcluster関数を使うと、何番のデータがどのクラスタに所属するのかを確認することができます。

パラメータcriterionに“maxclust”を指定し、クラスタ数が3の場合の、各データの所属データを確認します。

[Google Colaboratory]

1
2
3
4

from scipy.cluster.hierarchy import fcluster
clusters = fcluster(Z, t=3, criterion="maxclust")
for i, c in enumerate(clusters):
    print(i, c)

[実行結果]

樹形図の通りに３グループに分かれていることが確認できました。

距離を指定

次は、パラメータcriterionに“distance”を指定し、距離で閾値を指定してみます。

樹形図に横線を引いてその位置で分けるようなイメージになります。

樹形図の縦軸に表示されている数値が各点の距離になります。

距離（縦軸に対応）に1.6を指定して実行してみます。

[Google Colaboratory]

1
2
3

clusters1 = fcluster(Z, 1.6, criterion="distance")
for i, c in enumerate(clusters1):
    print(i, c)

[実行結果]

今度はクラスタ数4に分かれることが確認できました。

樹形図をみると、縦軸の8から横線を引くとクラスタ数が2に、縦軸の3から横線を引くとクラスタ数が3に、縦軸の1.6から横線を引くとクラスタ数が4になることが一目瞭然に確認できますね。

階層クラスタリングは、これまでの非階層クラスタリングとは違い最も距離が近くて似ている（類似度が高い）組み合わせからまとめていく方法です。

結果として出力される樹形図から、分類の過程でできるクラスタがどのように結合されていくかをひとつずつ確認できるため、最適なクラスタ数を後から決めることができます。

樹形図で視覚的に解釈がしやすいというメリットがある反面、非階層クラスタリングよりも計算量が多くなる傾向があるので結果がでるまでに時間がかかることがあります。

データセットの準備

階層クラスタリング用のデータセットを準備します。

まずscikit-learnのアヤメのデータセットを読み込みます。

解釈がしやすくなるようにサンプル数を１０分の１に減らし、４つある説明変数も２つに減らします。(4行目)

選択したデータを最後に散布図に表示します。(10行目)

[Google Colaboratory]

1
2
3
4
5
6
7
8
9
10

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

[実行結果]

散布図に表示されている番号は、樹形図に出てくる番号に対応します。

（樹形図は次回表示します。）

階層クラスタリングの実行

scipyを使って階層クラスタリングを行います。

4行目のmethodには、ウォード法(ward)を設定しています。

ウォード法(ward)は、あるクラスタ同士が結合すると仮定したとき、結合後の全てのクラスタにおいて、クラスタの重心とクラスタ内の各点の距離の２乗和の合計が最小となるようにクラスタを結合させていく方法です。

同じく4行目のmetricには元のデータの点と点の距離の定義で、今回はユークリッド距離(euclidean)を指定しています。

ユークリッド距離とは、2点間の直線距離のことです。

[Google Colaboratory]

1
2
3
4
5

import pandas as pd
from scipy.cluster.hierarchy import linkage

Z = linkage(X, method="ward", metric="euclidean")
pd.DataFrame(Z)

[実行結果]

結果に出力したデータフレームZの意味は次の通りです。

• １～２列目
  結合されたクラスタの番号
• ３列目
  クラスタ間の距離
• ４列目
  結合後に新しくできたクラスタの中に入っている元のデータの数

次回は、今回の結果を樹形図に表示します。