scikit-learn機械学習㉚k-means実践編

前回は、k-meansの仕組みを紹介しました。

今回は、scikit-learnのKMeansを使った実験を行います。

このようのデータを生成し、k-meansでクラスタリングを行います。

また、最適なクラスタ数$${k}$$を決めるために、エルボー法やシルエット係数を試してみます。

さらに、k-means が苦手なケースとして複雑な境界線を持つデータを試します。

そして、k-means に対する外れ値の影響も調査します。

どんな結果になるでしょうか。さっそく始めましょう。

Python環境の準備

Pythonの仮想環境を作ってscikit-learnとJupyterなど必要なライブラリをインストールします。

# 作業フォルダを作って、そこへ移動する
mkdir kmeans_test
cd kmeans_test

# 仮想環境を作り、アクティベートする
python3 -m venv venv
source venv/bin/activate

# pip をアップグレードしておく
pip install --upgrade pip

# 必要なライブラリをインストール
pip install numpy scikit-learn jupyter matplotlib 

Jupyterノートブックを立ち上げてPython３のノートブックを作成してください。Jupyterノートブックに関しては、こちらを参照してください。

なお、ノートブックの閲覧編集にはVSCode（Visual Studio Code）も使えます。私は、どちらかというとVSCodeをよく使います。これについてもこちらで簡単に解説しています。

データの準備

まずは、クラスタリングをするためのデータを生成します。

import numpy as np
from sklearn.datasets import make_blobs

# データをランダムに生成
X, y = make_blobs(n_samples=300, centers=4, cluster_std=0.60, random_state=0)

• 300のデータポイントを生成（n_samples=300）

• クラスタの中心となる点は４つ（centers=4）

• クラスタ内の分散（標準偏差）は0.6（cluster_std=0.60）

• 乱数シードを固定（random_state=0）

生成したデータを表示します。

import matplotlib.pyplot as plt

# データをプロット
plt.figure(figsize=(8, 6))
plt.scatter(X[:, 0], X[:, 1], s=50)
plt.title("Generated Data")
plt.xlabel("Feature 1")
plt.ylabel("Feature 2")
plt.grid(True)
plt.show()

k-meansの初期化

from sklearn.cluster import KMeans

# k-means クラスタリング
kmeans = KMeans(init='k-means++', n_clusters=4, random_state=0)

init='k-means++' では、「k-means++」による初期化ロジックを指定しています。これはデフォルトなので指定する必要はありませんが、ここでは明示的に指定しました。

なお、n_clusters=4 となっているのが、クラスタ数$${k}$$を指定しているところです。

クラスタリングの実行

クラスタリングを実行します。

# クラスタリングの実行
y_kmeans = kmeans.fit_predict(X)

kmeans.fit_predict(X) でデータを４つのクラスタに分け、それぞれのデータポイントのクラスタ番号をリスト（NumPyの配列）にして返します。

y_kmeans

なお、kmeans.fit(X)とすると、クラスタの計算だけ行います。その場合、クラスタ番号のリストを取得するには、kmeans.labels_ を参照します。

クラスターを色分け表示

クラスタリングの結果に従ってデータポイントを色分けします。また、セントロイドの位置を表示します。

import matplotlib.pyplot as plt

# クラスタリングのセントロイド
centers = kmeans.cluster_centers_

plt.figure(figsize=(8, 6))

# データポイントをクラスタごとに色分け
plt.scatter(X[:, 0], X[:, 1], c=y_kmeans, s=50, cmap='viridis')

# セントロイドをプロット
plt.scatter(centers[:, 0], centers[:, 1], 
            c='red', s=200, alpha=0.75, marker='x')

plt.title("K-Means Clustering Results")
plt.xlabel("Feature 1")
plt.ylabel("Feature 2")
plt.grid(True)
plt.show()

クラスタ数 k をエルボー法で選ぶ

前回紹介したエルボー法を使って、クラスタ数$${k}$$の適切な値を選んでみましょう。

まず、クラスタ数を 2 から 10 へと変化させながら、WCSS（Within-Cluster Sum of Squares）を取得します。

# kのリスト
k_list = range(2, 11)

# WCSSを計算
wcss = []
for k in k_list:
    # フィットする
    kmeans = KMeans(n_clusters=k, random_state=0)
    kmeans.fit(X)
    
    # 計算されてあるWCSSを取得
    wcss.append(kmeans.inertia_)  

WCSSは、サンプルから最も近いクラスター中心（セントロイド）までの距離の二乗の合計。inertia_ という属性から取得することができます。

結果をグラフとしてプロットします。

# クラスタ数kとWCSSのプロット
plt.figure(figsize=(12, 8))
plt.plot(k_list, wcss, marker='o')
plt.title("k vs WCSS")
plt.xlabel("k")
plt.ylabel("WCSS")
plt.grid(True)
plt.tight_layout()
plt.show()

$${k=4}$$以降は、WCSSの下がり具合が鈍くなっています。よって、$${k=4}$$がエルボー（肘ひじ）です。

シルエット係数（スコア）を使う

今度は、シルエット係数（シルエットスコア、Silhouette Score）をプロットして最適な$${k}$$を見極めます。

シルエット係数は、クラスタリングの質を測る指標で、その値の範囲は$${-1}$$から$${+1}$$です。大きいほど質が良いことになります。

まずは、シルエット係数をクラスタ数ごとに計算します。

from sklearn.metrics import silhouette_score

# シルエットスコアを格納するリスト
silhouette_scores = []
for k in k_list:
    # フィットする
    kmeans = KMeans(n_clusters=k, random_state=0)
    kmeans.fit(X)

    # シルエットスコアを計算
    silhouette_scores.append(silhouette_score(X, kmeans.labels_))

結果をグラフとしてプロットします。

# クラスタ数kとシルエットスコアのプロット
plt.figure(figsize=(12, 8))
plt.plot(k_list, silhouette_scores, marker='o', color='orange')
plt.title("k vs Silhouette Score")
plt.xlabel("k")
plt.ylabel("Silhouette Score")
plt.grid(True)
plt.tight_layout()
plt.show()

$${k=4}$$のところが山の頂点となって逆エルボーになっています。

k-meansが苦手なケース

複雑な境界

k-meansは、クラスタが球状で均一な分布を持つ場合に強みを発揮しますが、複雑な境界を扱うのが苦手です。これを確認します。

まず、月型のデータをランダムに生成します。

from sklearn.datasets import make_moons

# 月型のデータをランダムに生成
X_moons, y_moons = make_moons(n_samples=300, noise=0.1, random_state=0)

結果をプロットします。

import matplotlib.pyplot as plt

# 生成されたデータをプロット
plt.figure(figsize=(8, 6))
plt.scatter(X_moons[:, 0], X_moons[:, 1], s=50)
plt.title("Generated Data")
plt.xlabel("Feature 1")
plt.ylabel("Feature 2")
plt.grid(True)
plt.show()

これは２つのグループなので、$${k=2}$$でフィットしてみます。

from sklearn.cluster import KMeans

# k-means クラスタリング（k=2）
kmeans_moons = KMeans(n_clusters=2, random_state=0)

y_kmeans_moons = kmeans_moons.fit_predict(X_moons)

この結果をプロットします。

import matplotlib.pyplot as plt

# クラスタリングのセントロイド
centers_moons = kmeans_moons.cluster_centers_

plt.figure(figsize=(8, 6))

# データポイントをクラスタごとに色分け
plt.scatter(X_moons[:, 0], X_moons[:, 1], 
            c=y_kmeans_moons, s=50, cmap='viridis')

# セントロイドをプロット
plt.scatter(centers_moons[:, 0], centers_moons[:, 1], 
            c='red', s=200, alpha=0.75, marker='x')

plt.title("K-Means Clustering Results")
plt.xlabel("Feature 1")
plt.ylabel("Feature 2")
plt.grid(True)
plt.show()

これはちょっとまずいですね。k-meansクラスタリングは適用するデータを間違えると役に立ちません。

外れ値の影響

外れ値の有無によって k-meansクラスタリングの結果がどう変わるのかを確認します。

まずは、外れ値がないデータを生成します。

import numpy as np
from sklearn.datasets import make_blobs

# 正常なデータを生成
X, y = make_blobs(n_samples=70, centers=2, cluster_std=0.60, random_state=0)

次に、外れ値を追加したものを用意します。

# 外れ値を追加
outliers = np.array([ [11, 6], [10, 7] ])

X_with_outliers = np.vstack([X, outliers])

これらのデータをプロットして確認します。

#　データのプロット
fig, axes = plt.subplots(2, 1, figsize=(8, 8))

# 外れ値なしのプロット
axes[0].scatter(X[:, 0], X[:, 1], s=50)
axes[0].set_title("Generated Data (No Outliers)")
axes[0].set_xlim(-1, 12)
axes[0].set_ylim(-1, 8)
axes[0].grid(True)

# 外れ値ありのプロット
centers_outliers = kmeans_with_outliers.cluster_centers_

axes[1].scatter(X_with_outliers[:, 0], X_with_outliers[:, 1], s=50)
axes[1].set_title("Generated Data (with Outliers)")
axes[1].set_xlim(-1, 12)
axes[1].set_ylim(-1, 8)
axes[1].grid(True)

この二つのデータセットに対して別々にk-meansクラスタリングを実行します。

from sklearn.cluster import KMeans

# 正常なデータでk-means クラスタリングを実行
kmeans_normal = KMeans(n_clusters=2, random_state=0)
y_kmeans_normal = kmeans_normal.fit_predict(X)

# 外れ値を含むデータでk-means クラスタリングを実行
kmeans_with_outliers = KMeans(n_clusters=2, random_state=0)
y_kmeans_with_outliers = kmeans_with_outliers.fit_predict(X_with_outliers)

この結果を並べてプロットします。

# 結果のプロット
fig, axes = plt.subplots(2, 1, figsize=(8, 8))

# 外れ値なしのプロット
centers_normal = kmeans_normal.cluster_centers_

axes[0].scatter(X[:, 0], X[:, 1], 
                c=y_kmeans_normal, cmap='viridis', s=50)
axes[0].scatter(centers_normal[:, 0], centers_normal[:, 1], 
                c='red', s=200, alpha=0.75, marker='x')
axes[0].set_title("K-Means Clustering (No Outliers)")
axes[0].set_xlim(-1, 12)
axes[0].set_ylim(-1, 8)
axes[0].grid(True)

# 外れ値ありのプロット
centers_outliers = kmeans_with_outliers.cluster_centers_

axes[1].scatter(X_with_outliers[:, 0], X_with_outliers[:, 1], 
                c=y_kmeans_with_outliers, cmap='viridis', s=50)
axes[1].scatter(centers_outliers[:, 0], centers_outliers[:, 1], 
                c='red', s=200, alpha=0.75, marker='x')
axes[1].set_title("K-Means Clustering (with Outliers)")
axes[1].set_xlim(-1, 12)
axes[1].set_ylim(-1, 8)
axes[1].grid(True)

外れ値を含むと判断されたクラスタの中心点が右にずれているのが確認できます。

また、あのハズレ値は、どちらのクラスタにも属さないようにも見えます。$${k=3}$$だとどうなるでしょうか。

外れ値がないデータは、もともと２つのクラスタに対して$${k=3}$$を指定しているので、余計なクラスタが発生しました。一方、外れ値がある方は３つに分かれており、$${k=3}$$が最善に見えます。

エルボー法で確認します。

# kのリスト
k_list = range(2, 11)

# WCSSを計算
wcss = []
for k in k_list:
    kmeans = KMeans(n_clusters=k, random_state=0)
    kmeans.fit(X_with_outliers)
    wcss.append(kmeans.inertia_)  # クラスターごとにセントロイドからの距離の二乗の和


# クラスタ数kとWCSSのプロット
plt.figure(figsize=(12, 8))
plt.plot(k_list, wcss, marker='o')
plt.title("k vs WCSS")
plt.xlabel("k")
plt.ylabel("WCSS")
plt.grid(True)
plt.tight_layout()
plt.show()

やはり、この場合は$${k=3}$$を選択するのがベストです。

以上より、k-meansを使う時には、あらかじめ外れ値を取り除くべきかどうかはよく検討して決める必要があります。何が外れ値なのかは扱う問題によるので一意には決められないからです。

次回予告

k-meansは単純明快で計算量も比較的少ない便利なアルゴリズムです。しかし、$${k}$$を選択する必要があったり、境界の形状や外れ値などによってはクラスタリングの結果が悪くなる弱点があります。

そこで次回は、あらかじめクラスタ数を指定する必要がなく、データの分布に応じて自動的にクラスタが決定するDBSCANを紹介します。複雑な形状のクラスタや外れ値に対してもk-meansより柔軟な対応が可能です。

お楽しみに！