原理簡介
模糊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\) 為樣本數，\(C\) 為聚類中心數，\(c_j\) 表示第 \(j\) 個聚類中心，和樣本特征維數相同，\(x_i\) 表示第 \(i\) 個樣本，\(u_{ij}\) 表示樣本 \(x_i\) 對聚類中心 \(c_j\) 的隸屬度(通俗的說就是 \(x_i\) 屬於 \(c_j\) 的概率)，顯然滿足
\[\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}\) 和 \(c_j\) 來迭代地優化目標函數Eq. (1)：
\[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\) ，其中 \(t\) 是迭代步數，\(\varepsilon\) 是一個很小的常數表示誤差閾值。也就是說迭代地更新 \(u_{ij}\) 和 \(c_j\) 直到前后兩次隸屬度最大變化值不超過誤差閾值。這個過程最終收斂於 \(J_m\) 的局部極小值點或鞍點。
算法步驟
可以將模糊c均值聚類的過程歸納為以下幾步：
- 初始化隸屬度矩陣 \(U^{(0)}\)，若有 \(N\)個樣本，指定類別數為 \(C\)，則隸屬度矩陣應當是 \(N*C\) 的矩陣；
- 根據式(5)更新聚類中心 \(c_j, j=1,...,C\)；
- 根據式(4)更新 \(U^{(t)}, U^{(t+1)}\)；
- 若滿足終止條件 \(\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實現