MATLAB数学建模——k-means++及可视化_matlab 绘制kmeans迭代图

% data set;随机初始化一个数据集
clc;clear;
% Sigma = [1, 0; 0, 1];
% mu1 = [1, -1];
% x1 = mvnrnd(mu1, Sigma, 1000);
% mu2 = [5, -4];
% x2 = mvnrnd(mu2, Sigma, 1000);
% mu3 = [1, 4];
% x3 = mvnrnd(mu3, Sigma, 1000);
% mu4 = [6, 4];
% x4 = mvnrnd(mu4, Sigma, 1000);
% mu5 = [7, 0.0];
% x5 = mvnrnd(mu5, Sigma, 1000);
 
 
x1 = randn(1000,6);
x2 = randn(1000,6)+1;
x3 = randn(1000,6)+2;
x4 = randn(1000,6)+3;
x5 = randn(1000,6)+4;
 
X = [x1;x2;x3;x4;x5];
func_zbf(X,8);
[idx,mu,d] = func_k_means_pro(X,5,1); % 使用自定义的K-means++函数

function func_zbf(X,max_k)
%%% ZBF:肘部法，作出不同K值聚类偏差图
% X:数据矩阵
% max_k:最大聚类数
 
dk = zeros(max_k,2); % 初始化聚类数及对应的偏差矩阵
for i = 1 : max_k % 最小聚类数设置为1
    [L,~,D] = func_k_means_pro(X,i,0);
%       [L,~,~,D] =kmeans(zscore(X),i);  
    for j = 1 : i
        d = D(:,j);
        tD1 = L < j; 
        tD2 = L > j;
        tD = tD1 + tD2;
        tD = all(tD,2);
        d(tD,:) = [];
        dk(i,2) = dk(i,2) + sum(d.^2);
    end
    dk(i,1) = i;
end
 
figure
plot(dk(2:end,1),dk(2:end,2),'-','LineWidth',1.2,'Color','b')
hold on
grid off
plot(dk(2:end,1),dk(2:end,2),'or','markerfacecolor', 'c',"LineWidth",1.1);
% title('不同K值聚类偏差图') 
set(gca, 'linewidth', 1.1, 'fontsize', 14, 'fontname', '宋体')
xlabel('分类数(K值)') 
ylabel('簇内误差平方和') 

function [indexs, mu, D] = func_k_means_pro(X,k,trigger)
%%% K_MEANS_PRO: k-means++聚类
% X: 数据集，每一行为一条数据
% k: 聚类数
% mu: 最终的聚类中心
% D: 数据点与各聚类中心的距离
% indexs: 数据所对应的类别
% trigger: 选择是否可视化
 
%% 数据标准化(可选)
X = zscore(X);
% X = (mapminmax(X',0,1))';
 
%% 获取数据总数
[n,length] = size(X);
 
%% 随机选择一个数据点作为第一个聚类中心
init_mu = X(randperm(size(X, 1), 1), :);
% 选择剩余的聚类中心
for i = 2 : k
    % 计算各数据点与已有聚类中心的平方欧式距离
    Dis = pdist2(X, init_mu, 'squaredeuclidean'); 
    Dis = min(Dis, [], 2); % 找出各数据点与最近聚类中心的平方欧式距离
    Dis = Dis / sum(Dis); % 计算各数据点被选为下一个聚类中心的概率
    %{ 
    如果累积概率大于随机数，则选取该点为新的聚类中心，确保初始聚类中心间的
    距离足够大
    %}
%     init_mu(i, :) = X(find(rand < cumsum(Dis), 1), :); 
    init_mu(i, :) = X(find(Dis == max(Dis), 1), :); 
end
 
%% 计算迭代，更新聚类中心
indexs = zeros(n, 1);
mu = init_mu;
new_mu = mu;
eps = 1e-6; % 设置精度
delta = 1;
while (delta > eps) % 若聚类中心不再变化则终止迭代
    mu = new_mu;
    for i =1:n
        y = repmat(X(i, :), k, 1);
        dist = y - mu;
        d = sum(dist.*dist,2); % 计算欧式距离
        j = find(d==min(d)); % 寻找最小距离对应的聚类中心
        indexs(i) = j;
    end
    for j = 1 : k
        order = indexs == j;
        new_mu(j, :) = mean(X(order, :), 1); % 以均值作为新的聚类中心
    end
 
    % 计算新旧聚类中心的偏移距离
    delta = sqrt(sum(sum((mu-new_mu).*(mu-new_mu))));
end
 
%% 根据最终聚类中心确定样本聚类
indexs = zeros(n, 1);
for i = 1 : n
    R = repmat(X(i,:),k,1) - mu;
    Res = sum(R.*R,2);
    j = find(Res == min(Res));
    indexs(i) = j;
end
 
%% 计算数据点与各个聚类中心的距离
D = zeros(n,k);
for i = 1 : n
    for j = 1 : k
        d = X(i,:)-mu(j,:);
        D (i,j) = sqrt(sum(d.^2,2));
    end
end
 
%% 将数据投影到二维平面，画出聚类结果图
if trigger == 1
    new_mu = mu;
    new_X = X;
    if length > 2 % 若为高维数据，则通过主成分分析投影到二维平面
        [~, score, ~, ~, ~, ~]=pca(X);
        new_X=score(:,1:2);   % 取转换后的矩阵score的前k列为PCA降维后特征
        new_mu = new_mu * (new_X\X)';
    end
    colours = ['r','g','b','m','k','y','c','w'];
    figure
    hold on
    grid off
    set(gca, 'linewidth', 1.1, 'fontsize', 14, 'fontname', '宋体')
    for i = 1 : k
        for j = 1 : n
            if indexs(j) == i
                plot(new_X(j,1),new_X(j,2),'.','Color',colours(i),'LineWidth',1)
            end
        end
        plot(new_mu(i,1),new_mu(i,2),'o','Color','#A500EB','markerfacecolor', colours(i),'LineWidth',2)
    end
    if length > 2
        xlabel("投影指标1")
        ylabel("投影指标2")
    else
        xlabel("指标1")
        ylabel("指标2")    
    end
end
end

结果如下图所示：