Python 機械学習 Scikit-learnによるクラスタリング分析(k-means法)の実践

Pythonでk-means法を用いてクラスタリングを実行します。とりあえず動かせるように、サンプルコードを用意しています。

※Pythonの導入方法(http://costudyinfodatabase.nagoya/2017/01/05/%e3%81%82-2/

今回は、irisのデータセットから三次元データ(3つの特徴量を有するデータ)をインポートし、それに対して、scikit-learnにあるk-means法を適用してクラスタリング分析を実行します。可視化のために3次元データを使っていますが、原理上、次元に制約はありません。
また、クラスタ数の最適化の可視化手段であるエルボー法を用いてみます。

 

◆クラスタリング・k-means法とは?

クラスタリングとは教師なし学習法の一つです。データの集合体から似た者同士をグルーピングしていきます。教師あり学習法と異なり、クラスタリングされた結果の解釈は後から行うことになります。

k-means法は、クラスタリングの手法の一つで、計算効率が良く、産業界・学術界で広く用いられています。映画・音楽等のリコメンデーションエンジンなどに応用されています。

クラスタリング手法は大枠以下のように分類できます。

  1. プロトタイプベース クラスタリング
  2. 階層的クラスタリング
  3. 密度ベースクラスタリング

k-means法はこのうち1のプロトタイプベース クラスタリングにあたります。非階層的クラスタリングとも言います。以下に載せている手法では、クラスタ毎にセントロイドという代表点を設けて、サンプル各点のセントロイドに対する近さでグルーピングを行います。具体的には

  • Step1 クラスタ数kを決定し、その代表点であるk個のセントロイドを全サンプルからランダムに抽出。
  • Step2 各サンプルを最も距離が近いセントロイドに割り当てる。(各サンプルが所属するクラスタを仮決め)
  • Step3 各クラスタのサンプルの重心にセントロイドを移す。

セントロイドの位置が動かなくなるまで2-3を繰り返し、セントロイドの位置、およびクラスタを決定するという流れです。

 

◆サンプルコード

実際に動かしてみましょう。

まずはscikit-learn標準のIrisのデータを分類することを念頭に、データをインポートし可視化します。

Irisのデータセットは「Iris setosa」「 Iris versicolor」「 Iris virginica 」の三種のアヤメの花について、sepal length(がくの長さ), sepal width(がくの幅), petal length (花弁の長さ)そして petal width(花弁の幅) の測定値が特徴量として格納されている150個 のデータセットです。

今回は3つの特徴量、sepal width(がくの幅), petal length (花弁の長さ)そして petal width(花弁の幅)のみ抽出します。

## data import
from distutils.version import LooseVersion as Version
from sklearn import __version__ as sklearn_version
from sklearn import datasets
import numpy as np

iris = datasets.load_iris()
names=iris.target_names
data1=1
data2=2
data3=3

X = iris.data[:, [data1,data2,data3]]#sepal width/petal length/petal length
y = iris.target


print('Class labels:', np.unique(y))

if Version(sklearn_version) < '0.18':
    from sklearn.cross_validation import train_test_split
else:
    from sklearn.model_selection import train_test_split

##Split the data into trainning data and test data at a certain ratio
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=0)

from sklearn.preprocessing import StandardScaler

sc = StandardScaler()
sc.fit(X_train)
X_train_std = sc.transform(X_train)
X_test_std = sc.transform(X_test)

#==========================================================
#data plot
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

fig = plt.figure()
ax = Axes3D(fig)
ax1=ax.scatter3D(np.ravel(X[y==0,0]),np.ravel(X[y==0,1]),np.ravel(X[y==0,2]),c='b', marker='o')
ax2=ax.scatter3D(np.ravel(X[y==1,0]),np.ravel(X[y==1,1]),np.ravel(X[y==1,2]),c='r', marker='o')
ax3=ax.scatter3D(np.ravel(X[y==2,0]),np.ravel(X[y==2,1]),np.ravel(X[y==2,2]),c='y', marker='o')
ax.set_xlabel(iris.feature_names[data1])
ax.set_ylabel(iris.feature_names[data2])
ax.set_zlabel(iris.feature_names[data3])
plt.legend((ax1, ax2, ax2),
           (names[0], names[1], names[2]),
           scatterpoints=1,
           loc='upper left',
           ncol=3,
           fontsize=8)
plt.show()

#========================================================

・実行結果

次に、k-means法を用いてクラスタリングし、可視化します。

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)

#Predicted data plot
fig = plt.figure()
ax = Axes3D(fig)
#ax.plot_wireframe(X,Y,Z)
#==============================================================================
# #original data 
#ax1=ax.scatter3D(np.ravel(X[y==0,0]),np.ravel(X[y==0,1]),np.ravel(X[y==0,2]),c='b', marker='o')
#ax2=ax.scatter3D(np.ravel(X[y==1,0]),np.ravel(X[y==1,1]),np.ravel(X[y==1,2]),c='r', marker='o')
#ax3=ax.scatter3D(np.ravel(X[y==2,0]),np.ravel(X[y==2,1]),np.ravel(X[y==2,2]),c='y', marker='o')
 #ax.scatter3D(np.ravel(X1),np.ravel(X2),np.ravel(X3),c='b', marker='x',label='1')
#==============================================================================
ax_km1=ax.scatter3D(np.ravel(X[y_km==0,0]),np.ravel(X[y_km==0,1]),np.ravel(X[y_km==0,2]),c='y', marker='x')
ax_km2=ax.scatter3D(np.ravel(X[y_km==1,0]),np.ravel(X[y_km==1,1]),np.ravel(X[y_km==1,2]),c='r', marker='x')
ax_km3=ax.scatter3D(np.ravel(X[y_km==2,0]),np.ravel(X[y_km==2,1]),np.ravel(X[y_km==2,2]),c='b', marker='x')
ax_km4=ax.scatter3D(km.cluster_centers_[:, 0],km.cluster_centers_[:, 1],km.cluster_centers_[:, 2],c='lightgreen', marker='s')
ax.set_xlabel(iris.feature_names[data1])
ax.set_ylabel(iris.feature_names[data2])
ax.set_zlabel(iris.feature_names[data3])
plt.legend((ax_km1, ax_km2, ax_km3,ax_km4),
           (names[0], names[1], names[2],"Centroid"),
           scatterpoints=1,
           loc='upper left',
           ncol=3,
           fontsize=8)
plt.show()

 

・実行結果

以下のように、ほぼ元データと同様の分類ができています。

また、黄緑色でセントロイドもプロットしています。

(誤分類率を算出してみましょう)

 

以上の方法ではクラスタ数は自分で設定することになりますが、最適なクラスタ数はどう評価するでしょうか?

簡単な手法であるエルボー法で評価してみましょう。振り分けるクラスタ数を横軸、縦軸にクラスタ内誤差平方和(SSE)をプロットし、SSEが急降下/飽和する点が最適なクラスタ数というわけです。

クラスタ内誤差平方和とは、クラスタ内のサンプルとセントロイドとの距離の総和であり、「クラスタ全体のまとまり度合」のようなものです。

 

##Elbow method###
#Elbow method is to find the saturation point ("Elbow") with changing the cluster numbers
print("Elbow method")
print('Distortion: %.2f' % km.inertia_)
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()


・実行結果

3個のクラスタで十分であることがわかります。