原理簡介

  模糊c均值聚類(Fuzzy C-Means)是引入了模糊理論的一種聚類演算法，通過隸屬度來表示樣本屬於某一類的概率，原因在於在很多情況下多個類別之間的界限並不是絕對的明確。顯然，相比於k-means的硬聚類，模糊c均值聚類得到的聚類結果更靈活。

  模糊c均值聚類通過最小化一下目標函式來得到聚類中心：

\[J_{m}=\sum_{i=1}^{N} \sum_{j=1}^{C} u_{i j}^{m}\left\|x_{i}-c_{j}\right\|^{2} \quad, \quad 1 \leq m<\infty \tag{1} \]

其中，\(m>1\) 為模糊係數(fuzzy coefficient)，\(N\)

\[\sum_{j=1}^{C} u_{i j}=1 \tag{2} \]

\(||*||\) 可以是任意度量資料相似性(距離)的範數，最常見的就是歐幾里得範數（又稱歐氏範數，L2範數，歐氏距離）：

\[d=\|x\|_2=\sqrt{\sum_i {x_i^2}} \tag{3} \]

  模糊c均值聚類通過更新 \(u_{ij}\)

\[u_{i j}=\frac{1}{\sum_{k=1}^{C}\left(\frac{\left\|x_{i}-c_{j}\right\|}{\left\|x_{i}-c_{k}\right\|}\right)^{\frac{2}{m-1}}} \tag{4} \]

\[c_{j}=\frac{\sum_{i=1}^{N} u_{i j}^{m} \cdot x_{i}}{\sum_{i=1}^{N} u_{i j}^{m}} \tag{5} \]

迭代的終止條件為 \(\max _{ij}\left\{\left|u_{ij}^{(t+1)}-u_{ij}^{(t)}\right|\right\}<\varepsilon\)

演算法步驟

可以將模糊c均值聚類的過程歸納為以下幾步：

1. 初始化隸屬度矩陣 \(U^{(0)}\)，若有 \(N\)個樣本，指定類別數為 \(C\)，則隸屬度矩陣應當是 \(N*C\) 的矩陣；
2. 根據式(5)更新聚類中心 \(c_j, j=1,...,C\)；
3. 根據式(4)更新 \(U^{(t)}, U^{(t+1)}\)；
4. 若滿足終止條件 \(\max _{ij}\left\{\left|u_{ij}^{(t+1)}-u_{ij}^{(t)}\right|\right\}<\varepsilon\) 則停止迭代，否則返回步驟2。

程式實現

下面程式碼以Iris資料集為例實現了fuzzy c-means。

#!/usr/bin/python3
# -*- coding: utf-8 -*-

'''
@Date    : 2019/9/11
@Author  : Rezero
'''

import numpy as np
import pandas as pd

def loadData(datapath):
    data = pd.read_csv(datapath, sep=',', header=None)
    data = data.sample(frac=1.0)   # 打亂資料順序
    dataX = data.iloc[:, :-1].values # 特徵
    labels = data.iloc[:, -1].values # 標籤
    # 將標籤類別用 0, 1, 2表示
    labels[np.where(labels == "Iris-setosa")] = 0
    labels[np.where(labels == "Iris-versicolor")] = 1
    labels[np.where(labels == "Iris-virginica")] = 2

    return dataX, labels


def initialize_U(samples, classes):
    U = np.random.rand(samples, classes)  # 先生成隨機矩陣
    sumU = 1 / np.sum(U, axis=1)   # 求每行的和
    U = np.multiply(U.T, sumU)   # 使隸屬度矩陣每一行和為1

    return U.T

# 計算樣本和簇中心的距離，這裡使用歐氏距離
def distance(X, centroid):
    return np.sqrt(np.sum((X-centroid)**2, axis=1))


def computeU(X, centroids, m=2):
    sampleNumber = X.shape[0]  # 樣本數
    classes = len(centroids)
    U = np.zeros((sampleNumber, classes))
    # 更新隸屬度矩陣
    for i in range(classes):
        for k in range(classes):
            U[:, i] += (distance(X, centroids[i]) / distance(X, centroids[k])) ** (2 / (m - 1))
    U = 1 / U

    return U


def ajustCentroid(centroids, U, labels):
    newCentroids = [[], [], []]
    curr = np.argmax(U, axis=1)  # 當前中心順序得到的標籤
    for i in range(len(centroids)):
        index = np.where(curr == i)   # 建立中心和類別的對映
        trueLabel = list(labels[index])  # 獲取labels[index]出現次數最多的元素，就是真實類別
        trueLabel = max(set(trueLabel), key=trueLabel.count)
        newCentroids[trueLabel] = centroids[i]
    return newCentroids

def cluster(data, labels, m, classes, EPS):
    """
    :param data: 資料集
    :param m: 模糊係數(fuzziness coefficient)
    :param classes: 類別數
    :return: 聚類中心
    """
    sampleNumber = data.shape[0]  # 樣本數
    cNumber = data.shape[1]       # 特徵數
    U = initialize_U(sampleNumber, classes)   # 初始化隸屬度矩陣
    U_old = np.zeros((sampleNumber, classes))

    while True:
        centroids = []
        # 更新簇中心
        for i in range(classes):
            centroid = np.dot(U[:, i]**m, data) / (np.sum(U[:, i]**m))
            centroids.append(centroid)

        U_old = U.copy()
        U = computeU(data, centroids, m)  # 計算新的隸屬度矩陣

        if np.max(np.abs(U - U_old)) < EPS:
            # 這裡的類別和資料標籤並不是一一對應的, 調整使得第i箇中心表示第i類
            centroids = ajustCentroid(centroids, U, labels)
            return centroids, U


# 預測所屬的類別
def predict(X, centroids):
    labels = np.zeros(X.shape[0])
    U = computeU(X, centroids)  # 計算隸屬度矩陣
    labels = np.argmax(U, axis=1)  # 找到隸屬度矩陣中每行的最大值，即該樣本最大可能所屬類別

    return labels


def main():
    datapath = "iris.data"
    dataX, labels = loadData(datapath)  # 讀取資料

    # 劃分訓練集和測試集
    ratio = 0.6  # 訓練集的比例
    trainLength = int(dataX.shape[0] * ratio)  # 訓練集長度
    trainX = dataX[:trainLength, :]
    trainLabels = labels[:trainLength]
    testX = dataX[trainLength:, :]
    testLabels = labels[trainLength:]

    EPS = 1e-6   # 停止誤差條件
    m = 2        # 模糊因子
    classes = 3  # 類別數
    # 得到各類別的中心
    centroids, U = cluster(trainX, trainLabels, m, classes, EPS)

    trainLabels_prediction = predict(trainX, centroids)
    testLabels_prediction = predict(testX, centroids)


    train_error = 1 - np.sum(np.abs(trainLabels_prediction - trainLabels)) / trainLength
    test_error = 1 - np.sum(np.abs(testLabels_prediction - testLabels)) / (dataX.shape[0] - trainLength)
    print("Clustering on traintset is %.2f%%" % (train_error*100))
    print("Clustering on testset is %.2f%%" % (test_error*100))



if __name__ == "__main__":
    main()

參考資料

A Tutorial on Clustering Algorithms——Fuzzy C-Means Clustering
Fuzzy C-Means（模糊C均值聚類）演算法原理詳解與python實現