简介一份用于在MATLAB环境中实现K均值K-means聚类算法的完整资源包并附带了可加载的多维数据集面向机器学习初学者、数据分析人员以及正在完成相关实验课程的学生。压缩包共53个文件以.m脚本为主另有少量.c源文件、dll动态库、txt说明和.mat数据集整体容量仅58KB便于快速下载和二次修改。算法部分不仅包含kmeans.m、kmeansdemo.m等核心实现还提供了EM聚类、混合高斯模型、聚类统计评估与可视化脚本并配有test.mat和iris.txt等示例数据可直接运行验证聚类效果。目前已有528人浏览学习资源有助于理解K均值迭代原理、质心更新方式及常用评估指标同时也可作为拓展到其他聚类算法的入门参考。1. 从 iris.txt 到 kmeans.m一个亲手实现的聚类工具箱MATLAB 里调内置 kmeans 只需要一行但这行背后藏着初始化方式、距离度量、空簇处理、收敛容差和标签对齐五类问题任何一类没想清楚聚类结果都解释不通。解压这份源码主线非常清楚kmeans.m 是迭代主循环assign.m 负责样本分配dmean.m 负责质心更新dist1.m、manhattan.m、sqrDist.dll 组成距离度量层critsse.m、clusterstats.m 负责评估聚类质量mixtureEM.m、projectpca.m 又把方法延伸到软聚类和可视化。配合 iris.txt 与 test.txt 两个公开数据集能一次性跑通从底层实现到效果验证的完整链路。对把 K 均值当作业或实验模板的人来说这是少见的完整参考实现对长期只用官方 kmeans API 的人来说把这份代码读一遍相当于把官方文档里没写清楚的那些假设全部摊开重看一遍。2. 核心循环assign.m、dmean.m 与三种距离度量的取舍2.1 一次迭代的两个动作分配样本与更新质心源码包的主循环集中在 kmeans.m 里结构可以简化成下面这样。先随机选 k 个样本作为初始质心然后进入「分配 → 更新 → 判收敛」的循环。function [idx, C, sse] kmeans_demo(X, k, maxIter, distType) % X: NxD 特征矩阵, 每行一个样本 % k: 指定聚类数 % maxIter: 最大迭代次数 % distType: sqeuclidean 或 manhattan N size(X, 1); C X(randperm(N, k), :); % 初始质心从样本中随机抽取 prevDelta 0; for iter 1:maxIter idx assign(X, C, distType); % 每个样本分到最近质心 Cnew dmean(X, idx, k); % 每簇取均值作为新质心 delta sqrt(sum(sum((Cnew - C).^2, 2))); if abs(delta - prevDelta) 1e-8 % 震荡保护 break; end prevDelta delta; C Cnew; end sse critsse(X, idx, C); end这里assign.m内部不是简单写一个 min而是先构造 N×k 的距离矩阵再对每行取最小值下标。dmean.m则按 idx 分组求均值。代码里用abs(delta - prevDelta)而不是直接delta tol是因为某些数据集上质心会反复横跳但变化量本身趋于稳定这种震荡保护是手写 KMeans 时很容易漏的细节。critsse.m返回簇内平方和作为后续比较初始化好坏的标准。2.2 距离层设计dist1、sqrDist 与 manhattan这个包的距离函数分了三档dist1.m 是基础欧氏距离sqrDist.m 返回距离平方manhattan.m 实现 L1 距离。第一版 dist1.m 通常是双层循环function D dist1(X, C) % X: NxD, C: kxD, 返回 Nxk 距离矩阵 N size(X, 1); k size(C, 1); D zeros(N, k); for i 1:N for j 1:k D(i, j) sqrt(sum((X(i,:) - C(j,:)).^2)); end end end这个写法逻辑最直白但每轮迭代都要计算 O(NkD) 次数据量稍大就会拖慢整个实验。sqrDist.m 用矩阵展开消掉循环性能差距明显function D sqrDist(A, B) % A: NxD, B: MxD % 返回 NxM 的欧氏距离平方矩阵 na sum(A.^2, 2); nb sum(B.^2, 2); D -2 * A * B na nb; D max(D, 0); % 浮点误差可能让距离变负截断掉 end这里不开平方因为比较距离大小是单调操作不影响分配结果。sqrDist.dll是预编译的 C 版本处理大规模矩阵时比纯 MATLAB 循环再快一截。三种距离的适用场景差别很大距离函数适用特征质心更新方式典型问题dist1 / sqrDist连续稠密特征dmean 均值对异常值敏感manhattan高维稀疏、含离群值应改为中位数包内 dmean 默认不自动切换自定义距离函数业务距离如余弦按原型找最密集点需要自行扩展代码manhattan 距离对离群值的惩罚比欧氏温柔但它的几何中心不再是均值应该是每个维度上的中位数。包里的dmean.m默认不区分距离类型现场使用时要额外写一个切换分支否则 L1 距离配均值质心很容易漂移。2.3 收敛判断与空簇兜底很多新手版本的 kmenas 用isequal(Cnew, C)判断收敛在浮点情况下几乎不可能相等最后只能跑满 maxIter。更工程化的做法是看质心移动量的范数或者配合 clusterstats.m 观察每轮簇内方差。clusterstats.m 会输出每个簇的样本数、均值、方差如果某一轮某个簇的方差突然反向增大多半是空簇导致质心被某个远端点拉走。空簇处理在包里由move.m承担当某个簇在分配后没有样本就把该质心重新放到距离所有现有质心最远的样本上。这个策略比原地保留空质心更稳因为它给了空簇重新吸引样本的机会。顺带提一个细节压缩包里出现的kmeansdemo.asv是 MATLAB 自动保存文件不是作者故意放的残留但把.asv和正式的 demo 脚本对比往往能看到空簇处理是后来补上的这对理解迭代设计很有说服力。3. K 值怎么定、初始化怎么救critsse、Replicates 与稳定复现3.1 手肘法critsse.m 与簇内平方和的拐点critsse.m 计算的是所有样本到所属质心的距离平方和也就是 SSE。手肘法的思路是K 太小则簇内距离很大K 太大又失去压缩意义所以取 SSE 曲线拐点处的 K。下面的脚本在 iris 数据上自动扫 Kclear; data load(iris.txt); % 每行一个样本最后一列是真实标签 X data(:, 1:4); sseList zeros(1, 7); for k 2:8 totalSSE 0; for rep 1:10 [idx, C] kmeans(X, k, Distance, sqEuclidean); totalSSE totalSSE critsse(X, idx, C) / 10; end sseList(k - 1) totalSSE; end plot(2:8, sseList, o-); xlabel(K); ylabel(SSE);这里对每个 K 做 10 次平均是因为 KMeans 单次运行受初始化影响很大直接画原始曲线可能抖动严重。iris 数据集大约在 k3 处出现拐点再往上增加 KSSE 的下降速度会明显放缓。手肘法不是自动算法它需要人看图判断但如果曲线一直平滑没有明显拐点那就说明数据本身没有天然簇结构这时候再纠结 K 没有意义。3.2 多次运行取最优Replicates 与 restartEM.m随机初始质心容易落在局部最优单次 KMeans 的 SSE 可能比全局最优高不少。标准做法是重复跑多次保留 SSE 最小的一次。内置kmeans的Replicates参数就是这个逻辑包内虽然没有完全封装成一个参数但我们自己可以循环实现bestIdx []; bestC []; bestSSE Inf; for rep 1:20 rng(rep); [idx, C] kmeans(X, k, Distance, sqEuclidean); currentSSE critsse(X, idx, C); if currentSSE bestSSE bestSSE currentSSE; bestIdx idx; bestC C; end end配合rng设置种子可以让实验可复现。这个模式可以类比压缩包里的restartEM.m它服务于 mixtureEM.m作用同样是搜索多个初始点后保留最佳解。实际使用时Replicates 数量不要盲目加到几百数据量上去后每次都很贵建议先跑 5 到 10 次看 SSE 的波动范围如果多次结果一致就没有必要加。3.3 初始化不只有 randomkmeans 与 kmeans.m 的 Start 参数随机从数据里抽质心虽然简单但很容易把两个初始质心放进同一个簇。KMeans 的思路是第一个质心随机之后每个新质心离已有质心越远、选中概率越大。代码上可以在选质心时插入距离判断function C kmeansPPinit(X, k) % KMeans 初始化 N size(X, 1); C zeros(k, size(X, 2)); rng(1); C(1, :) X(randi(N), :); for j 2:k D zeros(N, 1); for m 1:j-1 d sqrDist(X, C(m, :)); D min(D, d); end p D / sum(D); cumP cumsum(p); t rand(); C(j, :) X(find(cumP t, 1), :); end end这段代码里sqrDist返回的是距离平方直接用平方作为权重可以避免开平方计算。如果包里的 kmeans.m 不支持Start, plus就可以把上面的 C 作为初始质心传入。初始化对最终结果的影响经常被低估尤其是簇大小不均匀的数据好的初始质心可以让迭代少走很多弯路。4. 从硬聚类到软聚类mixtureEM.m、projectpca.m 与标签对齐4.1 mixtureEM.m当分配不再是“非此即彼”KMeans 的硬分配只告诉样本属于哪个簇mixtureEM.m 则给出后验概率。EM 的经典流程分两步E 步计算每个样本属于每个高斯分量的概率M 步用加权样本重新估计均值、协方差和混合系数。核心代码可以简化成这样function [w, mu, Sigma] mixtureEM(X, k, maxIter) % X: NxD 数据, k: 高斯分量数 % w: 1xk 混合系数, mu: kxD, Sigma: DxDxk [N, D] size(X); rng(42); mu X(randperm(N, k), :); Sigma repmat(0.1 * eye(D), 1, 1, k); w ones(1, k) / k; for iter 1:maxIter % E-step: 后验概率矩阵 R, Nxk R zeros(N, k); for j 1:k R(:, j) w(j) * mvnpdf(X, mu(j, :), Sigma(:, :, j)); end R R ./ max(sum(R, 2), eps); % M-step: 加权重估计 Nk sum(R, 1); w Nk / N; for j 1:k mu(j, :) sum(X .* R(:, j), 1) / Nk(j); dX X - mu(j, :); Sigma(:, :, j) (dX .* R(:, j)) * dX / Nk(j) 1e-6 * eye(D); end end endE 步里mvnpdf对每个样本计算高斯密度然后把密度乘上混合权重再归一化就是后验概率M 步的加权均值就是软版本的dmean.m。当所有分量协方差都趋近 0 时后验概率会退化成 0/1 硬分配所以 GMM 可以被看成 KMeans 的软聚类扩展。mixtureSelect.m用来在多个分量数之间选模型通常配合 BIC 或 AIC 一起使用比直接看 LLH 更稳。4.2 projectpca.m 与 showpca3.m在 3D 里检验聚类边界高维聚类结果最难判断的就是边界是否合理projectpca.m 用 SVD 做主成分降维把数据投到前 m 个主成分上function [score, eigvec] projectpca(X, m) % X: NxD, m: 降到 m 维 % score: Nxm 投影坐标, eigvec: Dxm 主成分方向 Xc X - mean(X, 1); [U, S, V] svd(Xc); score U(:, 1:m) * S(1:m, 1:m); eigvec V(:, 1:m); endshowpca3.m会取前三个主成分做三维散点图然后根据 kmeans 输出的 idx 给点上色。这种做法有两个用途第一是在聚类前降维观察数据是不是真的分得开第二是聚类后把结果投影到 3D 空间直观看到哪些样本被分到了边界处。注意 PCA 只是线性投影如果簇是非线性结构散点图上重叠并不代表 KMeans 一定失败需要再结合轮廓系数判断。4.3 misclass.m簇编号与真实标签的强对齐KMeans 输出的簇编号是 0/1/2但编号顺序和真实标签没有对应关系不能直接拿 idx 和真实标签算准确率。misclass.m 和 majority1.m 处理的正是这个问题先对每个簇做多数投票把簇编号映射到出现次数最多的真实类别。function [mappedLabel, acc] majority1(idx, trueLabel, k) % idx: KMeans 输出 Nx1, 值为 1..k % trueLabel: 真实标签, 也是 1..C mappedLabel zeros(size(idx)); for c 1:k member (idx c); mappedLabel(member) mode(trueLabel(member)); end acc mean(mappedLabel trueLabel); endmode取簇内出现最多的真实标签作为映射目标。这步完成后才能计算误分率。压缩包里的misclass.m在内部调用majority1.m最后输出混淆矩阵或错误样本下标。如果映射后准确率还很低不要急着调 K先检查数据标准化再检查距离度量是否与特征语义匹配。5. 把工具包迁移到自己的数据集预处理、mex 编译和轮廓系数定 K5.1 特征数据怎样进 kmeansloadiris.m 的读入套路loadiris.m 的职责是把 iris.txt 读成特征矩阵 X 和标签向量 y。iris.txt 前四列是花萼和花瓣的长宽最后一列是类别编号。读文本用load即可但更通用的做法是写一个带文件名入参的读取函数function [X, y] loadiris(filename) % 默认读 iris.txt if nargin 1 filename iris.txt; end data load(filename); X data(:, 1:4); y data(:, 5); end接入自己的数据时只需要保证文本或 CSV 文件里每一行是一个样本列和列之间用空格或制表符分开。真正隐蔽的坑是量纲欧氏距离默认把所有特征当成同尺度如果一列是毫米、另一列是温度量纲大的特征会完全主导距离。所以读取数据后第一步永远是标准化s std(X, 1); X (X - mean(X, 1)) ./ (s eps);eps是为了避免某列方差为 0 时产生除零错误。标准化之后再做 PCA、KMeans 和距离计算结果才有可比性。压缩包里的 test.txt 也是同样的结构用loadtest.m读取后可以先跑一遍kmeansTestdemo.m确认函数调用链没有断裂再替换成自己的数据。5.2 sqrDist.dll 不兼容时的 mex 编译路径sqrDist.dll和dist1.dll是预编译产物但如果 MATLAB 版本换到 64 位或者换到 Linux/macOSdll 就会加载失败。遇到这种情况不需要慌包里的 .c 文件可以重新编译。执行下面的命令mex -setup c mex sqrDist.c mex dist1.c mex mygetfield.cmex -setup c会列出当前系统的 C 编译器如果提示没有编译器需要在 MATLAB 附加功能里安装 MinGW-w64。编译完成后调用exist(sqrDist, file) 3可以确认 mex 文件是否可用。更稳妥的做法是让代码自动降级if exist(sqrDist, file) 3 D sqrDist(X, C); else D dist1(X, C); end这样即使 dll 无法运行也能退回纯 MATLAB 的 dist1.m。线上实验的时候我会提前把编译状态写进脚本日志避免哪一天换机器后得到完全不同的聚类结果。5.3 最后一个技巧用轮廓系数给 K 打分手肘法看的是整体 SSE轮廓系数则能反映单个样本的分离度。每个样本的轮廓值接近 1说明它离自己的簇很紧、离相邻簇很远接近 0 说明落在两个簇边界的模糊地带负值则说明它更可能被分错了。MATLAB 里可以直接用自带 silhouette 函数也可以调用 evalclusters 自动扫 Kdata load(iris.txt); X data(:, 1:4); X (X - mean(X, 1)) ./ (std(X, 1) eps); eva evalclusters(X, kmeans, silhouette, KList, 2:6); fprintf(最佳 K %d\n, eva.OptimalK);evalclusters 会自动把每个 K 对应的平均轮廓值算出来。实际使用时我会把手肘法和轮廓系数交叉验证手肘法说 k3轮廓系数也接近峰值就放心用 3如果两个指标打架再把Replicates从 10 提到 30排除初始化带来的偶然性。换数据集时把数据路径、K 列表和 mex 编译状态三个配置项写在脚本开头后续实验就只需要改路径不需要再动算法逻辑。本文还有配套的精品资源点击获取 MATLAB手写KMeans聚类：从原理到工程实践
简介一份用于在MATLAB环境中实现K均值K-means聚类算法的完整资源包并附带了可加载的多维数据集面向机器学习初学者、数据分析人员以及正在完成相关实验课程的学生。压缩包共53个文件以.m脚本为主另有少量.c源文件、dll动态库、txt说明和.mat数据集整体容量仅58KB便于快速下载和二次修改。算法部分不仅包含kmeans.m、kmeansdemo.m等核心实现还提供了EM聚类、混合高斯模型、聚类统计评估与可视化脚本并配有test.mat和iris.txt等示例数据可直接运行验证聚类效果。目前已有528人浏览学习资源有助于理解K均值迭代原理、质心更新方式及常用评估指标同时也可作为拓展到其他聚类算法的入门参考。1. 从 iris.txt 到 kmeans.m一个亲手实现的聚类工具箱MATLAB 里调内置 kmeans 只需要一行但这行背后藏着初始化方式、距离度量、空簇处理、收敛容差和标签对齐五类问题任何一类没想清楚聚类结果都解释不通。解压这份源码主线非常清楚kmeans.m 是迭代主循环assign.m 负责样本分配dmean.m 负责质心更新dist1.m、manhattan.m、sqrDist.dll 组成距离度量层critsse.m、clusterstats.m 负责评估聚类质量mixtureEM.m、projectpca.m 又把方法延伸到软聚类和可视化。配合 iris.txt 与 test.txt 两个公开数据集能一次性跑通从底层实现到效果验证的完整链路。对把 K 均值当作业或实验模板的人来说这是少见的完整参考实现对长期只用官方 kmeans API 的人来说把这份代码读一遍相当于把官方文档里没写清楚的那些假设全部摊开重看一遍。2. 核心循环assign.m、dmean.m 与三种距离度量的取舍2.1 一次迭代的两个动作分配样本与更新质心源码包的主循环集中在 kmeans.m 里结构可以简化成下面这样。先随机选 k 个样本作为初始质心然后进入「分配 → 更新 → 判收敛」的循环。function [idx, C, sse] kmeans_demo(X, k, maxIter, distType) % X: NxD 特征矩阵, 每行一个样本 % k: 指定聚类数 % maxIter: 最大迭代次数 % distType: sqeuclidean 或 manhattan N size(X, 1); C X(randperm(N, k), :); % 初始质心从样本中随机抽取 prevDelta 0; for iter 1:maxIter idx assign(X, C, distType); % 每个样本分到最近质心 Cnew dmean(X, idx, k); % 每簇取均值作为新质心 delta sqrt(sum(sum((Cnew - C).^2, 2))); if abs(delta - prevDelta) 1e-8 % 震荡保护 break; end prevDelta delta; C Cnew; end sse critsse(X, idx, C); end这里assign.m内部不是简单写一个 min而是先构造 N×k 的距离矩阵再对每行取最小值下标。dmean.m则按 idx 分组求均值。代码里用abs(delta - prevDelta)而不是直接delta tol是因为某些数据集上质心会反复横跳但变化量本身趋于稳定这种震荡保护是手写 KMeans 时很容易漏的细节。critsse.m返回簇内平方和作为后续比较初始化好坏的标准。2.2 距离层设计dist1、sqrDist 与 manhattan这个包的距离函数分了三档dist1.m 是基础欧氏距离sqrDist.m 返回距离平方manhattan.m 实现 L1 距离。第一版 dist1.m 通常是双层循环function D dist1(X, C) % X: NxD, C: kxD, 返回 Nxk 距离矩阵 N size(X, 1); k size(C, 1); D zeros(N, k); for i 1:N for j 1:k D(i, j) sqrt(sum((X(i,:) - C(j,:)).^2)); end end end这个写法逻辑最直白但每轮迭代都要计算 O(NkD) 次数据量稍大就会拖慢整个实验。sqrDist.m 用矩阵展开消掉循环性能差距明显function D sqrDist(A, B) % A: NxD, B: MxD % 返回 NxM 的欧氏距离平方矩阵 na sum(A.^2, 2); nb sum(B.^2, 2); D -2 * A * B na nb; D max(D, 0); % 浮点误差可能让距离变负截断掉 end这里不开平方因为比较距离大小是单调操作不影响分配结果。sqrDist.dll是预编译的 C 版本处理大规模矩阵时比纯 MATLAB 循环再快一截。三种距离的适用场景差别很大距离函数适用特征质心更新方式典型问题dist1 / sqrDist连续稠密特征dmean 均值对异常值敏感manhattan高维稀疏、含离群值应改为中位数包内 dmean 默认不自动切换自定义距离函数业务距离如余弦按原型找最密集点需要自行扩展代码manhattan 距离对离群值的惩罚比欧氏温柔但它的几何中心不再是均值应该是每个维度上的中位数。包里的dmean.m默认不区分距离类型现场使用时要额外写一个切换分支否则 L1 距离配均值质心很容易漂移。2.3 收敛判断与空簇兜底很多新手版本的 kmenas 用isequal(Cnew, C)判断收敛在浮点情况下几乎不可能相等最后只能跑满 maxIter。更工程化的做法是看质心移动量的范数或者配合 clusterstats.m 观察每轮簇内方差。clusterstats.m 会输出每个簇的样本数、均值、方差如果某一轮某个簇的方差突然反向增大多半是空簇导致质心被某个远端点拉走。空簇处理在包里由move.m承担当某个簇在分配后没有样本就把该质心重新放到距离所有现有质心最远的样本上。这个策略比原地保留空质心更稳因为它给了空簇重新吸引样本的机会。顺带提一个细节压缩包里出现的kmeansdemo.asv是 MATLAB 自动保存文件不是作者故意放的残留但把.asv和正式的 demo 脚本对比往往能看到空簇处理是后来补上的这对理解迭代设计很有说服力。3. K 值怎么定、初始化怎么救critsse、Replicates 与稳定复现3.1 手肘法critsse.m 与簇内平方和的拐点critsse.m 计算的是所有样本到所属质心的距离平方和也就是 SSE。手肘法的思路是K 太小则簇内距离很大K 太大又失去压缩意义所以取 SSE 曲线拐点处的 K。下面的脚本在 iris 数据上自动扫 Kclear; data load(iris.txt); % 每行一个样本最后一列是真实标签 X data(:, 1:4); sseList zeros(1, 7); for k 2:8 totalSSE 0; for rep 1:10 [idx, C] kmeans(X, k, Distance, sqEuclidean); totalSSE totalSSE critsse(X, idx, C) / 10; end sseList(k - 1) totalSSE; end plot(2:8, sseList, o-); xlabel(K); ylabel(SSE);这里对每个 K 做 10 次平均是因为 KMeans 单次运行受初始化影响很大直接画原始曲线可能抖动严重。iris 数据集大约在 k3 处出现拐点再往上增加 KSSE 的下降速度会明显放缓。手肘法不是自动算法它需要人看图判断但如果曲线一直平滑没有明显拐点那就说明数据本身没有天然簇结构这时候再纠结 K 没有意义。3.2 多次运行取最优Replicates 与 restartEM.m随机初始质心容易落在局部最优单次 KMeans 的 SSE 可能比全局最优高不少。标准做法是重复跑多次保留 SSE 最小的一次。内置kmeans的Replicates参数就是这个逻辑包内虽然没有完全封装成一个参数但我们自己可以循环实现bestIdx []; bestC []; bestSSE Inf; for rep 1:20 rng(rep); [idx, C] kmeans(X, k, Distance, sqEuclidean); currentSSE critsse(X, idx, C); if currentSSE bestSSE bestSSE currentSSE; bestIdx idx; bestC C; end end配合rng设置种子可以让实验可复现。这个模式可以类比压缩包里的restartEM.m它服务于 mixtureEM.m作用同样是搜索多个初始点后保留最佳解。实际使用时Replicates 数量不要盲目加到几百数据量上去后每次都很贵建议先跑 5 到 10 次看 SSE 的波动范围如果多次结果一致就没有必要加。3.3 初始化不只有 randomkmeans 与 kmeans.m 的 Start 参数随机从数据里抽质心虽然简单但很容易把两个初始质心放进同一个簇。KMeans 的思路是第一个质心随机之后每个新质心离已有质心越远、选中概率越大。代码上可以在选质心时插入距离判断function C kmeansPPinit(X, k) % KMeans 初始化 N size(X, 1); C zeros(k, size(X, 2)); rng(1); C(1, :) X(randi(N), :); for j 2:k D zeros(N, 1); for m 1:j-1 d sqrDist(X, C(m, :)); D min(D, d); end p D / sum(D); cumP cumsum(p); t rand(); C(j, :) X(find(cumP t, 1), :); end end这段代码里sqrDist返回的是距离平方直接用平方作为权重可以避免开平方计算。如果包里的 kmeans.m 不支持Start, plus就可以把上面的 C 作为初始质心传入。初始化对最终结果的影响经常被低估尤其是簇大小不均匀的数据好的初始质心可以让迭代少走很多弯路。4. 从硬聚类到软聚类mixtureEM.m、projectpca.m 与标签对齐4.1 mixtureEM.m当分配不再是“非此即彼”KMeans 的硬分配只告诉样本属于哪个簇mixtureEM.m 则给出后验概率。EM 的经典流程分两步E 步计算每个样本属于每个高斯分量的概率M 步用加权样本重新估计均值、协方差和混合系数。核心代码可以简化成这样function [w, mu, Sigma] mixtureEM(X, k, maxIter) % X: NxD 数据, k: 高斯分量数 % w: 1xk 混合系数, mu: kxD, Sigma: DxDxk [N, D] size(X); rng(42); mu X(randperm(N, k), :); Sigma repmat(0.1 * eye(D), 1, 1, k); w ones(1, k) / k; for iter 1:maxIter % E-step: 后验概率矩阵 R, Nxk R zeros(N, k); for j 1:k R(:, j) w(j) * mvnpdf(X, mu(j, :), Sigma(:, :, j)); end R R ./ max(sum(R, 2), eps); % M-step: 加权重估计 Nk sum(R, 1); w Nk / N; for j 1:k mu(j, :) sum(X .* R(:, j), 1) / Nk(j); dX X - mu(j, :); Sigma(:, :, j) (dX .* R(:, j)) * dX / Nk(j) 1e-6 * eye(D); end end endE 步里mvnpdf对每个样本计算高斯密度然后把密度乘上混合权重再归一化就是后验概率M 步的加权均值就是软版本的dmean.m。当所有分量协方差都趋近 0 时后验概率会退化成 0/1 硬分配所以 GMM 可以被看成 KMeans 的软聚类扩展。mixtureSelect.m用来在多个分量数之间选模型通常配合 BIC 或 AIC 一起使用比直接看 LLH 更稳。4.2 projectpca.m 与 showpca3.m在 3D 里检验聚类边界高维聚类结果最难判断的就是边界是否合理projectpca.m 用 SVD 做主成分降维把数据投到前 m 个主成分上function [score, eigvec] projectpca(X, m) % X: NxD, m: 降到 m 维 % score: Nxm 投影坐标, eigvec: Dxm 主成分方向 Xc X - mean(X, 1); [U, S, V] svd(Xc); score U(:, 1:m) * S(1:m, 1:m); eigvec V(:, 1:m); endshowpca3.m会取前三个主成分做三维散点图然后根据 kmeans 输出的 idx 给点上色。这种做法有两个用途第一是在聚类前降维观察数据是不是真的分得开第二是聚类后把结果投影到 3D 空间直观看到哪些样本被分到了边界处。注意 PCA 只是线性投影如果簇是非线性结构散点图上重叠并不代表 KMeans 一定失败需要再结合轮廓系数判断。4.3 misclass.m簇编号与真实标签的强对齐KMeans 输出的簇编号是 0/1/2但编号顺序和真实标签没有对应关系不能直接拿 idx 和真实标签算准确率。misclass.m 和 majority1.m 处理的正是这个问题先对每个簇做多数投票把簇编号映射到出现次数最多的真实类别。function [mappedLabel, acc] majority1(idx, trueLabel, k) % idx: KMeans 输出 Nx1, 值为 1..k % trueLabel: 真实标签, 也是 1..C mappedLabel zeros(size(idx)); for c 1:k member (idx c); mappedLabel(member) mode(trueLabel(member)); end acc mean(mappedLabel trueLabel); endmode取簇内出现最多的真实标签作为映射目标。这步完成后才能计算误分率。压缩包里的misclass.m在内部调用majority1.m最后输出混淆矩阵或错误样本下标。如果映射后准确率还很低不要急着调 K先检查数据标准化再检查距离度量是否与特征语义匹配。5. 把工具包迁移到自己的数据集预处理、mex 编译和轮廓系数定 K5.1 特征数据怎样进 kmeansloadiris.m 的读入套路loadiris.m 的职责是把 iris.txt 读成特征矩阵 X 和标签向量 y。iris.txt 前四列是花萼和花瓣的长宽最后一列是类别编号。读文本用load即可但更通用的做法是写一个带文件名入参的读取函数function [X, y] loadiris(filename) % 默认读 iris.txt if nargin 1 filename iris.txt; end data load(filename); X data(:, 1:4); y data(:, 5); end接入自己的数据时只需要保证文本或 CSV 文件里每一行是一个样本列和列之间用空格或制表符分开。真正隐蔽的坑是量纲欧氏距离默认把所有特征当成同尺度如果一列是毫米、另一列是温度量纲大的特征会完全主导距离。所以读取数据后第一步永远是标准化s std(X, 1); X (X - mean(X, 1)) ./ (s eps);eps是为了避免某列方差为 0 时产生除零错误。标准化之后再做 PCA、KMeans 和距离计算结果才有可比性。压缩包里的 test.txt 也是同样的结构用loadtest.m读取后可以先跑一遍kmeansTestdemo.m确认函数调用链没有断裂再替换成自己的数据。5.2 sqrDist.dll 不兼容时的 mex 编译路径sqrDist.dll和dist1.dll是预编译产物但如果 MATLAB 版本换到 64 位或者换到 Linux/macOSdll 就会加载失败。遇到这种情况不需要慌包里的 .c 文件可以重新编译。执行下面的命令mex -setup c mex sqrDist.c mex dist1.c mex mygetfield.cmex -setup c会列出当前系统的 C 编译器如果提示没有编译器需要在 MATLAB 附加功能里安装 MinGW-w64。编译完成后调用exist(sqrDist, file) 3可以确认 mex 文件是否可用。更稳妥的做法是让代码自动降级if exist(sqrDist, file) 3 D sqrDist(X, C); else D dist1(X, C); end这样即使 dll 无法运行也能退回纯 MATLAB 的 dist1.m。线上实验的时候我会提前把编译状态写进脚本日志避免哪一天换机器后得到完全不同的聚类结果。5.3 最后一个技巧用轮廓系数给 K 打分手肘法看的是整体 SSE轮廓系数则能反映单个样本的分离度。每个样本的轮廓值接近 1说明它离自己的簇很紧、离相邻簇很远接近 0 说明落在两个簇边界的模糊地带负值则说明它更可能被分错了。MATLAB 里可以直接用自带 silhouette 函数也可以调用 evalclusters 自动扫 Kdata load(iris.txt); X data(:, 1:4); X (X - mean(X, 1)) ./ (std(X, 1) eps); eva evalclusters(X, kmeans, silhouette, KList, 2:6); fprintf(最佳 K %d\n, eva.OptimalK);evalclusters 会自动把每个 K 对应的平均轮廓值算出来。实际使用时我会把手肘法和轮廓系数交叉验证手肘法说 k3轮廓系数也接近峰值就放心用 3如果两个指标打架再把Replicates从 10 提到 30排除初始化带来的偶然性。换数据集时把数据路径、K 列表和 mex 编译状态三个配置项写在脚本开头后续实验就只需要改路径不需要再动算法逻辑。本文还有配套的精品资源点击获取