fcm聚類算法原理及應(yīng)用

　　模糊C均值（Fuzzy C-means）算法簡稱FCM算法，是一種基于目標(biāo)函數(shù)的模糊聚類算法，主要用于數(shù)據(jù)的聚類分析。理論成熟，應(yīng)用廣泛，是一種優(yōu)秀的聚類算法。本文關(guān)于FCM算法的一些原理推導(dǎo)部分介紹等參考下面視頻，加上自己的理解以文字的形式呈現(xiàn)出來，視頻參考如下，比較長，看不懂的可以再去看看：

　　FCM原理介紹

　　FCM分析1

　　FCM分析2

　　FCM分析3

　　首先介紹一下模糊這個(gè)概念，所謂模糊就是不確定，確定性的東西是什么那就是什么，而不確定性的東西就說很像什么。比如說把20歲作為年輕不年輕的標(biāo)準(zhǔn)，那么一個(gè)人21歲按照確定性的劃分就屬于不年輕，而我們印象中的觀念是21歲也很年輕，這個(gè)時(shí)候可以模糊一下，認(rèn)為21歲有0.9分像年輕，有0.1分像不年輕，這里0.9與0.1不是概率，而是一種相似的程度，把這種一個(gè)樣本屬于結(jié)果的這種相似的程度稱為樣本的隸屬度，一般用u表示，表示一個(gè)樣本相似于不同結(jié)果的一個(gè)程度指標(biāo)。

　　基于此，假定數(shù)據(jù)集為X，如果把這些數(shù)據(jù)劃分成c類的話，那么對(duì)應(yīng)的就有c個(gè)類中心為C，每個(gè)樣本j屬于某一類i的隸屬度為uij，那么定義一個(gè)FCM目標(biāo)函數(shù)（1）及其約束條件（2）如下所示：

　　看一下目標(biāo)函數(shù)（式1）而知，由相應(yīng)樣本的隸屬度與該樣本到各個(gè)類中心的距離相乘組成的，m是一個(gè)隸屬度的因子，個(gè)人理解為屬于樣本的輕緩程度，就像x2與x3這種一樣。式（2）為約束條件，也就是一個(gè)樣本屬于所有類的隸屬度之和要為1。觀察式（1）可以發(fā)現(xiàn)，其中的變量有uij、ci，并且還有約束條件，那么如何求這個(gè)目標(biāo)函數(shù)的極值呢？

　　這里首先采用拉格朗日乘數(shù)法將約束條件拿到目標(biāo)函數(shù)中去，前面加上系數(shù)，并把式（2）的所有j展開，那么式（1）變成下列所示：

　　現(xiàn)在要求該式的目標(biāo)函數(shù)極值，那么分別對(duì)其中的變量uij、ci求導(dǎo)數(shù)，首先對(duì)uij求導(dǎo)。

　　分析式（3），先對(duì)第一部分的兩級(jí)求和的uij求導(dǎo)，對(duì)求和形式下如果直接求導(dǎo)不熟悉，可以把求和展開如下：

　　再來看后面那個(gè)對(duì)uij求導(dǎo)，同樣把求和展開，再去除和uij不相關(guān)的（求導(dǎo)為0），那么只剩下這一項(xiàng)：λj（uij?1），它對(duì)uij求導(dǎo)就是λj了。

　　那么最終J對(duì)uij的求導(dǎo)結(jié)果并讓其等于0就是：

　　我們發(fā)現(xiàn)uij與ci是相互關(guān)聯(lián)的，彼此包含對(duì)方，有一個(gè)問題就是在fcm算法開始的時(shí)候既沒有uij也沒有ci，那要怎么求解呢？很簡單，程序開始的時(shí)候我們隨便賦值給uij或者ci其中的一個(gè)，只要數(shù)值滿足條件即可。然后就開始迭代，比如一般的都賦值給uij，那么有了uij就可以計(jì)算ci，然后有了ci又可以計(jì)算uij，反反復(fù)復(fù)，在這個(gè)過程中還有一個(gè)目標(biāo)函數(shù)J一直在變化，逐漸趨向穩(wěn)定值。那么當(dāng)J不在變化的時(shí)候就認(rèn)為算法收斂到一個(gè)比較好的結(jié)了?？梢钥吹絬ij和ci在目標(biāo)函數(shù)J下似乎構(gòu)成了一個(gè)正反饋一樣，這一點(diǎn)很像EM算法，先E在M，有了M在E，在M直至達(dá)到最優(yōu)。

　　公式（5），（6）是算法的關(guān)鍵?，F(xiàn)在來重新從宏觀的角度來整體看看這兩個(gè)公式，先看（5），在寫一遍

　　假設(shè)看樣本集中的樣本1到各個(gè)類中心的隸屬度，那么此時(shí)j=1，i從1到c類，此時(shí)上述式中分母里面求和中，分子就是這個(gè)點(diǎn)相對(duì)于某一類的類中心距離，而分母是這個(gè)點(diǎn)相對(duì)于所有類的類中心的距離求和，那么它們兩相除表示什么，是不是表示這個(gè)點(diǎn)到某個(gè)類中心在這個(gè)點(diǎn)到所有類中心的距離和的比重。當(dāng)求和里面的分子越小，是不是說就越接近于這個(gè)類，那么整體這個(gè)分?jǐn)?shù)就越大，也就是對(duì)應(yīng)的uij就越大，表示越屬于這個(gè)類，形象的圖如下：

　　再來宏觀看看公式（6），考慮當(dāng)類i確定后，式（6）的分母求和其實(shí)是一個(gè)常數(shù)，那么式（6）可以寫成：

　　這是類中心的更新法則。說這之前，首先讓我們想想kmeans的類中心是怎么更新的，一般最簡單的就是找到屬于某一類的所有樣本點(diǎn)，然后這一類的類中心就是這些樣本點(diǎn)的平均值。那么FCM類中心怎么樣了？看式子可以發(fā)現(xiàn)也是一個(gè)加權(quán)平均，類i確定后，首先將所有點(diǎn)到該類的隸屬度u求和，然后對(duì)每個(gè)點(diǎn)，隸屬度除以這個(gè)和就是所占的比重，乘以xj就是這個(gè)點(diǎn)對(duì)于這個(gè)類i的貢獻(xiàn)值了。畫個(gè)形象的圖如下：

　　由上述的宏觀分析可知，這兩個(gè)公式的迭代關(guān)系式是這樣的也是可以理解的。

　　（二）簡單程序?qū)崿F(xiàn)

　　下面我們?cè)?a href="http://srfitnesspt.com/tags/matlab/" target="_blank">matlab下用最基礎(chǔ)的循環(huán)實(shí)現(xiàn)上述的式（5）與式（6）的FCM過程。首先，我們需要產(chǎn)生可用于FCM的數(shù)據(jù)，為了可視化方便，我們產(chǎn)生一個(gè)二維數(shù)據(jù)便于在坐標(biāo)軸上顯示，也就是每個(gè)樣本由兩個(gè)特征（或者x坐標(biāo)與y坐標(biāo)構(gòu)成），生成100個(gè)這樣的點(diǎn)，當(dāng)然我們?cè)谌藶楦淖円幌拢屵@些點(diǎn)看起來至少屬于不同的類。生成的點(diǎn)畫出來如下：

　　那么我們說FCM算法的一般步驟為：

　?。?）確定分類數(shù)，指數(shù)m的值，確定迭代次數(shù)（這是結(jié)束的條件，當(dāng)然結(jié)束的條件可以有多種）。

　　（2）初始化一個(gè)隸屬度U（注意條件—和為1）；

　?。?）根據(jù)U計(jì)算聚類中心C；

　?。?）這個(gè)時(shí)候可以計(jì)算目標(biāo)函數(shù)J了

　?。?）根據(jù)C返回去計(jì)算U，回到步驟3，一直循環(huán)直到結(jié)束。

　　還需要說一點(diǎn)的是，當(dāng)程序結(jié)束后，怎么去判斷哪個(gè)點(diǎn)屬于哪個(gè)類呢？在結(jié)束后，肯定有最后一次計(jì)算的U吧，對(duì)于每一個(gè)點(diǎn)，它屬于各個(gè)類都會(huì)有一個(gè)u，那么找到其中的最大的u就認(rèn)為這個(gè)點(diǎn)就屬于這一類?；诖艘粋€(gè)基礎(chǔ)的程序如下：

　　clc

　　clear

　　close all

　　%% create samples：

　　for i=1:100

　　x1（i） = rand（）*5; %人為保證差異性

　　y1（i） = rand（）*5;

　　x2（i） = rand（）*5 + 3; %人為保證差異性

　　y2（i） = rand（）*5 + 3;

　　end

　　x = ［x1，x2］;

　　y = ［y1，y2］;

　　data = ［x;y］;

　　data = data‘;%一般數(shù)據(jù)每一行代表一個(gè)樣本

　　%plot（data（：，1），data（：，2），’*‘）; %畫出來

　　%%---

　　cluster_n = 2;%類別數(shù)

　　iter = 50;%迭代次數(shù)

　　m = 2;%指數(shù)

　　num_data = size（data，1）;%樣本個(gè)數(shù)

　　num_d = size（data，2）;%樣本維度

　　%--初始化隸屬度u，條件是每一列和為1

　　U = rand（cluster_n，num_data）;

　　col_sum = sum（U）;

　　U = U./col_sum（ones（cluster_n，1），：）;

　　%% 循環(huán)--規(guī)定迭代次數(shù)作為結(jié)束條件

　　for i = 1:iter

　　%更新c

　　for j = 1:cluster_n

　　u_ij_m = U（j，：）.^m;

　　sum_u_ij = sum（u_ij_m）;

　　sum_1d = u_ij_m./sum_u_ij;

　　c（j，：） = u_ij_m*data./sum_u_ij;

　　end

　　%-計(jì)算目標(biāo)函數(shù)J

　　temp1 = zeros（cluster_n，num_data）;

　　for j = 1:cluster_n

　　for k = 1:num_data

　　temp1（j，k） = U（j，k）^m*（norm（data（k，：）-c（j，：）））^2;

　　end

　　end

　　J（i） = sum（sum（temp1））;

　　%更新U

　　for j = 1:cluster_n

　　for k = 1:num_data

　　sum1 = 0;

　　for j1 = 1:cluster_n

　　temp = （norm（data（k，：）-c（j，：））/norm（data（k，：）-c（j1，：）））.^（2/（m-1））;

　　sum1 = sum1 + temp;

　　end

　　U（j，k） = 1./sum1;

　　end

　　end

　　end

　　figure;

　　subplot（1，2，1），plot（J）;

　　［~，label］ = max（U）; %找到所屬的類

　　subplot（1，2，2）;

　　gscatter（data（：，1），data（：，2），label）123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960

　　得到結(jié)果如下：

　　分成3類看看：

　　基于此，結(jié)果還算正確。但是不得不說的一個(gè)問題就是算法的效率問題。為了和公式計(jì)算方式吻合，便于理解，這個(gè)程序里面有很多的循環(huán)操作，當(dāng)分類數(shù)大一點(diǎn)，樣本多一點(diǎn)的時(shí)候，這么寫是很慢的，matlab號(hào)稱矩陣實(shí)驗(yàn)室，所以要盡量少的使用循環(huán)，直接矩陣操作，那么上述的操作很多地方是可以把循環(huán)改成矩陣計(jì)算的，這里來介紹下matlab自帶的fcm函數(shù)，就是使用矩陣運(yùn)算來的。

　　Matlab下help fcm既可以查閱相關(guān)用法們這里只是簡單介紹，fcm函數(shù)輸入需要2個(gè)或者3個(gè)參數(shù)，返回3個(gè)參數(shù)，如下：

　?。踓enter， U， obj_fcn］ = fcm（data， cluster_n， options）

　　對(duì)于輸入：data數(shù)據(jù)集（注意每一行代表一個(gè)樣本，列是樣本個(gè)數(shù)）

　　cluster_n為聚類數(shù)。

　　options是可選參數(shù)，完整的options包括：

　　OPTIONS（1）： U的指數(shù) （default： 2.0）

　　OPTIONS（2）： 最大迭代次數(shù) （default： 100）

　　OPTIONS（3）： 目標(biāo)函數(shù)的最小誤差 （default： 1e-5）

　　OPTIONS（4）： 是否顯示結(jié)果 （default： 1，顯示）

　　options都有默認(rèn)值，自帶的fcm結(jié)束的條件是OPTIONS（2）與OPTIONS（3）有一個(gè)滿足就結(jié)束。

　　對(duì)于輸出：center聚類中心，U隸屬度，obj_fcn目標(biāo)函數(shù)值，這個(gè)迭代過程的每一代的J都在這里面存著。

　　為了驗(yàn)證我們寫的算法是否正確，用它的fcm去試試我們的數(shù)據(jù)（前提是數(shù)據(jù)一樣），分成3類，畫出它們的obj_fcn看看如下：

　　可以看到，雖然迭代的中間過程不一樣，但是結(jié)果卻是一樣的。

　?。ㄈ┻M(jìn)階應(yīng)用

　　了解了fcm，再來看看它的幾個(gè)應(yīng)用。

　　3.1）基于fcm的圖像分割

　　我們知道fcm主要用于聚類，而圖像分割本身就是一個(gè)聚類的過程。所以可以用fcm去實(shí)現(xiàn)圖像分割。

　　這里以matlab下的灰度圖像為例?；叶葓D像一圖像的角度看是二維的，但是我們知道，決定圖像的無非是里面的灰度值。而灰度值就是一個(gè)值，所以當(dāng)我們把圖像變成1維，也就是拉成一行或者一列的時(shí)候，其實(shí)灰度圖像就是一個(gè)一維數(shù)據(jù)（上面那個(gè)例子生成的隨機(jī)點(diǎn)是二維的）。

　　那么我們就可以對(duì)這個(gè)一維數(shù)據(jù)進(jìn)行聚類，待得到了分類結(jié)果后，再把這個(gè)結(jié)果返回到二維圖像空間去顯示就可以了。

　　一個(gè)例子如下：

　　clc

　　clear

　　close all

　　img = double（imread（’lena.jpg‘））;

　　subplot（1，2，1），imshow（img，［］）;

　　data = img（：）;

　　%分成4類

　　［center，U，obj_fcn］ = fcm（data，4）;

　?。踾，label］ = max（U）; %找到所屬的類

　　%變化到圖像的大小

　　img_new = reshape（label，size（img））;

　　subplot（1，2，2），imshow（img_new，［］）;123456789101112

　　需要注意的是label出來的是標(biāo)簽類別（1-4），并非真實(shí)的灰度，這里不過是把它顯示出來就行了。

　　3.2）實(shí)際數(shù)據(jù)的分類

　　這里介紹一個(gè)常用于機(jī)器學(xué)習(xí)、模式劃分的數(shù)據(jù)下載網(wǎng)站：

　　http://archive.ics.uci.edu/ml/datasets.html

　　這里面包含眾多的數(shù)據(jù)庫可用分類劃分等。這里我們選擇其中一個(gè)數(shù)據(jù)庫：

　　http://archive.ics.uci.edu/ml/datasets/seeds#

　　這個(gè)數(shù)據(jù)庫看介紹好像是關(guān)于種子分類的，里面共包含3類種子，每類種子通過什么x射線技術(shù)等等采集他們的特征，反正最后每個(gè)種子共有7個(gè)特征值來表示它（也就是說在數(shù)據(jù)里面相當(dāng)于7維），每類種子又有70個(gè)樣本，那么整個(gè)數(shù)據(jù)集就是210*7的樣本集。從上面那個(gè)地方下載完樣本集存為txt文件，并放到matlab工作目錄下就可以使用了（注意看看下下來的數(shù)據(jù)有沒有串位的，有的話要手動(dòng)調(diào)整回去）。因?yàn)閙atlab只能顯示低于3維的數(shù)據(jù)，這里有7維，我們現(xiàn)在在二維下顯示其中的兩維以及正確的分類結(jié)果看看什么情況：

　　clc

　　clear

　　close all

　　data = importdata（’data.txt‘）;

　　%data中還有第8列，正確的標(biāo)簽列

　　subplot（2，2，1）;

　　gscatter（data（：，1），data（：，6），data（：，8）），title（’choose:1，6 列‘）

　　subplot（2，2，2）;

　　gscatter（data（：，2），data（：，4），data（：，8）），title（’choose:2，4 列‘）

　　subplot（2，2，3）;

　　gscatter（data（：，3），data（：，5），data（：，8）），title（’choose:3，5 列‘）

　　subplot（2，2，4）;

　　gscatter（data（：，4），data（：，7），data（：，8）），title（’choose:4，7 列‘）12345678910111213

　　組合有限，隨便組合幾種看看，發(fā)現(xiàn)似乎任意兩個(gè)特征都可以把他們分開，當(dāng)然還是有一些分不開的，其中最后一個(gè)選擇特征4，7似乎很好的分開了。

　　Ok，看過之后我們來試試fcm算法對(duì)其進(jìn)行分類，并計(jì)算一下準(zhǔn)確率，我們先把7個(gè)特征都用上看看：

　　clc

　　clear

　　close all

　　data = importdata（’data.txt‘）;

　　%data中還有第8列，正確的標(biāo)簽列

　?。踓enter，U，obj_fcn］ = fcm（data（：，1:7），3）;

　　［~，label］ = max（U）; %找到所屬的類

　　subplot（1，2，1）;

　　gscatter（data（：，4），data（：，7），data（：，8）），title（’choose:4，7列，理論結(jié)果‘）

　　% cal accuracy

　　a_1 = size（find（label（1:70）==1），2）;

　　a_2 = size（find（label（1:70）==2），2）;

　　a_3 = size（find（label（1:70）==3），2）;

　　a = max（［a_1，a_2，a_3］）;

　　b_1 = size（find（label（71:140）==1），2）;

　　b_2 = size（find（label（71:140）==2），2）;

　　b_3 = size（find（label（71:140）==3），2）;

　　b = max（［b_1，b_2，b_3］）;

　　c_1 = size（find（label（141:210）==1），2）;

　　c_2 = size（find（label（141:210）==2），2）;

　　c_3 = size（find（label（141:210）==3），2）;

　　c = max（［c_1，c_2，c_3］）;

　　accuracy = （a+b+c）/210;

　　% plot answer

　　subplot（1，2，2）;

　　gscatter（data（：，4），data（：，7），label），title（［’實(shí)際結(jié)果，accuracy=‘，num2str（accuracy）］）1234567891011121314151617181920212223242526

　　這里選擇以第1與6維的數(shù)據(jù)來可視化這個(gè)結(jié)果?？梢钥吹綔?zhǔn)確率為0.89524。

　　這就是用了所有特征來實(shí)驗(yàn)的，這與用哪個(gè)特征能到達(dá)更好的結(jié)果、怎么樣吧特征進(jìn)行處理下能達(dá)到更好的結(jié)果，這都是機(jī)器學(xué)習(xí)與分類領(lǐng)域在研究的事情。上面我們感覺特征4，7不錯(cuò)，那么當(dāng)我們只用特征4與7去進(jìn)行fcm會(huì)怎樣呢？

　　好像并不是很好，想想只用特征4與7結(jié)果本來就是這樣的，不好就對(duì)了，fcm是根據(jù)數(shù)據(jù)距離劃分來的，所以結(jié)果就是這樣。

　　試了很多組特征，都沒有超過0.89524的，那就所有特征都用上吧。其實(shí)這個(gè)準(zhǔn)確率是可以提高的，我們看到這7個(gè)特征似乎有點(diǎn)重復(fù)有沒有，如果我們把這7個(gè)特征采用pca降維到3，4個(gè)特征了再去fcm實(shí)驗(yàn)?zāi)兀靠梢匀ピ囋?，有待?shí)驗(yàn)……