k-meansクラスタリング実装メモ

k-meansクラスタリングを実装してみたのでメモ．

gist:

k-means clustering · GitHub

k-meansクラスタリングのアルゴリズム

データのクラスター分けに使われる．（SVMのような分類問題ではない）

データをランダムにクラスター1～kに分ける．
各クラスターの重心（平均位置）を求める．（n次元データであれば重心もn次元になる）
データと各クラスターの重心との距離を計算する．よく使うのは2二乗距離．
データに割り振ったクラスターをもっとも距離の近いクラスターに振り直す．d_c[i] = cluster[j] = argmin(|d_i - cluster[j].centroid|^2)
2～４をクラスターが変わらなくなるまで，または，上限回数まで繰り返す．

初期クラスターに依存する．

実装・コード

仕様：

クラスター数3の場合，クラスターのラベル[0, 1, 2]に割り振る．
クラスターの更新でクラスターが消えてしまった場合は消えたままにする．
クラスターが消えた場合，クラスターラベルを連番になるよう振り直す．ex. k=3で[0, 2]に減った場合，[0, 1]に振り直す．

code:

import scipy as sp
from scipy import stats

def k_means(data, cluster_n = 2, try_num = 1000):
    def calc_cluster_centroid(cluster_list, data_c):
        # update centroid
        '''
        途中で消える場合を想定
        '''
        cluster_centroid = sp.empty((len(cluster_list), data_c.shape[1]-1))
        c_bool = []
        for c in cluster_list:
            one_cluster_bool = (data_c[:,-1]==c)
            is_cluster = sp.any(one_cluster_bool)
            c_bool.append(is_cluster)
            if not is_cluster:
                continue

            idx = sp.where(cluster_list == c)
            cluster_centroid[idx] = sp.mean(data_c[one_cluster_bool, :-1], axis=0)

        return cluster_list[c_bool], cluster_centroid[c_bool]

    def rename_cluster(old_name, data_c, new_name):
        if sp.all(old_name==new_name):
            return data_c

        for idx in sp.arange(len(old_name)):
            data_c[data_c[:, -1]==old_name[idx], -1] = new_name[idx]
        return data_c

    cluster_list = sp.arange(0, cluster_n)
    data = sp.array(data)
    data_n = data.shape[0]

    # add random label
    c_l = None
    for i in sp.arange(try_num):
        random_cluster = sp.random.choice(cluster_list, data_n).reshape(data_n, 1)

        # check
        bool_arr = [sp.any(c==cluster_list) for c in random_cluster]
        if sp.all(bool_arr):
            c_l = random_cluster
            break
    if sp.any(c_l == None):
        print("can't select initial random cluster.")
        return 0

    data_c = sp.hstack((data, c_l))
    old_d = data_c
    for i in sp.arange(try_num):
        new_d = old_d.copy()

        # update centroid
        cluster_list, cluster_centroid = calc_cluster_centroid(cluster_list, new_d)

        # update label
        for d in new_d:
            distance = sp.sum((d[:-1] - cluster_centroid)**2, axis=1)
            min_idx = sp.argmin(distance)
            d[-1] = cluster_list[min_idx]

        # end check
        if sp.all(new_d[:, -1]==old_d[:, -1]):
            return len(cluster_list), sp.arange(len(cluster_list)), rename_cluster(cluster_list, new_d, sp.arange(len(cluster_list)))

        # update for next loop
        old_d = new_d
    return len(cluster_list), sp.arange(len(cluster_list)), rename_cluster(cluster_list, new_d, sp.arange(len(cluster_list)))

実行

d = stats.multivariate_normal(mean=(0,0), cov=[[1,0],[0,1]]).rvs(10)
means = k_means(d)
means

test

2次元標準正規分布をクラスタリング

2次元標準正規分布からサンプリングした値をデータとして， k=2, 3, 4, 5でクラスタリングしてみる．

描画は，gridspecを使って2x2で描く．

%matplotlib inline

import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

gs = gridspec.GridSpec(2, 2)
fig = plt.figure(figsize=(10, 10))

d = stats.multivariate_normal(mean=(0,0), cov=[[1,0],[0,1]]).rvs(10000)

cluster_list = [2, 3, 4, 5]
for idx in sp.arange(len(cluster_list)):
    n, cluster_names, c_d = k_means(d, cluster_n=cluster_list[idx])

    ax = fig.add_subplot(gs[int(idx/2), int(idx%2)])
    for c in cluster_names:
        ax.scatter(c_d[c_d[:,-1]==c][:,0], c_d[c_d[:,-1]==c][:,1], label=f"cluster{c}")

    ax.set_aspect("equal")
    ax.set_xlim(-6, 6)
    ax.set_ylim(-6, 6)
    ax.legend()

plt.tight_layout()
plt.show()