matlab實作kmeans聚類算法

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(dis(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);
        end
end
           

下面一段代碼用于測試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);
toc
M=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