matlab練習程序（馬爾可夫聚類MCL）

本文主要參考：

https://wenku.baidu.com/view/b7907665caaedd3383c4d31b.html

http://www.javashuo.com/article/p-ofbnxekt-ct.html

雖然參考上面兩篇文章，不過這裏和他給出的算法步驟不徹底一致。

由於上面文章是針對Graphs的，矩陣中有邊爲0的數據，而個人數據是Vector的，邊是根據距離計算出來的，應該仍是有點區別的，並且我還用了kmeans。

首先給出他的算法步驟：

　　1.輸入一個無向圖，Expansion的冪e和Inflation的參數r。

　　2.建立鄰接矩陣。

　　3.對每一個結點添加自循環（可選的）。

　　4.標準化矩陣（每一個元素除以所在列的全部元素之和）。

　　5.計算矩陣的第e次冪。

　　6.用參數r對求得的矩陣進行Inflation處理。

　　7.重複第5步和第6步，直到狀態穩定不變（收斂）。

　　8.把最終結果矩陣轉換成聚簇。

而後是這裏的算法步驟：

　　1.輸入數據，建立鄰接矩陣。

　　2.標準化矩陣（每一個元素除以所在列的全部元素之和）。

　　3.對標準化後的矩陣進行馬爾可夫狀態轉移。

　　4.重複第2步和第3步，直到狀態穩定不變（收斂）。

　　5.使用kmeans把最終結果矩陣轉換成聚簇。（不明白kmeans的能夠看看這篇文章）

代碼以下：

clear all;
close all;
clc;

K=3;
theta=0:0.01:2*pi;
p1=[3*cos(theta) + rand(1,length(theta))/2;3*sin(theta)+ rand(1,length(theta))/2];
p2=[2*cos(theta) + rand(1,length(theta))/2;2*sin(theta)+ rand(1,length(theta))/2];
p3=[cos(theta) + rand(1,length(theta))/2;sin(theta)+ rand(1,length(theta))/2];
p=[p1 p2 p3]';

randIndex = randperm(length(p))';       %打亂數據順序
p=p(randIndex,:);
plot(p(:,1),p(:,2),'.')

for i = 1:length(p)
    for j =1:length(p)
        W(i,j) = sqrt(sum((p(i,:)-p(j,:)).^2));　　%根據距離初始化無向圖的邊
    end
end

preW=W;
while 1
    x=repmat(sum(W),length(p),1);
    W=W./x;
    W=W*W;                      %馬爾科夫狀態轉移
    
    if sum(sum(preW-W))<1e-15
        break;
    end
    
    preW=W;
end

[idx,ctrs] = kmeans(W(:,1),K);          %用kmeans將收斂矩陣轉換爲聚簇
figure;
plot(p(idx==1,1),p(idx==1,2),'r.') 
hold on;
plot(p(idx==2,1),p(idx==2,2),'g.')    
plot(p(idx==3,1),p(idx==3,2),'b.') 

原始數據：

聚類後：