ARTICLE DETAIL

资讯详情

深耕网站视觉设计与运营推广的一线实战洞察。

MATLAB多维矩阵K-means聚类实战:从数据标准化到PCA可视化与结果验证

MATLAB多维矩阵K-means聚类实战:从数据标准化到PCA可视化与结果验证 简介本资源面向MATLAB环境下从事数据分析与机器学习的学习者聚焦K-means聚类在多维矩阵上的实现与可视化。包内共7个文件以5个m脚本和2个mat数据文件为主脚本涵盖质心初始化、样本分配、质心更新与主流程调用mat文件提供足球与空手道等测试数据集压缩包约62KB结构紧凑便于直接运行调试。已有2270人学习下载说明其在聚类入门与实战练习中具有一定参考价值。读者可借助完整代码掌握kmeans函数调用、类簇索引与质心矩阵的获取方式理解迭代收敛过程并学习用散点图、三维散点图及颜色编码展示聚类结果同时接触高维数据可能面临的维度灾难与降维思路适合课程实验、项目原型与自学练手。1. 多维矩阵聚类到底在聚什么从一摞 Excel 到能看懂的簇手里拿到一份 200 行 × 30 列的数据矩阵每一行是一个样本每一列是一个特征老板要你「把相似的样本归归类画个图看看」。这时候 K-means 聚类是最常见也最省事的起点。它做的事很朴素在 30 维空间里找 K 个中心点让每个样本归到离自己最近的那个中心反复调整中心位置直到稳定。多维矩阵聚类的难点不在算法本身而在于维度一高距离就变得不直观可视化也不能直接画 30 维散点。这篇讲的就是用 MATLAB 把这条链路走通读入多维矩阵、标准化、选 K、跑 K-means、把高维结果压到二维或三维画出来再判断聚类到底靠不靠谱。适合手头有数值矩阵、想做探索性分组又不想上深度学习的人。MATLAB 的统计与机器学习工具箱把 kmeans、pca、silhouette 这些函数都备齐了几十行就能跑通但参数设错、标准化漏做、K 选歪结果会直接翻车。下面按「先立住原理、再动手复现、最后避坑」的顺序拆开讲。2. K-means 在多维矩阵上的原理与选型理由2.1 目标函数与迭代逻辑为什么它只适合「球状簇」K-means 的优化目标是最小化簇内平方误差和写成公式就是每个样本到所属簇中心的欧氏距离平方之和。算法分两步循环分配步把每个样本分给最近的中心更新步把每个中心移到本簇样本的均值位置。这两步单调降低目标函数所以一定收敛但收敛到的可能是局部最优跟初始中心位置强相关。多维矩阵下这个逻辑有个隐含前提它默认簇是各向同性的球状、大小差不多。如果数据是长条形、环形或者密度差异极大K-means 会强行切出错误的边界。这也是为什么很多人跑完发现「图上看明明两团它却分错」——不是代码错是算法假设和数据形状不匹配。选型时先想清楚你的数据簇是不是近似球形、量纲是否统一这两条不满足就该考虑 DBSCAN 或层次聚类。2.2 距离度量与标准化量纲不统一等于白跑多维矩阵里各列量纲常常天差地别比如一列是 0 到 1 的比例另一列是上万的数量级。欧氏距离会被大数量级的列主导小量级的列几乎不起作用。所以跑 K-means 之前必须做标准化常见做法是 z-score让每列均值 0、标准差 1。% 假设原始数据矩阵 X 为 n×pn 个样本 p 个特征 X readmatrix(data.xlsx); % 读入多维矩阵 X rmmissing(X); % 去掉含缺失值的行K-means 不接受 NaN Xz zscore(X); % 按列标准化消除量纲影响readmatrix直接吃 Excel 和 csvrmmissing处理缺失zscore按列做标准化。注意zscore默认对每列独立处理这正是我们要的。如果某列方差为 0全是同一个值标准化会出 NaN需要提前删掉这种常量列。标准化后的矩阵才是喂给 kmeans 的输入别拿原始矩阵直接跑。2.3 K 值怎么定肘部法加轮廓系数双验证K 是 K-means 唯一需要人给的关键参数给错整个结果没意义。常用两招肘部法看簇内误差随 K 下降的拐点轮廓系数看每个样本跟自己簇的紧密度和跟最近邻簇的分离度。两个指标一起看比单看一个稳。maxK 10; wss zeros(maxK,1); % 簇内误差和 sil zeros(maxK,1); % 平均轮廓系数 for k 2:maxK [idx, C, sumd] kmeans(Xz, k, Replicates, 5); wss(k) sum(sumd); sil(k) mean(silhouette(Xz, idx)); end plot(2:maxK, wss(2:maxK), -o); xlabel(K); ylabel(簇内误差和); figure; plot(2:maxK, sil(2:maxK), -o); xlabel(K); ylabel(平均轮廓系数);Replicates设 5 表示用 5 组不同初始中心各跑一遍取最优能明显降低陷入局部最优的概率代价是耗时翻几倍。sumd是每个簇内样本到中心的距离平方和加总就是 WSS。肘部法找 WSS 曲线由陡变缓的拐点轮廓系数找峰值两者指向的 K 若一致就基本可信。轮廓系数越接近 1 越好低于 0.2 说明簇结构很弱这时候要回头怀疑数据本身有没有可分性。3. 从矩阵到聚类结果MATLAB 完整实现步骤3.1 数据准备与异常值处理真实矩阵里常有异常值一个离群点能把某个中心拽偏。跑之前先看分布用箱线图或马氏距离筛一遍。下面这段把读入、去缺失、去常量列、标准化串起来是每次跑聚类前的固定动作。X readmatrix(data.xlsx); X rmmissing(X); % 删除方差为 0 的常量列 keep std(X) 0; X X(:, keep); % 用马氏距离标记异常值阈值取卡方分位数 d mahal(X, mean(X)); outlier d chi2inv(0.975, size(X,2)); fprintf(检测到 %d 个异常样本\n, sum(outlier)); Xz zscore(X);mahal算马氏距离考虑了特征间相关性比单纯欧氏距离更适合多维。chi2inv(0.975, p)是 p 维卡方分布的 97.5% 分位点超过就判为异常。异常样本是删是留取决于业务如果是录入错误就删如果是真实极端情况就保留但心里有数。这一步不做后面聚类结果里冒出一个孤立小簇多半就是它在作怪。3.2 跑 K-means 并提取聚类标签与中心确定 K 之后正式跑一次把标签、中心、簇内距离都拿出来。中心是标准化空间里的坐标要还原回原始量纲才能解释业务含义。K 4; % 由肘部法和轮廓系数确定 rng(42); % 固定随机种子保证结果可复现 [idx, C, sumd, D] kmeans(Xz, K, Replicates, 10, Distance, sqeuclidean); % 把标准化空间的中心还原到原始量纲 mu mean(X); sg std(X); C_orig C .* sg mu; disp(各簇原始量纲中心); disp(C_orig); % 统计每簇样本数 tabulate(idx);rng(42)固定种子很关键否则每次跑初始中心不同结果会飘做汇报时说不清。Distance默认就是平方欧氏距离写出来更明确。C_orig C .* sg mu是标准化的逆变换因为 z-score 是 (x-mu)/sg反解就是乘 sg 加 mu。tabulate看每簇样本数如果某簇只有一两个样本基本是异常值或 K 给大了。3.3 用 PCA 把多维结果压到二维可视化30 维没法直接画最常见做法是用主成分分析降到 2 维或 3 维再按聚类标签上色。要注意 PCA 降维只是为了画图聚类本身是在原始标准化空间做的别搞反顺序。[coeff, score, latent] pca(Xz); explained cumsum(latent) / sum(latent) * 100; fprintf(前两个主成分累计解释方差%.1f%%\n, explained(2)); figure; gscatter(score(:,1), score(:,2), idx, lines(K), ., 12); xlabel(sprintf(PC1 (%.1f%%), explained(1))); ylabel(sprintf(PC2 (%.1f%%), explained(2))); title(K-means 聚类结果PCA 二维投影); % 叠加聚类中心投影 hold on; C_score (C - mean(Xz)) * coeff(:,1:2); plot(C_score(:,1), C_score(:,2), kx, MarkerSize, 14, LineWidth, 2); hold off;pca返回的score是样本在主成分上的坐标latent是各主成分方差explained算累计解释比例。如果前两个主成分加起来不到 50%说明二维投影丢了很多信息图上的重叠不代表原始空间里也重叠这点要跟看图的同事讲清楚。gscatter按idx自动分色。中心投影用(C - mean(Xz)) * coeff(:,1:2)因为 pca 内部对数据做了中心化投影中心时也要减均值。3.4 三维投影与轮廓图辅助判断二维不够时可以画三维MATLAB 的 scatter3 直接支持。再配一张轮廓系数图能看出哪些样本归类不踏实。figure; scatter3(score(:,1), score(:,2), score(:,3), 30, idx, filled); xlabel(PC1); ylabel(PC2); zlabel(PC3); title(三维 PCA 投影下的聚类); colorbar; figure; silhouette(Xz, idx); title(轮廓系数分布);scatter3第四个参数是点大小第五个是颜色映射filled填充实心。轮廓图里每簇一条横条条越长越整齐出现负值说明那些样本可能被分错了簇。三维图旋转一下能从不同角度看簇的分离情况比二维更直观但投影到三维仍可能失真判断聚类质量还是以轮廓系数和业务解释为主。4. 聚类结果怎么验证与调参别只看图好看4.1 内部指标轮廓系数、Calinski-Harabasz 与 Davies-Bouldin图好看不等于聚类对。内部指标不需要真实标签靠数据本身的紧密度和分离度打分。轮廓系数前面用过Calinski-Harabasz 越大越好Davies-Bouldin 越小越好。三个一起看能交叉验证。eva evalclusters(Xz, kmeans, CalinskiHarabasz, KList, 2:10); fprintf(CH 指标最优 K %d\n, eva.OptimalK); eva2 evalclusters(Xz, kmeans, DaviesBouldin, KList, 2:10); fprintf(DB 指标最优 K %d\n, eva2.OptimalK);evalclusters是 MATLAB 自带的聚类评估函数直接传数据、算法名和指标名它会自动遍历 K 列表给出最优值。CH 指标对簇数多的情形有偏好DB 指标对球形簇敏感两者结论不一致时以轮廓系数和业务可解释性为准。别迷信单一指标它们各有偏向。4.2 稳定性检验换种子、换子集看标签是否一致K-means 对初始值敏感一个靠谱的聚类应该在多次随机初始化下结果稳定。做法是跑多次用调整兰德指数比较标签一致性。nRun 20; labels zeros(size(Xz,1), nRun); for r 1:nRun labels(:,r) kmeans(Xz, K, Replicates, 1); end ari zeros(nRun, nRun); for i 1:nRun for j 1:nRun ari(i,j) rand_index(labels(:,i), labels(:,j)); end end fprintf(平均成对兰德指数%.3f\n, mean(ari(triu(true(nRun),1))));rand_index需要自己实现或用第三方函数核心是统计两次聚类中「同簇且同簇」「异簇且异簇」的比例。平均 ARI 高于 0.8 说明结果很稳低于 0.5 说明聚类结构本身模糊换种子就变这种结果拿去指导决策风险很大。这一步很多人省掉但恰恰是判断「能不能信」的关键。4.3 业务解释把簇中心翻译成人话统计指标过关后最后一步是把每个簇的中心还原到原始量纲看它在各特征上的高低给簇起个业务名字。比如客户分群里某簇「消费频次高、客单价低」就能叫「高频低客单群体」。这一步没有代码能替你完成得结合领域知识逐列对比。for k 1:K fprintf(簇 %d%d 个样本特征均值\n, k, sum(idxk)); disp(array2table(C_orig(k,:), VariableNames, featureNames)); endfeatureNames是你原始矩阵的列名提前存好。逐簇打印中心跟全局均值比高于全局的标出来重点看。簇中心差异越明显聚类越有业务价值如果几个簇中心几乎一样说明 K 给多了或者数据本来就没结构。5. 避坑与排查多维 K-means 最容易翻车的五个地方5.1 没标准化导致某几列主导距离现象聚类结果几乎只按某一列分组其他列看不出影响。原因那列数值量级远大于其他列欧氏距离被它垄断。解决跑之前统一 z-score 标准化跑完把中心还原回原始量纲解释。这是最高频的翻车点血泪经验是拿到数据先std看一眼各列量级。5.2 K 值拍脑袋定结果无法解释现象K 设成 5 跑出来某簇只有 3 个样本或者两个簇中心几乎重合。原因K 没经过肘部法和轮廓系数验证纯凭感觉。解决用evalclusters或手写循环遍历 K2 到 10结合业务可解释性定 K。簇太小往往是异常值先回去做异常检测。5.3 缺失值没处理直接报错或结果失真现象kmeans报 NaN 错误或者某些样本距离算出来是 NaN 被随机分配。原因矩阵里有缺失值。解决rmmissing删行或用均值、中位数填补。删行会损失样本填补会引入偏差样本量大就删量小就填并记录填补比例。5.4 PCA 投影图重叠就以为聚类失败现象二维图上几簇糊在一起以为算法不行。原因前两个主成分解释方差太低二维投影丢信息。解决看explained累计比例低于 60% 就换三维或改用 t-SNE 投影同时以轮廓系数而非肉眼判断聚类质量。5.5 每次跑结果不一样汇报时说不清现象同一份数据两次运行标签完全不同。原因初始中心随机没固定种子也没开Replicates。解决rng固定种子Replicates设 10 以上取最优再用兰德指数验证稳定性。可复现是工程底线别让结果变成玄学。6. 进阶技巧用轮廓系数热力图定位「骑墙」样本跑通基础流程后真正拉开差距的是对边界样本的处理。轮廓系数只给一个平均值会掩盖问题把每个样本的轮廓值画成热力图能一眼看出哪些样本归类不踏实。这些「骑墙」样本往往是业务上最值得单独关注的群体比如客户分群里行为跨两类的过渡人群。s silhouette(Xz, idx); % 每个样本的轮廓值 [~, order] sort(idx); % 按簇排序让同簇样本聚在一起 figure; imagesc(s(order)); colorbar; caxis([-1 1]); xlabel(样本按簇排序); ylabel(样本序号); title(样本轮廓系数热力图); % 找出轮廓值低于 0 的样本 bad find(s 0); fprintf(轮廓值为负的样本数%d\n, numel(bad));silhouette不带输出参数时直接画图带输出时返回每个样本的值。imagesc把一维轮廓值铺成色带负值区域偏蓝就是可能分错的样本。caxis([-1 1])固定色标范围方便多次对比。找出负值样本后回到原始矩阵看它们的特征往往能发现是异常值、录入错误或者确实处于两类之间的过渡态。我自己的习惯是任何一次聚类汇报前先跑这张热力图负值样本超过 10% 就不急着下结论先回头查数据质量。聚类从来不是跑完就完事验证和解释占的功夫比调参多得多。把标准化、K 值验证、稳定性检验、边界样本排查这四步做成固定流程多维矩阵聚类才真正能拿来支撑决策而不是画一张好看的图交差。希望帮到你。本文还有配套的精品资源点击获取
返回列表