K-meansクラスタリング - 人工知能に関する断創録

9章のK-meansをPythonで実装してみます。データx_nをあらかじめ指定したK個のクラスタにわけることを考えます。各クラスタの重心をμ_kとします。K個のデータ平均（means）=重心があるからK-meansですね。さらに、2値指示変数r_nkを用意します。これは、データ点x_nがk番目のクラスタに含まれるとき1、それ以外は0になります。各データx_nはただ1つのクラスタに属するという仮定があるためr[n]のK次元ベクトルは1つだけ1であとは0になります。このとき、最小化したい目的関数Jは（9.1）で与えられます。

$(9.1) \hspace{50pt} \displaystyle J = \sum_{n=1}^N \sum_{k=1}^K r_{nk} ||\bf{x_n} - \bf{\mu_k}||^2$

上の式は各データ点x_nからその点が所属するクラスタの重心までの距離を最小化することを意味しています。ここでの目的は、Jを最小にするr_nkとμ_kを求めることです。r_nkとμ_kの2つともわからないので、r_nkとμ_kをそれぞれ最適化する2つのステップを交互に何度も繰り返す手続きで求めようというのがミソのようです。

(1) μ_kの初期値を適当に選ぶ
while (収束するまで）:
    (2) 現在のμ_kの値を用いてJを最小化するr_nkを求める（Eステップに対応）
    (3) 現在のr_nkの値を用いてJを最小化するμ_kを求める（Mステップに対応）

r_nkとμ_kは上のループを回すとJを最小かする最適解に収束していく。r_nkは各データ点がどのクラスタに属するかを表す潜在変数、μ_kが推定すべきパラメータと考えるEMアルゴリズムに対応するのだろうと理解しました。次の混合モデルでも各データ点がどのガウス分布から出たかを表すzを潜在変数としたモデルになっています。ただ、K-meansでは、r_nkが事後分布ではないし、Jも尤度関数ではないので厳密にはEMアルゴリズムとは言わないのかも。

(2) を解くと（9.2）が得られます。

$(9.2) \hspace{50pt} \displaystyle r_{nk} = 1 \hspace{50pt} \arg \min_{j} ||\bf{x_n} - \bf{mu_j}||^2$
$(9.2) \hspace{50pt} \displaystyle r_{nk} = 0 \hspace{50pt} otherwise$

また、(3)を解くと（9.4）が得られます。

$(9.4) \hspace{50pt} \displaystyle \mu_k = \frac{\sum_{n} r_{nk} \bf{x_n}}{\sum_n r_{nk}}$

では、Pythonでプログラミングしてみます。PRMLの原著サポートページのfaithful.txtというデータを使うので同じフォルダにおいてください。

#coding:utf-8

# K-meansアルゴリズム

import numpy as np
from pylab import *

K = 2  # クラスターの数（固定）

def scale(X):
    """データ行列Xを属性ごとに標準化したデータを返す"""
    # 属性の数（=列の数）
    col = X.shape[1]
    
    # 属性ごとに平均値と標準偏差を計算
    mu = np.mean(X, axis=0)
    sigma = np.std(X, axis=0)
    
    # 属性ごとデータを標準化
    for i in range(col):
        X[:,i] = (X[:,i] - mu[i]) / sigma[i]
    
    return X

def J(X, mean, r):
    """K-meansの目的関数（最小化を目指す）"""
    sum = 0.0
    for n in range(len(X)):
        temp = 0.0
        for k in range(K):
            temp += r[n, k] * np.linalg.norm(X[n] - mean[k]) ** 2
        sum += temp
    return sum

if __name__ == "__main__":
    # 訓練データをロード
    data = np.genfromtxt("faithful.txt")
    X = data[:, 0:2]
    X = scale(X)
    N = len(X)
    
    # 平均を初期化
    mean = np.random.rand(K, 2)
    
    # クラスタ割当変数の配列を用意
    r = np.zeros( (N, K) )
    r[:, 0] = 1
    
    # 目的関数の初期値を計算
    target = J(X, mean, r)
    
    turn = 0
    while True:
        print turn, target
        
        # E-step : 現在のパラメータmeanを使って、クラスタ割当変数rを計算 (9.2)
        for n in range(N):
            idx = -1
            min = sys.maxint
            for k in range(K):
                temp = np.linalg.norm(X[n] - mean[k]) ** 2
                if temp < min:
                    idx = k
                    min = temp
            for k in range(K):
                r[n, k] = 1 if k == idx else 0
        
        # M-step : 現在のクラスタ割当変数rを用いてパラメータmeanを更新 (9.4)
        for k in range(K):
            numerator = 0.0
            denominator = 0.0
            for n in range(N):
                numerator += r[n, k] * X[n]
                denominator += r[n, k]
            mean[k] = numerator / denominator
        
        # 収束判定
        new_target = J(X, mean, r)
        diff = target - new_target
        if diff < 0.01:
            break
        target = new_target
        turn += 1
    
    # 訓練データを描画
    color = ['r', 'b', 'g', 'c', 'm', 'y', 'k', 'w']  # K < 8と想定
    for n in range(N):
        for k in range(K):
            if r[n, k] == 1:
                c = color[k]
        scatter(X[n,0], X[n,1], c=c, marker='o')
    
    # クラスタの平均を描画
    for k in range(K):
        scatter(mean[k, 0], mean[k, 1], s=120, c='y', marker='s')
    
    xlim(-2.5, 2.5)
    ylim(-2.5, 2.5)
    show()

収束判定は、目的関数Jの変化が0.01より小さくなったらにしています。目的関数Jの値を出力してみると

0 742.330759235
1 269.478511661
2 55.271004262
3 50.7190592081
4 48.7905321676
5 47.8034529317
6 47.1466162636
7 46.7667915463
8 46.5271172318
9 46.4756638583

のようにループが進むごとに小さくなっていき収束するのがわかります。クラスタ数K=2、K=3、K=4で実行すると下のようになります。ただ、クラスタ重心のμの初期値によって違う結果になります。

f:id:aidiary:20100515104056p:plain
f:id:aidiary:20100515104057p:plain
f:id:aidiary:20100515104058p:plain