K均值聚類算法的MATLAB實現

1.K-均值聚類法的概述

   之前在參加數學建模的過程中用到過這種聚類方法，但是當時只是簡單知道了在matlab中如何調用工具箱進行聚類，并不是特別清楚它的原理。最近因為在學模式識別，又重新接觸了這種聚類算法，所以便仔細地研究了一下它的原理。弄懂了之后就自己手工用matlab編程實現了，最后的結果還不錯，嘿嘿~~~

  簡單來說，K-均值聚類就是在給定了一組樣本(x1, x2, ...xn) （xi, i = 1, 2, ... n均是向量) 之后，假設要將其聚為 m(<n) 類，可以按照如下的步驟實現:

  Step 1: 從 (x1, x2, ...xn) 中隨機選擇 m 個向量(y1,y2,...ym) 作為初始的聚類中心(可以隨意指定,不在n個向量中選擇也可以);

  Step 2: 計算 (x1, x2, ...xn) 到這 m 個聚類中心的距離(嚴格來說為 2階范數);

  Step 3: 對于每一個 xi(i = 1,2,...n)比較其到 (y1,y2,...ym) 距離,找出其中的最小值,若到 yj 的距離最小，則將 xi 歸為第j類;

  Step 4: m 類分好之后, 計算每一類的均值向量作為每一類新的聚類中心;

  Step 5: 比較新的聚類中心與老的聚類中心之間的距離,若大于設定的閾值，則跳到 Step2; 否則輸出分類結果和聚類中心，算法結束。

  OK，廢話不多說，直接上Matlab代碼。

% 利用K-均值聚類的原理，實現對一組數據的分類。這里以一組二維的點為例。
N = 40; % 點的個數
X = 10*rand(1,N);
Y = 10*rand(1,N); % 隨機生成一組橫縱坐標取值均在（0,10）之間的點，X Y 分別代表橫縱坐標
plot(X, Y, 'r*'); % 繪出原始的數據點
xlabel('X');
ylabel('Y');
title('聚類之前的數據點');
n = 2; %將所有的數據點分為兩類
m = 1; %迭代次數
eps = 1e-7; % 迭代結束的閾值
u1 = [X(1),Y(1)]; %初始化第一個聚類中心
u2 = [X(2),Y(2)]; %初始化第二個聚類中心
U1 = zeros(2,100);
U2 = zeros(2,100); %U1,U2 用于存放各次迭代兩個聚類中心的橫縱坐標
U1(:,2) = u1;
U2(:,2) = u2;
D = zeros(2,N); %初始化數據點與聚類中心的距離while(abs(U1(1,m) - U1(1,m+1)) > eps || abs(U1(2,m) - U1(2,m+1) > eps || abs(U2(1,m) - U2(1,m+1)) > eps || abs(U2(2,m) - U2(2,m+1)) > eps))
    m = m +1;    % 計算所有點到兩個聚類中心的距離for i = 1 : N
    D(1,i) = sqrt((X(i) - U1(1,m))^2 + (Y(i) - U1(2,m))^2);
endfor i = 1 : N
    D(2,i) = sqrt((X(i) - U2(1,m))^2 + (Y(i) - U2(2,m))^2);
end
A = zeros(2,N); % A用于存放第一類的數據點
B = zeros(2,N); % B用于存放第二類的數據點for k = 1: N
    [MIN,index] = min(D(:,k)); 
    if index == 1  % 點屬于第一個聚類中心
        A(1,k) = X(k);
        A(2,k) = Y(k);    else           % 點屬于第二個聚類中心
        B(1,k) = X(k);
        B(2,k) = Y(k);
    end
end
indexA = find(A(1,:) ~= 0); % 找出第一類中的點
indexB = find(B(1,:) ~= 0); % 找出第二類中的點
U1(1,m+1) = mean(A(1,indexA));
U1(2,m+1) = mean(A(2,indexA));
U2(1,m+1) = mean(B(1,indexB));
U2(2,m+1) = mean(B(2,indexB)); % 更新兩個聚類中心
end
figure;
plot(A(1,indexA) , A(2,indexA), '*b'); % 作出第一類點的圖形
hold on
plot(B(1,indexB) , B(2,indexB), 'oy'); %作出第二類點的圖形
hold on
centerx = [U1(1,m) U2(1,m)];
centery = [U1(2,m) U2(2,m)];
plot(centerx , centery, '+g'); % 畫出兩個聚類中心點
xlabel('X');
ylabel('Y');
title('聚類之后的數據點');
disp(['迭代的次數為:',num2str(m)]);

得到的分類結果如下:

 50個隨機生成的點分為兩類迭代只需要4步，從上圖來看，分類的效果還是不錯的。但是每次運行可能分類的結果會不一樣，這是因為這些點是隨機生成的,而且也沒有明確的分類標準的緣故。