function [cx,cost] = kMedoids(K,data,num)

%   生成将data聚成K类的最佳聚类

%   K为聚类数目，data为数据集,num为随机初始化次数

    [cx,cost] = kMedoids1(K,data);

    for i = 2:num

        [cx1,min] = kMedoids1(K,data);

        if min<cost

            cost = min;

            cx = cx1;

        end

    end

end

function [cx,cost] = kMedoids1(K,data)

%   把分类数据集data聚成K类

%   [cx,cost] = kmeans(K,data)

%   K为聚类数目，data为数据集

%   cx为样本所属聚类，cost为此聚类的代价值

% 选择需要聚类的数目

% 随机选择聚类中心

    centroids = data(randperm(size(data,1),K),:);

% 迭代聚类

    centroids_temp = zeros(size(centroids));

    num = 0;

    while (~isequal(centroids_temp,centroids)&&num<20)

        centroids_temp = centroids;

        [cx,cost] = findClosest(data,centroids,K);

        centroids = compueCentroids(data,cx,K);

        num = num+1;

    end

%     cost = cost/size(data,1);

end

function [cx,cost] = findClosest(data,centroids,K)

% 将样本划分到最近的聚类中心

    cost = 0;

    n = size(data,1);

    cx = zeros(n,1);

    for i = 1:n

%       曼哈顿距离

        [M,I] = min(sum(abs((data(i,:)-centroids))'));

        cx(i) = I;

        cost = cost+M;

    end

end

function centroids = compueCentroids(data,cx,K)

% 计算新的聚类中心

    centroids = zeros(K,size(data,2));

    for i = 1:K

%       寻找代价值最小的当前聚类中心

        temp = data((cx==i),:);

        [~,I] = min(sum(squareform(pdist(temp))));

        centroids(i,:) = temp(I,:);

    end

end

% 主函数

% 生成符合高斯分布的数据

mu = [5,5];

sigma = [16,0;0,16];

sigma1 = [0.5,0;0,0.5];

data =  gaussianSample(8,50,mu,sigma,sigma1);

% 聚类

K = 6;

[cx,cost] = kMedoids(K,data,10);

plotMedoids(data,cx,K);