ITKeyword,专注技术干货聚合推荐

注册 | 登录

Gap Statistic算法 Matlab代码

hyanglu1573 分享于 2013-10-19

推荐:matlab版K均值聚类

先看一篇关于K均值聚类的论文。 链接http://www.cnblogs.com/CBDoctor/archive/2011/10/24/2222358.html K均值聚类作用于灰度图上相当于根据灰度级别进行图像的

2020腾讯云8月秒杀活动,优惠非常大!(领取2860元代金券),
地址https://cloud.tencent.com/act/cps/redirect?redirect=1040

2020阿里云最低价产品入口,含代金券(新老用户有优惠),
地址https://www.aliyun.com/minisite/goods

总共三个文件

gap_stat.m

clear all;close all;
SampleNum=30;                                %样本观测点的数量
B=1000;                                      %参考数据集的数量
MaxK=10;                                     %最大的聚类数
u=[0 0;4 4;4 -4];                            %各类样本的均值向量
sigma=[1 0;0 1];                             %协方差矩阵


SampleSet=[];
[uM,~]=size(u);
for i=1:uM                                   %生成样本数据
    SampleSet=[SampleSet;mvnrnd(u(i,:),sigma,fix(SampleNum/uM))];
end
figure(1);
plot(SampleSet(:,1),SampleSet(:,2),'.');     %画二维样本分布图
title('样本数据');


Wk=log(CompuWk(SampleSet,MaxK));             %计算样本的Wk


interval=minmax(SampleSet');                 %样本各维度的取值区间
[SampleM,SampleN]=size(SampleSet);
RefSet=zeros(SampleM,SampleN,B);
for b=1:B                                    %生成参考数据集
    for i=1:SampleN
        TempSet(i,:)=unifrnd(interval(i,1),interval(i,2),1,SampleM);
    end
    RefSet(:,:,b)=TempSet';
end


for b=1:B                                   %计算参考数据集的Wkb
    Wkb(:,b)=log(CompuWk(RefSet(:,:,b),MaxK))';
end


for k=1:MaxK                                %计算Gap(k)
    Gap(k)=sum(Wkb(k,:))/B-Wk(k);
    l(k)=sum(Wkb(k,:))/B;
    sdk(k)=norm(Wkb(k,:)-l(k)*ones(1,B))/sqrt(B);
    sk(k)=sdk(k)*sqrt(1+1/B);
end


OptimusK=1;
for k=2:MaxK                                %计算最优的聚类数k
    if Gap(k)<=Gap(k-1)+sk(k)&&OptimusK==1
        OptimusK=k-1;
    end
end


figure(2);                                  %画Gap值和最优聚类数
plot(1:MaxK,Gap,'.-',OptimusK,Gap(OptimusK),'dr');
legend('Gap值','最优类别数','Location','NorthEast');
xlabel('类别数k');ylabel('Gap值');title('Gap值');




CompuWk.m

function Wk=CompuWk(SampleSet,MaxK)
%Compute Wk, the within cluster dispersion.
[SampleM,SampleN]=size(SampleSet);
%为了兼容Matlab自带的kmeans算法,这里先单独计算聚类数为1时的Wk,
%因为Matlab的kmeans函数不能计算聚类数为1的情况。
Dr=0;                                              %计算Wk(1)
for i=1:SampleM
    for j=i:SampleM
        Dr=Dr+norm(SampleSet(i,:)-SampleSet(j,:))^2;
    end
end
Wk(1)=0.5*Dr/SampleM;
for k=2:MaxK                                       %计算Wk(2:MaxK)
    %idx=kmeans(SampleSet,k,'emptyaction','drop'); %使用Matlab自带kmeans算法
    idx=cmeans(SampleSet,k);                       %使用自己编写的cmeans算法
    Dr=zeros(1,k);
    for i=1:k
        Cr=[];
        for j=1:SampleM
            if idx(j)==i
                Cr=[Cr;SampleSet(j,:)];
            end
        end
        [CrM,~]=size(Cr);
        if CrM~=0
            for m=1:CrM
                for n=m:CrM
                    Dr(i)=Dr(i)+norm(Cr(m,:)-Cr(n,:))^2;
                end
            end
            Dr(i)=0.5*Dr(i)/CrM;
        end
    end
    Wk(k)=sum(Dr);
end
end





cmeans.m

function idx=cmeans(y,c)
% C-means Algorithm.
% idx=cmeans(y,c).
% y: M*N Matrix,Sample Data.
% c: Cluster Number.
% idx: M*1 Matrix,contains cluster information.
[SampleM,SampleN] = size(y);
idx = zeros(SampleM,1);
CenterPre = y(1:c,:);
CenterNow = zeros(c,SampleN);
dis = zeros(SampleM,c);
w = zeros(SampleM,c);
while 1
    for i = 1:SampleM                    %计算类内离差距离
        for j = 1:c
            for k = 1:SampleN
                dis(i,j) = dis(i,j) + (y(i,k) - CenterPre(j,k))^2;
            end
        end   
    end 
    for i = 1:SampleM                    %更新分类情况
        k=1;
        for j = 1:c-1
            if dis(i,k) > dis(i,j+1)
                k = j+1;
            end
        end
        w(i,k) = 1;
        for l = 1:SampleN
            CenterNow(k,l) = CenterNow(k,l) + y(i,l);
        end
    end
    for i = 1:c                         %计算新的类心
        k = 0;
        for j = 1:SampleM
            k = k + w(j,i);
        end
        for j = 1:SampleN
            if k ~= 0
                CenterNow(i,j) = CenterNow(i,j)/k;
            end
        end
    end
    if CenterPre == CenterNow           %两次聚类后的类心相等,退出聚类
        break;
    end
    CenterPre = CenterNow;              %更新临时变量
    CenterNow = zeros(c,SampleN);
    dis = zeros(SampleM,c);
    w = zeros(SampleM,c);
end
for i=1:SampleM                         %idx记录样本分类情况
    for j=1:c
        if w(i,j)==1
            idx(i)=j;
        end
    end
end

总共三个文件 gap_stat.m clear all;close all; SampleNum=30;                                %样本观测点的数量 B=1000;                                      %参考数据集的数量 MaxK=10;

相关阅读排行


相关内容推荐

最新文章

×

×

请激活账号

为了能正常使用评论、编辑功能及以后陆续为用户提供的其他产品,请激活账号。

您的注册邮箱: 修改

重新发送激活邮件 进入我的邮箱

如果您没有收到激活邮件,请注意检查垃圾箱。