尧图网页设计 
简介面向数据分析和机器学习初学者的Matlab聚类分析源代码资源围绕K-Means、层次聚类、DBSCAN、谱聚类与Fuzzy C-Means等常用算法提供可直接运行的m脚本和配套示例数据适合教学演示、算法对比与二次开发。压缩包共9个文件含5个m源代码和4个xls数据表整体约18KB体量紧凑但覆盖算法实现、有效性指标、数据预处理与结果可视化等关键环节。已有1891人学习下载是经典教材配套例程的典型整理。通过运行这些代码可观察K-Means质心迭代、层次聚类树生成、DBSCAN密度簇划分、谱分解与模糊聚类隶属度计算等完整过程还能基于Silhouette等指标评估聚类效果并修改样本或参数快速拓展实验。1. 拿到 MATLAB 聚类分析源代码先搞懂这份代码到底在算哪一步聚类不是单个函数很多人从网上下载一份 matlab 聚类分析源代码文件名写着 kmeans.m 或者 cluster_analysis.m拿来直接 load 数据就跑结果往往跟 README 里贴的图对不上聚类数不对、散点图里颜色乱跳、有时还直接报错说矩阵维度不对。这不是源代码写得差而是聚类分析在 MATLAB 里从来不是一个函数的事。它是一串由数据矩阵、距离度量、聚类数 k、评价指标组成的流程源代码只是这条流程的一个快照。我一般会把这类标题下的代码分成两类一类是教学演示型重点在算法本身另一类是项目落地型包含数据预处理和结果导出。你要做的第一件事是分清手上这份属于哪类然后重新核对数据矩阵的方向和参数范围。这篇笔记就按这个思路从读代码一路讲到改代码、避坑和验证结果。2. 聚类源代码的四个关键件数据矩阵、距离度量、聚类数 k 与评价指标2.1 数据矩阵怎么喂行是样本、列是特征别把维度搞反聚类分析源代码里几乎第一个可执行语句就是[idx, C] kmeans(X, k)这里的 X 必须是 n 行 p 列每一行是一个观测样本每一列是一个特征。很多人拿着 Excel 表直接复制粘贴把特征放到了行上跑出来的距离矩阵是 p 乘 p 而不是 n 乘 n结果要么内存爆掉要么每个样本变成一个特征而失去聚类意义。检查维度最笨也最可靠的方法是 size 和 isa% 检查数据矩阵方向行样本列特征 X readmatrix(your_data.csv); % 读取数据 fprintf(样本数 %d, 特征数 %d\n, size(X,1), size(X,2)); if size(X,1) size(X,2) warning(行数小于列数确认是否应该转置X X;); end这段代码的逻辑很简单先用 size 拿到形状如果样本数比特征数还少十有八九是转置了。注意 readmatrix 是 R2019a 之后才有的函数老版本可以用 csvread 或 importdata。另一个容易忽略的问题是数据里有缺失值或文本列kmeans 遇到 NaN 会直接报错所以读进来以后还要做一句if any(ismissing(X)), error(数据含缺失值请先填充或删除); end。矩阵方向确定后再想算法的事。2.2 距离度量与算法选型欧氏、相关系数与马氏距离聚类分析源代码的第二个关键件是距离度量。K-Means 默认用欧氏距离因为它每轮迭代都在做点到质心的更新层次聚类可以支持更多度量而 DBSCAN 则依赖密度对距离的定义更敏感。选度量的原则是如果特征是连续数值且量纲一致直接用欧氏如果特征是量纲不同的混合类型先归一化再欧氏如果特征之间存在明显共线性可以试试马氏距离。在 MATLAB 里距离矩阵可以用 pdist 或 pdist2 生成% 生成距离矩阵的两个函数pdist 与 pdist2 D_vec pdist(X, euclidean); % 返回行向量长度 n*(n-1)/2 D_mat squareform(D_vec); % 转成 n x n 对称矩阵 % 或者直接用 pdist2 得到矩阵 D_mat2 pdist2(X, X, squaredeuclidean);这里 pdist 返回的是压缩距离向量省内存如果你要交给 linkage 用直接用 pdist 的输出即可。squaredeuclidean是平方欧氏距离但很多源码里直接用 pdist2 算矩阵再喂给自定义聚类算法这时要注意它和后端算法期望的度量是不是一致。相关系数距离也有场景比如基因表达谱聚类MATLAB 里可以写成pdist2(X, X, correlation)它把每个样本当作一个向量计算 1 减去皮尔逊相关系数。我见过不少源代码在层次聚类里默认用cosine结果对向量的模长敏感的项目就会翻车所以选度量前先看样本向量的长度是否有意义。如果样本量不大但特征之间相关性很高马氏距离是更合理的选择。MATLAB 里没有直接的马氏距离 pdist 选项但可以用mahal(X, X)或者先做 PCA 白化再算欧氏距离。实际项目中我很少直接上马氏距离因为协方差矩阵估计不稳定样本量少于特征数时尤其危险。对比一下欧氏距离简单、快、对异常值敏感相关系数距离适合形状相似但幅值不同的样本马氏距离能消除特征间相关性但代价是计算复杂度高。源代码里的距离度量一般写在算法开头你改一行就能切换。2.3 聚类数 k 怎么定肘部法、轮廓系数与 gap statistic源代码里最常见的硬编码是k 3。这个值如果拍脑袋定后面所有可视化都是自欺欺人。我常用的定 k 流程是先跑肘部法缩小范围再对候选 k 算轮廓系数最后用 gap statistic 交叉验证。肘部法看的是每个 k 对应的组内离差平方和WSS下降速度% 肘部法对 k1:8 计算总组内平方和 rng(1); % 固定随机种子保证可复现 X zscore(X); % 先标准化 wss zeros(8,1); for k 1:8 [~, ~, sumd] kmeans(X, k, Replicates, 5, MaxIter, 300); wss(k) sum(sumd); % sumd 是每个样本到所属质心的距离平方和 end plot(1:8, wss, o-); xlabel(k); ylabel(WSS);这段代码里 kmeans 的第三个输出 sumd 是 k 乘 1 的向量每个元素是该类所有样本到质心的距离平方和所有类加起来就是 WSS。Replicates设为 5是因为 K-Means 每次初始中心不同多次重复取最优可避免局部最优。肘部法的判读是找曲线拐点但真实数据往往没有清晰拐点所以要用轮廓系数再确认% 对候选 k 计算平均轮廓系数 avg_s zeros(1, 5); for k 2:6 idx kmeans(X, k, Replicates, 10, MaxIter, 500); s silhouette(X, idx); avg_s(k-1) mean(s); % 保存平均轮廓系数 end [best, loc] max(avg_s); fprintf(平均轮廓系数最高的 k %d, 值为 %.3f\n, loc1, best);轮廓系数取值范围 -1 到 1越接近 1 表示样本离自己的簇而远离其他簇。注意 silhouette 函数返回的是每个样本的轮廓值平均之后才能比较 k。gap statistic 在 MATLAB 里没有直接函数但evalclusters(X, (X,k) kmeans(X,k), gap)可以一行替代它比较的是真实数据与随机均匀数据的 WSS 差距计算结果最稳但耗时也最久。一般小数据集少于 5000 条特征 50 以内可以直接用。2.4 评价指标内部指标与外部指标源代码里至少要有一个聚类分析源代码如果只输出idx不输出评价指标那这份代码只能算完成了一半。评价指标分两类内部指标只看数据和聚类结果比如轮廓系数、Davies-Bouldin 指数、Calinski-Harabasz 指数外部指标需要真实标签比如调整兰德指数ARI、归一化互信息NMI。在 MATLAB 里内部指标用evalclusters最省事% 用 evalclusters 同时计算多种内部指标 eva evalclusters(X, (X,k) kmeans(X,k,Replicates,5), CalinskiHarabasz, KList, 2:6); plot(eva); % 画出 CH 指数随 k 的变化 bestK eva.OptimalK;这里第二个参数是一个函数句柄evalclusters 会在每个候选 k 上调用它。Calinski-Harabasz 指数越大表明聚类越紧致、类间越分散它比轮廓系数计算快很多适合先粗筛。如果你有真实标签 y计算 ARI 可以用文件交换上的randindexMATLAB 自带的confusionmat只能做标签对齐后的精确度不能直接当 ARI 用。外部指标还有一个容易踩的坑聚类返回的类号是 1、2、3真实标签可能是 0、1、2先把类号统一成连续正整数再算指标否则结果会低到离谱。我的建议是源代码里至少保留一个内部指标这样没有真实标签时也能向别人证明这个聚类不是随便分的。3. 用 MATLAB 跑通三种最常用聚类源代码K-Means、层次聚类与 DBSCAN 的最小可复现命令3.1 K-Means 最小可执行代码从 load 到 silhouetteK-Means 是绝大多数聚类分析源代码的默认算法因为它快、参数少、结果好解释。下面这一段是可以在 MATLAB 里直接复制运行的最小案例用的是自带数据集 fisheriris省去造数环节% K-Means 最小可复现案例 load fisheriris; % 加载内置鸢尾花数据 X meas; % 150 x 4 的特征矩阵 rng(1); % 固定随机种子保证结果可复现 [idx, C, sumd] kmeans(X, 3, Replicates, 10, MaxIter, 500); figure; gscatter(X(:,1), X(:,2), idx); % 用前两列画散点按聚类着色 hold on; plot(C(:,1), C(:,2), kx, MarkerSize, 15, LineWidth, 2); % 标出质心 title(K-Means 聚类结果 (k3));这里Replicates是重复次数设为 10 表示从 10 组不同初始中心里挑 WSS 最小的一次MaxIter是单次迭代上限500 次在 150 个样本上完全够用。C 是最终的 k 乘 p 质心矩阵sumd 是每个类内距离平方和后续算总 WSS 或做肘部图都要用到它。跑完这段代码你应该看到三类样本大致把山鸢尾和另外两种分开但 versicolor 和 virginica 会有部分重叠这说明单纯用 K-Means 不加核技巧线性边界很难完美区分重叠类别。如果样本量很大几十万行把Replicates降为 1MaxIter按收敛情况调。K-Means 的复杂度接近 O(nkp)n 是样本数k 是类别数p 是特征数。在 MATLAB 里它默认使用 kmeans 初始化比纯随机初始化稳定得多但依然要固定 rng否则你每次跑出的类标签顺序会漂移。3.2 层次聚类源代码怎么改linkage dendrogram cluster层次聚类适合你不知道聚类数、想看层级关系的场景源代码里一般就是 linkage、dendrogram、cluster 三连。它不要求预置 k而是先生成谱系图再按距离阈值或类别数截断。下面这段是常见写法% 层次聚类最小案例Ward 连接 谱系图 按类别数截断 load fisheriris; X meas; Z linkage(X, ward, euclidean); % 生成聚类树 figure; dendrogram(Z, 0); % 画出完整谱系图0 表示不限制叶子数 idx cluster(Z, MaxClust, 3); % 按最多 3 类截断linkage 的第一个输入可以是原始数据矩阵也可以是 pdist 输出的距离向量。ward方法要求距离度量是欧氏距离所以第三个参数写euclidean如果你用average或complete可以搭配pdist2得到任意距离矩阵再喂给 linkage但注意 linkage 第二个参数指定的是连接准则第三个参数才是距离度量。dendrogram 的 0 表示把所有叶子都画出来如果样本超过 300 个画出来会密密麻麻这时建议只画前 30 个叶子dendrogram(Z, 30)并按相干系数排序。层次聚类最容易被忽略的是cluster的截断方式。MaxClust是告诉 MATLAB 找到一种截断使类别数不超过指定值Cutoff是根据不一致系数截断对类间距离阈值更直观。如果你发现截断出来的类数比期望少先检查 linkage 是否用的是wardWard 法倾向于生成球形簇对异常值敏感换成average往往能缓解。谱系图的纵轴是连接距离距离越大表示两类越不像你可以在图上画一条横线看应该切哪里这比硬调 MaxClust 更符合直觉。3.3 DBSCAN 源代码与参数敏感性epsilon 和 minpts 的选法DBSCAN 在 MATLAB 里不是老版本标配R2019a 之后才有内置dbscan函数。如果你手上的源代码是早年下载的可能是一段自实现版本但用法大同小异。它的核心是两个参数邻域半径 epsilon 和最小邻域点数 minpts。下面用合成数据演示两个密度不同的簇加上一圈噪声正好暴露 DBSCAN 的优势% DBSCAN 最小案例合成双簇 噪声 rng(2); X [randn(100, 2) * 0.4; randn(80, 2) 2; rand(30, 2) * 3 - 1]; % 簇1、簇2、噪声 [idx, corepts] dbscan(X, 0.5, 10); % epsilon0.5, minpts10 gscatter(X(:,1), X(:,2), idx); title(DBSCAN 聚类结果噪声标为0);运行后你会看到 idx 中大部分样本被分为 1 和 2 两类但靠左下角那摊噪声会被标记为 0。dbscan 将 0 视为噪声点不归属任何簇。corepts是逻辑索引标出每个点是核心点还是边界点。epsilon 的取值是最大的坑设太小一个簇会被拆成好几块设太大所有点都变成一类噪声也吞进去。我常用的定参方法是 k-距离图% 用 k-距离图选 epsilon找曲线拐点 k 10; % 一般取 minpts [~, distK] knnsearch(X, X, K, k); % 每个样本到第 k 近邻的距离 distK sort(distK(:, end)); % 按升序排序 plot(distK); ylabel(第 k 近邻距离); xlabel(样本排序);这里的思路是对每个样本计算它到第 k 个近邻的距离按距离从小到大排序后画曲线。曲线会出现一个膝盖位置对应密度分布的分界那个位置的距离就是合适的 epsilon。minpts 的经验值是取特征维数的两倍以上二维数据取 10 可以高维数据要更大。DBSCAN 不能像 K-Means 那样直接指定聚类数如果你非要固定簇数得循环试参这也是很多聚类分析源代码里 DBSCAN 写得最潦草的原因。4. 聚类分析源代码避坑指南五个让结果翻车的隐藏问题4.1 现象聚类结果每次跑都不一样今天和明天对不上原因K-Means 和 DBSCAN 的初始状态依赖随机数。K-Means 在 MATLAB 里虽然用 kmeans 初始化但随机数种子不固定时每次调用得到的质心初始位置不同虽然大多数情况下会收敛到同一组簇但在数据分布重叠较大时可能收敛到不同的局部最优类标签的数字顺序也会变。DBSCAN 本身是确定性的但如果你用了随机抽样或打乱数据顺序结果也会漂移。解决在源代码开头写一行rng(2024)或者在调用 kmeans 时设置Options statset(UseSubstreams, true)。更稳妥的做法是把 rng 种子的值保存在结果结构体里每次跑记录 seed这样别人复现时能拿到完全一致的输出。我在自己的工程脚本里一律先rng(default)再调算法避免上次能跑出好结果这次死活复现不了的窘境。4.2 现象聚类结果里异常点被单独分成一类甚至把整个簇拽偏原因没有做特征归一化或离群点处理。K-Means 基于欧氏距离如果某个特征的量纲是 1000 而其他是 0.1距离几乎完全由那个大数特征决定离群点则会把质心拉向自己导致几个正常点被划成一类。解决先对特征做 zscore 标准化或 min-max 缩放再用 medfilt 或马氏距离剔除离群点。zscore 在 MATLAB 里一句话X zscore(X);这个操作应该放在聚类之前而不是之后。如果数据里存在极端长尾用robustcov配合马氏距离更稳。注意 zscore 的均值和标准差要保存下来等新数据进来时用同一组参数变换别直接重新算全局的 zscore。4.3 现象层次聚类画出的谱系图与 cluster 截断结果对不上原因dendrogram 的叶子顺序与 cluster 输出的标签顺序没有直接映射关系尤其是样本量超过 100 时谱系图会压缩显示。另外用Cutoff截断时MATLAB 默认用不一致系数而不是直接按距离阈值很多人没读文档就乱填数字。解决不要靠眼睛在图上数类数直接用cophenet(Z, D)检查聚类树与原始距离的相干系数数值越接近 1 表示谱系图保留的距离信息越好。要按距离截断明确写Cutoff, 1.5, Depth, 3其中 Depth 指定计算不一致系数时考虑的层级深度。最省心的方法是像我一样用cluster(Z, MaxClust, k)只负责定类数不碰阈值。4.4 现象DBSCAN 把所有点都归为一类或者全部标成噪声原因epsilon 选得太大所有点的邻域都彼此连通自然只剩一个簇epsilon 选得太小每个点都找不到足够邻居全成了噪声点。还有一种隐蔽情况输入数据没有标准化二维数据横纵轴量纲不同圆形邻域变成椭圆密度判断完全失真。解决先用 k-距离图定 epsilon再对每个候选 minpts 画一下聚类结果看稳定性。我这里有一个通用参数起点dbscan(X, 0.5, 10)对二维标准化数据通常可用但业务数据不要直接抄。更实际的办法是写一个 3 行循环把 epsilon 按几何间隔扫 10 个值记录每次的非噪声点数选一个非噪声点数突然上升后又稳定的区域。4.5 现象聚类结果与论文或同事的源代码对不上怀疑代码被改坏原因数据版本不同、特征顺序不同、随机种子不同这三样是高频元凶。聚类的类标签本身没有绝对顺序K-Means 里第 1 类可能是别人代码里的第 3 类直接对比标签数字完全没意义。解决不要比较 idx 的数字要比较聚类结果与已知标签的互信息或 ARI。先把两份代码的类标签用confusionmat做映射再看对角线占比。排查时依次检查读入文件的列顺序、是否做了同样的标准化、rng 种子是否一致。我自己的做法是在源代码里加一行注释记录数据文件的 MD5 校验和从源头避免数据不对的扯皮。5. 把源代码封装成自己的聚类工具函数自动搜 k、自动出图、落盘结果5.1 一个函数收编三种算法参数接口这样设计下载的源代码往往面向单次运行参数写死在脚本里。要应对换数据、换算法、换 k的日常最好封装成一个函数。下面是我常用的模板支持三种算法统一返回 idx 和最优 k并且自动固定随机种子function result clusterPipeline(X, varargin) %CLUSTERPIPELINE 聚类分析一站式封装 % result clusterPipeline(X, Method, kmeans, KRange, 2:8) % 返回值 result 包含 idx/bestK/silhouette 等字段 % 参数解析 p inputParser; addParameter(p, Method, kmeans, (x) ismember(x, {kmeans,hierarchical,dbscan})); addParameter(p, KRange, 2:8, (x) isnumeric(x) isvector(x)); addParameter(p, Norm, true, islogical); addParameter(p, Seed, 1, isscalar); parse(p, varargin{:}); opts p.Results; % 固定随机种子 归一化 rng(opts.Seed); if opts.Norm X zscore(X); end bestS -inf; bestK 2; switch opts.Method case kmeans for k opts.KRange idx kmeans(X, k, Replicates, 10, MaxIter, 500); s mean(silhouette(X, idx)); if s bestS bestS s; bestK k; bestIdx idx; end end case hierarchical D pdist(X, euclidean); Z linkage(D, ward); figure; dendrogram(Z, 0); title(层次聚类谱系图); for k opts.KRange idx cluster(Z, MaxClust, k); if numel(unique(idx)) 2, continue; end s mean(silhouette(X, idx)); if s bestS bestS s; bestK k; bestIdx idx; end end case dbscan % DBSCAN 不适合自动搜 k这里用固定参数演示 [bestIdx, core] dbscan(X, 0.5, 10); bestK numel(unique(bestIdx(bestIdx 0))); bestS mean(silhouette(X, bestIdx(bestIdx 0))); % 仅对非噪声计算 end % 打包结果 result struct(); result.idx bestIdx; result.bestK bestK; result.silhouette bestS; result.Method opts.Method; result.Seed opts.Seed; end参数设计上Method决定走哪条分支KRange只在 kmeans 和 hierarchical 分支下生效Norm控制是否标准化Seed固定随机性。注意 dbscan 分支里我把silhouette只用在非噪声样本上因为噪声点的类号为 0轮廓系数对 0 类会报错。这个函数的好处是你换数据时只需要改 X换算法时改一行 Method不用每处脚本都翻找参数。5.2 输出结构体与落盘结果、图、日志一次到位这段函数返回的 result 结构体在 MATLAB 里可以直接被后续脚本使用。我习惯在调用后追加一段落盘逻辑% 调用示例K-Means 聚类搜索 k2:6并保存全部结果 res clusterPipeline(X, Method, kmeans, KRange, 2:6, Seed, 42); save(cluster_result.mat, res); writematrix(res.idx, cluster_label.csv); % 把每行的类别标签导出writematrix 把 idx 写成 CSV适合后续做业务分析save 保存的 mat 文件则留着复现。如果你想连聚类中心一起存在函数返回里再加一个 C 字段。日志方面我一般用fprintf输出一行摘要包含 seed、最优 k 和轮廓系数这样跑批量实验时能在命令窗口直接看到对比不必打开 mat 文件。5.3 调用示例与边界数据量超过几万行怎么办上面的函数适合几千行、几十维以内的数据。如果样本量到十万行有几处必须改silhouette在十万样本上要算 n 乘 n 的距离矩阵内存直接爆炸。这时不要用轮廓系数选 k改用 Calinski-Harabasz 或干脆只跑肘部法。DBSCAN 的 knnsearch 在十万样本上也慢建议先抽样估计 epsilon再在全量上跑。我实际处理过二十万行的用户分群做法是先随机抽 20% 样本定参数再用那个参数跑全量时间从一小时降到五分钟。聚类分析源代码从教学到工程差别就在这些规模边界上。6. 验证聚类结果不是玄学一个置换检验的 MATLAB 实现费了大力气跑出聚类结果下一步最怕的是被别人问一句你怎么知道这个聚类不是把随机数据硬分出来的这个问题我吃过亏以前直接拿轮廓系数当证据结果审稿人让我证明簇的显著性。我后来用置换检验解决对每个特征独立打乱顺序让原始数据的结构被破坏再重新聚类看轮廓系数下降多少。如果真实数据的轮廓系数显著高于随机打乱后的分布说明聚类结果不是偶然。% 置换检验打乱特征验证聚类显著性 rng(1); res_obs clusterPipeline(X, Method, kmeans, KRange, 2:6, Seed, 1); obs_s res_obs.silhouette; null_s zeros(100, 1); for i 1:100 X_perm X(:, randperm(size(X,2))); % 按列打乱特征 % 注意randperm 只打乱列的顺序不破坏行内关系属于轻量置换 res_null clusterPipeline(X_perm, Method, kmeans, KRange, 2:6, Seed, i); null_s(i) res_null.silhouette; end p_value sum(null_s obs_s) / 100; fprintf(观测轮廓系数 %.3f, 置换检验 p %.3f\n, obs_s, p_value);这段代码跑了 100 次置换每次用不同的种子。p 值小于 0.05 就说明聚类结构显著强于随机打乱后的结果。更严格的做法是对每一列做随机重排而不是整列顺序不变但那样会破坏列间相关性检验力度过强实际项目中我常用列置换作为快速版本。这个技巧对 K-Means 和层次聚类都适用DBSCAN 也可以把轮廓系数换成非噪声点的平均数。跑完不要只贴 p 值把 null_s 的直方图画出来和 obs_s 的位置一起展示这是最有说服力的证据。我现在的习惯是任何聚类分析源代码交付前至少要跑一次这个置换检验然后把它写进注释里。固定种子、保留原始数据、记录参数这三件事看似不起眼却能让你的结果经得起反复问怎么复现。聚类分析在 MATLAB 里不难难的是让结果站得住脚。希望这篇笔记能帮你把源代码跑通也帮你把结果讲清楚。本文还有配套的精品资源点击获取 MATLAB聚类分析源代码从跑通到验证：K-Means、层次聚类与DBSCAN
简介面向数据分析和机器学习初学者的Matlab聚类分析源代码资源围绕K-Means、层次聚类、DBSCAN、谱聚类与Fuzzy C-Means等常用算法提供可直接运行的m脚本和配套示例数据适合教学演示、算法对比与二次开发。压缩包共9个文件含5个m源代码和4个xls数据表整体约18KB体量紧凑但覆盖算法实现、有效性指标、数据预处理与结果可视化等关键环节。已有1891人学习下载是经典教材配套例程的典型整理。通过运行这些代码可观察K-Means质心迭代、层次聚类树生成、DBSCAN密度簇划分、谱分解与模糊聚类隶属度计算等完整过程还能基于Silhouette等指标评估聚类效果并修改样本或参数快速拓展实验。1. 拿到 MATLAB 聚类分析源代码先搞懂这份代码到底在算哪一步聚类不是单个函数很多人从网上下载一份 matlab 聚类分析源代码文件名写着 kmeans.m 或者 cluster_analysis.m拿来直接 load 数据就跑结果往往跟 README 里贴的图对不上聚类数不对、散点图里颜色乱跳、有时还直接报错说矩阵维度不对。这不是源代码写得差而是聚类分析在 MATLAB 里从来不是一个函数的事。它是一串由数据矩阵、距离度量、聚类数 k、评价指标组成的流程源代码只是这条流程的一个快照。我一般会把这类标题下的代码分成两类一类是教学演示型重点在算法本身另一类是项目落地型包含数据预处理和结果导出。你要做的第一件事是分清手上这份属于哪类然后重新核对数据矩阵的方向和参数范围。这篇笔记就按这个思路从读代码一路讲到改代码、避坑和验证结果。2. 聚类源代码的四个关键件数据矩阵、距离度量、聚类数 k 与评价指标2.1 数据矩阵怎么喂行是样本、列是特征别把维度搞反聚类分析源代码里几乎第一个可执行语句就是[idx, C] kmeans(X, k)这里的 X 必须是 n 行 p 列每一行是一个观测样本每一列是一个特征。很多人拿着 Excel 表直接复制粘贴把特征放到了行上跑出来的距离矩阵是 p 乘 p 而不是 n 乘 n结果要么内存爆掉要么每个样本变成一个特征而失去聚类意义。检查维度最笨也最可靠的方法是 size 和 isa% 检查数据矩阵方向行样本列特征 X readmatrix(your_data.csv); % 读取数据 fprintf(样本数 %d, 特征数 %d\n, size(X,1), size(X,2)); if size(X,1) size(X,2) warning(行数小于列数确认是否应该转置X X;); end这段代码的逻辑很简单先用 size 拿到形状如果样本数比特征数还少十有八九是转置了。注意 readmatrix 是 R2019a 之后才有的函数老版本可以用 csvread 或 importdata。另一个容易忽略的问题是数据里有缺失值或文本列kmeans 遇到 NaN 会直接报错所以读进来以后还要做一句if any(ismissing(X)), error(数据含缺失值请先填充或删除); end。矩阵方向确定后再想算法的事。2.2 距离度量与算法选型欧氏、相关系数与马氏距离聚类分析源代码的第二个关键件是距离度量。K-Means 默认用欧氏距离因为它每轮迭代都在做点到质心的更新层次聚类可以支持更多度量而 DBSCAN 则依赖密度对距离的定义更敏感。选度量的原则是如果特征是连续数值且量纲一致直接用欧氏如果特征是量纲不同的混合类型先归一化再欧氏如果特征之间存在明显共线性可以试试马氏距离。在 MATLAB 里距离矩阵可以用 pdist 或 pdist2 生成% 生成距离矩阵的两个函数pdist 与 pdist2 D_vec pdist(X, euclidean); % 返回行向量长度 n*(n-1)/2 D_mat squareform(D_vec); % 转成 n x n 对称矩阵 % 或者直接用 pdist2 得到矩阵 D_mat2 pdist2(X, X, squaredeuclidean);这里 pdist 返回的是压缩距离向量省内存如果你要交给 linkage 用直接用 pdist 的输出即可。squaredeuclidean是平方欧氏距离但很多源码里直接用 pdist2 算矩阵再喂给自定义聚类算法这时要注意它和后端算法期望的度量是不是一致。相关系数距离也有场景比如基因表达谱聚类MATLAB 里可以写成pdist2(X, X, correlation)它把每个样本当作一个向量计算 1 减去皮尔逊相关系数。我见过不少源代码在层次聚类里默认用cosine结果对向量的模长敏感的项目就会翻车所以选度量前先看样本向量的长度是否有意义。如果样本量不大但特征之间相关性很高马氏距离是更合理的选择。MATLAB 里没有直接的马氏距离 pdist 选项但可以用mahal(X, X)或者先做 PCA 白化再算欧氏距离。实际项目中我很少直接上马氏距离因为协方差矩阵估计不稳定样本量少于特征数时尤其危险。对比一下欧氏距离简单、快、对异常值敏感相关系数距离适合形状相似但幅值不同的样本马氏距离能消除特征间相关性但代价是计算复杂度高。源代码里的距离度量一般写在算法开头你改一行就能切换。2.3 聚类数 k 怎么定肘部法、轮廓系数与 gap statistic源代码里最常见的硬编码是k 3。这个值如果拍脑袋定后面所有可视化都是自欺欺人。我常用的定 k 流程是先跑肘部法缩小范围再对候选 k 算轮廓系数最后用 gap statistic 交叉验证。肘部法看的是每个 k 对应的组内离差平方和WSS下降速度% 肘部法对 k1:8 计算总组内平方和 rng(1); % 固定随机种子保证可复现 X zscore(X); % 先标准化 wss zeros(8,1); for k 1:8 [~, ~, sumd] kmeans(X, k, Replicates, 5, MaxIter, 300); wss(k) sum(sumd); % sumd 是每个样本到所属质心的距离平方和 end plot(1:8, wss, o-); xlabel(k); ylabel(WSS);这段代码里 kmeans 的第三个输出 sumd 是 k 乘 1 的向量每个元素是该类所有样本到质心的距离平方和所有类加起来就是 WSS。Replicates设为 5是因为 K-Means 每次初始中心不同多次重复取最优可避免局部最优。肘部法的判读是找曲线拐点但真实数据往往没有清晰拐点所以要用轮廓系数再确认% 对候选 k 计算平均轮廓系数 avg_s zeros(1, 5); for k 2:6 idx kmeans(X, k, Replicates, 10, MaxIter, 500); s silhouette(X, idx); avg_s(k-1) mean(s); % 保存平均轮廓系数 end [best, loc] max(avg_s); fprintf(平均轮廓系数最高的 k %d, 值为 %.3f\n, loc1, best);轮廓系数取值范围 -1 到 1越接近 1 表示样本离自己的簇而远离其他簇。注意 silhouette 函数返回的是每个样本的轮廓值平均之后才能比较 k。gap statistic 在 MATLAB 里没有直接函数但evalclusters(X, (X,k) kmeans(X,k), gap)可以一行替代它比较的是真实数据与随机均匀数据的 WSS 差距计算结果最稳但耗时也最久。一般小数据集少于 5000 条特征 50 以内可以直接用。2.4 评价指标内部指标与外部指标源代码里至少要有一个聚类分析源代码如果只输出idx不输出评价指标那这份代码只能算完成了一半。评价指标分两类内部指标只看数据和聚类结果比如轮廓系数、Davies-Bouldin 指数、Calinski-Harabasz 指数外部指标需要真实标签比如调整兰德指数ARI、归一化互信息NMI。在 MATLAB 里内部指标用evalclusters最省事% 用 evalclusters 同时计算多种内部指标 eva evalclusters(X, (X,k) kmeans(X,k,Replicates,5), CalinskiHarabasz, KList, 2:6); plot(eva); % 画出 CH 指数随 k 的变化 bestK eva.OptimalK;这里第二个参数是一个函数句柄evalclusters 会在每个候选 k 上调用它。Calinski-Harabasz 指数越大表明聚类越紧致、类间越分散它比轮廓系数计算快很多适合先粗筛。如果你有真实标签 y计算 ARI 可以用文件交换上的randindexMATLAB 自带的confusionmat只能做标签对齐后的精确度不能直接当 ARI 用。外部指标还有一个容易踩的坑聚类返回的类号是 1、2、3真实标签可能是 0、1、2先把类号统一成连续正整数再算指标否则结果会低到离谱。我的建议是源代码里至少保留一个内部指标这样没有真实标签时也能向别人证明这个聚类不是随便分的。3. 用 MATLAB 跑通三种最常用聚类源代码K-Means、层次聚类与 DBSCAN 的最小可复现命令3.1 K-Means 最小可执行代码从 load 到 silhouetteK-Means 是绝大多数聚类分析源代码的默认算法因为它快、参数少、结果好解释。下面这一段是可以在 MATLAB 里直接复制运行的最小案例用的是自带数据集 fisheriris省去造数环节% K-Means 最小可复现案例 load fisheriris; % 加载内置鸢尾花数据 X meas; % 150 x 4 的特征矩阵 rng(1); % 固定随机种子保证结果可复现 [idx, C, sumd] kmeans(X, 3, Replicates, 10, MaxIter, 500); figure; gscatter(X(:,1), X(:,2), idx); % 用前两列画散点按聚类着色 hold on; plot(C(:,1), C(:,2), kx, MarkerSize, 15, LineWidth, 2); % 标出质心 title(K-Means 聚类结果 (k3));这里Replicates是重复次数设为 10 表示从 10 组不同初始中心里挑 WSS 最小的一次MaxIter是单次迭代上限500 次在 150 个样本上完全够用。C 是最终的 k 乘 p 质心矩阵sumd 是每个类内距离平方和后续算总 WSS 或做肘部图都要用到它。跑完这段代码你应该看到三类样本大致把山鸢尾和另外两种分开但 versicolor 和 virginica 会有部分重叠这说明单纯用 K-Means 不加核技巧线性边界很难完美区分重叠类别。如果样本量很大几十万行把Replicates降为 1MaxIter按收敛情况调。K-Means 的复杂度接近 O(nkp)n 是样本数k 是类别数p 是特征数。在 MATLAB 里它默认使用 kmeans 初始化比纯随机初始化稳定得多但依然要固定 rng否则你每次跑出的类标签顺序会漂移。3.2 层次聚类源代码怎么改linkage dendrogram cluster层次聚类适合你不知道聚类数、想看层级关系的场景源代码里一般就是 linkage、dendrogram、cluster 三连。它不要求预置 k而是先生成谱系图再按距离阈值或类别数截断。下面这段是常见写法% 层次聚类最小案例Ward 连接 谱系图 按类别数截断 load fisheriris; X meas; Z linkage(X, ward, euclidean); % 生成聚类树 figure; dendrogram(Z, 0); % 画出完整谱系图0 表示不限制叶子数 idx cluster(Z, MaxClust, 3); % 按最多 3 类截断linkage 的第一个输入可以是原始数据矩阵也可以是 pdist 输出的距离向量。ward方法要求距离度量是欧氏距离所以第三个参数写euclidean如果你用average或complete可以搭配pdist2得到任意距离矩阵再喂给 linkage但注意 linkage 第二个参数指定的是连接准则第三个参数才是距离度量。dendrogram 的 0 表示把所有叶子都画出来如果样本超过 300 个画出来会密密麻麻这时建议只画前 30 个叶子dendrogram(Z, 30)并按相干系数排序。层次聚类最容易被忽略的是cluster的截断方式。MaxClust是告诉 MATLAB 找到一种截断使类别数不超过指定值Cutoff是根据不一致系数截断对类间距离阈值更直观。如果你发现截断出来的类数比期望少先检查 linkage 是否用的是wardWard 法倾向于生成球形簇对异常值敏感换成average往往能缓解。谱系图的纵轴是连接距离距离越大表示两类越不像你可以在图上画一条横线看应该切哪里这比硬调 MaxClust 更符合直觉。3.3 DBSCAN 源代码与参数敏感性epsilon 和 minpts 的选法DBSCAN 在 MATLAB 里不是老版本标配R2019a 之后才有内置dbscan函数。如果你手上的源代码是早年下载的可能是一段自实现版本但用法大同小异。它的核心是两个参数邻域半径 epsilon 和最小邻域点数 minpts。下面用合成数据演示两个密度不同的簇加上一圈噪声正好暴露 DBSCAN 的优势% DBSCAN 最小案例合成双簇 噪声 rng(2); X [randn(100, 2) * 0.4; randn(80, 2) 2; rand(30, 2) * 3 - 1]; % 簇1、簇2、噪声 [idx, corepts] dbscan(X, 0.5, 10); % epsilon0.5, minpts10 gscatter(X(:,1), X(:,2), idx); title(DBSCAN 聚类结果噪声标为0);运行后你会看到 idx 中大部分样本被分为 1 和 2 两类但靠左下角那摊噪声会被标记为 0。dbscan 将 0 视为噪声点不归属任何簇。corepts是逻辑索引标出每个点是核心点还是边界点。epsilon 的取值是最大的坑设太小一个簇会被拆成好几块设太大所有点都变成一类噪声也吞进去。我常用的定参方法是 k-距离图% 用 k-距离图选 epsilon找曲线拐点 k 10; % 一般取 minpts [~, distK] knnsearch(X, X, K, k); % 每个样本到第 k 近邻的距离 distK sort(distK(:, end)); % 按升序排序 plot(distK); ylabel(第 k 近邻距离); xlabel(样本排序);这里的思路是对每个样本计算它到第 k 个近邻的距离按距离从小到大排序后画曲线。曲线会出现一个膝盖位置对应密度分布的分界那个位置的距离就是合适的 epsilon。minpts 的经验值是取特征维数的两倍以上二维数据取 10 可以高维数据要更大。DBSCAN 不能像 K-Means 那样直接指定聚类数如果你非要固定簇数得循环试参这也是很多聚类分析源代码里 DBSCAN 写得最潦草的原因。4. 聚类分析源代码避坑指南五个让结果翻车的隐藏问题4.1 现象聚类结果每次跑都不一样今天和明天对不上原因K-Means 和 DBSCAN 的初始状态依赖随机数。K-Means 在 MATLAB 里虽然用 kmeans 初始化但随机数种子不固定时每次调用得到的质心初始位置不同虽然大多数情况下会收敛到同一组簇但在数据分布重叠较大时可能收敛到不同的局部最优类标签的数字顺序也会变。DBSCAN 本身是确定性的但如果你用了随机抽样或打乱数据顺序结果也会漂移。解决在源代码开头写一行rng(2024)或者在调用 kmeans 时设置Options statset(UseSubstreams, true)。更稳妥的做法是把 rng 种子的值保存在结果结构体里每次跑记录 seed这样别人复现时能拿到完全一致的输出。我在自己的工程脚本里一律先rng(default)再调算法避免上次能跑出好结果这次死活复现不了的窘境。4.2 现象聚类结果里异常点被单独分成一类甚至把整个簇拽偏原因没有做特征归一化或离群点处理。K-Means 基于欧氏距离如果某个特征的量纲是 1000 而其他是 0.1距离几乎完全由那个大数特征决定离群点则会把质心拉向自己导致几个正常点被划成一类。解决先对特征做 zscore 标准化或 min-max 缩放再用 medfilt 或马氏距离剔除离群点。zscore 在 MATLAB 里一句话X zscore(X);这个操作应该放在聚类之前而不是之后。如果数据里存在极端长尾用robustcov配合马氏距离更稳。注意 zscore 的均值和标准差要保存下来等新数据进来时用同一组参数变换别直接重新算全局的 zscore。4.3 现象层次聚类画出的谱系图与 cluster 截断结果对不上原因dendrogram 的叶子顺序与 cluster 输出的标签顺序没有直接映射关系尤其是样本量超过 100 时谱系图会压缩显示。另外用Cutoff截断时MATLAB 默认用不一致系数而不是直接按距离阈值很多人没读文档就乱填数字。解决不要靠眼睛在图上数类数直接用cophenet(Z, D)检查聚类树与原始距离的相干系数数值越接近 1 表示谱系图保留的距离信息越好。要按距离截断明确写Cutoff, 1.5, Depth, 3其中 Depth 指定计算不一致系数时考虑的层级深度。最省心的方法是像我一样用cluster(Z, MaxClust, k)只负责定类数不碰阈值。4.4 现象DBSCAN 把所有点都归为一类或者全部标成噪声原因epsilon 选得太大所有点的邻域都彼此连通自然只剩一个簇epsilon 选得太小每个点都找不到足够邻居全成了噪声点。还有一种隐蔽情况输入数据没有标准化二维数据横纵轴量纲不同圆形邻域变成椭圆密度判断完全失真。解决先用 k-距离图定 epsilon再对每个候选 minpts 画一下聚类结果看稳定性。我这里有一个通用参数起点dbscan(X, 0.5, 10)对二维标准化数据通常可用但业务数据不要直接抄。更实际的办法是写一个 3 行循环把 epsilon 按几何间隔扫 10 个值记录每次的非噪声点数选一个非噪声点数突然上升后又稳定的区域。4.5 现象聚类结果与论文或同事的源代码对不上怀疑代码被改坏原因数据版本不同、特征顺序不同、随机种子不同这三样是高频元凶。聚类的类标签本身没有绝对顺序K-Means 里第 1 类可能是别人代码里的第 3 类直接对比标签数字完全没意义。解决不要比较 idx 的数字要比较聚类结果与已知标签的互信息或 ARI。先把两份代码的类标签用confusionmat做映射再看对角线占比。排查时依次检查读入文件的列顺序、是否做了同样的标准化、rng 种子是否一致。我自己的做法是在源代码里加一行注释记录数据文件的 MD5 校验和从源头避免数据不对的扯皮。5. 把源代码封装成自己的聚类工具函数自动搜 k、自动出图、落盘结果5.1 一个函数收编三种算法参数接口这样设计下载的源代码往往面向单次运行参数写死在脚本里。要应对换数据、换算法、换 k的日常最好封装成一个函数。下面是我常用的模板支持三种算法统一返回 idx 和最优 k并且自动固定随机种子function result clusterPipeline(X, varargin) %CLUSTERPIPELINE 聚类分析一站式封装 % result clusterPipeline(X, Method, kmeans, KRange, 2:8) % 返回值 result 包含 idx/bestK/silhouette 等字段 % 参数解析 p inputParser; addParameter(p, Method, kmeans, (x) ismember(x, {kmeans,hierarchical,dbscan})); addParameter(p, KRange, 2:8, (x) isnumeric(x) isvector(x)); addParameter(p, Norm, true, islogical); addParameter(p, Seed, 1, isscalar); parse(p, varargin{:}); opts p.Results; % 固定随机种子 归一化 rng(opts.Seed); if opts.Norm X zscore(X); end bestS -inf; bestK 2; switch opts.Method case kmeans for k opts.KRange idx kmeans(X, k, Replicates, 10, MaxIter, 500); s mean(silhouette(X, idx)); if s bestS bestS s; bestK k; bestIdx idx; end end case hierarchical D pdist(X, euclidean); Z linkage(D, ward); figure; dendrogram(Z, 0); title(层次聚类谱系图); for k opts.KRange idx cluster(Z, MaxClust, k); if numel(unique(idx)) 2, continue; end s mean(silhouette(X, idx)); if s bestS bestS s; bestK k; bestIdx idx; end end case dbscan % DBSCAN 不适合自动搜 k这里用固定参数演示 [bestIdx, core] dbscan(X, 0.5, 10); bestK numel(unique(bestIdx(bestIdx 0))); bestS mean(silhouette(X, bestIdx(bestIdx 0))); % 仅对非噪声计算 end % 打包结果 result struct(); result.idx bestIdx; result.bestK bestK; result.silhouette bestS; result.Method opts.Method; result.Seed opts.Seed; end参数设计上Method决定走哪条分支KRange只在 kmeans 和 hierarchical 分支下生效Norm控制是否标准化Seed固定随机性。注意 dbscan 分支里我把silhouette只用在非噪声样本上因为噪声点的类号为 0轮廓系数对 0 类会报错。这个函数的好处是你换数据时只需要改 X换算法时改一行 Method不用每处脚本都翻找参数。5.2 输出结构体与落盘结果、图、日志一次到位这段函数返回的 result 结构体在 MATLAB 里可以直接被后续脚本使用。我习惯在调用后追加一段落盘逻辑% 调用示例K-Means 聚类搜索 k2:6并保存全部结果 res clusterPipeline(X, Method, kmeans, KRange, 2:6, Seed, 42); save(cluster_result.mat, res); writematrix(res.idx, cluster_label.csv); % 把每行的类别标签导出writematrix 把 idx 写成 CSV适合后续做业务分析save 保存的 mat 文件则留着复现。如果你想连聚类中心一起存在函数返回里再加一个 C 字段。日志方面我一般用fprintf输出一行摘要包含 seed、最优 k 和轮廓系数这样跑批量实验时能在命令窗口直接看到对比不必打开 mat 文件。5.3 调用示例与边界数据量超过几万行怎么办上面的函数适合几千行、几十维以内的数据。如果样本量到十万行有几处必须改silhouette在十万样本上要算 n 乘 n 的距离矩阵内存直接爆炸。这时不要用轮廓系数选 k改用 Calinski-Harabasz 或干脆只跑肘部法。DBSCAN 的 knnsearch 在十万样本上也慢建议先抽样估计 epsilon再在全量上跑。我实际处理过二十万行的用户分群做法是先随机抽 20% 样本定参数再用那个参数跑全量时间从一小时降到五分钟。聚类分析源代码从教学到工程差别就在这些规模边界上。6. 验证聚类结果不是玄学一个置换检验的 MATLAB 实现费了大力气跑出聚类结果下一步最怕的是被别人问一句你怎么知道这个聚类不是把随机数据硬分出来的这个问题我吃过亏以前直接拿轮廓系数当证据结果审稿人让我证明簇的显著性。我后来用置换检验解决对每个特征独立打乱顺序让原始数据的结构被破坏再重新聚类看轮廓系数下降多少。如果真实数据的轮廓系数显著高于随机打乱后的分布说明聚类结果不是偶然。% 置换检验打乱特征验证聚类显著性 rng(1); res_obs clusterPipeline(X, Method, kmeans, KRange, 2:6, Seed, 1); obs_s res_obs.silhouette; null_s zeros(100, 1); for i 1:100 X_perm X(:, randperm(size(X,2))); % 按列打乱特征 % 注意randperm 只打乱列的顺序不破坏行内关系属于轻量置换 res_null clusterPipeline(X_perm, Method, kmeans, KRange, 2:6, Seed, i); null_s(i) res_null.silhouette; end p_value sum(null_s obs_s) / 100; fprintf(观测轮廓系数 %.3f, 置换检验 p %.3f\n, obs_s, p_value);这段代码跑了 100 次置换每次用不同的种子。p 值小于 0.05 就说明聚类结构显著强于随机打乱后的结果。更严格的做法是对每一列做随机重排而不是整列顺序不变但那样会破坏列间相关性检验力度过强实际项目中我常用列置换作为快速版本。这个技巧对 K-Means 和层次聚类都适用DBSCAN 也可以把轮廓系数换成非噪声点的平均数。跑完不要只贴 p 值把 null_s 的直方图画出来和 obs_s 的位置一起展示这是最有说服力的证据。我现在的习惯是任何聚类分析源代码交付前至少要跑一次这个置换检验然后把它写进注释里。固定种子、保留原始数据、记录参数这三件事看似不起眼却能让你的结果经得起反复问怎么复现。聚类分析在 MATLAB 里不难难的是让结果站得住脚。希望这篇笔记能帮你把源代码跑通也帮你把结果讲清楚。本文还有配套的精品资源点击获取