Matlab實現K-Means聚類算法

      人生如戲！！！！

1、理論準備

      聚類算法，不是分類算法。分類算法是給一個數據，而後判斷這個數據屬於已分好的類中的具體哪一類。聚類算法是給一大堆原始數據，而後經過算法將其中具備類似特徵的數據聚爲一類。

      K-Means算法的基本思想是初始隨機給定K個簇中心，按照最鄰近原則把待分類樣本點分到各個簇。而後按平均法從新計算各個簇的質心，從而肯定新的簇心。一直迭代，直到簇心的移動距離小於某個給定的值。

      算法大體思路：
      一、從給定樣本中任選幾個點做爲初始中心（我取k=2）
      二、計算其他點分別和初始中心點的距離，跟哪一個初始中心近就跟那個中心點歸爲一類(歐式距離公式)，直到各自爲「派別」
      三、在分好類的基礎上按平均值的方法從新計算聚類中心點，再重複第二步...以此類推
      四、直到最後算法收斂（能夠理解爲中心點再也不變更）則結束。

      小知識點：

    （1）s=size(A),
         當只有一個輸出參數時，返回一個行向量，該行向量的第一個元素時矩陣的行數，第二個元素是矩陣的列數。
   （2）[r,c]=size(A),
         當有兩個輸出參數時，size函數將矩陣的行數返回到第一個輸出變量r，將矩陣的列數返回到第二個輸出變量c。
   （3）size(A,n)若是在size函數的輸入參數中再添加一項n，並用1或2爲n賦值，則 size將返回矩陣的行數或列數。其中r=size(A,1)該語句返回的時矩陣A的行數， c=size(A,2) 該語句返回的時矩陣A的列數。
另外，length()=max(size())。

 

2、算法實現

  1: %K-means算法主程序 

  2: k=4; 

  3: x =[ 1.2126   2.1338   0.5115   0.2044 

  4:   -0.9316   0.7634   0.0125   -0.2752 

  5:   -2.9593   0.1813   -0.8833   0.8505 

  6:   3.1104   -2.5393   -0.0588   0.1808 

  7:   -3.1141   -0.1244   -0.6811   0.9891 

  8:   -3.2008   0.0024   -1.2901   0.9748 

  9:   -1.0777   1.1438   0.1996   0.0139 

 10:   -2.7213   -0.1909   0.1184   0.1013 

 11:   -1.1467   1.3820   0.1427   -0.2239 

 12:   1.1497   1.9414   -0.3035   0.3464 

 13:   2.6993   -2.2556   0.1637   -0.0139 

 14:   -3.0311   0.1417   0.0888   0.1791 

 15:   -2.8403   -0.1809   -0.0965   0.0817 

 16:   1.0118   2.0372   0.1638   -0.0349 

 17:   -0.8968   1.0260   -0.1013   0.2369 

 18:   1.1112   1.8802   -0.0291   -0.1506 

 19:   1.1907   2.2041   -0.1060   0.2167 

 20:   -1.0114   0.8029   -0.1317   0.0153 

 21:   -3.1715   0.1041   -0.3338   0.0321 

 22:   0.9718   1.9634   0.0305   -0.3259 

 23:   -1.0377   0.8889   -0.2834   0.2301 

 24:   -0.8989   1.0185   -0.0289   0.0213 

 25:   -2.9815   -0.4798   0.2245   0.3085 

 26:   -0.8576   0.9231   -0.2752   -0.0091 

 27:   -3.1356   0.0026   -1.2138   0.7733 

 28:   3.4470   -2.2418   0.2014   -0.1556 

 29:   2.9143   -1.7951   0.1992   -0.2146 

 30:   3.4961   -2.4969   -0.0121   0.1315 

 31:   -2.9341   -0.1071   -0.7712   0.8911 

 32:   -2.8105   -0.0884   -0.0287   -0.1279 

 33:   3.1006   -2.0677   -0.2002   -0.1303 

 34:   0.8209   2.1724   0.1548   0.3516 

 35:   -2.8500   0.3196   0.1359   -0.1179 

 36:   -2.8679   0.1365   -0.5702   0.7626 

 37:   -2.8245   -0.1312   0.0881   -0.1305 

 38:   -0.8322   1.3014   -0.3837   0.2400 

 39:   -2.6063   0.1431   0.1880   0.0487 

 40:   -3.1341   -0.0854   -0.0359   -0.2080 

 41:   0.6893   2.0854   -0.3250   -0.1007 

 42:   1.0894   1.7271   -0.0176   0.6553 

 43:   -2.9851   -0.0113   0.0666   -0.0802 

 44:   1.0371   2.2724   0.1044   0.3982 

 45:   -2.8032   -0.2737   -0.7391   1.0277 

 46:   -2.6856   0.0619   -1.1066   1.0485 

 47:   -2.9445   -0.1602   -0.0019   0.0093 

 48:   1.2004   2.1302   -0.1650   0.3413 

 49:   3.2505   -1.9279   0.4462   -0.2405 

 50:   -1.2080   0.8222   0.1671   0.1576 

 51:   -2.8274   0.1515   -0.9636   1.0675 

 52:   2.8190   -1.8626   0.2702   0.0026 

 53:   1.0507   1.7776   -0.1421   0.0999 

 54:   -2.8946   0.1446   -0.1645   0.3071 

 55:   -1.0105   1.0973   0.0241   0.1628 

 56:   -2.9138   -0.3404   0.0627   0.1286 

 57:   -3.0646   -0.0008   0.3819   -0.1541 

 58:   1.2531   1.9830   -0.0774   0.2413 

 59:   1.1486   2.0440   -0.0582   -0.0650 

 60:   -3.1401   -0.1447   -0.6580   0.9562 

 61:   -2.9591   0.1598   -0.6581   1.1937 

 62:   -2.9219   -0.3637   -0.1538   -0.2085 

 63:   2.8948   -2.2745   0.2332   -0.0312 

 64:   -3.2972   -0.0219   -0.0288   -0.1436 

 65:   -1.2737   0.7648   0.0643   0.0858 

 66:   -1.0690   0.8108   -0.2723   0.3231 

 67:   -0.5908   0.7508   -0.5456   0.0190 

 68:   0.5808   2.0573   -0.1658   0.1709 

 69:   2.8227   -2.2461   0.2255   -0.3684 

 70:   0.6174   1.7654   -0.3999   0.4125 

 71:   3.2587   -1.9310   0.2021   0.0800 

 72:   1.0999   1.8852   -0.0475   -0.0585 

 73:   -2.7395   0.2585   -0.8441   0.9987 

 74:   -1.2223   1.0542   -0.2480   -0.2795 

 75:   -2.9212   -0.0605   -0.0259   0.2591 

 76:   3.1598   -2.2631   0.1746   0.1485 

 77:   0.8476   1.8760   -0.2894   -0.0354 

 78:   2.9205   -2.2418   0.4137   -0.2499 

 79:   2.7656   -2.1768   0.0719   -0.1848 

 80:   -0.8698   1.0249   -0.2084   -0.0008 

 81:   -1.1444   0.7787   -0.4958   0.3676 

 82:   -1.0711   1.0450   -0.0477   -0.4030 

 83:   0.5350   1.8110   -0.0377   0.1622 

 84:   0.9076   1.8845   -0.1121   0.5700 

 85:   -2.7887   -0.2119   0.0566   0.0120 

 86:   -1.2567   0.9274   0.1104   0.1581 

 87:   -2.9946   -0.2086   -0.8169   0.6662 

 88:   1.0536   1.9818   -0.0631   0.2581 

 89:   -2.8465   -0.2222   0.2745   0.1997 

 90:   -2.8516   0.1649   -0.7566   0.8616 

 91:   -3.2470   0.0770   0.1173   -0.1092 

 92:   -2.9322   -0.0631   -0.0062   -0.0511 

 93:   -2.7919   0.0438   -0.1935   -0.5023 

 94:   0.9894   1.9475   -0.0146   -0.0390 

 95:   -2.9659   -0.1300   0.1144   0.3410 

 96:   -2.7322   -0.0427   -1.0758   0.9718 

 97:   -1.4852   0.8592   -0.0503   -0.1373 

 98:   2.8845   -2.1465   -0.0533   -0.1044 

 99:   -3.1470   0.0536   0.1073   0.3323 

100:   2.9423   -2.1572   0.0505   0.1180 

101:   -3.0683   0.3434   -0.6563   0.8960 

102:   1.3215   2.0951   -0.1557   0.3994 

103:   -0.7681   1.2075   -0.2781   0.2372 

104:   -0.6964   1.2360   -0.3342   0.1662 

105:   -0.6382   0.8204   -0.2587   0.3344 

106:   -3.0233   -0.1496   -0.2607   -0.0400 

107:   -0.8952   0.9872   0.0019   0.3138 

108:   -0.8172   0.6814   -0.0691   0.1009 

109:   -3.3032   0.0571   -0.0243   -0.1405 

110:   0.7810   1.9013   -0.3996   0.7374 

111:   -0.9030   0.8646   -0.1498   0.1112 

112:   -0.8461   0.9261   -0.1295   -0.0727 

113:   2.8182   -2.0818   -0.1430   -0.0547 

114:   2.9295   -2.3846   -0.0244   -0.1400 

115:   1.0587   2.2227   -0.1250   0.0957 

116:   3.0755   -1.7365   -0.0511   0.1500 

117:   -1.3076   0.8791   -0.3720   0.0331 

118:   -2.8252   -0.0366   -0.6790   0.7374 

119:   -2.6551   -0.1875   0.3222   0.0483 

120:   -2.9659   -0.1585   0.4013   -0.1402 

121:   -3.2859   -0.1546   0.0104   -0.1781 

122:   -0.6679   1.1999   0.1396   -0.3195 

123:   -1.0205   1.2226   0.1850   0.0050 

124:   -3.0091   -0.0186   -0.9111   0.9663 

125:   -3.0339   0.1377   -0.9662   1.0664 

126:   0.8952   1.9594   -0.3221   0.3579 

127:   -2.8481   0.1963   -0.1428   0.0382 

128:   1.0796   2.1353   -0.0792   0.6491 

129:   -0.8732   0.8985   -0.0049   0.0068 

130:   1.0620   2.1478   -0.1275   0.3553 

131:   3.4509   -1.9975   0.1285   -0.1575 

132:   -3.2280   -0.0640   -1.1513   0.8235 

133:   -0.6654   0.9402   0.0577   -0.0175 

134:   -3.2100   0.2762   -0.1053   0.0626 

135:   3.0793   -2.0043   0.2948   0.0411 

136:   1.3596   1.9481   -0.0167   0.3958 

137:   -3.1267   0.1801   0.2228   0.1179 

138:   -0.7979   0.9892   -0.2673   0.4734 

139:   2.5580   -1.7623   -0.1049   -0.0521 

140:   -0.9172   1.0621   -0.0826   0.1501 

141:   -0.7817   1.1658   0.1922   0.0803 

142:   3.1747   -2.1442   0.1472   -0.3411 

143:   2.8476   -1.8056   -0.0680   0.1536 

144:   -0.6175   1.4349   -0.1970   -0.1085 

145:   0.7308   1.9656   0.2602   0.2801 

146:   -1.0310   1.0553   -0.2928   -0.1647 

147:   -2.9251   -0.2095   0.0582   -0.1813 

148:   -0.9827   1.2720   -0.2225   0.2563 

149:   -1.0830   1.1158   -0.0405   -0.1181 

150:   -2.8744   0.0195   -0.3811   0.1455 

151:   3.1663   -1.9241   0.0455   0.1684 

152:   -1.0734   0.7681   -0.4725   -0.1976]; 

153: [n,d] = size(x); 

154: bn=round(n/k*rand);%第一個隨機數在前1/K的範圍內 

155: %；表示按列顯示，都好表示按行顯示

156: nc=[x(bn,:);x(2*bn,:);x(3*bn,:);x(4*bn,:)];%初始聚類中心 

157: %x(bn,:) 選擇某一行數據做爲聚類中心，其列值爲所有

158: 

159: %x數據源，k聚類數目，nc表示k個初始化聚類中心

160: %cid表示每一個數據屬於哪一類，nr表示每一類的個數，centers表示聚類中心

161: [cid,nr,centers] = kmeans(x,k,nc)%調用kmeans函數 

162: %認爲不應是150，或者說不應是個肯定值，該是size(x,1)就是x行數

163: for i=1:150

164:   if cid(i)==1, 

165:      plot(x(i,1),x(i,2),'r*') % 顯示第一類

166:     %plot(x(i,2),'r*') % 顯示第一類

167:     hold on 

168:   else 

169:       if cid(i)==2, 

170:         plot(x(i,1),x(i,2),'b*') %顯示第二類 

171:         %  plot(x(i,2),'b*') % 顯示第一類

172:           hold on 

173:       else 

174:               if cid(i)==3, 

175:                     plot(x(i,1),x(i,2),'g*') %顯示第三類 

176:                    % plot(x(i,2),'g*') % 顯示第一類

177:                     hold on 

178:                   else 

179:                       if cid(i)==4, 

180:                        plot(x(i,1),x(i,2),'k*') %顯示第四類

181:                         % plot(x(i,2),'k*') % 顯示第一類

182:                         hold on 

183:                       end 

184:               end 

185:         end 

186:   end 

187: end 

188: strt=['紅色*爲第一類；藍色*爲第二類；綠色*爲第三類；黑色*爲第四類' ]; 

189: text(-4,-3.6,strt); 

 

  1: %BasicKMeans.m主類 

  2: %x數據源，k聚類數目，nc表示k個初始化聚類中心

  3: %cid表示每一個數據屬於哪一類，nr表示每一類的個數，centers表示聚類中心

  4: function [cid,nr,centers] = kmeans(x,k,nc) 

  5:     [n,d] = size(x); 

  6:     % 設置cid爲分類結果顯示矩陣 

  7:     cid = zeros(1,n); 

  8:     % Make this different to get the loop started. 

  9:     oldcid = ones(1,n); 

 10:     % The number in each cluster. 

 11:     nr = zeros(1,k); 

 12:     % Set up maximum number of iterations. 

 13:     maxgn= 100; 

 14:     iter = 1; 

 15:     %計算每一個數據到聚類中心的距離 ,選擇最小的值得位置到cid

 16:     %我記得是聚類中心近乎再也不變化迭代中止

 17:     while iter < maxgn  

 18:         for i = 1:n 

 19:             %repmat 即 Replicate Matrix ，複製和平鋪矩陣，是 MATLAB 裏面的一個函數。

 20:             %B = repmat(A,m,n)將矩陣 A 複製 m×n 塊，即把 A 做爲 B 的元素，B 由 m×n 個 A 平鋪而成。B 的維數是 [size(A,1)*m, size(A,2)*n] 。

 21:             %點乘方a.^b,矩陣a中每一個元素按b中對應元素乘方或者b是常數

 22:             %sum(x,2)表示矩陣x的橫向相加，求每行的和，結果是列向量。 而缺省的sum(x)就是豎向相加，求每列的和，結果是行向量。

 23:             dist = sum((repmat(x(i,:),k,1)-nc).^2,2); 

 24:             [m,ind] = min(dist); % 將當前聚類結果存入cid中 

 25:             cid(i) = ind; 

 26:         end 

 27:          %找到每一類的全部數據，計算他們的平均值，做爲下次計算的聚類中心 

 28:         for i = 1:k 

 29:             %find（A>m,4）返回矩陣A中前四個數值大於m的元素所在位置

 30:             ind = find(cid==i); 

 31:             %mean(a,1)=mean(a)縱向;mean(a,2)橫向

 32:             nc(i,:) = mean(x(ind,:)); 

 33:             %統計每一類的數據個數 

 34:             nr(i) = length(ind); 

 35:         end 

 36:         iter = iter + 1; 

 37:     end 

 38: 

 39:     % Now check each observation to see if the error can be minimized some more. 

 40:     % Loop through all points. 

 41:     

 42:     maxiter = 2;

 43:     iter = 1; 

 44:     move = 1; 

 45:     %j~=k這是一個邏輯表達式，j不等於k，若是j不等於k，返回值爲1，不然爲0

 46:     while iter < maxiter & move ~= 0 

 47:         move = 0; 

 48:         %對全部的數據進行再次判斷，尋求最佳聚類結果 

 49:         for i = 1:n 

 50:             dist = sum((repmat(x(i,:),k,1)-nc).^2,2); 

 51:             r = cid(i); % 將當前數據屬於的類給r 

 52:             %點除，a.\b表示矩陣b的每一個元素除以a中對應元素或者除以常數a，a./b表示常數a除以矩陣b中每一個元素或者矩陣a除以矩陣b對應元素或者常數b

 53:             %nr是沒一類的的個數

 54:             

 55:             

 56:             

 57:             %這個調整看不懂

 58:             %點乘(對應元素相乘),必須同維或者其中一個是標量，a.*b

 59:             dadj = nr./(nr+1).*dist'; % 計算調整後的距離 

 60:             

 61:             [m,ind] = min(dadj); % 找到該數據距哪一個聚類中心最近 

 62:             if ind ~= r % 若是不等則聚類中心移動 

 63:                   cid(i) = ind;%將新的聚類結果送給cid 

 64:                   ic = find(cid == ind);%從新計算調整當前類別的聚類中心 

 65:                   nc(ind,:) = mean(x(ic,:)); 

 66:                   move = 1; 

 67:             end 

 68:         end 

 69:         iter = iter+1; 

 70:     end 

 71:     centers = nc; 

 72:     if move == 0 

 73:         disp('No points were moved after the initial clustering procedure.') 

 74:     else 

 75:         disp('Some points were moved after the initial clustering procedure.') 

 76:     end

 

3、算法結果

      控制檯自動輸出的結果以下，我很奇怪怎麼本身輸出了。

  1: >> main

  2: No points were moved after the initial clustering procedure.

  3: cid =

  4:   Columns 1 through 22

  5:      2     3     1     4     1     1     3     1     3     2     4     1     1     2     3     2     2     3     1     2     3     3

  6:   Columns 23 through 44

  7:      1     3     1     4     4     4     1     1     4     2     1     1     1     3     1     1     2     2     1     2     1     1

  8:   Columns 45 through 66

  9:      1     2     4     3     1     4     2     1     3     1     1     2     2     1     1     1     4     1     3     3     3     2

 10:   Columns 67 through 88

 11:      4     2     4     2     1     3     1     4     2     4     4     3     3     3     2     2     1     3     1     2     1     1

 12:   Columns 89 through 110

 13:      1     1     1     2     1     1     3     4     1     4     1     2     3     3     3     1     3     3     1     2     3     3

 14:   Columns 111 through 132

 15:      4     4     2     4     3     1     1     1     1     3     3     1     1     2     1     2     3     2     4     1     3     1

 16:   Columns 133 through 150

 17:      4     2     1     3     4     3     3     4     4     3     2     3     1     3     3     1     4     3

 18: nr =

 19:     55    30    40    25

 20: centers =

 21:   -2.962918181818183  -0.023009090909091  -0.297021818181818   0.341136363636364

 22:    0.995233333333333   1.997873333333334  -0.078486666666667   0.229650000000000

 23:   -0.956882500000000   0.997800000000000  -0.123667500000000   0.049320000000000

 24:    3.023444000000000  -2.098592000000001   0.102096000000000  -0.050580000000000

      每次運行結果獲得的圖片都不同，最奇怪的是第二個圖片居然重疊類別。

                

4、結果分析

      不適合處理離散型屬性，可是對於連續型具備較好的聚類效果。

      對於不一樣的初始值，可能會致使不一樣結果：多設置一些不一樣的初值，但比較耗時和浪費資源。

      分類數目K不肯定：經過類的自動合併和分裂，獲得較爲合理的類型數目K，例如ISODATA算法。相同點：聚類中心都是經過樣本均值的迭代運算來決定的；不一樣點：主要是在選代過程當中可將一類一分爲二，亦可能二類合二爲一，即「自組織」，這種算法具備啓發式的特色。因爲算法有自我調整的能力，於是須要設置若干個控制用參數，如聚類數指望值K、最小類內樣本數、類間中心距離參數、每次迭代容許合併的最大聚類對數L及容許迭代次數I等。

      參考了老王的課件。