% 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
结果如下图所示: