function [data] = gaussianSample(n,m,mu,sigma,sigma1)
%   生成n个符合多元高斯分布的样本
%   data = gaussianSample(n,m,mu,sigma,mu1,sigma1)
%   生成一个符合高斯分布的数据集data
%   n表示生成的数量，m表示每个簇的点的数量
%   先通过mu与sigma生成簇中心的分布情况
%   在通过簇中心分布mu1与sigma1生成每个簇的的分布情况
%   mu与mu1为均值，sigma与sigma1为协方差矩阵。


% % 默认值生成2维数据
% m = 100;
% mu = [5,5];
% sigma = [16 0;0 16];
% sigma1 = [0.5,0;0,0.5];

% 生成中心点分布
mu1 =  mvnrnd(mu,sigma,n);
% 生成数据
data = [];
for i = 1:n
    temp = mvnrnd(mu1(i,:),sigma1,m);
    data = [data;temp];
end

% % 可视化样本,以二维数据为例
% hold on;
% plot(data(:,1),data(:,2),\'ko\',\'MarkerFaceColor\',\'y\');
% plot(mu1(:,1),mu1(:,2),\'r+\',\'LineWidth\',2,\'MarkerSize\',7);
% hold off;

function [cx,cost] = myKmeans(K,data,num)
% 生成将data聚成K类的最佳聚类
    [cx,cost] = kmeans1(K,data);
    for i = 2:num
        [cx1,min] = kmeans1(K,data);
        if min<cost
            cost = min;
            cx = cx1;
        end
    end
%     plotMeans(data,cx,K); 
end

function [cx,cost] = kmeans1(K,data)
%KMEANS 把数据集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((centroids-data(i,:))\'.^2));
        cx(i) = I;
        cost = cost+M;
    end
  
end



function centroids = compueCentroids(data,cx,K)
% 计算新的聚类中心
    centroids = zeros(K,size(data,2));
    for i = 1:K
        centroids(i,:) = mean(data(cx==i,:));
    end
end

function [num,cx] = bestKmeans(data,m,n)
%   返回数据集的最佳聚类数目与聚类结果
%   data为数据集，m为寻找的最大聚类数量，n为每次聚类寻找次数
%   返回num最佳聚类数量，cx聚类结果。
costs = zeros(m,1)\';
for i = 1:m
    [~,cost] = myKmeans(i,data,n);
    costs(i) = cost;
    fprintf(\'最小代价值为：%.4f\n\',cost);
end
 costs = costs./costs(2)*m;
 % 绘制拐点图
X = (1:m)\';
X = [X,costs\'];

plot(X(2:m,1),X(2:m,2),\'ko\',\'MarkerFaceColor\',\'y\',\'MarkerSize\',2);
title(\'Inflection point map\');
% 寻找最佳聚类
min = -1;
for i = 3:m
    x1 = X(2,:)-X(i,:);
    x2 = X(m,:)-X(i,:);
    k = x1*x2\'/((x1*x1\')*(x2*x2\'));
    if k>min&&k<0
        num = i;
        min = k;
    end
end
% fprintf(\'聚类结果为：%d\n\',num);
[cx,~] = myKmeans(num,data,n);
end

function [] = plotMeans(data,cx,K)
% 可视化数据聚类效果
    figure;
    color=\'cbygmkr\';
    hold on;
    for i = 1:K
        px = find(cx==i);
        plot(data(px,1),data(px,2),\'ko\',\'MarkerFaceColor\',color(i),\'MarkerSize\',5);
    end
    title(\'K-means\')
    hold off;
end

% 主函数

% 生成符合高斯分布的数据
mu = [5,5];
sigma = [16,0;0,16];
sigma1 = [0.5,0;0,0.5];
data =  gaussianSample(4,50,mu,sigma,sigma1);

% 计算最佳聚类数量与结果
[num,cx] = bestKmeans(data,10,10);
fprintf(\'聚类结果为：%d\n\',num);
plotMeans(data,cx,num);