我的K均值演算法的matlab實現
這是我的第一篇部落格;
K-Means演算法過程,略;
這是一次課程的任務2333,是利用所學K-means聚類分析方法,對iris資料集進行聚類分析,並利用已知的樣本類別標
籤進行聚類分析評價;
我的K均值演算法以iris.data為例(附在文末);
資料集:Iris資料集
(http://archive.ics.uci.edu/ml/datasets/Iris)
資料描述:Iris資料集包含150個鳶尾花模式樣本,其中 每個模式樣本採用5維的特徵描述
X = (x1,x2,x3,x4,w);
x1: 萼片長度(釐米)
x2: 萼片寬度(釐米)
x3: 花瓣長度(釐米)
x4:花瓣寬度(釐米)
w(類別屬性 ): 山鳶尾 (Iris setosa),變色鳶尾(Iris versicolor)和維吉尼亞鳶尾(Iris virginica)
先貼上我的函式結構:
函式結構
—— FindCluster(~) 聚類演算法主函式
|
MyKmeans —— MyPlot2(~),MyPlot3(~) 畫圖
|
—— Accuracy(~) 聚類精度評價
MyKmean.m,程式執行的主函式
% 作者 : YG1501 LQY % 日期 : 2018.01.13 星期六 % 函式功能 : 實現對iris.data資料集的分類,並根據分類結果進行精度評價 clear;clc;close all; %手動選取檔案,選用iris.data [filename,fpath] = uigetfile(... {... '*.m;*.txt;*.data',... 'Data(*.m;*.txt;*.data)';... '*.*','All Files(*.*)'... },'Select data'); % filename = 'iris.data'; %嫌麻煩可以用這個 [X1,X2,X3,X4,X5] = ... textread(filename,'%f%f%f%f%s','delimiter',','); %Get data clear filename fpath X = [X1 X2 X3 X4]; [m,n] = size(X); DataLabel = zeros(m,1); for i = 1:size(X5,1) if (strcmp(X5(i),'Iris-setosa')) DataLabel(i) = 1 ; elseif(strcmp(X5(i),'Iris-versicolor')) DataLabel(i) = 2 ; else DataLabel(i) = 3 ; end end clear m n i %二維結果 % [MyCenter1,ClusterLabel12] = FindCluster([X1 X2],3,DataLabel); % [MyCenter2,ClusterLabel13] = FindCluster([X1 X3],3,DataLabel); % [MyCenter3,ClusterLabel14] = FindCluster([X1 X4],3,DataLabel); % [MyCenter4,ClusterLabel23] = FindCluster([X2 X3],3,DataLabel); % [MyCenter5,ClusterLabel24] = FindCluster([X2 X4],3,DataLabel); % [MyCenter6,ClusterLabel34] = FindCluster([X3 X4],3,DataLabel); % hold on; % figure(1),title('2D'); % subplot(231) % MyPlot2(X1,X2,DataLabel,MyCenter1),xlabel('X1'),ylabel('X2') % subplot(232) % MyPlot2(X1,X3,DataLabel,MyCenter2),xlabel('X1'),ylabel('X3') % subplot(233) % MyPlot2(X1,X4,DataLabel,MyCenter3),xlabel('X1'),ylabel('X4') % subplot(234) % MyPlot2(X2,X3,DataLabel,MyCenter4),xlabel('X2'),ylabel('X3') % subplot(235) % MyPlot2(X2,X4,DataLabel,MyCenter5),xlabel('X2'),ylabel('X4') % subplot(236) % MyPlot2(X3,X4,DataLabel,MyCenter6),xlabel('X3'),ylabel('X4') %三維結果 [MyCenter7,ClusterLabel123] = FindCluster([X1,X2,X3],3,DataLabel); [MyCenter8,ClusterLabel124] = FindCluster([X1,X2,X4],3,DataLabel); [MyCenter9,ClusterLabel134] = FindCluster([X1,X3,X4],3,DataLabel); [MyCenter10,ClusterLabel234] = FindCluster([X2,X3,X4],3,DataLabel); hold on; figure(2),title('3D'); subplot(221) MyPlot3(X1,X2,X3,DataLabel,MyCenter7) xlabel('X1'),ylabel('X2'),zlabel('X3'); subplot(222) MyPlot3(X1,X2,X4,DataLabel,MyCenter8) xlabel('X1'),ylabel('X2'),zlabel('X4'); subplot(223) MyPlot3(X1,X3,X4,DataLabel,MyCenter9) xlabel('X1'),ylabel('X3'),zlabel('X4'); subplot(224) MyPlot3(X2,X3,X4,DataLabel,MyCenter10) xlabel('X2'),ylabel('X3'),zlabel('X4'); %聚類精度評價 %二維結果 % ClusterLabel = [ClusterLabel12,ClusterLabel13,ClusterLabel14... % ClusterLabel23,ClusterLabel24,ClusterLabel34] % ClusterAccuracy = Accuracy(DataLabel,ClusterLabel) %三維結果 ClusterLabel = [ClusterLabel123,ClusterLabel124,ClusterLabel134,ClusterLabel234]; ClusterAccuracy = Accuracy(DataLabel,ClusterLabel)
FindCluster.m
%函式功能 : 輸入資料集、聚類中心個數與樣本標籤 % 得到聚類中心與聚類樣本標籤 function [ClusterCenter,ClusterLabel] = FindCluster(MyData,ClusterCounts,DataLabel) [m,n] = size(MyData); ClusterLabel = zeros(m,1); %用於儲存聚類標籤 % MyLabel = unique(DataLabel,'rows'); % for i = 1:size(MyLabel,2); % LabelIndex(1,i) = i; %為資料標籤建立索引 % end %已知資料集的每個樣本的中心 OriginCenter = zeros(ClusterCounts,n); for j = 1:ClusterCounts DataCounts = 0; for i = 1:m %按照資料標籤,計算樣本中心 if DataLabel(i) == j OriginCenter(j,:) = OriginCenter(j,:) + MyData(i,:); DataCounts = DataCounts + 1; end end OriginCenter(j,:) = OriginCenter(j,:) ./ DataCounts; end %按照第一列對樣本中心排序 SortCenter1 = sortrows(OriginCenter,1); FalseTimes = 0; CalcuateTimes = 0; %此迴圈用於糾正分類錯誤的情況 while (CalcuateTimes < 15) %這裡注意如果顯示 @It;, 把它改成小於號 ClusterCenter = zeros(ClusterCounts,n); %初始化聚類中心 for i = 1:ClusterCounts ClusterCenter(i,:) = MyData( randi(m,1),:); %隨機選取一個點作為中心 end %此迴圈用於尋找聚類中心 %目前還未解決該迴圈陷入死迴圈的問題,所以設定一個引數來終止迴圈 kk = 0; while (kk < 15) %這裡注意如果顯示 @It;, 把它改成小於號 Distance = zeros(1,ClusterCounts); %儲存單個樣本到每個聚類中心的距離 DataCounts = zeros(1,ClusterCounts); %記錄每個聚類的樣本數目 NewCenter = zeros(ClusterCounts,n); for i = 1:m for j = 1:ClusterCounts Distance(j) = norm(MyData(i,:) - ClusterCenter(j,:)); end %index返回最小距離的索引,亦即聚類中心的標號 [~,index] = min(Distance); ClusterLabel(i) = index; end k = 0; for j = 1:ClusterCounts for i = 1:m %按照標記,對樣本進行分類,並計算聚類中心 if ClusterLabel(i) == j NewCenter(j,:) = NewCenter(j,:) + MyData(i,:); DataCounts(j) = DataCounts(j) + 1; end end NewCenter(j,:) = NewCenter(j,:) ./ DataCounts(j); %若新的聚類中心與上一個聚類中心相同,則認為演算法收斂 if norm(NewCenter(j,:) - ClusterCenter(j,:)) < 0.1 %這裡注意如果顯示 @It;, 把它改成小於號 k = k + 1; end end ClusterCenter = NewCenter; %判斷是否全部收斂 if k == ClusterCounts break; end kk = kk + 1 ; end %再次判斷每個聚類是否分類正確,若分類錯誤則進行懲罰 t = ClusterCounts; SortCenter2 = sortrows(ClusterCenter,1); for i = 1:ClusterCounts if norm(SortCenter1(i,:) - SortCenter2(i,:)) > 0.5 %這裡注意如果顯示 @gt;, 把它改成大於號 t = t - 1; FalseTimes = FalseTimes + 1; break; end end if t == ClusterCounts break; end CalcuateTimes = CalcuateTimes + 1; end % FalseTimes % CalcuateTimes % t % kk % DataCounts % OriginCenter % NewCenter %理論上每個聚類的標籤應是123排列的,但實際上,由於每個聚類中心都是隨機選取的, %實際分類的順序可能是213,132等,所以需要對分類標籤進行糾正,這對之後的精度評 %價帶來了方便。如果不需要進行精度評價,這一段可以不要 %對分類標籤進行糾正: %演算法原理:從第一個已知的樣本中心開始,尋找離其最近的聚類中心,然後將歸類於該 % 聚類中心的樣本的聚類標籤更換為i for i = 1:ClusterCounts for j = 1:ClusterCounts if norm(OriginCenter(i,:) - ClusterCenter(j,:)) < 0.6 %這裡注意如果顯示 @It;, 把它改成小於號 for p = 1:m if ClusterLabel(p) == j ClusterLabel(p) = 2 * ClusterCounts + i; end end %break; end end end ClusterLabel = ClusterLabel - 2 * ClusterCounts; %Temp = [MyData(:,:),ClusterLabel]
至此已經完成了聚類中心的計算。該方法基本解決了因隨機選取初始中心而導致最後聚類中心明顯錯誤的情況,但缺點在於迴圈太多,時間複雜度O(?),而且偶爾會陷入死迴圈,原因在於有兩個或者以上的聚類中心被選到了同一點。若判明進入了死迴圈,還是重新執行下程式吧
畫圖的函式:
MyPlot2.m
%函式功能 : 輸入樣本,樣本標籤及求出的聚類中心,顯示二維影象,實現資料視覺化
function MyPlot2(X1,X2,DataLabel,ClusterCenter)
[m,~] = size(X1);
hold on;
p = find(DataLabels == 1);
plot(X1(p),X2(p),'*r')
p = find(DataLabels == 2);
plot(X1(p),X2(p),'*g')
p = find(DataLabels == 3);
plot(X1(p),X2(p),'*b')
% xlabel(who('X1'));
% ylabel(who('X2'));
%PS : 我想在座標軸中根據輸入的形參的矩陣的名字,轉換成字串,來動態輸出
% 座標軸名稱,不知道該怎麼做?用上面註釋的語句不行。。
[n,~] = size(ClusterCenter);
for i = 1:n
plot(ClusterCenter(i,1),ClusterCenter(i,2),'ok')
end
grid on;
MyPlot3.m
function MyPlot3(X1,X2,X3,DataLabels,ClusterCenter)
[m,~] = size(X1);
hold on;
p = find(DataLabels == 1);
plot3(X1(p),X2(p),X3(p),'.r')
p = find(DataLabels == 2);
plot3(X1(p),X2(p),X3(p),'.g')
p = find(DataLabels == 3);
plot3(X1(p),X2(p),X3(p),'.b')
[n,~] = size(ClusterCenter);
plot3(ClusterCenter(1:1:n,1),ClusterCenter(1:1:n,2),ClusterCenter(1:1:n,3),'ok')
view([1 1 1]);
grid on;
精度評價Accuracy.m
%函式功能:根據聚類結果進行精度評價
%精度評價,返回每一種分類的精度值(正確率)
function ClusterAccuracy = Accuracy(DataLable,CLusterLabel)
[m,n] = size(CLusterLabel);
ClusterAccuracy = zeros(1,n);
%理論上的聚類標籤應為 1,2,3,但實際上可能變成了 213 , 132等,導致計算失誤
%因此需要對分類標籤進行糾正,而這一步驟已經在FindCluster函式中完成了
for i = 1:n
%原理:假設某樣本在已知資料集中屬於第一類,而其聚類後的也同樣被分到了第一類,那麼它們的標籤
%都是1,這樣相減後結果就為0,表明已經分類正確,否則不為0,分類錯誤
Temp(:,i) = DataLable - CLusterLabel(:,i);
end
for j = 1:n
for i = 1:m
if Temp(i,j) == 0
ClusterAccuracy(1,j) = ClusterAccuracy(1,j) + 1;
end
end
end
ClusterAccuracy = ClusterAccuracy ./ m;
執行結果展示(二維結果):
可以看到,在二維的情況下,比較X3: 花瓣長度(釐米)和X4:花瓣寬度(釐米)的精度更高(94.67%),也就是說只比較兩種特徵時,比較花瓣長度和花瓣寬度,區分三種花的效果更好,分類結果更可靠。
三維分類結果:
可以看到,在三維的情況下,比較X2:萼片寬度(釐米),X3: 花瓣長度(釐米)和X4:花瓣寬度(釐米)的精度更高(93.33%),也就是說比較三種特徵時,取X2,X3,X4,區分三種花的效果更好,分類結果更可靠。
最後:本程式目前還存在諸多不足,比如時間複雜度高,效率較低;目前對Matlab語言還不是很熟,寫的程式也比較C-Style,各位看官若是有改進的建議,歡迎留言,或者直接聯絡我(聯絡方式附在文最末),大家一起探討:D
附iris.data資料集:
5.1,3.5,1.4,0.2,Iris-setosa
4.9,3.0,1.4,0.2,Iris-setosa
4.7,3.2,1.3,0.2,Iris-setosa
4.6,3.1,1.5,0.2,Iris-setosa
5.0,3.6,1.4,0.2,Iris-setosa
5.4,3.9,1.7,0.4,Iris-setosa
4.6,3.4,1.4,0.3,Iris-setosa
5.0,3.4,1.5,0.2,Iris-setosa
4.4,2.9,1.4,0.2,Iris-setosa
4.9,3.1,1.5,0.1,Iris-setosa
5.4,3.7,1.5,0.2,Iris-setosa
4.8,3.4,1.6,0.2,Iris-setosa
4.8,3.0,1.4,0.1,Iris-setosa
4.3,3.0,1.1,0.1,Iris-setosa
5.8,4.0,1.2,0.2,Iris-setosa
5.7,4.4,1.5,0.4,Iris-setosa
5.4,3.9,1.3,0.4,Iris-setosa
5.1,3.5,1.4,0.3,Iris-setosa
5.7,3.8,1.7,0.3,Iris-setosa
5.1,3.8,1.5,0.3,Iris-setosa
5.4,3.4,1.7,0.2,Iris-setosa
5.1,3.7,1.5,0.4,Iris-setosa
4.6,3.6,1.0,0.2,Iris-setosa
5.1,3.3,1.7,0.5,Iris-setosa
4.8,3.4,1.9,0.2,Iris-setosa
5.0,3.0,1.6,0.2,Iris-setosa
5.0,3.4,1.6,0.4,Iris-setosa
5.2,3.5,1.5,0.2,Iris-setosa
5.2,3.4,1.4,0.2,Iris-setosa
4.7,3.2,1.6,0.2,Iris-setosa
4.8,3.1,1.6,0.2,Iris-setosa
5.4,3.4,1.5,0.4,Iris-setosa
5.2,4.1,1.5,0.1,Iris-setosa
5.5,4.2,1.4,0.2,Iris-setosa
4.9,3.1,1.5,0.1,Iris-setosa
5.0,3.2,1.2,0.2,Iris-setosa
5.5,3.5,1.3,0.2,Iris-setosa
4.9,3.1,1.5,0.1,Iris-setosa
4.4,3.0,1.3,0.2,Iris-setosa
5.1,3.4,1.5,0.2,Iris-setosa
5.0,3.5,1.3,0.3,Iris-setosa
4.5,2.3,1.3,0.3,Iris-setosa
4.4,3.2,1.3,0.2,Iris-setosa
5.0,3.5,1.6,0.6,Iris-setosa
5.1,3.8,1.9,0.4,Iris-setosa
4.8,3.0,1.4,0.3,Iris-setosa
5.1,3.8,1.6,0.2,Iris-setosa
4.6,3.2,1.4,0.2,Iris-setosa
5.3,3.7,1.5,0.2,Iris-setosa
5.0,3.3,1.4,0.2,Iris-setosa
7.0,3.2,4.7,1.4,Iris-versicolor
6.4,3.2,4.5,1.5,Iris-versicolor
6.9,3.1,4.9,1.5,Iris-versicolor
5.5,2.3,4.0,1.3,Iris-versicolor
6.5,2.8,4.6,1.5,Iris-versicolor
5.7,2.8,4.5,1.3,Iris-versicolor
6.3,3.3,4.7,1.6,Iris-versicolor
4.9,2.4,3.3,1.0,Iris-versicolor
6.6,2.9,4.6,1.3,Iris-versicolor
5.2,2.7,3.9,1.4,Iris-versicolor
5.0,2.0,3.5,1.0,Iris-versicolor
5.9,3.0,4.2,1.5,Iris-versicolor
6.0,2.2,4.0,1.0,Iris-versicolor
6.1,2.9,4.7,1.4,Iris-versicolor
5.6,2.9,3.6,1.3,Iris-versicolor
6.7,3.1,4.4,1.4,Iris-versicolor
5.6,3.0,4.5,1.5,Iris-versicolor
5.8,2.7,4.1,1.0,Iris-versicolor
6.2,2.2,4.5,1.5,Iris-versicolor
5.6,2.5,3.9,1.1,Iris-versicolor
5.9,3.2,4.8,1.8,Iris-versicolor
6.1,2.8,4.0,1.3,Iris-versicolor
6.3,2.5,4.9,1.5,Iris-versicolor
6.1,2.8,4.7,1.2,Iris-versicolor
6.4,2.9,4.3,1.3,Iris-versicolor
6.6,3.0,4.4,1.4,Iris-versicolor
6.8,2.8,4.8,1.4,Iris-versicolor
6.7,3.0,5.0,1.7,Iris-versicolor
6.0,2.9,4.5,1.5,Iris-versicolor
5.7,2.6,3.5,1.0,Iris-versicolor
5.5,2.4,3.8,1.1,Iris-versicolor
5.5,2.4,3.7,1.0,Iris-versicolor
5.8,2.7,3.9,1.2,Iris-versicolor
6.0,2.7,5.1,1.6,Iris-versicolor
5.4,3.0,4.5,1.5,Iris-versicolor
6.0,3.4,4.5,1.6,Iris-versicolor
6.7,3.1,4.7,1.5,Iris-versicolor
6.3,2.3,4.4,1.3,Iris-versicolor
5.6,3.0,4.1,1.3,Iris-versicolor
5.5,2.5,4.0,1.3,Iris-versicolor
5.5,2.6,4.4,1.2,Iris-versicolor
6.1,3.0,4.6,1.4,Iris-versicolor
5.8,2.6,4.0,1.2,Iris-versicolor
5.0,2.3,3.3,1.0,Iris-versicolor
5.6,2.7,4.2,1.3,Iris-versicolor
5.7,3.0,4.2,1.2,Iris-versicolor
5.7,2.9,4.2,1.3,Iris-versicolor
6.2,2.9,4.3,1.3,Iris-versicolor
5.1,2.5,3.0,1.1,Iris-versicolor
5.7,2.8,4.1,1.3,Iris-versicolor
6.3,3.3,6.0,2.5,Iris-virginica
5.8,2.7,5.1,1.9,Iris-virginica
7.1,3.0,5.9,2.1,Iris-virginica
6.3,2.9,5.6,1.8,Iris-virginica
6.5,3.0,5.8,2.2,Iris-virginica
7.6,3.0,6.6,2.1,Iris-virginica
4.9,2.5,4.5,1.7,Iris-virginica
7.3,2.9,6.3,1.8,Iris-virginica
6.7,2.5,5.8,1.8,Iris-virginica
7.2,3.6,6.1,2.5,Iris-virginica
6.5,3.2,5.1,2.0,Iris-virginica
6.4,2.7,5.3,1.9,Iris-virginica
6.8,3.0,5.5,2.1,Iris-virginica
5.7,2.5,5.0,2.0,Iris-virginica
5.8,2.8,5.1,2.4,Iris-virginica
6.4,3.2,5.3,2.3,Iris-virginica
6.5,3.0,5.5,1.8,Iris-virginica
7.7,3.8,6.7,2.2,Iris-virginica
7.7,2.6,6.9,2.3,Iris-virginica
6.0,2.2,5.0,1.5,Iris-virginica
6.9,3.2,5.7,2.3,Iris-virginica
5.6,2.8,4.9,2.0,Iris-virginica
7.7,2.8,6.7,2.0,Iris-virginica
6.3,2.7,4.9,1.8,Iris-virginica
6.7,3.3,5.7,2.1,Iris-virginica
7.2,3.2,6.0,1.8,Iris-virginica
6.2,2.8,4.8,1.8,Iris-virginica
6.1,3.0,4.9,1.8,Iris-virginica
6.4,2.8,5.6,2.1,Iris-virginica
7.2,3.0,5.8,1.6,Iris-virginica
7.4,2.8,6.1,1.9,Iris-virginica
7.9,3.8,6.4,2.0,Iris-virginica
6.4,2.8,5.6,2.2,Iris-virginica
6.3,2.8,5.1,1.5,Iris-virginica
6.1,2.6,5.6,1.4,Iris-virginica
7.7,3.0,6.1,2.3,Iris-virginica
6.3,3.4,5.6,2.4,Iris-virginica
6.4,3.1,5.5,1.8,Iris-virginica
6.0,3.0,4.8,1.8,Iris-virginica
6.9,3.1,5.4,2.1,Iris-virginica
6.7,3.1,5.6,2.4,Iris-virginica
6.9,3.1,5.1,2.3,Iris-virginica
5.8,2.7,5.1,1.9,Iris-virginica
6.8,3.2,5.9,2.3,Iris-virginica
6.7,3.3,5.7,2.5,Iris-virginica
6.7,3.0,5.2,2.3,Iris-virginica
6.3,2.5,5.0,1.9,Iris-virginica
6.5,3.0,5.2,2.0,Iris-virginica
6.2,3.4,5.4,2.3,Iris-virginica
5.9,3.0,5.1,1.8,Iris-virginica
程式碼與資料可直接下載,
連結:https://pan.baidu.com/s/1xEGpUIDA2-ayPtxReKrd-w
提取碼:yh65
博主聯絡方式 : QQ:1765928683
2018 / 01 / 27 補充:
我想到了一些提高執行效率的方法,那就是儘量減少for迴圈,因為matlab對for迴圈的處理效率是很低的!
可拱參考的優化思路:
1) 使用parfor;
2)使用find;
3)向量化,即使用向量操作代替迴圈;
4)bsxfun(),sum(), ‘./’ '.*'
5)待補充
2018 / 05 / 03 更新:
程式之所以跑得慢,是因為畫圖函式用了for迴圈,把for迴圈去掉直接用plot,會快很多。
2018/07/12更新:
優化了畫圖函式,現在可以很快看到結果了;
直接從網頁把程式碼複製過去的話,函式執行會出問題,修改方法如下:看到FindCluster函式裡,找到所有while迴圈語句和 if 判斷語句,把 @It; 改成 < , 把 @gt; 改成 > 就行了(這個問題我解決不了,修改的時候能正確顯示,但儲存後就出問題了.....)