k-means+matlab 之辣雞學算法

首先，把 偷來 借鑒來的博文貼上，然後就開始漫漫求生路。

matlab練習程式（k-means聚類）

研究的過程中發現不少人引以為用，足以見其效用啦。

k-means算法分析

• 一、了解算法含義
• 二、引用例子了解
• 三、算法代碼“超”詳解
• main.m
• • 第一類資料
  • 顯示資料
  • k-means聚類
  • 顯示聚類後的資料
• KMeans.m
• • 部分一：
  • zeros（）函數
  • rand（）函數
  • 部分二：
  • 關于花括号{}
  • 部分三：
  • norm()函數
  • min()函數
  • 部分四：
  • 點乘和叉乘
  • sum()函數
  • 部分五：
  • 回頭看看main函數的最後一部分
  • plot3（）函數
• 四、顯示結果

一、了解算法含義

1. 首先準備原始資料{x1，x2，…xn}，是不帶标記的。

   在現實生活中可能是你待聚類的大資料。

2. 初始化k個随機資料{u1，u2，…uk}，即初始配置設定的聚類中心。

在1，2中 x ’ 空格 ’ 和 u ’ 空格 ’ 都是向量。其中，概念“向量”表示1×n 為行向量、n×1 為列向量。如果要表示一個人的身高體重，用向量（high，weight）。是以這裡的向量相當于被描述對象的屬性。

根據兩公式求出最終所有的u，u是最終所有類的聚類中心。

公式一：

求出所有原始資料和初始聚類中心的距離，然後找出離每個初始聚類中心最近的原始資料。

變量：

1. arg j min 函數用法

例如：arg x min f(x)

如果隻有一個值使函數 f(x) 取最小值，則arg為該值，x。其中下标 x 的位置寫後面函數的變量，如果有取值範圍則寫出，如：x屬于[0，4Π]或x屬于R。

如果有多個值使函數取最小值，則arg為集合。

如果無法取最小值，則無定義。

如果min換為max，同理。

在這裡j是變量，代表先确定原始資料，之後不斷地改變聚類中心。也就是拿聚類中心去試原始資料，看原始離哪個聚類中心最近。

1. || x(i) - uj ||下标2

   例如 ||x||下标2， 表示2-範數（歐式距離） = 各分量 x(i) 的平方和再開方。

   是以這裡是x(i) - uj，原始資料的各個分量減去聚類中心的各個分量再開方。

感覺公式裡c、x、|| ||外的上标該改成下标。

公式一總結：求出原始資料與聚類中心的距離，并把距離該原始資料最近的聚類中心的值賦給 c(i) 。

問題：如果算出距離相等怎麼辦？

答：随機劃分。反正要是下次換聚類中心了距離可能就不一樣了。

公式二：

恕我直言，這個公式沒看懂。

解釋長這樣：

求出所有和這個初始聚類中心最近的原始資料的距離的均值，做聚類中心

然後不斷疊代這兩個公式，直到所有的 u 都不怎麼變就完成了。

二、引用例子了解

1. 大哥法（原始資料一維，也就是大家都隻是一個數）

   K-means算法之門派風雲

2. 大哥法（原始資料二維，也就是個坐标）

   K-Means算法之經典大哥法

   害怕連結丢了，我還是偷一部分貼上吧。

炒雞有助于了解，一定要結合前面的公式弄懂了。

K-Means 聚類算法的大緻意思就是“物以類聚，人以群分”：

1、首先輸入 k 的值，即我們指定希望通過聚類得到 k 個分組；

2、從資料集中随機選取 k 個資料點作為初始大佬（質心）；

3、對集合中每一個小弟，計算與每一個大佬的距離，離哪個大佬距離近，就跟定哪個大佬。

4、這時每一個大佬手下都聚集了一票小弟，這時候召開選舉大會，每一群選出新的大佬（即通過算法選出新的質心）。

5、如果新大佬和老大佬之間的距離小于某一個設定的門檻值（表示重新計算的質心的位置變化不大，趨于穩定，或者說收斂），可以認為我們進行的聚類已經達到期望的結果，算法終止。

6、如果新大佬和老大佬距離變化很大，需要疊代3~5步驟。

說了這麼多，估計還是有點糊塗，下面舉個非常形象簡單的例子：

有6個點，從圖上看應該可以分成兩堆，前三個點一堆，後三個點另一堆。現在我手工地把 k-means 計算過程示範一下，同時檢驗是不是和預期一緻：

1.設定 k 值為2

2.選擇初始大佬（就選 P1 和 P2）

3.計算小弟與大佬的距離：

從上圖可以看出，所有的小弟都離 P2 更近，是以次站隊的結果是：

A 組：P1

B 組：P2、P3、P4、P5、P6

4.召開選舉大會：

A 組沒什麼可選的，大佬就是自己

B 組有5個人，需要重新選大佬，這裡要注意選大佬的方法是每個人 X 坐标的平均值和 Y 坐标的平均值組成的新的點，為新大佬，也就是說這個大佬是“虛拟的”。是以，B 組選出新大哥的坐标為：P 哥（（1+3+8+9+10）/5，（2+1+8+10+7）/5）=（6.2，5.6）。

綜合兩組，新大哥為 P1（0，0），P哥（6.2，5.6），而P2-P6重新成為小弟。

5.再次計算小弟到大佬的距離：

這時可以看到P2、P3離P1更近，P4、P5、P6離P哥更近，是以第二次站隊的結果是：

A 組：P1、P2、P3

B 組：P4、P5、P6（虛拟大哥這時候消失）

6.第二屆選舉大會：

同樣的方法選出新的虛拟大佬：P哥1（1.33，1），P哥2（9，8.33），P1-P6都成為小弟。

7.第三次計算小弟到大佬的距離：

這時可以看到 P1、P2、P3 離 P哥1 更近，P4、P5、P6離 P哥2 更近，是以第二次站隊的結果是：

A 組：P1、P2、P3

B 組：P4、P5、P6

我們可以發現，這次站隊的結果和上次沒有任何變化了，說明已經收斂，聚類結束，聚類結果和我們最開始設想的結果完全一緻。

三、算法代碼“超”詳解

main.m

clear all;%清除工作空間的所有變量，函數，和MEX檔案
close all;%關閉所有的Figure視窗 
clc;%清除指令視窗的内容，對工作環境中的全部變量無任何影響 

%第一類資料
mu1=[0 0 0];  %均值，見詳解
S1=[0.3 0 0;0 0.35 0;0 0 0.3];  %協方差，見詳解
data1=mvnrnd(mu1,S1,100);   %産生高斯分布資料，見詳解

%%第二類資料
mu2=[1.25 1.25 1.25];
S2=[0.3 0 0;0 0.35 0;0 0 0.3];
data2=mvnrnd(mu2,S2,100);

%第三類資料
mu3=[-1.25 1.25 -1.25];
S3=[0.3 0 0;0 0.35 0;0 0 0.3];
data3=mvnrnd(mu3,S3,100);


%顯示資料
plot3(data1(:,1),data1(:,2),data1(:,3),'+');
hold on;
plot3(data2(:,1),data2(:,2),data2(:,3),'r+');
plot3(data3(:,1),data3(:,2),data3(:,3),'g+');
grid on;//顯示坐标軸網格線


%三類資料合成一個不帶标号的資料類
data=[data1;data2;data3];   %把三組資料看成整體再劃分

%k-means聚類
[re u]=KMeans(data,3);  %最後産生帶标号的資料，标号在所有資料的最後，資料再加一次元
[m n]=size(re);%表示re矩陣的大小m行n列

%最後顯示聚類後的資料
figure;%打開一個視窗，浮于目前界面之上，直到新的視窗被喚醒
hold on;
for i=1:m 
    if re(i,4)==1   
         plot3(re(i,1),re(i,2),re(i,3),'ro'); 
    elseif re(i,4)==2
         plot3(re(i,1),re(i,2),re(i,3),'go'); 
    else 
         plot3(re(i,1),re(i,2),re(i,3),'bo'); 
    end
end
plot3(u(:,1),u(:,2),u(:,3),'kx','MarkerSize',14,'LineWidth',4);
grid on;
           

第一類資料

mu1=[0 0 0]; %均值

如果内部3（2）個數，生成三個（兩個）數組成的行向量

S1=[0.3 0 0;0 0.35 0;0 0 0.3]; %協方差

生成一個3×3（2×2）的矩陣

data1=mvnrnd(mu1,S1,100); %産生高斯分布資料

生成100行，3（2）列的資料矩陣，其中三列（兩列）分别為100個值的100個坐标（x，y，z）在三維空間（（x，y）在二維空間）

mvnrnd()：生成多元（正态分布 = 高斯分布）資料，n維d列的均值為mu，協方差為SIGMA（S1）的矩陣。

mu1：表示需要生成資料的均值。

S1：需要生成的資料的自相關矩陣（相關系數矩陣）

産生高斯分布資料，需要均值為0，即mu1=0，标準差S1可取任意正數。S1越小，随機變量數值越集中在mu1附近，平均值越接近0。

想産生幾類資料就設幾遍 mu S data，最後顯示數值時若兩類則有兩組，三類則三組資料。

顯示資料

plot3(data1(:,1),data1(:,2),data1(:,3),’+’);

表示取data1的第一列，第二列，第三列，以加号表示畫圖。

plot3（）函數可繪制3D圖形

plot函數用法：plot(X，Y，LineSpec)

X由所有輸入點坐标的x組成

Y由X中包含的x對應的y所組成

LineSpec是使用者指定的繪圖樣式

• +…加号
• o…圓
• *…星号
• ·…點
• g…綠色
• r…紅色
• r+…紅色加号
• g+…綠色加号

grid on：顯示坐标軸網格線

grid off：隐藏坐标軸網格線

hold on：在目前坐标系中畫了一幅圖後再畫一副，原來的與新圖共存。

hold off：在目前坐标系中畫了一幅圖後再畫一副時之前的圖被替代。

k-means聚類

[re u]=KMeans(data,3);

data：MN的矩陣，待聚類的原始資料

3 = K：是data的聚類個數，與u相關

re：N1的向量，存儲的是每組原始資料的聚類标号，類似于label

u：K*N的矩陣，存儲的是K個聚類中心的位置

Kmeans（）函數表示把data矩陣分3類，之後傳回每個點的聚類标号給u，最終聚類中心的位置給re

還有其他形式的Kmeans（）函數

• Index = Kmeans（A，K）
• [Index，C，sumD] = Kmeans（A，K）
• […] = Kmeans（…，‘Param1’，‘Val1’，‘Param2’，‘Val2’，…）

有興趣可自行查閱

顯示聚類後的資料

for i=1:m 
if re(i,4)==1   
%表示矩陣re中的第i行（從第一組資料到最後一組）第4個資料，也就是聚類後新加的一維标簽列，判斷是否為1，也就是是否屬于第一組聚類？
     plot3(re(i,1),re(i,2),re(i,3),'ro'); 
     %  如果是，則根據該組資料的坐标畫出聚類後的圖。
elseif re(i,4)==2
     plot3(re(i,1),re(i,2),re(i,3),'go'); 
else 
     plot3(re(i,1),re(i,2),re(i,3),'bo'); 
end
end
           

從main函數調用kmeans函數開始就涉及到kmeans函數的内容了，是以分析所調用的kmeans函數會有助于了解。

首先說明：

由代碼裡的function文字識别出，Kmeans函數是一個function函數檔案，不是腳本檔案。

把main函數和kmeans函數單獨寫也是為了友善函數之間互相調用。

KMeans.m

%N是資料一共分多少類
%data是輸入的不帶分類标号的資料
%u是每一類的中心
%re是傳回的帶分類标号的資料
%部分一：這一部分屬于以設零矩陣的方式提前限定聚類中心的取值範圍
function [re u]=KMeans(data,N)   
    [m n]=size(data);   %m是資料個數，n是資料維數
    ma=zeros(n);        %每一維最大的數
    mi=zeros(n);        %每一維最小的數
    u=zeros(N,n);       %随機初始化，在用u(i,j)之前都要這麼做。最終疊代到每一類的中心位置
    for i=1:n          
       ma(i)=max(data(:,i));    %每一維最大的數
       mi(i)=min(data(:,i));    %每一維最小的數
       for j=1:N
            u(j,i)=ma(i)+(mi(i)-ma(i))*rand();  %随機初始化，不過最好在每一維[min max]中初始化
       end      
    end
%部分二：求出中心與原始資料的距離   
    while 1
        pre_u=u;            %首先把上一次求得的中心位置放到pre_u中，為到末尾進行pre_u與u的比較做準備
        for i=1:N           %固定聚類中心個數
            tmp{i}=[];      % 公式一中的x(i)-uj,為公式一實作做準備
            for j=1:m       %周遊原始資料行數
                tmp{i}=[tmp{i};data(j,:)-u(i,:)];
            end
        end
%部分三：實作公式一        
        quan=zeros(m,N);
        for i=1:m        %公式一的實作
            c=[];
            for j=1:N
                c=[c norm(tmp{j}(i,:))];
            end
            [junk index]=min(c);
            quan(i,index)=1;              
        end
        
%部分四：實作公式二，設定打破循環條件
        for i=1:N            %公式二的實作
           for j=1:n
                u(i,j)=sum(quan(:,i).*data(:,j))/sum(quan(:,i));
           end           
        end
        
        if norm(pre_u-u)<0.1  %不斷疊代直到位置不再變化
            break;
        end
    end
%部分五：    定義要傳回的re值
    re=[];
    for i=1:m
        tmp=[];
        for j=1:N
            tmp=[tmp norm(data(i,:)-u(j,:))];
        end
        [junk index]=min(tmp);
        re=[re;data(i,:) index];
    end
    
end
           

部分一：

zeros（）函數

ma=zeros(n); %每一維最大的數

其中zeros()函數表示生成一個n×n的零矩陣

同理

u=zeros(N,n);

表示生成N×n的零矩陣

rand（）函數

rand () 函數産生在（0，1）之間均勻分布的随機數組成的數組。

不帶參數的情況rand () 表示産生均值為0，方差為1的随機矩陣，也是标準正态分布的随機矩陣。

為什麼要生成随機矩陣？

• 初始化時如果b = np.zeros((2,1))，b全賦0值，會産生對稱性。為避免該現象，随機生成一些資料指派則不用擔心對稱性。

為什麼用ma(i)+(mi(i)-ma(i))*rand();？

• 是為了保證初始化的中心值在最大值和最小值之間。

為什麼出現u(j,i)這麼别扭的數？

• 因為想先以i循環，後以j循環。同時又必須先以維數為準固定住，再周遊聚類中心行數。是以令i = n，j = N。随之出現先固定維數，行數改變的u（j，i）

部分二：

從 while 1開始進入大循環。由于1永遠為真，是以while永遠循環，直到找到break打破循環

關于花括号{}

tmp{i}=[];       
           

%賦空矩陣給名為tmp的元胞數組。

花括号 {} 表示一個元胞數組，用于元胞型的數組的配置設定或引用。元胞數組的每個單元中可存放不同的資料類型，如數值、字元、矩陣、數組。

小括号（）引用數組的元素

X（3）代表X的第三個元素。

中括号 [] 建構向量或矩陣

X（[1 2 3]）代表X的頭三個元素。

tmp{i}=[tmp{i};data(j,:)-u(i,:)];   % 計算x(i)-uj,為公式一實作做準備
           

關于 tmp{i}=[tmp{i};…];且tmp{i}=[];

解釋為先弄一個空矩陣，之後把後面（…）輸出的資料放到空矩陣中。

如：

Data = [];

Data = [Data;a];

表示把a作為一列加入到Data中，a是列向量就豎着一列一列加，a是行向量就一行一行加。

由此可知，tmp{i}中放的是先固定聚類中心，得到每一行原始資料與所有聚類中心的距離。

部分三：

quan=zeros(m,N);

把m行N列零矩陣賦給quan，初始化。

c=[c norm(tmp{j}(i,:))];

就是公式一，還是往空矩陣中加資料

[]中的空格表示一行一行加資料，

如果是；則是一列一列加資料

如下圖：

i——[1…m];外循環周遊

j——[1…N];内循環周遊

norm()函數

其中norm（）函數表示求範數。

如果norm（A）中的A是矩陣，則求矩陣的最大奇異值。

如果norm（A）中的A是向量，則求向量二範數。

如果norm(A,p)，則求向量A的p範數。

lp範數（p範數）：

二範數（p=2)也叫歐幾裡得範數，也是k-means算法所用的計算距離公式。

norm()函數内部的變量：

tmp{j}(i,:)

如a{i}(j)中i表示通路第i個元胞，j表示通路第i個元胞中的第j個數。

由此可知，要找tmp 第j個元胞的第i個數。由取值範圍可知，j取值1-N，與部分二中的i取值相同，是以部分三中的tmp{j}=部分二中的tmp{i}，而tmp{i}中裝着原始資料與中心的距離。求第二部分tmp第i個元胞的第j個數，表示先拿第一個原始資料和第一個中心的距離（坐标），用歐幾裡得距離（norm()函數）計算成好比較的數，并不斷指派到名為c的空矩陣中。之後周遊聚類中心，表示保持第一個原始資料不動拿出距所有聚類中心的距離。之後依次換原始資料。

min()函數

[junk index]=min（c);

如[m,j] = min(e)表示傳回矩陣e中最小值的位置

m記錄e的每列最小值，j記錄每列最小值行号。

運作時提示junk改為~，因為在matlab運作之前由于沒有使用，是以替換

~：

①表示“與”或“非”中的“非”；

②表示參數預設（比如size()函數會傳回[row,col]，如果我們隻要row，可以寫成[row,~]）

c矩陣如下圖所示，是以找到的每列最小值就是離每個原始資料最近的聚類中心的距離，并得到聚類中心的序号，也就聚類出所有資料屬于哪些聚類中心。

quan(i,index)=1;

表示把quan零矩陣中i行（原始資料）對應的自身聚類中心序号的位置指派為1，這樣選中為1，未選中為0，聚類中心顯而易見。

如圖，假如所有原始資料屬于的聚類中心長這樣

部分四：

u(i,j)=sum(quan(:,i).*data(:,j))/sum(quan(:,i));

公式二是以每個類為機關，自己算自己的新的聚類中心。

方法為：求每個類中所有原始資料與中心的距離的均值。

quan(:,i)中的i範圍為1:N，是部分三中的quan(:,index)，是以要從quan(:,i)開始，在quan的第一列找為1的位置，也就是通過相乘的方式找序号為1的聚類中心包括哪些原始資料。有1可以計算，剩餘位置為0。

點乘和叉乘

點乘：對應位置的數相乘

叉乘：矩陣相乘

如：

a = [1 2 3]

b = [4;5;6]

c = a*b

輸出：c = [1 2 3]·[4;5;6] = 32

d = a.*b

輸出：d =[1 2 3][4;5;6]

輸出錯誤

e = [7 8 9]

d = a.*e

輸出：d = [1 2 3][7 8 9] = [7 16 27]

sum()函數

若X為向量，sum(X)為向量元素之和

若X為矩陣，sum(X)為一個行向量（由每一列元素之群組成）

data(:,j)表示原始資料的第1-n列

由此可知：

sum(quan(:,i).*data(:,j))如圖

是以sum中的值為向量，需要把向量的值相加，就是把所有屬于第一個聚類中心的原始資料都相加。再除以屬于第一個聚類中心的原始資料的個數，平均後的值放入u(1,1)，就完成了公式二

當data(:,j)中j從1-n周遊時表示抽取構成原始資料的所有坐标，依次去乘判斷是否屬于聚類中心的[0-1]矩陣，平均後再指派給u(i,j)

if norm(pre_u-u)<0.1

表示如果前一個儲存好的聚類中心與新産生的聚類中心之間的歐幾裡得距離小于0.1，應該就是聚類中心不變了，是以打破while循環

部分五：

tmp=[tmp norm(data(i,:)-u(j,:))];

把所有原始資料與最終确定的中心的距離算出，并放入tmp中

[junk index]=min(tmp);

把tmp矩陣中每列最小值輸出為junk，每列最小值行号輸出為index

re=[re;data(i,:) index];

index表示加一列，記錄行号

把矩陣data(i,:)和index縱向拼接

至此，kmeans徹底結束。

回頭看看main函數的最後一部分

打開視窗，周遊原始資料從1-m

if re(i,4)==1

plot3(re(i,1),re(i,2),re(i,3),‘ro’);

表示如果re矩陣周遊到的某個原始資料的最後一列，聚類中心标号列為1，則，根據該資料的（x,y,z）輸出該值的位置，并以不同的顔色區分聚類中心。

同理，輸出中心标号為2，3的值。

plot3（）函數

plot3(u(:,1),u(:,2),u(:,3),‘kx’,‘MarkerSize’,14,‘LineWidth’,4);

取出聚類中心矩陣u的第一列（y1,y2,y3)、第二列(k1,k2,k3)、第三列(r1,r2,r3)資料，整理出給了三個點的坐标，（y1,k1,r1),(y2,k2,r2), (y3,k3,r3)。kx中k是黑色，x是x線。‘Markersize’表示标記尺寸為14，也就是黑色x線的大小為14，'LineWidth’表示線寬4号。

如plot3([x,x,x+1000,x+1000,x],[y,y-1000,y-1000,y,y],[z,z,z,z,z]);

給了5個點的坐标，第一個點是（x，y，z），第三個點是（x+1000，y-1000，z）；預設情況下是把點連成線，可以更改參數，如畫點而不連線。

==當遇到不了解的指令時，在指令行輸入 'doc+你不了解的指令’即可 ==

四、顯示結果

東抄抄西學學，希望接下來的研究順順利利。

下一個不出意外的話是DTW算法代碼研究。

然後看能不能銅鼓到一塊去再改進。