kmeans聚類算法是一種簡單實用的聚類算法，matlab自帶函數kmeans可直接對數據進行kmeans聚類。為了方便更好地掌握kmeans聚類算法，今天我們自己來實現一個弱化的版本mykmeans。mykmeans輸入包含三項，分別為聚類所使用的數據data，data每一行代表一個樣本，每一列代表一個特徵；聚類中心數量numclass；第三項為所使用的距離的定義，默認情況下為歐式距離。
function [cluster,clusterhistory]=mykmeans(data,numclass,varargin)% kmeans聚類。% 聚類過程動畫顯示% INPUTS:% data：每一行代表一個樣本，每一列代表一個特徵% numclass：聚類中心的數量% varargin{1}: 距離定義。支持euclidean(默認)、hamming% OUTPUT:% cluster: numsample行的列向量，代表每個樣本的分類歸屬% clusterhistory{i}.cluster第i次迭代的聚類% clusterhistory{i}.core第i次迭代的聚類中心%% 公眾號【數學建模公會】，HCLO4，20190902
if nargin==2 method='euclidean';else method=varargin{1};end% 初始化聚類中心坐標，在數據點中隨機選擇numclass個點作為初始聚類中心。numsample=size(data,1);temp=randperm(numsample);core=data(temp(1:numclass),:);dis=caldis(data,core,method); %存儲每個樣本到當前聚類中心的距離
% 執行迭代過程maxIter=20; % 最大迭代次數numiter=1;clusterhistory=cell(1,maxIter);while 1 newcore=zeros(size(core)); % 迭代聚類中心 [~,ind]=min(dis,[],2); %計算每個sample歸屬於哪個聚類中心，如果某個聚類中心沒有一個點？ for i=1:numclass newcore(i,:)=mean(data(ind==i,:)); end clusterhistory{numiter}.cluster=ind; clusterhistory{numiter}.core=core; if all(newcore(:)==core(:))||numiter>=maxIter % 迭代終止條件，聚類中心不再改變 cluster=ind; break end core=newcore; dis=caldis(data,core,method); numiter=numiter+1; end
clusterhistory=clusterhistory(1:numiter);
end % mykmeans
function dis=caldis(data,core,method)%計算每個樣本到當前聚類中心的距離numsample=size(data,1);numclass=size(core,1);dis=zeros(numsample,numclass);switch method case 'euclidean' for i=1:numclass dis(:,i)=sqrt(sum((data-repmat(core(i,:),numsample,1)).^2,2)); end case 'hamming' for i=1:numclass dis(:,i)=mean(data~=repmat(core(i,:),numsample,1),2); endend下面一段代碼用於測試mykmeans和kmeans的運行結果。生成兩組服從二維高斯分布的數據，作為兩個數據類別，分別使用mykmeans和matlab自帶的kmeans函數進行聚類分析。%%%%%%%%%%%%%%%%%%%%%%% 測試mykmeans函數 %%%%%%%%%%%%%%%%%%%%%%%%% 公眾號【數學建模公會】，HCLO4，20190902
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 生成模擬數據,平面上的兩組點，均服從二維高斯分布mu1=[2,8];sigma1=[2,2];mu2=[8,2];sigma2=[3,3];
numsample1=100;numsample2=200;data1=zeros(numsample1,2);data1(:,1)=normrnd(mu1(1),sigma1(1),numsample1,1);data1(:,2)=normrnd(mu1(2),sigma1(2),numsample1,1);
data2=zeros(numsample2,2);data2(:,1)=normrnd(mu2(1),sigma2(1),numsample2,1);data2(:,2)=normrnd(mu2(2),sigma2(2),numsample2,1);
data=[data1;data2];
tic,[cluster1,clusterhistory]=mykmeans(data,2);tocM=getFrame_Kmeans(data,clusterhistory); % 生成聚類動圖
tic,cluster2=kmeans(data,2);toc下面的函數實現mykmeans聚類過程的可視化，僅針對二維數據和三維數據類型。function M=getFrame_Kmeans(data,clusterhistory)% kmeans聚類過程可視化程序。只能對二維和三維數據可視化！%% INPUTS: % data: 聚類用的數據，每行代表一個樣本，每列代表一個特徵% clusterhistory: mykmeans函數輸出的聚類過程數據%% OUTPUT:% kmeans聚類過程的動畫。
dimension=size(data,2);if dimension>3 % 若三維以上，只能通過pca降維處理 error('無法實現高於三維的數據的聚類可視化')end
colorset=cell(1,1000);colorset(1:8)={'r','g','b','k','c','m','y','k'};for i=9:1000 %如果聚類中心數量大於8，就使用隨機的顏色。 colorset{i}=rand(1,3);end
numcore=length(unique(clusterhistory{1}.cluster));numiter=length(clusterhistory);
for k=1:numiter figure hold on switch dimension case 2 for i=1:numcore ind=clusterhistory{k}.cluster==i; scatter(data(ind,1),data(ind,2),6,colorset{i}) core=clusterhistory{k}.core; plot(core(i,1),core(i,2),[colorset{i},'.'],'MarkerSize',20) end case 3 for i=1:numcore ind=clusterhistory{k}.cluster==i; scatter3(data(ind,1),data(ind,2),data(ind,3),6,colorset{i}) core=clusterhistory{k}.core; plot3(core(i,1),core(i,2),core(i,3),colorset{i},'MarkerSize',6) end end M(k)=getframe(gcf); frame = getframe(gcf); im = frame2im(frame); %將影片動畫轉換為編址圖像,因為圖像必須是index索引圖像 imshow(im); [I,map] = rgb2ind(im,20); %將真彩色圖像轉化為索引圖像 if k==1 imwrite(I,map,'kmeans.gif','gif','Loopcount',inf,'DelayTime',0.3); %Loopcount只是在i==1的時候才有用 else imwrite(I,map,'kmeans.gif','gif','WriteMode','append','DelayTime',1);%DelayTime:幀與幀之間的時間間隔 end close(gcf); end