薬剤師のプログラミング学習日記

プログラミングやコンピュータに関する記事を書いていきます

k-meansによるクラスタリング

k-meansはデータを自動的にクラスタリング(グループ化)する手法で、k平均法ともいわれます。同じクラスタ内のデータは類似するように、別のクラスタのデータは異なるようにデータを分類します。ナイーブベイズなどでは正解を与えて「教師あり」で分類を行いましたが、k-meansクラスタリングは正解ラベルのないデータで行う「教師なし」学習のアルゴリズムです。


k-meansのアルゴリズム

  1. ランダムにクラスタ重心を決める
  2. クラスタ重心からの距離が近い方にデータを分類する
  3. 分類されたデータの平均にクラスタ重心を設定する
  4. 新たに分類の更新がなくなるまで2と3を繰り返す

2と3のステップを繰り返す、単純なアルゴリズムです。
2次元データでアルゴリズムの動きを可視化すると次の図のようになります。
f:id:enokisaute:20200907002237g:plain
この例では2次元データですが、軸(次元)が複数あっても同じようにk-meansを適用することができます。

注意点としては、1のステップの最初にランダムに決めた点によって結果が変わることがある、ということです。このような事態を避けるために、複数回のランダム初期化を行って最良のものを選択する方法や、最初の点の決め方そのものを工夫するといった方法があります。

人間がパッとみてわからないような、きれいにクラスタに分かれていないようなデータに対しても最初にクラスタ数を与えることで複数のクラスタにデータを分けることができます。その反面、データをいくつに分けるかのクラスタ数の選択は人間が行わなくてはいけません。このへんは記事の後半でみていこうと思います。

k-meansアルゴリズムをPythonで実装する

まずは2次元のデータセットを使ったコードを書いていきます。重心の移動とクラスタリングの過程を可視化する上図のものです。
なお、コードはcourseraの課題を参考にしました。

import sys
import random
import scipy.io
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
 
# k-meansクラスタリング
 
# 各データ点から最も距離が近い重心のインデックスを割り当てる
def find_closest_centroids(X, centroids):
    K = len(centroids)
    m = len(X)
    idx = np.zeros((m, 1))
    for i in range(len(X)):
        min = sys.maxsize
        for j in range(K):
            # データiと重心jの距離を計算
            distance = np.linalg.norm(X[i, :] - centroids[j, :])
            if distance < min:
                min = distance
                idx[i, :] = j
    return idx
 
# クラスタ毎のデータ点の平均を計算し重心を求める
def compute_centroids(X, idx, K):
    n = X.shape[1]
    centroids = np.zeros((K, n))
    for i in range(K):
        c = np.where(idx == i)
        centroids[i, :] = np.mean(X[c[0], :], axis=0)
    return centroids
 
# 各データ点を色分けしてプロット
def plot_data_points(X, idx, K):
    x = X[:, 0]
    y = X[:, 1]
    cmap = plt.get_cmap("hsv")
    for i in range(K):
        # 重心iのインデックスの要素番号を抽出
        c = np.where(idx == i)
        plt.scatter(x[c[0]], y[c[0]], marker='o', color=cmap(float(i) / K), alpha=0.2)
 
# 重心の前回点からの移動(線分)をプロット
def draw_line(p1, p2):
    plt.plot([p1[0], p2[0]], [p1[1], p2[1]], color='black', linewidth=1)
 
# クラスタと重心の変化をプロットする. 描画は単純に重ねていくだけ.
def plot_progress_kmeans(X, centroids, previous, idx, K, i):
    plot_data_points(X, idx, K)
    plt.scatter(centroids[:, 0], centroids[:, 1], marker='x', c='black', linewidth=2)
    for j in range(K):
        draw_line(centroids[j, :], previous[j, :])
    plt.title('Iteration number ' + str(i + 1))
    plt.pause(0.5)  # 0.5秒sleep
 
def run_kmeans(X, init_centroids, max_iters, K, plot_progress):
    centroids = init_centroids
    prev_centroids = init_centroids
    for i in range(max_iters):
        idx = find_closest_centroids(X, centroids)
        if plot_progress:
            plot_progress_kmeans(X, centroids, prev_centroids, idx, K, i)
        prev_centroids = centroids
        centroids = compute_centroids(X, idx, K)
    if plot_progress:
        plt.show()
    return centroids
 
# 重心の初期化. XからランダムにK個抽出する
def kmeans_init_centroids(X, K):
    randidx = np.array(random.sample(range(len(X)), K))
    centroids = X[randidx, :]
    return centroids
 
 
if __name__ == '__main__':
    load_data = scipy.io.loadmat('./save/ex7data2.mat')
    data = load_data['X']
 
    K = 3             # クラスタ数
    max_iters = 10    # ループの実行回数
 
    # init_centroids = np.array([[3, 3], [6, 2], [8, 5]])   # デモの重心の初期値
    init_centroids = kmeans_init_centroids(data, K)
    centroids = run_kmeans(data, init_centroids, max_iters, K, True)

    print(centroids)    # 重心の最終地点

プログラムは実質アルゴリズムの2と3のステップが実装できればほぼ完成で、残りは描画部分のコードです。上で示したアルゴリズムに沿って書いていきます。

ステップ1 ーランダムにクラスタ重心を決めるー

kmeans_init_centroids()関数で、与えられたデータの中からランダムにクラスタ数(K個)分の点を選んで返すようにしています。

ステップ2 ークラスタの割り付けー

このステップでは各データに対してクラスタの割り付けを行います。データを各クラスタ重心との距離が最小であるクラスタに割り振っていきます。
コード中では、find_closest_centroids関数でデータ点毎に各重心centroidsとの距離を計算し、最も小さい距離である重心のインデックスを割り振っていきます。

ステップ3 ー重心の移動ー

ここでは重心の移動を行います。同じ重心に割り付けされた全てのデータの平均を計算して、それを新たな重心とします。
コードではcompute_centroids関数がこれを行います。np.where()を使って同じ重心を持つデータ点を抽出し、そのデータ点の平均を取っています。

この2と3のステップをrun_kmeans内で一定回数ループで繰り返しています。

使用するデータについて

適当に作成してもかまいませんが、上のプログラムと同じものを使用したい場合は「ex7data2.mat」などでネットで検索していただくと見つかると思います。
ロード後の変数load_data['X']には次のような2次元データが入っています。

[[ 1.84207953  4.6075716 ]
 [ 5.65858312  4.79996405]
 [ 6.35257892  3.2908545 ]
・・・


何度かプログラムを実行すると、下図のように重心の初期値によって結果が変わることも確認できます。
f:id:enokisaute:20200906161404p:plain

k-meansで画像を圧縮する

これもcourseraにあった課題です。画像ファイルの各ピクセルに対してk-meansを行い、何千とある色の数をクラスタ数のK個に減らす(圧縮する)というものです。
各ピクセルには次のように赤、緑、青の強度値を表す0~255の数値が入っています。今回はWindowsのpaintを使って編集したため、4列目にアルファ値(透明度)が入っていますが、減色には関係ないのでこの列はk-means実行前に削除しておきます。

[[[ 18  22  17 255]
  [ 13  21  18 255]
  [ 11  23  20 255]
  ...


実行結果がこちらです。それぞれ、元画像と16色と8色へ減色した画像です。
f:id:enokisaute:20200906224032p:plain

コードはこちら。関数はすべて上のものを使用しています。ここではPNG画像で行っていますが、JPEG画像でも可能です。ただし、データサイズが大きければそれだけプログラムの処理時間が長くなります。

if __name__ == '__main__':
 
    # 色に対してk-meansを適用して画像を圧縮する.
    img_org = Image.open('./save/fish.png')
    img_array = np.array(img_org)
    # print(img_array.shape)
    # print(img_array)
 
    K = 16            # クラスタ数
    max_iters = 10    # ループの実行回数
 
    # windowsのpaintで編集した場合は4チャンネル(アルファ情報あり)なので, 4列目の軸を削除しておく. jpegの場合は必要ない
    img_array = np.delete(img_array, [3], axis=2)
    w = img_array.shape[0]
    h = img_array.shape[1]
 
    img_array = img_array.reshape(w * h, 3)     # ピクセル数×3に変換
    img_array = img_array / 255     # 0~1の範囲になるよう255で割る
    # k-means実行
    init_centroids = kmeans_init_centroids(img_array, K)
    centroids = run_kmeans(img_array, init_centroids, max_iters, K, False)
    idx = find_closest_centroids(img_array, centroids)
 
    # クラスタリングされたK個のRGB値を元のピクセルに戻す
    for i in range(K):
        c = np.where(idx == i)
        img_array[c, :] = centroids[i, :]
    img_recovered = img_array.reshape(w, h, 3)
 
    # 減色(圧縮)した画像を保存
    img_rec = Image.fromarray(np.uint8(img_recovered))
    img_rec.save('./save/fish_recoverd.png')
 
    # 元画像とk-means適用後の画像の表示
    plt.subplot(1, 2, 1)
    plt.imshow(img_org)
    plt.title('Original')
    plt.subplot(1, 2, 2)
    plt.imshow(img_recovered)
    plt.title('Compressed, with ' + str(K) + ' colors')
    plt.show()


多次元データにk-meansを適用する

多次元の特徴量を持つデータ例として、ここではsklearnのワインデータセットを用います。
データ数は178、次のような13の特徴量を持つデータセットです。

>>> from sklearn.datasets import load_wine
>>> wine = load_wine()
>>> wine.data.shape
(178, 13)
>>> wine.feature_names
['alcohol', 'malic_acid', 'ash', 'alcalinity_of_ash', 'magnesium', 'total_phenols', 'flavanoids', 'nonflavanoid_
phenols', 'proanthocyanins', 'color_intensity', 'hue', 'od280/od315_of_diluted_wines', 'proline']
>>> wine.target
array([0, 0, 0, ... 1, 1, 1, ...2, 2, 2])

特徴量の名称としては”アルコール度数”や"色の強さ"、”色合い”などがあり、これらを数値化したものが入っています。また、正解ラベル(target)は0, 1, 2の3クラスのデータセットとなります。
詳しくはこちらをご参照ください。
sklearn.datasets.load_wine — scikit-learn 0.24.0 documentation
このデータセットにはもともとラベルがあり3クラスということはわかっていますが、ここでは多次元のデータに対してもk-meansがうまく機能するということをみてみたいと思います。

if __name__ == '__main__':

    from sklearn.datasets import load_wine
    from sklearn.preprocessing import StandardScaler

    # ワインデータセットの読み込み
    wine = load_wine()

    # 各特徴量のスケールを合わせておく
    scaler = StandardScaler()
    data = scaler.fit_transform(wine.data)
 
    K = 3              # クラスタ数
    max_iters = 20    # ループの実行回数
 
    # 13の特徴量すべてを使ってk-meansを実行する
    init_centroids = kmeans_init_centroids(data, K)
    centroids = run_kmeans(data, init_centroids, max_iters, K, True)
    idx = find_closest_centroids(data, centroids)
 
    # ラベルで色分けして元データをプロットする
    cmap = plt.get_cmap("hsv")
    for i in set(wine.target):
        cls = np.where(wine.target == i)
        plt.scatter(data[cls, 0], data[cls, 1], marker='o', color=cmap(float(i) / K), label='class_' + str(i))
    plt.legend()
    plt.show()

だいたいどの程度正解ラベルと一致しているかを見るために2つの図を並べてみました。
左は13の特徴量すべてを使って上のプログラムk-meansを行い、最初と2番目の軸(alcoholとmalic acid)だけを使ってプロットしたもの、右図は同じ軸を使いデータを正解ラベルで色分けしてプロットしたものです。*1(クラスタ重心の×印は13の特徴量を使って更新されるので、色分けと合わないところもあるかと思います。)
f:id:enokisaute:20200906182828p:plain
ざっくりとですが、K=3を与えた場合は正解ラベルに近い形でクラスタリングできているのがわかります。
通常、特徴量の多いデータを可視化する場合は主成分分析などで次元を削減することが多いと思いますが、今回はすべての特徴量を減らすことなく使ったのでこのようにしました。

どのようにしてクラスタ数を選択するか?

最初の方でデータをいくつに分けるかのクラスタ数の選択は人間が行わなくてはいけないと書きましたが、そもそもデータにいくつのクラスタがあるのかは曖昧な場合も多く、明確な正解があるとも限りません。
ここでは最適なクラスタ数を決める際の指標となる手法のひとつであるエルボー法について書きます。
k-meansもデータ点と重心との2乗距離を最小化するという、最適化のための目的関数を持っています。エルボー法は、各クラスタ数におけるこの値(SSE:sum of squared errors of prediction)をプロットすることで最適と思われるクラスタ数を見つけよう、という手法です。
sklearnのKMeansクラスでは、inertia_というアトリビュートでこのSSEを取得することができます。

ここでは、「正しい」クラスタの数がわかっているデータに対して、エルボー法でうまくクラスタ数を見つけられるか試してみます。
まずはK=4のクラスタから。右側の図をみると横軸が4のところで肘(ひじ)のように折れ曲がっている点があります。エルボー法では、この点のクラスター数である4を最適なものとして選びます。
f:id:enokisaute:20200906223224p:plain
次ははっきりしたクラスタを持たないデータを想定しました。右図では折れ曲がったエルボーはありません。
f:id:enokisaute:20200906222336p:plain
続いてはクラスタ数2を想定したデータです。左図では一見はっきり2つのクラスタに分かれているように見えますが、SSEが落ち切っておらずK=4のようにわかりやすいエルボーにはなっていません。右図からだけでクラスタ数が2というのは微妙です。
f:id:enokisaute:20200907001403p:plain
これはk-meansのアルゴリズムがすべての方向の距離に対して同じように働くためです。要するに、k-meansでは丸くない、歪な形のクラスタをうまく識別することは難しいのです。
縦横のスケールを同じにしてこのデータを表示すると、下のようになります。
f:id:enokisaute:20200913122015p:plain
実際にK=2でsklearnのk-meansを試してみたところ、データを2つに分けることはできますが右図のような結果となりました。

courseraの講義によると、試してみる価値はあるが、上の例のようにはっきりとしたエルボーがない場合も多いようです。エルボー法の他にもシルエット分析という方法もあるようですが、結論としては、決定的に優れた方法というものはないとのこと。
このKをいくつにするかの選択というのは、どのようなモデルなのかを考えたり、データを可視化するなどしながら最終的には人間の洞察によって手動で決めるしかないようです。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
 
# 最適なクラスタ数をエルボー法で確認する
 
if __name__ == '__main__':
    # K=4, s1=0.1, s2=0.1
    # K=2, s1=0.05, s2=3
    # エルボーなし s1=2, s2=2
    s1 = 0.1
    s2 = 0.1
    # 平均と分散共分散行列
    mean1 = np.array([6, 6])
    cov1 = np.array([[s1, 0], [0, s2]])
 
    mean2 = np.array([6, 8])
    cov2 = np.array([[s1, 0], [0, s2]])
 
    mean3 = np.array([8, 6])
    cov3 = np.array([[s1, 0], [0, s2]])
 
    mean4 = np.array([8, 8])
    cov4 = np.array([[s1, 0], [0, s2]])
 
    data_1 = np.random.multivariate_normal(mean1, cov1, size=100)
    data_2 = np.random.multivariate_normal(mean2, cov2, size=100)
    data_3 = np.random.multivariate_normal(mean3, cov3, size=100)
    data_4 = np.random.multivariate_normal(mean4, cov4, size=100)
 
    # データを連結する
    data = np.vstack((data_1, data_2, data_3, data_4))
    
    plt.figure(1)
    plt.scatter(data[:,0], data[:, 1])
    # 縦横比を同じスケールにする場合
    # plt.axes().set_aspect('equal', 'datalim')
    # plt.xticks(np.arange(2, 14, 2))
 
    sse = []
 
    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(data)
        sse.append(km.inertia_)
 
    plt.figure(2)
    plt.plot(range(1, 11), sse, marker='o')
    plt.xlabel('num of clusters')
    plt.ylabel('SSE')
    plt.show()


参考

k-meansの最適なクラスター数を調べる方法 - Qiita
・Pythonではじめる機械学習 ―scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎

*1:クラスタの番号がtargetと一致するとは限らないので、何度かプログラムを実行する必要があります