企业官网建设流程全解析

简介：这份《k-means数据聚类实验报告（MATLAB详细代码）》文档面向学习聚类分析、数据挖掘课程实验及需要MATLAB实现参考的学生与初级算法实践者，系统梳理k-means的无监督划分思路与实验过程。文档围绕算法原理、流程步骤、Iris数据集聚类实验展开，给出k=3与k=4的聚类效果对比，并讨论初始类中心、数据输入顺序、k值选择对迭代次数和结果稳定性的影响。包内共1个doc文件，约126KB，以实验报告正文和MATLAB代码附录为主，便于直接对照阅读算法描述与实现细节。已有2743人学习，适合用于课程实验报告撰写、算法复现、参数敏感性分析与聚类方法入门。读者可从中获取完整的实验目的、步骤、结果分析、结论心得及可参考的MATLAB主程序代码，快速理解k-means的收敛条件、平方误差准则和实际应用局限。

1. 一份 k-means 实验报告被退回：能跑通和能复现差在哪

某高校实验室交上去的一份 k-means 数据聚类实验报告被打回来两次，理由不在结果对错，而在导师的三个追问：换台机器重跑，散点图的颜色还会一样吗；K 取 4，除了肘部图上的那个拐点，还有没有别的依据；表里那一列质心坐标，能不能指着代码说清它是哪一步算出来的。

这三个问题指向同一件事：k-means 在 MATLAB 里跑通只要十几行，但要让别人照着复现、让结论经得起问，得把随机种子、初始化方式、K 的选择准则、结果回读四件事写进代码，而不是留在口头。下面的顺序是：先把目标函数写成看得懂的矩阵运算，再把数据管线搭起来，然后讨论 K 怎么定，最后把踩过的坑和固化脚本的做法摆出来。适合正在写课程实验报告、要给业务方出一版分群结果，或者打算把 MATLAB 原型交付出去的人。

2. 把 k-means 的目标函数写进 MATLAB：向量化距离、质心更新与空簇处理

k-means 的全部内容就一句话：反复做“把样本分给最近的质心”和“把质心挪到簇内均值”两件事，直到质心不再动。最小化的目标是簇内平方误差和 $J=\sum_{i}|x_i-c_{l_i}|^2$。它收敛是因为每一步都不会让 $J$ 变大，但它只保证收敛到局部最优——这一点决定了后面所有关于初始化和多次重启的讨论。

要写成能交差的实验报告，光调kmeans函数不够。工具箱函数把初始化、重启、空簇处理都收进了黑匣子，报告里说不清哪一步影响了结果，也没法在缺工具箱的机器上跑。自己写一份核心循环，二十来行，可控性完全不一样。

2.1 距离矩阵用展开式算，别拼三维数组

新手最容易写成三层循环或者n×k×d的三维数组，样本上到五万行、特征几十维，内存直接顶到天花板。平方欧氏距离有现成的展开式：$|x-c|^2=|x|^2-2x^\top c+|c|^2$，右边三项都是矩阵乘法，交给 BLAS 处理。

function [labels, centroids, Jhist] = kmeans_core(X, k, seed, maxIter, tol) % 输入 % X n x d 数值特征矩阵，一行一个样本，不能含 NaN % k 簇数 % seed 随机种子，实验报告里必须写死，否则每次图都不一样 % maxIter 最大迭代轮数，默认 300 % tol 质心最大位移容差，小于它就认为收敛 % 输出 % labels n x 1 簇编号 % centroids k x d 质心坐标 % Jhist 每轮的簇内平方和，用来画收敛曲线 if nargin < 4 || isempty(maxIter), maxIter = 300; end if nargin < 5 || isempty(tol), tol = 1e-6; end rng(seed); % 固定随机流，复现的第一步 n = size(X, 1); if k > n, error('簇数不能大于样本数'); end centroids = X(randperm(n, k), :); % 从真实样本里挑 k 个做初值 labels = zeros(n, 1); Jhist = zeros(maxIter, 1); for it = 1:maxIter % 平方欧氏距离展开：n x k，一次矩阵乘法出结果 D = sum(X.^2, 2) + sum(centroids.^2, 2)' - 2 * (X * centroids'); D = max(D, 0); % 浮点误差可能给出 -1e-14，截掉 [~, labels] = min(D, [], 2); % 第二个参数 2 表示沿列方向取最小 counts = accumarray(labels, 1, [k 1]); % 每簇样本数 sums = accumarray(labels, X, [k size(X, 2)]); % 每簇特征求和 newC = sums ./ max(counts, 1); % 空簇先占位，下面补 % 空簇兜底：把离自己质心最远的样本提上来当新质心 for e = find(counts == 0)' own = D(sub2ind(size(D), (1:n)', labels)); [~, farIdx] = max(own); newC(e, :) = X(farIdx, :); labels(farIdx) = e; end Jhist(it) = sum(min(D, [], 2)); % 当前分配的 J shift = max(sqrt(sum((newC - centroids).^2, 2))); % 质心最大位移 centroids = newC; if shift < tol, break; end % 位移够小就停 end Jhist = Jhist(1:it); end

X必须是double，整型矩阵做.^2会溢出；D = max(D,0)这行看着多余，但在维度高、数值量级大的数据上不加会出负距离，紧接着min选出的归属就会莫名其妙。min(D,[],2)里的2是维度参数，漏了它会把整个矩阵压成一个标量，症状是所有样本被分到同一簇。

这里没用pdist2，因为它是统计工具箱的函数，换个机房、换台笔记本就可能报未定义。展开式只依赖基础矩阵运算，代价是数值稳定性略差，靠max(D,0)兜住就够了。距离度量也不一定非欧氏不可：特征稀疏、维度高（比如文本向量）时余弦距离更合适，把D换成1 - X*centroids' ./ (vecnorm(X,2,2) * vecnorm(centroids,2,2)')即可；特征是有量纲的计数且差异体现在绝对偏移上时，曼哈顿距离对异常值更稳。选哪种要在报告里写一句理由，别默认欧氏。

2.2 质心更新、空簇与收敛判据的三个细节

accumarray是这段代码里最值得记住的函数：第一个参数是分组下标（正整数），第二个是待累加的值，第三个是输出尺寸。accumarray(labels, X, [k size(X,2)])直接得到 k 行、每行是该簇所有样本的特征和，比for j = 1:k, sum(X(labels==j,:))快一个量级，因为后者每轮都要做 k 次全表比较。max(counts,1)是为了防零除——如果某簇在这一轮变空，sums那一行是全零，除以 0 会得到NaN，然后NaN会沿着后面的距离计算污染整个矩阵，最终所有标签都变成 1。

空簇在真实数据上并不罕见，尤其是初始质心挨得太近、或者某一簇本身只有几个离群点。两种处理方式：一是代码里的重新播种，把离自己质心最远的样本点拉过来当新质心，好处是几乎不增加计算量；二是合并到最近的簇，然后 k 减一。前者适合报告里演示，后者适合 k 有业务含义、不允许变的场景。我一般用前者，但在日志里打一行警告，因为空簇频繁出现说明 k 给大了或者数据本身有高度重复的点。

tol控制的是质心最大欧氏位移，不是 J 的相对变化。1e-6 在标准化后的数据上通常 10 到 30 轮收敛，如果特征量级在 1e6 以上，这个阈值等于没有，得按特征尺度调到 1e-3 之类。Jhist别看它只是个数列，画出来能立刻分辨两种失败：一种是曲线一路平到底、第 2 轮就停——初始质心选得太好或者 k 太小；另一种是曲线反复小幅抖动几十轮——数据里有大量近似重复点，标签在两组等距质心之间来回跳。

2.3 先拿合成数据验一遍，再碰真实数据

真实数据出问题时，你分不清是代码错还是数据怪。先用两组结构已知的合成数据把管线跑通，是省时间的做法。

rng(7); % 场景一：三个各向同性高斯簇，k-means 的舒适区 mu = [0 0; 6 0; 3 6]; Xg = []; for i = 1:3 Xg = [Xg; mu(i,:) + 0.9 * randn(150, 2)]; end trueG = repelem((1:3)', 150); % 真实标签，用于算纯度 % 场景二：两个同心半圆，非凸结构，用来暴露 k-means 的失效边界 t = linspace(0, pi, 200)'; Xm = [[cos(t), sin(t)]; [2*cos(t), 2*sin(t)]] + 0.05 * randn(400, 2); trueM = repelem((1:2)', 200); [Lg, Cg, Jg] = kmeans_core(Xg, 3, 1); [Lm, Cm, Jm] = kmeans_core(Xm, 2, 1); fprintf('高斯数据 纯度 = %.3f，迭代 %d 轮，J = %.2f\n', ... purity(trueG, Lg), numel(Jg), Jg(end)); fprintf('半圆数据 纯度 = %.3f，迭代 %d 轮，J = %.2f\n', ... purity(trueM, Lm), numel(Jm), Jm(end)); function p = purity(trueLab, predLab) % 纯度：每个预测簇取占比最高的真实类别，按簇大小加权 u = unique(predLab); hit = 0; for i = 1:numel(u) hit = hit + max(accumarray(trueLab(predLab == u(i)), 1)); end p = hit / numel(trueLab); end

跑完你会看到，高斯数据上纯度接近 1.0，半圆数据上纯度只有 0.5 出头，两个半圆各被劈成一半再拼起来。这不是代码写错了，是 k-means 的目标函数本身就假设簇是凸的、各向同性的、大小相近的——它只会用直线（高维下是超平面）去切空间。报告里把这两组对照放进去，比只贴一张漂亮的三色散点图有说服力得多，因为它同时说明了方法的适用边界。

纯度只看“每个簇里最多的那类占多少”，簇数给多了它会虚高，所以还得配一个对簇数不敏感的指标，这点放到第 4 章讲。

3. 数据管线决定上限：缺失值、量纲与 PCA 降维前后的聚类差异

聚类结果不好，八成问题出在进kmeans_core之前。算法本身没有“理解”能力，它只认数值和距离，字段里有一个量级差三位的金额列，整个空间就被这一维主宰了。

3.1 缺失值、异常值与量纲的处理顺序

顺序很重要：先删列，再补缺，再压异常值，最后标准化。反过来做，先标准化会把缺失值算成NaN传染整列，先补缺再删列等于白补。

T = readtable('features.csv', 'VariableNamingRule', 'preserve'); % 只把数值列拿进特征矩阵，ID 和文本列留档 isNum = varfun(@isnumeric, T, 'OutputFormat', 'uniform'); ids = T(:, ~isNum); Xraw = table2array(T(:, isNum)); % 一、按缺失率删列 missRate = mean(ismissing(Xraw), 1); Xraw = Xraw(:, missRate < 0.30); % 二、补齐：列内线性插值，两端用最近值兜底 Xraw = fillmissing(Xraw, 'linear', 'EndValues', 'nearest'); % 三、MAD 稳健 z 分数压离群点，先标 NaN 再补 med = median(Xraw, 1); madz = 1.4826 * median(abs(Xraw - med), 1); % 1.4826 让 MAD 对齐正态标准差 Z = (Xraw - med) ./ max(madz, eps); % eps 防常数列为 0 Xraw(abs(Z) > 3.5) = NaN; Xraw = fillmissing(Xraw, 'nearest'); % 四、标准化：让每一维在距离里权重相同 X = zscore(Xraw);

1.4826这个常数不是拍脑袋来的，正态分布下 MAD 乘上它才等于标准差，所以 3.5 这个阈值和常见的 3σ 大致可比。max(madz, eps)是给常数列留活路——某一列全表一个值，MAD 是 0，不兜底会出Inf，后面全变NaN。

fillmissing的'EndValues','nearest'只在'linear'方法下有效，管首尾两端的NaN。数据量大、时间有序时，插值比填均值合理；如果行与行之间没有顺序关系（比如问卷），直接填中位数反而更老实。

3.2 标准化与 PCA 之后，聚类结果差多少

很多人以为 PCA 只是降维提速，实际上它会改变簇的形状。标准化把每一维拉到同一尺度，PCA 再旋转坐标系、按方差排序，前几个主成分往往集中了绝大部分方差——这时候欧氏距离又被前两三个主成分主导了，等于换了个方式重复量纲问题。所以 PCA 之后要不要再对主成分做白化，是个得实测的选择。

% 标准化后做 PCA，看累计方差 [coeff, score, latent] = pca(X); cumVar = cumsum(latent) / sum(latent); p = find(cumVar >= 0.85, 1); % 取到 85% 方差所需的主成分数 Xp = score(:, 1:p); % 三种输入下的聚类结果对比 [L1, ~] = kmeans_core(X, 3, 1); % 标准化原空间 [L2, ~] = kmeans_core(Xp, 3, 1); % PCA 前 p 维 Xw = zscore(score(:, 1:p)); % 主成分白化 [L3, ~] = kmeans_core(Xw, 3, 1); fprintf('原空间 vs PCA ARI = %.3f\n', ari(L1, L2)); fprintf('原空间 vs 白化PCA ARI = %.3f\n', ari(L1, L3)); fprintf('前三主成分解释 %.1f%% 方差\n', 100 * cumVar(p)); function a = ari(l1, l2) % 调整兰德指数：对标签编号置换不敏感，衡量两次划分的一致程度 n = numel(l1); [~, ~, i1] = unique(l1); [~, ~, i2] = unique(l2); M = accumarray([i1 i2], 1); nij = sum(M(:) .* (M(:) - 1) / 2); ai = sum(sum(M,2) .* (sum(M,2) - 1) / 2); bj = sum(sum(M,1)' .* (sum(M,1)' - 1) / 2); tot = n * (n - 1) / 2; expct = ai * bj / tot; mx = (ai + bj) / 2; a = (nij - expct) / (mx - expct + eps); end

ARI 取值大致落在 0 到 1 之间，1 表示两次划分完全一致（标签编号不同也算一致），0 附近表示跟随机分组差不多。这个指标的好处是不用真实标签也能用，特别适合写报告时比较“不同预处理方案下结论稳不稳”。如果你的报告里要对比三种方案，一张 ARI 矩阵比三张散点图直观。

白化这一步值得单独说：zscore(score(:,1:p))把每个主成分也拉到单位方差，等于告诉算法“第 1 主成分和第 5 主成分一样重要”。原始特征维度之间本来没有可比性，这么做通常合理；但如果某个主成分确实承载了核心业务信号，白化会把它稀释掉。我的习惯是两种都跑，ARI 差距小于 0.1 就用白化版，因为它的结果对随机种子更不敏感。

3.3 结果回写与出图：报告里该出现的三样东西

散点图只画前两个特征，但聚类用的是全部维度，图上看不出的东西得靠数字补。报告里至少要有：带质心的散点图、轮廓值分布图、以及一张把标签写回原始表的结果表。

% 距离矩阵用于取每个样本到自身质心的距离 D = sum(X.^2, 2) + sum(C1.^2, 2)' - 2 * (X * C1'); D = max(D, 0); [~, L1] = min(D, [], 2); figure('Position', [100 100 960 400]); subplot(1,2,1); gscatter(X(:,1), X(:,2), L1, lines(3), '.', 8); hold on; plot(C1(:,1), C1(:,2), 'kp', 'MarkerSize', 14, 'MarkerFaceColor', 'y'); title('前两维散点与质心'); hold off; subplot(1,2,2); silhouette(X, L1, 'Euclidean'); % 每个点的轮廓值，越接近 1 越贴合本簇 title('轮廓值分布'); % 标签和“到质心距离”一起写回，方便后续按簇下钻 T.cluster = L1; T.dist2centroid = sqrt(min(D, [], 2)); writetable(T, 'clustered_result.csv');

dist2centroid这一列常被忽略，实际很有用：同一个簇里距离 0.1 的点和距离 4.0 的点，业务含义完全不同。把它一起导出，后续做分群画像时可以直接取每簇最核心的样本看典型特征，也可以把距离最大的那批单独拉出来看是不是异常样本。

subplot里两张图并排，是因为散点图会让你误判——看起来分得很开的两个簇，在完整维度上可能轮廓值只有 0.2。轮廓值分布图才是判断这次聚类到底靠不靠谱的那张。

4. 簇数 K 怎么定：肘部法、轮廓系数与多次重启的实测对照

K 是 k-means 唯一的外生参数，也是最容易在报告里被追问的地方。常见的三种做法——肘部法、轮廓系数、Gap 统计量——各有各的失效场景，指望某一个自动给出正确答案，基本是玄学。稳妥的做法是三种都算，看它们是否指向同一个区间。

4.1 肘部法的读数陷阱与二阶差分

肘部法的逻辑是：K 增大时簇内平方和 J 单调下降，但下降速度会在真实簇数附近突然放缓，画出来像个手肘。问题在于“突然”这个判断太主观，曲线平滑的时候，十个人能看出五个肘部。

ks = 2:10; sse = zeros(size(ks)); for i = 1:numel(ks) [~, ~, Jh] = kmeans_core(X, ks(i), 42); % 固定种子，保证可比较 sse(i) = Jh(end); end d1 = diff(sse); % 一阶差分：J 的下降量 d2 = diff(sse, 2); % 二阶差分：下降量的变化率 [~, idx] = max(d2); % 曲率最大处就是肘部 kElbow = ks(idx + 1); fprintf('肘部法给出 K = %d\n', kElbow); plot(ks, sse, '-o'); xlabel('K'); ylabel('簇内平方和 J');

diff(sse,2)里第二个参数是差分阶数，意思是“下降速度本身的变化”。曲线上肘部明显时，二阶差分会在那里出现一个尖峰；曲线整体平滑时，这个尖峰只有零点几个单位，和噪声没区别，这时max挑出来的位置基本随机。所以我会在报告里同时贴原始曲线和二阶差分，如果尖峰不突出，就明说“肘部法在这里不可靠，改用轮廓系数”。

另外两个坑：一是ks从 2 开始而不是 1，K=1 时 J 必然最大，放进去会把整个曲线的纵向尺度拉爆，后面的变化被压平；二是每轮必须用同一个种子，否则 J 的抖动来自随机初始化而不是 K 的变化，曲线会呈锯齿状，看着像有一堆肘部。

4.2 轮廓系数：能给出绝对好坏，但对非凸簇天生苛刻

轮廓系数衡量的是“点离自己簇有多近、离最近的其他簇有多远”，单点取值在 -1 到 1 之间，取所有点的均值作为该 K 的得分。相比肘部法，它有个明显优势：给的是绝对评价，0.6 就是比 0.3 好，不依赖眼睛。

ks = 2:10; s = zeros(size(ks)); for i = 1:numel(ks) [L, ~] = kmeans_core(X, ks(i), 42); if numel(X) > 20000 % 大数据抽样，silhouette 是 O(n^2) idx = randperm(numel(L), 20000); s(i) = mean(silhouette(X(idx,:), L(idx))); else s(i) = mean(silhouette(X, L)); end end [bestS, bi] = max(s); fprintf('轮廓系数最优 K = %d，得分 %.3f\n', ks(bi), bestS);

silhouette内部要算点对之间的距离，样本量上万时内存和时间都会吃紧，抽样是常规处理，但要在报告里写清抽了多少。得分解读上，0.5 以上说明结构清晰，0.25 到 0.5 之间算“有结构但不强”，0.25 以下基本可以认为这组特征上不存在明显的簇结构——这时候不要硬凑 K，先回头检查特征选得对不对。

前面半圆数据那个例子，轮廓系数会在 K=2 时给出一个相当不错的分数，但聚类结果肉眼可见是错的。原因是轮廓系数同样基于欧氏距离，它衡量的是“球形簇”的紧致程度。遇到这种矛盾（轮廓系数高但业务上看不合理），说明数据分布本身不满足 k-means 的假设，该换密度聚类或者谱聚类，而不是继续调 K。

4.3 多次重启加 k-means++：把随机性压到可接受范围

同一个 K、不同的随机初值，J 能差出十几个百分点，这是 k-means 最让人不放心的地方。两个手段组合起来能把它压住。

function C = kmeanspp_init(X, k, seed) % k-means++：第一个质心随机取，之后每个质心按 D^2 的概率挑 rng(seed); n = size(X, 1); C = zeros(k, size(X, 2)); C(1, :) = X(randi(n), :); d2 = sum((X - C(1,:)).^2, 2); for j = 2:k p = d2 / sum(d2); idx = find(rand <= cumsum(p), 1, 'first'); % 不用 randsample，免工具箱 C(j, :) = X(idx, :); d2 = min(d2, sum((X - C(j,:)).^2, 2)); % 只保留更近的距离 end end % 多次重启，留最优 bestJ = inf; bestL = []; bestC = []; for r = 1:20 C0 = kmeanspp_init(X, 4, r); % 用 k-means++ 播种 [L, C, Jh] = kmeans_core_from(X, 4, C0, 300, 1e-6); if Jh(end) < bestJ bestJ = Jh(end); bestL = L; bestC = C; end end fprintf('20 次重启：最优 J = %.2f\n', bestJ);

kmeanspp_init里的关键在于p = d2 / sum(d2)：已经被选中质心覆盖过的区域，距离小、被选中的概率也小，新质心倾向于落在还没被代表的地方。这比纯随机播种好多少？在合成高斯数据上，纯随机 20 次重启的 J 极差通常在 5% 上下，k-means++ 能压到 1% 以内，迭代轮数也少一半。

kmeans_core_from是我把kmeans_core拆出来的一个变体——接受外部传入的初始质心，其余逻辑一致，代码这里省略。报告里值得记录的一个数字是：20 次重启中 J 的最大值、最小值和标准差。标准差大说明这次聚类本身不稳定，即使拿到最优解，结论也要打折扣。

还有一点，K 的扫描和重启次数是乘法关系。K 从 2 扫到 10、每次重启 20 遍，就是 180 次完整聚类。小数据没问题，几十万行时建议先抽 10% 做 K 的粗扫，选定 K 之后再用全量数据跑最后一次，并把这次的种子写进报告。

5. 避坑与排查：五条让 k-means 结果翻车的记录

下面这些是真实踩过、并且在别人的报告里反复见到的。每条按“现象、原因、解决”写，症状描述尽量具体到能对上号。

5.1 结果层面的三个坑

坑一：同一份数据跑两次，簇的颜色全换了。

现象是散点图的形状一模一样，但簇 1 变成了簇 2，质心坐标表对不上上一版。原因有两层：一是没调rng固定种子，初始化不同导致收敛到不同局部最优；二是即便 J 一样，k-means 输出的标签编号本身没有语义，第一次的第 1 簇可能就是第二次的第 3 簇。解决上，种子必须写进代码且记进报告；标签不要直接比对编号，用质心坐标排序后再对齐——按第一主成分或者按簇大小从大到小重排，或者在报告里直接贴质心坐标表，别人一看就懂。真要严谨，用第 3 章的 ARI 比对两次划分。

坑二：K 从 2 加到 8，J 一路下降，报告里写“K=8 最优”。

这是最典型的误读。J 对 K 单调不增，K 等于样本数时 J 恰好是 0。光看 J 挑 K，永远会挑到最大的那个。解决是加约束：轮廓系数取最大、肘部法看拐点、或者直接按业务可解释性定上限。我一般把 K 的搜索范围限制在 2 到 $\sqrt{n/2}$ 之间，超过这个量级，簇的平均样本数已经少到做画像没有统计意义了。

坑三：某一列金额把距离全吃掉了，结果只按金额分群。

现象是分出来的簇在其他特征上几乎没差别，唯独金额从低到高排得整整齐齐。原因是该列方差比其他列大三个数量级，平方之后在距离里占了 99% 的权重。解决就是标准化，但要注意标准化必须用同一套均值方差——训练集算出来的均值方差应用到新数据上，不能每次重算。如果该列本身是长尾分布（多数人在低位、极少数极高），zscore之后仍会被极值拉扯，先取对数再标准化更稳。

5.2 实现层面的两个坑

坑四：样本量一大就内存溢出，任务被系统杀掉。

现象是十万行、五十维的数据，跑到距离计算那一步 MATLAB 直接卡死或者报内存不足。原因是有人为了“算清楚”，先构造了n×n的距离矩阵（十万乘十万就是八十亿个 double）。k-means 全程只需要n×k的距离，k 通常不超过 20，差了五个数量级。解决是坚持用第 2 章的展开式，任何时候都不要在 k-means 里出现pdist生成的全距离矩阵；如果确实要用pdist2算到质心的距离，也要显式传质心矩阵而不是整个数据集。

坑五：结果表写回之后行序错乱，标签和样本对不上。

现象是导出 CSV 打开一看，某个 ID 对应的簇号跟图上的位置明显不符。原因是在预处理阶段用Xraw(completeIdx,:)删过缺失行，但写回时用的是原始表T，两者行数不同，MATLAB 按位置赋值，错位就发生了。解决的铁律是全程保留一个唯一 ID 列，所有过滤、抽样都只操作索引，最后用 ID 做join合并标签，永远不要依赖行号对齐。这条在跨步骤的脚本里尤其重要，一个中间步骤改了行序，后面所有分析都是错的。

6. 把实验报告固化成脚本：种子、日志与一键复现的验证清单

到这一步，代码能跑、图能出、K 有依据了，但离“交出去别人能复现”还差一层。我的习惯是把整份实验拆成四个文件，用主脚本串起来，任何一个环节改了都能重跑。

分开的好处很直接：导师问“换个标准化方式会怎样”，你只改prep_data.m重跑，其余三个文件一行不动，对比实验十分钟出结果。全塞在一个脚本里，改一处就得从头跑，跑到一半忘了哪步改过。

验证清单我固定查四项。第一，连跑两次，用ari(L_first, L_second)比对，结果必须是 1.000，不是就说明某个环节的种子没固定住。第二，把select_k.m的输出存成k_selection.mat，里面除了 K 还要存扫描曲线和轮廓得分，重跑时先比对这两条曲线是否逐点相同。第三，日志里记下数据文件的行数、列名、每列的缺失率和标准化后的均值方差，这些数字比“数据已预处理”这句话有用得多。第四，最终结果表的行数必须等于原始表的行数，用height(T) == height(Tout)断言一次，能挡住第 5 章坑五那类错位。

顺带说一个提速的技巧：K 扫描阶段用抽样，但不要把抽样的结论直接当最终结果。我的做法是先抽 10% 跑一遍拿到 K 的候选区间，再在全量数据上只跑区间内的两三个 K，最后只对选定的 K 做完整的 20 次重启。这样总耗时约为全量扫描的三分之一，结论基本一致。

血泪经验是：报告里最值钱的不是那张漂亮的散点图，而是能让人把你的脚本下载下来、改个数据路径、跑出同样数字的那套东西。我早期交报告只贴图和结论，被问过一次“你这个 K 是试出来的还是算出来的”，答不上来，后来才养成把种子、日志、断言都写进代码的习惯。这份后悔药不难吃，只是要在一开始就想到。希望帮到你。

本文还有配套的精品资源，点击获取