Fuzzy C Means 算法及其 Python 实现

1. K-Means 算法向 FCM 算法的扩展

K-Means 算法中,如果要将数据集合 X=\{ X_1,X_2,X_3,\dots,X_n \} 划分为 k\,(1 \le k \le n) 个类,使得任意数据对象 X_i 必须属于并且仅属于一个类,同时每一个类至少包含一个数据对象,那么可以用一个 k\times n 的矩阵 U 来表示,矩阵中的任意一个元素 u_{ij} 可以表示为:

    \[{u_{ij}} = \left{ \begin{cases} 1& {X_i \in G_j} \\ 0& {X_i \notin G_j} \end{cases} \right.\]

其中 {G_j}\left( {j = 1,2, \ldots ,k} \right) 表示第 j 个类。并且 U 需要满足如下条件 (1 \le i \le k,\,1 \le j \le n)

    \[\left{ \begin{cases} {u_{ij}} \in \{0,1\} \\ \sum\limits_{i = 1}^k {{u_{ij}}} = 1 \\ \sum\limits_{j = 1}^n {{u_{ij}}} > 0 \end{cases} \right.\]

如果上述矩阵 U 中的元素 u_{ij} 的取值范围不仅仅是 0 或者 1,那么就可以推广到模糊集合上的划分,U 就变成了模糊判定矩阵。此时 {u_{ij}} 需满足:

(1)   \begin{equation*}  \left{ \begin{cases} {u_{ij}} \in [ 0,1 ]} \\ \sum\limits_{i = 1}^k {{u_{ij}}} = 1 \\ \sum\limits_{j = 1}^n {{u_{ij}}} > 0 \end{cases} \right. \end{equation*}

2. 目标函数与聚类中心

K-Means 算法在度量数据对象的非相似性(或者说距离)时一般使用欧几里得距离,要求每个类的聚类中心与数据对象的距离平方之和最小,目标函数可以表示为:

    \[J = \sum\limits_{i = 1}^k {\sum\limits_{j = 1}^n {s_{ij}^2} }\]

    \[{s_{ij}} = Eculid({C_i},{X_j})\]

其中 C_i 表示任意聚类中心,而聚类中心一般取类内所有对象在各属性上的平均值,因此可以表示为:

    \[{C_i} = \frac{{\sum\limits_{j,{X_j} \in {G_i}} {{X_j}} }}{{\sum\limits_{j = 1}^n {{u_{ij}}} }}\]

{G_i}{\kern 1pt} \left( {1 \le i \le k} \right) 表示任意一个类。

将算法推广到模糊集后,Dunn 对样本与类中心之间的距离采用隶属度的平方来加权,Bezdek 则进一步引入了隶属度的加权指数 m 从而得到了新的目标函数:

(2)   \begin{equation*}  J = \sum\limits_{i = 1}^k {\sum\limits_{j = 1}^n {{{\left( {{u_{ij}}} \right)}^m}s_{ij}^2} } \end{equation*}

要使得 (2) 式达到最小值则要求聚类中心 C_i 和隶属度 u_{ij} 满足如下条件:

(3)   \begin{equation*}  {C_i} = \frac{{\sum\limits_{j = 1}^n {u_{ij}^m{X_j}} }}{{\sum\limits_{j = 1}^n {u_{ij}^m} }} \end{equation*}

(4)   \begin{equation*}  {u_{ij}} = \frac{1}{{\sum\limits_{l = 1}^k {{{\left( {\frac{{{u_{ij}}}}{{{u_{lj}}}}} \right)}^{2/\left( {m - 1} \right)}}} }} \end{equation*}

3. FCM 算法计算过程

由于聚类中心和隶属度两个必要条件的存在,FCM 算法的可以表示为一个简单的迭代过程。将 n 个拥有 m 个属性的对象构建成一个 n\times m 的矩阵记为 DFCM 的算法可以描述为:

输入:隶属度加权指数 m,数据矩阵 D,聚类的个数 k,最大迭代次数 iterMax,阈值 threshold
输出:隶属度矩阵 U
算法过程:
Step1: 用区间 [0,1] 内的满足 (1) 式的随机数初始化隶属度矩阵,记为 D^{(0)}。初始化迭代计数器 iterCounter
Step2: 根据 (3) 式和 U^{(iterCounter)} 计算 k 个聚类中心,记做 C^{(iterCounter)}
Step3: 根据 (2) 式计算目标函数的值。如果本次迭代的目标函数值与上一次计算的目标函数值之差小于 threshold 或者迭代计数器 iterCounter 大于等于 iterMax 那么迭代结束,输出结果,算法停止。
Step4: 依据 C^{(iterCounter)} 和 (4) 式计算 U^{(iterCounter+1)}iterCounter++ ,返回Step2。

4. Python 实现代码

import numpy as np

def loadDataFromTxt(fileName):
    dataSet=np.mat(np.loadtxt(fileName,delimiter='\t'))
    return dataSet

def normalization(dataSet):
    colNum=np.shape(dataSet)[1]
    for index in range(colNum):
        colMax=np.max(dataSet[:,index])
        dataSet[:,index]=dataSet[:,index]/colMax
    return dataSet

def initWithFuzzyMat(n,k):
    fuzzyMat=np.mat(np.zeros((k,n)))
    for colIndex in range(n):
        memDegreeSum=0
        randoms=np.random.rand(k-1,1)
        for rowIndex in range(k-1):
            fuzzyMat[rowIndex,colIndex]=randoms[rowIndex,0]*(1-memDegreeSum)
            memDegreeSum+=fuzzyMat[rowIndex,colIndex]
        fuzzyMat[-1,colIndex]=1-memDegreeSum
    return fuzzyMat

def eculidDistance(vectA,vectB):
    return np.sqrt(np.sum(np.power(vectA-vectB,2)))

def calCentWithFuzzyMat(dataSet,fuzzyMat,p):
    n,m=dataSet.shape
    k=fuzzyMat.shape[0]
    centroids=np.mat(np.zeros((k,m)))
    for rowIndex in range(k):
        degExpArray=np.power(fuzzyMat[rowIndex,:],p)
        denominator=np.sum(degExpArray)
        numerator=np.array(np.zeros((1,m)))
        for colIndex in range(n):
            numerator+=dataSet[colIndex]*degExpArray[0,colIndex]
        centroids[rowIndex,:]=numerator/denominator
    return centroids

def calFuzzyMatWithCent(dataSet,centroids,p):
    n,m=dataSet.shape
    c=centroids.shape[0]
    fuzzyMat=np.mat(np.zeros((c,n)))
    for rowIndex in range(c):
        for colIndex in range(n):
            d_ij=eculidDistance(centroids[rowIndex,:],dataSet[colIndex,:])
            fuzzyMat[rowIndex,colIndex]=1/np.sum([np.power(d_ij/eculidDistance(centroid,dataSet[colIndex,:]),2/(p-1)) for centroid in centroids])
    return fuzzyMat

def calTargetFunc(dataSet,fuzzyMat,centroids,k,p):
    n,m=dataSet.shape
    c=fuzzyMat.shape[0]
    targetFunc=0
    for rowIndex in range(c):
        for colIndex in range(n):
            targetFunc+=eculidDistance(centroids[rowIndex,:],dataSet[colIndex,:])**2*np.power(fuzzyMat[rowIndex,colIndex],p)
    return targetFunc

def fuzzyCMean(dataSet,k,p,initMethod=initWithFuzzyMat):
    n,m=dataSet.shape
    fuzzyMat=initWithFuzzyMat(n,k)
    centroids=calCentWithFuzzyMat(dataSet,fuzzyMat,p)
    lastTargetFunc=calTargetFunc(dataSet,fuzzyMat,centroids,k,p)

    fuzzyMat=calFuzzyMatWithCent(dataSet,centroids,p)
    centroids=calCentWithFuzzyMat(dataSet,fuzzyMat,p)
    targetFunc=calTargetFunc(dataSet,fuzzyMat,centroids,k,p)
    while lastTargetFunc*0.99>targetFunc:
        lastTargetFunc=targetFunc
        fuzzyMat=calFuzzyMatWithCent(dataSet,centroids,p)
        centroids=calCentWithFuzzyMat(dataSet,fuzzyMat,p)
        targetFunc=calTargetFunc(dataSet,fuzzyMat,centroids,k,p)
    return fuzzyMat,centroids


if __name__=='__main__':
    dataSet=loadDataFromTxt('data.txt')
    dataSet=normalization(dataSet)
    fuzzyMat,centroids=fuzzyCMean(dataSet,3,2)
    print 'fuzzyMat=\n',fuzzyMat
    print np.sum(fuzzyMat,axis=0)
    print 'centroids=\n',centroids

使用Matplotlib绘制散点图

由于要做算法结果的可视化输出,需要绘制大量的散点图,所以打算使用 matplotlib,以下给出一个绘制散点图例子。

效果图

scatter

源代码

import matplotlib.pyplot as plt
import numpy as np

n = 100

for color in ['red','blue','green']:
    x,y=np.random.rand(2,n)
    scale=100*np.random.rand(n)
    plt.scatter(x,y,c=color,s=scale,label=color,alpha=0.6,edgecolors='white')

plt.title('Scatter')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.show()

函数说明

numpy.random.rand(d0,d1,d2…)

返回维度为 (d0,d1,d2…) 的位于区间 [0,1) 满足均匀分布的随机数,返回值结果类型为 ndarray。

matplotlib.pyplot.scatter(x,y,c=color,s=scale,label=color,alpha=0.6,edgecolors=’white’)

scatter() 函数用来绘制散点图, x 和 y 为输入数据,形如 shape (n, )。c 表示散点的颜色,指定一个颜色或者色序。s 表示散点的大小,形如 shape (n, )。label 表示显示在图例中的标注。alpha 是 RGBA 颜色的透明分量。edgecolors 指定三点圆周的颜色。

其他函数

函数名 作用
title 图标的标题
xlabel x轴的名称
ylabel y轴的名称
legend 显示右上角的图例
grid 显示网格
show 显示图像

使用Numpy实现K-Means算法

1. K-Means算法简介

1967年,MacQueen 首次提出了 K 均值聚类算法( K-Means算法 )。迄今为止,很多聚类任务都选择该经典算法。K-Means算法的优点是快速高效,其算法复杂度为 O(tKmn),其中 t 表示算法迭代的次数,K 表示聚类的数目,m 表示每个对象拥有的属性个数,n 表示待聚类的对象的个数。

该算法的核心思想是找出 k 个聚类中心,使得每一个数据点 X_i 和与其最近的聚类中心 C_l 的平方距离和被最小化( 该平方距离和被称为代价函数值 E )。

    \[D(X_i,C_l) = \sum\limits_{j = 1}^{ m } {Euclid^2(X_{ij},C_{lj})}\]

    \[E = \sum\limits_{l = 1}^k {\sum\limits_{i = 1}^{{n}} {D({X_i},{C_l})} }\]

K-Means算法的 mean 为类簇的中心,是一个类簇中所有对象在所有属性上的平均。在聚类的初始阶段,我们一般随机指定待聚类对象中的 K 个对象作为 mean

但是该方法也存在着一些缺陷:(1)K-Means 算法往往只能收敛到一个局部最优值。(2)聚类结果对聚类数目 K 和聚类开始时选取的 K 个初始对象非常敏感。

2. K-Means算法描述

输入:n 个带聚类对象,聚类中心个数 K
输出:K 个类簇(包括类簇中心和类簇中的对象)
算法过程:

  1. n 个对象的数据进行标准化。
  2. 随机选取 K 个对象作为初始的聚类中心。
  3. 计算每个对象到 K 个聚类中心的距离,将对象加入距离最近的类簇中。
  4. 根据类簇内的对象更新类簇的中心。
  5. 计算代价函数值 E
  6. 迭代第3步到第5步直至 E收敛。

3. Python代码

下文的代码实现使用了Python的Numpy库,关于Numpy的安装和使用请参照这里

import numpy as np

def loadDataFromTxt(fileName):
    dataMat=[]
    fr=open(fileName)
    for line in fr.readlines():
        curLine=line.strip().split('\t')
        fltLine=map(float,curLine)
        dataMat.append(fltLine)
    dataMat=np.mat(dataMat)
    return dataMat

def autoNormal(dataSet):
    dataShape=np.shape(dataSet)
    n=dataShape[0]
    m=dataShape[1]
    for j in range(m):
        colMax=np.max(dataSet[:,j])
        colMin=np.min(dataSet[:,j])
        colRange=colMax-colMin
        dataSet[:,j]=(dataSet[:,j]-colMin)/colRange
    return dataSet

def randomInitCentroids(dataSet,k):
    m=np.shape(dataSet)[1]
    centroids=np.mat(np.zeros((k,m)))
    for j in range(m):
        colMin=np.min(dataSet[:,j])
        colMax=np.max(dataSet[:,j])
        rangeJ=float(colMax-colMin)
        centroids[:,j]=colMin+np.random.rand(k,1)*rangeJ
    return centroids

def euclidDistance(vectA,vectB):
    return np.sqrt(np.sum(np.power(vectA-vectB,2)))

def KMean(dataSet,k,normalMethod=autoNormal,intialMethod=randomInitCentroids,distMethod=euclidDistance):
    n=dataSet.shape[0]
    dataSet=normalMethod(dataSet)
    centroids=intialMethod(dataSet,k)
    clusterAssignment=np.mat(np.zeros((n,2)))
    clusterChanged=True
    while clusterChanged:
        clusterChanged=False
        for i in range(n):
            minIndex=-1;minDist=np.inf
            for j in range(k):
                distance=distMethod(dataSet[i,:],centroids[j,:])
                if distance<minDist:
                    minIndex=j
                    minDist=distance
            if clusterAssignment[i,0] != minIndex :
                clusterChanged=True
                clusterAssignment[i,:]=minIndex,minDist        
        for cent in range(k):
            ptsInCluster=dataSet[np.where(clusterAssignment[:,0].A==cent)[0]]
            centroids[cent,:]=np.mean(ptsInCluster,axis=0)
    return centroids,clusterAssignment

dataSet=loadDataFromTxt('data.txt')
centroids,clusterAssignment=KMean(dataSet,2)
print clusterAssignment
print centroids

Numpy和matplotlib的安装和使用

最近要写几个数据挖掘的算法,可是由于Java和C#没有比较成熟和知名的可视化类库,所以想到了用Python的Numpy计算,用matplotlib来做结果的可视化。

Numpy与matplotlib

Numpy
NumPy系统是Python的一种开源的数字扩展。这种工具可用来存储和处理大型矩阵,比Python自身的嵌套列表(nested list structure)结构要高效的多(该结构也可以用来表示矩阵matrix)。我用它来进行矩阵的计算,非常便捷。
matplotlib
Matplotlib 可能是 Python 2D- 绘图领域使用最广泛的套件。它能让使用者很轻松地将数据图形化,并且提供多样化的输出格式。

在Windows上安装软件包

在Windows上安装Numpy和matplotlib是一个和痛苦的过程。总体来说,安装这两个包的官方安装程序是远远不够的。总会提示你缺这个缺那个。我的解决办法接单粗暴,提示缺少什么了我就谷歌然后安装这个包。所有的包都可以在UCI的网站上找到。

最后为了安装Nump和matplotlib,我一共安装了6个包:

  • matplotlib-1.4.2.win-amd64-py2.7
  • numpy-MKL-1.9.1.win-amd64-py2.7
  • pyparsing-2.0.3.win-amd64-py2.7
  • python-dateutil-2.2.win-amd64-py2.7
  • scipy-0.15.0b1.win-amd64-py2.7
  • six-1.8.0.win-amd64-py2.7

在Ubuntu Linux上安装

在Ubuntu上安装可就简单多了,软件源里就包括这两个包,一句话搞定:

sudo apt-get install python-numpy python-matplotlib

所有依赖的包都自动安装上了,给个赞!
真心说一句,Ubuntu除了网络问题更新软件不给力之外,装软件真心是很方便、很赞的!

通用方法

其实只要用Python都可以用Python包管理工具pip,两句命令即可搞定:

pip install numpy
pip install matplotlib

而安装pip依赖于setuptools,可以点击setuptoolspip来进行下载安装。

测试代码

安装完成后,可以使用引自此页面的代码进行测试。

import numpy as np
import matplotlib.pyplot as plt

N = 5
menMeans = (20, 35, 30, 35, 27)
menStd =   (2, 3, 4, 1, 2)

ind = np.arange(N)  # the x locations for the groups
width = 0.35       # the width of the bars

fig, ax = plt.subplots()
rects1 = ax.bar(ind, menMeans, width, color='r', yerr=menStd)

womenMeans = (25, 32, 34, 20, 25)
womenStd =   (3, 5, 2, 3, 3)
rects2 = ax.bar(ind+width, womenMeans, width, color='y', yerr=womenStd)

# add some
ax.set_ylabel('Scores')
ax.set_title('Scores by group and gender')
ax.set_xticks(ind+width)
ax.set_xticklabels( ('G1', 'G2', 'G3', 'G4', 'G5') )

ax.legend( (rects1[0], rects2[0]), ('Men', 'Women') )

def autolabel(rects):
    # attach some text labels
    for rect in rects:
        height = rect.get_height()
        ax.text(rect.get_x()+rect.get_width()/2., 1.05*height, '%d'%int(height),
                ha='center', va='bottom')

autolabel(rects1)
autolabel(rects2)

plt.show()

参考文献

查阅到了几篇不错的博客和网站可以用大家参考一下: