k-meansクラスタリングをnumpyで実装してみる

k-means

k-meansとは

k-meansとは非階層型クラスタリングの手法です。
クラスタ重心との距離情報を使うことで、あらかじめ決めたクラスタ数に分けていきます。
よく使われる手法ですが、アルゴリズム自体は結構単純です。

1. ランダムにクラスタ重心を決める
   あらかじめ決めたクラスタ数分、ランダムに重心となる座標を決定します。
   

   

2. 各データのラベルを更新
   すべてのデータに対して各クラスタ重心との距離を計算し、最も近いクラスタのラベルを割り当てます。 距離は何を使っても大丈夫ですが、一般的にはユークリッド距離が使われることが多いみたいです。
   

3. クラスタ重心を更新
   各クラスタの重心を、それぞれのラベルが割り当てられたデータの平均に更新します。
   

4. 2, 3を収束するまで繰り返す
   ラベル更新とクラスタ重心更新を何度も繰り返すことで近いデータに同じラベルが割り当てられます。 このとき、どれくらい繰り返し処理を行うかは実装次第ですが、一般的には「あらかじめ決められた繰り返し回数処理を行った」や「ラベルの更新が行われなくなった」あたりを条件にすることが多いようです。
   

numpyでの実装

scikit-learnとかを使えば簡単に使うこともできるのですが、勉強も兼ねてnumpyを使って自分で実装してみます。
実装は「機械学習のエッセンス」を参考にさせてもらいました。（特にデータとクラスタ重心の距離をまとめて計算する方法は参考になりました。）

import numpy as np
import matplotlib.pyplot as plt

class KMeans:
    def __init__(self, n_clusters, max_iter=100, random_state=0):
        self.n_clusters = n_clusters
        self.max_iter = max_iter
        self.cluster_centers = None
        np.random.seed(random_state)
    
    def fit_predict(self, X):
        # 1. クラスタ中心をランダムに決定する
        n = X.shape[0]
        pr = np.repeat(1/n,n)
        self.cluster_centers = X[np.random.choice(np.arange(n), self.n_clusters, p=pr, replace=False), :]

        labels = np.zeros(X.shape[0])
        pre_labels = np.ones(X.shape[0])
        iter_counter = 0
        while iter_counter < self.max_iter:
            if (labels==pre_labels).all():
                break

            # 2. 各データのラベルを最も近いクラスタのラベルに変更する
            dists = ((X[:, :, np.newaxis]-self.cluster_centers.T[np.newaxis, :, :])**2).sum(axis=1)
            pre_labels = labels
            labels = np.argmin(dists, axis=1)

            # 3. クラスタ重心を更新する
            for c in range(self.n_clusters):
                cluster_X = X[labels==c, :]
                self.cluster_centers[c, :] = cluster_X.mean(axis=0) 
            
            iter_counter +=  1
        return labels


#データの生成
np.random.seed(123)
x1 = np.r_[np.random.normal(size=20,loc=1,scale=2),np.random.normal(size=20,loc=8,scale=2),np.random.normal(size=20,loc=15,scale=2),np.random.normal(size=20,loc=25,scale=2)]
x2 = np.r_[np.random.normal(size=20,loc=15,scale=2),np.random.normal(size=20,loc=1,scale=2),np.random.normal(size=20,loc=20,scale=2),np.random.normal(size=20,loc=0,scale=2)]
X = np.c_[x1,x2]

# kmeansでのクラスタリング
kmeans = KMeans(n_clusters=4)
labels = kmeans.fit_predict(X)

#可視化
plt.scatter(X[:,0],X[:,1],cmap="jet", c=labels, s=10,alpha=0.5)
plt.show()

今回は4つの塊になるようなデータを人工的に生成して試してみました。
結果はこちらの通りです。

ちなみにクラスタ重心の初期値は、全データにランダムにラベルを割り振ってその重心を初期値にするのが本来のやり方のようなのですが、 その方法で初期化をすると更新時にどのデータにも割り当てられないラベルが発生して上手く動作しないことがありました。 なので今回はランダムにデータを選択して、それをクラスタ重心の初期値にする方法を取っています。（実装方法が悪いのかもしれませんが...）

実装のとき少し躓いた点としては、初期クラスタ重心を選択するときに同じデータが初期重心として選ばれて上手く動作しないことがありました。
np.random.choiceの引数にreplace=Falseを追加して重複を許さないようにすることで対策しています。