コードで理解するk-means

はじめに

今回はk-meansを自前実装することで、理解を深めてみようと思います。

解説

k-meansの動作

概ね、以下のような流れで実装します。

1. K個の重心を決める(ランダムで置く)
2. 各点と各重心の距離を計算する
3. 前回の距離と今回の距離が変化していなかったら更新を終了する
4. 各点から見て、計算した各重心との距離の中から最も近い重心を求め、その点の所属するクラスタの重心をそれとして扱う
5. 2.に戻る

実装

ライブラリのインポート

今回はこのあたりを使っていきます

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn import datasets
from sklearn.decomposition import PCA

データの用意

データはあやめデータを使いますが、どの軸を使うか悩むので、PCAで2次元に落とし込んだものを利用することにします。
PCAで落とし込んだ結果でDataFrameを作成し、追加で本当のlabel列を加えています。

dataset = datasets.load_iris()
pca = PCA(n_components=2)
X = pca.fit_transform(dataset.data)
df = pd.DataFrame(X, columns=["PCA1", "PCA2"])
df["label"] = dataset.target
sns.scatterplot(x="PCA1", y="PCA2", hue="label", data=df)

本体の実装

# 再現性のため、乱数シードを固定
np.random.seed(5)
def kMeans(data, K=3):
    # 重心の初期値を決定
    center = []
    ndim = data.shape[1]
    max_vals = np.max(data, axis=0)
    min_vals = np.min(data, axis=0)
    std_vals = np.std(data, axis=0)
    avr_vals = (max_vals + min_vals) / 2
    for k in range(K):
        c = []
        center.append(c)
        for i in range(ndim):
            c.append(np.random.normal(loc=avr_vals[i], scale=std_vals[i], size=1)[0])
        
    # 更新処理
    distance = [[np.inf for i in range(data.shape[0])] for j in range(K)]
    while True:
        changed = False
        for k in range(K):
            d = []
            for i in range(data.shape[0]):
                d.append(np.linalg.norm(center[k] -data[i,:]))
            if d != distance[k]:
                changed = True
                distance[k] = d
        if not changed:
            break
        cluster = np.argmin(distance, axis=0)
        for k in range(K):
            d = data[cluster == k]
            center[k] = np.sum(d, axis=0) / d.shape[0]
    return (center, cluster)

やってみる

df_result = df.loc[:]
center, df_result["cluster"] = kMeans(df.drop(columns=["label"]).values)
ax = sns.scatterplot(x="PCA1", y="PCA2", hue="cluster", data=df_result)
sns.scatterplot(x="x", y="y", s=100, marker="x", color="blue", ax=ax, data=pd.DataFrame(center, columns=["x", "y"]))

結果は以下のようになります。青いバツがそれぞれの重心です。

アルゴリズムの特性を理解していると、それらしい形になったな、という感覚が持てます。

さらにやってみる（エルボー法）

k-meansで難しいのはkをいくつにするか、というところです。
調べてみたところ、エルボー法という方法を使うことによってある程度いい感じのkを選択できるということがわかりましたので実装してみました。
エルボー法はクラスタ内の重心とクラスタに属するデータの残差の和([tex: sse = \sum{\sum{p - c_i}})の値を見ていき、それが人間の肘のように「カクン」といったところが一番ちょうどいいkであるという感じで選択する方法です。

sses = []
for k in range(1,6):
    center, distance, cluster = kMeans(df.drop(columns=["label"]).values, K=k)
    sse = 0
    for j in range(k):
        sse += np.sum((np.array(distance[j])[cluster==j]) ** 2)
    sses.append(sse)
plt.plot(range(1,len(sses)+1), sses)

結果です。だいたい3くらいが丁度いいのかなぁという印象です。

まとめ

普段何気なしに使ってしまうアルゴリズムも実装してみることで本来の姿が見えてきます。
私としては、（あくまで2次元のデータでですが）こういうロジックだからこういうクラスタリングになる、ということがわかったことでk-meansを少しだけ理解できたように思います。
これからも少しずつこういう形で勉強していこうと思います。