k平均法（k-means）アルゴリズムをわかりやすく解説【Pythonコード有】

Contents

k平均法とは
- k平均法のアルゴリズム
- k平均法の利点、欠点
k平均法の実装
- sklearnによるデータセットの作成
エルボー法による最適なkの決定
- エルボー法とはどんな手法なのか
- エルボー法の実装
まとめ

k平均法とは

k平均法(k-means)は機械学習のアルゴリズムの一つです。簡単に説明すると、k-means法とは互いに近いデータ同士は同じクラスタに属するという考えに基づき、データ群をk個に分類するクラスタリング(分類)手法です。これだけではいまいちイメージがわかないと思いますので、より詳しく見ていきましょう。

k平均法のアルゴリズム

では、k-meansのアルゴリズムについて見ていきましょう。以下のようなクラスタリングされていないデータを分類したいとします。

ここで、まずそれぞれの点をランダムにクラス分けします。

次に、各クラスについて、重心を計算します。

そして、各点について、計算した重心からの距離を計算し、距離が一番近いクラスタに割り当てなおします。

重心の計算と割り当てなおしのプロセスを、各点の割り当てられるクラスタが変化しなくなるまで繰り返します。こうすることで、最終的には完全に分類された状態となります。

これがk-meansです。また、k-meansはクラスタリングの手法ですので、教師なし学習に分類されます。

k平均法の利点、欠点

k-meansでは重心からの距離のみを計算すればよいため、計算量が少なく済むという利点があります。一方、最初の重心の位置によっては計算量が変動するため、初期値に依存するというデメリットがあります。また、いくつのクラスタに分類するのか、ということはあらかじめ決めておかなければいけません。

k平均法の実装

sklearnによるデータセットの作成

今回はscikit-learnのmake_blobsを用いて作成したデータを用いることにします。make_blobsを用いることでガウス分布に従って生成されるクラスタリング用のデータ群を作成することができます。

ここではガウス分布中心の数を3つに設定しましたが、実際に分布はだいたい3つにクラス分けできそうであることが見て取れると思います。また、yは正解ラベルですが今回は教師なし学習ですので学習には用いません。

データの形状も確認してみましょう。

データの個数は150個、その特徴量の数は2つに設定してあります。

では、データの準備ができたので早速k-meansを実装していきましょう。

まず、各点のクラスを適当に初期化する処理を実装します。

n_clusters = k #クラスタ数

def init(X):
    #0からクラスタ数-1までの整数からランダムにラベルを振る
    labels_ = np.random.randint(0, n_clusters, X.shape[0]) 
 
    #ラベルの変動確認用の配列
    old_labels = np.ones(X.shape[0]) * -1 

    cluster_centers_ = np.zeros((n_clusters, X.shape[1])) #クラスタの重心を格納する配列
    dist = np.zeros((X.shape[0], n_clusters)) #距離を格納する配列

ランダムなラベルの配列、重心格納用の配列、距離の格納用の配列を初期化します。ここで、old_labelsは次のループ処理においてループを抜ける条件のための配列です。

次に、重心の計算処理、ラベルの更新処理を実装しましょう。ループはラベルの変動がなくなるまで繰り返します。距離は各重心について全ての点との距離を計算するということに注意してください。

def predict(X):
    #ラベルの変動がなくなるまで繰り返す
    while (not(labels_ == old_labels).all() ):
        for i in range(n_clusters):
            X_clusterd = X[labels_ == i,:] #ラベルに応じたXを抜き出す

            #個数方向の平均として抜き出したXの重心を計算
            cluster_centers_[i,:] = X_clusterd.mean(axis = 0) 

        #全ての点について、重心からの距離を計算
        for i, x in enumerate(X):
            for j, center in enumerate(cluster_centers_):
                dist[i][j] = math.sqrt(distance_2(x, center))

        #更新前のラベルを保存      
        old_labels = labels_ 

        #最も重心との距離が近かったクラスにラベルを更新
        labels_ = dist.argmin(axis = 1)

    return labels_

#ユークリッド距離の計算
def distance_2(self, p, q):
    return (p[0] - q[0])**2 + (p[1] - q[1])**2

ここで、空のクラスができてしまったときのことを考慮する必要があります。上のコードでは要素数が0のクラスができたとき、重心が計算されません。ここで場合分けをする必要があります。これらを踏まえて、k-meansモデルの実装は以下のようになります。

import numpy as np
import math

class KMeans2D():  
    def __init__(self, k=3):
        self.n_clusters = k
        #0からクラスタ数-1までの整数からランダムにラベルを振る
        self.labels_ = np.random.randint(0, self.n_clusters, X.shape[0]) 
        
        #ラベルの変動確認用の配列
        self.old_labels = np.ones(X.shape[0]) * -1 

        self.cluster_centers_ = np.zeros((self.n_clusters, X.shape[1])) #クラスタの重心を格納する配列
        self.dist = np.zeros((X.shape[0], self.n_clusters)) #距離を格納する配列
        self.dist_2 =np.zeros(X.shape[0]) #空のクラスタができたときの対策用
        
    def predict(self, X):  
        #ラベルの変動がなくなるまで繰り返す
        while (not(self.labels_ == self.old_labels).all() ):

            for i in range(self.n_clusters):
                X_clusterd = X[self.labels_ == i,:] #ラベルに応じたXを抜き出す

                #個数方向の平均として抜き出したXの重心を計算
                if X_clusterd.size != 0: #空のクラスタじゃないときはそのまま重心計算
                    self.cluster_centers_[i,:] = X_clusterd.mean(axis = 0)
                else :
                    #空のクラスタが存在する場合はそのときの重心から最も遠い点を新たな重心に設定
                    for n, x in enumerate(X):
                        self.dist_2[n] = self.distance_2(x,self.cluster_centers_[i,:])    

                    self.cluster_centers_[i,:] = X[self.dist_2.argmax(axis = 0)]

            #全ての点について、重心からの距離を計算
            for i, x in enumerate(X):
                for j, center in enumerate(self.cluster_centers_):
                   self.dist[i][j] = math.sqrt(self.distance_2(x, center))

            #更新前のラベルを保存      
            self.old_labels = self.labels_ 

            #最も重心との距離が近かったクラスにラベルを更新
            self.labels_ = self.dist.argmin(axis = 1)

        return self.labels_

    def distance_2(self, p, q):
        return (p[0] - q[0])**2 + (p[1] - q[1])**2

では、実際にクラスタリング処理を行ってみましょう。

上手く分類することが確認できたのではないかと思います。上記のプログラムは、点の描画については3クラスまでしか対応していませんが、分類自体はクラスタ数が増えても対応できます。また、3クラスまでは対応しているので、k=2やk=1のときを試してみてください。

エルボー法による最適なkの決定

ここまでの実装ではkの数を事前に決める必要がありました。この最適なkの値はなんとなくの目視や試行錯誤で決定するのもよいですが、エルボー法という手法でも決定することができます。

エルボー法とはどんな手法なのか

エルボー法では、まず以下に示す式に従い、残差平方和(SSEもしくはRSS)という値を算出します。残差平方和は簡単に言うと、データと予測値との差異を評価する尺度です。ということは、低ければ低いほどいいということですよね。残差平方和は以下の式で与えられます。ここでyiはデータ、f(xi)は予測値、nはデータ総数です。