Fuzzy C-means法

概要

  • クラスタ中心との距離の逆数に応じた所属確率を用いる
  • K-means法では最近傍のクラスタのみに所属させていた

イメージ

自宅が吉祥寺駅三鷹駅の間にあるとき、吉祥寺駅との距離をk、三鷹駅との距離をmとすると k < m のとき、吉祥寺駅の近くに住んでいると自慢してよい度合いが増す。逆に k = m に近づけば近づくほど、(住所というものを考慮しなければ) 同時に三鷹の住人でもある度合いが増す (自慢しにくくなる)。

K-means法の場合は、単に最寄り駅がどこかだけを考慮する (どれだけ同じような距離だけ、二つの駅から自宅が離れていたとしても)。

一方、C-meansでは7割吉祥寺の住人で、かつ3割程度三鷹の住人である、というような表現をすることになる。(あくまで駅からの距離で考える場合、) 全員が吉祥寺に住んでいるし三鷹の住人でもあるのだが、その割合が人により異なる

ところで、その駅の住人度というのは駅からの相対的な距離で求めたのであるから、当然その割合が大きい住人の自宅近くが駅である可能性が高い (吉祥寺な度合いが9割なら吉祥寺駅に近いはず)。逆に度合いの低い住人の住所は、(度合いが低い程度に応じて) 考慮に入れなければよい (※実際にはクラスタ中心である駅が、更新毎に位置を変えることになるが..)。

はっきりしないはfuzzyだが (恥をかかない程度に) 役に立つby 市民s

実装してみた (Python、※sse()の実装だけ自信がない)

# fcm.py
import numpy as np
from numpy.linalg import norm


class FCMeans(object):
    
    def __init__(self, n_clusters, n_iter=5, m=2, save_hist=False):
        self.n_clusters = n_clusters    # クラスタ数
        self.n_iter = n_iter    # 中心点の更新回数
        self.m = m    # ファジー係数
        self.fuz = 2.0 / (m - 1.0)    # ファジー性 (曖昧さ)
        self.save_hist = save_hist    # 中心点の移動履歴を保存するかどうか

        self.X = None
        self._cents = None    # 中心点
        self._cents_history = []    # 中心点の移動履歴
        self._sse_history = []    # クラスタ内誤差平方和の履歴
        self._memsp = None    # 各サンプルの、クラスタへの所属確率

    @staticmethod
    def _fspan(xi):
        """特徴値の値の範囲を求める"""
        fmin, fmax = np.min(xi), np.max(xi)
        return fmax - fmin

    def _rand_cent(self):
        """特徴値の範囲に応じて、ランダムな中心点を決定する"""
        ndim = self.X.shape[1]
        cents = np.random.rand(self.n_clusters, ndim)
        for di in range(ndim):
            fspan = FCMeans._fspan(self.X[:, di])
            cents[:, di] *= fspan
        self._cents = cents
        return cents

    @staticmethod
    def _dist_ratio(x, cent_a, cent_b):
        """サンプル点の、それぞれのクラスタ中心との距離の比を求める"""
        da = norm(x - cent_a)
        db = norm(x - cent_b)
        return da / db

    def _update_memsp(self):
        """各サンプルの所属確率を更新する"""
        memsp = []
        for x in self.X:
            probs = []

            # ある中心点と、他の中心点との距離の比を算出
            for target_cent in self._cents:
                dist_ratios = [
                    FCMeans._dist_ratio(x, target_cent, cent)
                        for cent in self._cents
                ]
                dist_ratios_fuz = np.array(dist_ratios) ** self.fuz
                ttl = np.sum(dist_ratios_fuz)
                ttl_inv = 1.0 / ttl
                probs.append(ttl_inv)

            memsp.append(probs)
        memsp = np.array(memsp)
        self._memsp = memsp

    def _update_cents(self):
        """クラスタ中心を更新する"""
        ujs = []
        for j in range(self.n_clusters):
            wj = self._memsp[:, j].reshape(-1, 1)
            wjm = wj ** self.m
            X_wjm = self.X * wjm
            X_wjm_sum = np.sum(X_wjm, axis=0)
            wjm_sum = np.sum(wjm, axis=0)
            uj = X_wjm_sum / wjm_sum
            ujs.append(uj)
        new_cents = np.array(ujs)

        # 中心点の更新履歴を保存
        if self.save_hist:
            self._cents_history.append(new_cents.copy())

        self._cents = new_cents

    def _iter(self):
        """繰り返し部分を実行"""

        # 各クラスタへの所属確率を求める
        self._update_memsp()

        # クラスタ中心を更新
        self._update_cents()

    def fit(self, X):
        """学習させる"""
        self.X = X

        # ランダムな中心点を決定
        self._rand_cent()

        # 繰り返し部分を実行する (所属確率の算出→中心点の更新)
        for i in range(self.n_iter):
            self._iter()
            if self.save_hist:
                new_sse = self.sse()
                self._sse_history.append(new_sse)

    def predict(self, X):
        """サンプルの所属先クラスタを予測する"""

        # 所属先クラスタを更新
        self._update_memsp()

        # 所属確率が高いほうを返す
        return np.argmax(self._memsp, axis=1)

    def sse(self):
        """クラスタ内誤差平方和 (SSE) を算出"""

        # 各クラスタ中心との距離を算出
        d = []
        for x in self.X:
            ds = [norm(x - cent) ** 2 for cent in self._cents]
            d.append(ds)
        d = np.array(d)
        
        return np.sum(self._memsp * d)

使用例

import numpy as np
from matplotlib import pyplot as plt


def make_two_clusters(n_samples, mu1, std1, mu2, std2):
    """
    指定された平均、標準偏差に基づく正規分布から
    クラスターを2つ生成する(次元数=2)
    """
    ndim = 2    # 次元数
    a = np.random.randn(n_samples, ndim) * std1 + mu1
    b = np.random.randn(n_samples, ndim) * std2 + mu2
    X = np.vstack((a, b))
    return X, a, b


def main():
    np.random.seed(123)

    # サンプルデータを生成
    n_samples = 500
    X, a, b = make_two_clusters(n_samples,
        mu1=3.8, std1=0.5,
        mu2=6.0, std2=1.05)

    # 学習させる
    n_clusters = 2
    n_iter = 7
    fcm = FCMeans(n_clusters, n_iter=n_iter, save_hist=True)
    fcm.fit(X)
    y_pred = fcm.predict(X)

    # 分類結果をプロットする
    plot_result = True
    if plot_result:
        colors = ("darkcyan", "orange", "crimson")
        for i, label in enumerate(np.unique(y_pred)):
            mask = y_pred == label
            plt.scatter(X[mask, 0], X[mask, 1],
                color=colors[i], edgecolor="black", alpha=0.7,
                label="cluster %d" % (i + 1,))

        # 中心点の移動履歴をプロット
        cents_hist = np.array(fcm._cents_history)
        for i in range(n_clusters):
            plt.plot(cents_hist[:, i, 0], cents_hist[:, i, 1],
                color="red", label="centroid %d" % (i + 1,))

        plt.legend(loc="best")
        plt.tight_layout()
        plt.show()

    # クラスタ内誤差平方和 (SSE) の推移をプロット
    plot_sse = True
    if plot_sse:
        plt.title("Transition of sum squared error")
        plt.plot(fcm._sse_history, color="orange")
        plt.xlabel("number of iteration")

        ax = plt.gca()
        labels = [str(i) for i in range(n_iter + 1)]
        ax.set_xticklabels(labels)

        plt.grid(True, alpha=0.3)
        plt.show()

以下のようになる。
f:id:zdassen:20170616155216p:plain
走査回数がK-meansより少なくなる性質がある*1
f:id:zdassen:20170616155357p:plain

体感すべき部分

以下のように、サンプルxとクラスタ中心 (現時点での候補) の位置は変化しない状況でファジー係数を大きくしていくと、サンプルxの所属確率はそれぞれ50%に近づく (曖昧な状態)。プロットすると以下のようになる。パッと見た感じではサンプルxは左下のクラスタ中心に所属しているように見える。
f:id:zdassen:20170616174433p:plain
f:id:zdassen:20170616174030p:plain

以下はプロット用のコード。

def dist_ratio(x, cent_a, cent_b):
    """中心点との距離の比を求める"""
    da = norm(x - cent_a)
    db = norm(x - cent_b)
    return da / db


def main():
    """ファジー係数による所属確率の変化の度合いを調べる"""

    # 座標 (0, 0)、(1, 1) に中心点の候補が存在する場合
    cents = np.array([[0., 0.], [1., 1.]])

    # あるサンプル点 ((0, 0) のほうが近そう)
    x = np.array([0.2, 0.2]).reshape(1, -1)
    
    # 現在の状況をプロット
    plot_situation = False
    if plot_situation:
        colors = ("darkcyan", "red")
        labels = ("sample x", "centroids")
        for i, p in enumerate((x, cents)):
            plt.scatter(p[:, 0], p[:, 1],
                color=colors[i], edgecolor="black",
                alpha=0.7, label=labels[i])
        plt.legend(loc="upper left")
        plt.grid(True, alpha=0.3)
        plt.show()

    # ファジー係数を変化させながら
    # サンプルxの所属確率の変化度合いを調べる
    m_max = 10
    memsp = []
    for m in range(2, m_max + 1):
        
        # クラスタへの所属確率
        probs = []

        for target_cent in cents:

            # ファジー性 (曖昧さ)
            fuz = 2.0 / (m - 1.0)

            # (サンプルxの)ある中心点との距離 &
            # (サンプルXの)他の中心点との距離、の比
            dist_ratios = np.array([
                dist_ratio(x, target_cent, cent)
                    for cent in cents
            ])

            # ファジー性 (曖昧さ) を考慮する
            dist_ratios_fuz = dist_ratios ** fuz

            # クラスタへの所属確率
            drf_sum = np.sum(dist_ratios_fuz)
            drf_sum_inv = 1.0 / drf_sum
            probs.append(drf_sum_inv)

        memsp.append(probs)
    memsp = np.array(memsp)

    # 所属確率の変化度合いをプロット
    plot_memsp = True
    if plot_memsp:
        n_clusters = cents.shape[0]
        colors = ("darkcyan", "orange")
        labels = tuple(["ownership of cluster %d" % (i + 1)
            for i in range(n_clusters)])
        for i in range(n_clusters):
            plt.plot(memsp[:, i], color=colors[i], label=labels[i])

        x_labels = [str(i + 1) for i in range(m_max)]
        ax = plt.gca()
        ax.set_xticklabels(x_labels)
        
        plt.xlabel("m")
        plt.ylabel("ownership")
        plt.title("Transition of ownership")
        plt.grid(True, alpha=0.3)
        plt.legend(loc="best")
        plt.show()


if __name__ == '__main__':
    main()

*1:Python機械学習プログラミング 達人データサイエンティストによる理論と実践』 p305