简介一份面向全国大学生数学建模竞赛备赛者的完整解题文档聚焦2022年C题古代玻璃制品的成分分析与鉴别。内容系统覆盖数据预处理、缺失值与异常值处理、特征工程、主成分分析降维、因子分析构建转换矩阵、敏感性分析验证鲁棒性以及灰色关联度分析比较类别差异完整呈现四个子问题的建模流程、公式推导和结果解释可作为赛前训练、方法对比和论文写作的实用参考尤其适合需要从数据探索到模型评价全流程参考的参赛者。压缩包内唯一文件为PDF格式整体大小约3.55MB包含摘要、问题重述、数据探索、模型构建与评价等章节结构清晰层级分明便于按模块研读。目前已有927人学习适合高校数学建模参赛队伍及对玻璃文物成分分析感兴趣的研究者。1. 这道C题在考什么从数据到结论的完整链路2022年国赛C题“古代玻璃制品的成分分析与鉴别”表面上是分类题按化学成分把未知玻璃样品分到高钾或铅钡亚类。但真正让队伍翻车的不是分类精度而是“风化”这两个字——风化样品的成分已经失真直接拿去聚类欧氏距离会把风化差异误当成类间差异亚类全串位。这道题值不值得投入取决于你的目标如果只求一篇交差的论文抄个分类模型就够了如果想让聚类、差异检验、判别模型这条链路成为以后参赛都能复用的模板按这篇笔记走会省下大量试错时间。适合正在备赛、熟悉 pandas 基本操作、想用 sklearn 解决真实题目的同学。2. 先把玻璃成分数据洗干净缺失值、百分化和那口必须填的坑2.1 理解题给数据的结构两个表单、两套指标、三种坑2022年C题的原始数据常见习惯分成两份表单一份是样品信息表记录样品编号、纹饰类型、颜色、是否风化另一份是化学成分表每一列是一种氧化物一行是一个样品含量以百分比形式给出。两份表单靠“样品编号”关联这个关联字段通常是字符串不是默认的数字索引第一步就该确认它有没有重复、有没有空值。实际数据分析里成分数据有三个非常固定的坑。第一个坑是闭合效应每行百分比加总应该是100但题给数据常常加总成98、102甚至更离谱这是因为检测手段测不到所有组分空白值被记成了空。第二个坑是全空列个别氧化物在本批样品里一个都没测出来这一列直接是NaN留着会污染标准化。第三个坑是半定量数据有些低含量组分比如SnO2、P2O5大量为0或空直接log变换会得到负无穷。处理顺序上我一般会先看形状和缺失分布再决定“删列、填值、变换”的具体操作而不是一上来就StandardScaler。标准化解决的是量纲解决不了闭合效应这一步顺序颠倒后面聚类和判别都会出问题。2.2 成分数据为什么要做变换从欧氏距离说到中心对数比成分数据本质上是“单纯形”上的点各组分总和固定为100%一个组分升高其他组分必然被压缩。在这种约束下原始百分比空间的欧氏距离并没有几何意义。最典型的现象是SiO2在很多玻璃样品里占比60%以上PbO、BaO在铅钡玻璃里各占十几到三十个百分点原始空间里两点之间的距离几乎被SiO2一个变量支配其他元素的差异全部淹没。解决这个问题最常用的变换是中心对数比变换centered log-ratioCLR。做法是先对每个样品的所有组分求几何平均再把每个组分除以这个几何平均后取对数。变换之后的数据不再受“总和为100”的约束变量之间依然保留成分的相对关系而且正值负值都有可以直接丢进K-Means、PCA、LDA这些依赖欧氏距离或协方差的算法。CLR有个前提是组分必须严格大于0log(0)没有定义。所以做CLR之前必须填充质控样本里出现0的组分常见做法是用这一列非零最小值的一半去替代整列全空的组分直接删掉。注意填充值是“假的”对低含量组分的CLR影响很大后面做敏感性分析时要把填充策略换一种再跑一遍确认结论没变这一步很多人会漏。注意CLR变换必须在填充之后进行几何平均是逐行算的一行里只要有一个0整行都会变成NaN。2.3 数据清洗的最小可复现代码读表、缺失值、CLR三步走下面这段代码覆盖读入、合并、删列、填充、CLR变换五个动作是这套题最稳定的起点。import pandas as pd import numpy as np from scipy.stats import gmean # 读入两份表单sheet名和列名以实际附件为准 meta pd.read_excel(附件1_样品信息.xlsx) raw pd.read_excel(附件2_化学成分.xlsx) print(样品编号重复数, raw[样品编号].duplicated().sum()) print(空值分布) print(raw.isna().sum().sort_values(ascendingFalse).head(10)) # 剔除全空列 comp_cols [c for c in raw.columns if c ! 样品编号] comp_cols [c for c in comp_cols if raw[c].notna().sum() 0] print(保留的成分列数, len(comp_cols)) # 用非零最小值的一半填充空值和0 for c in comp_cols: col raw[c].replace(0, np.nan) fill col.min() / 2 raw[c] raw[c].fillna(fill) # 合并样品信息成分列此时已是填充后的干净数据 df raw.merge( meta[[样品编号, 纹饰, 颜色, 风化与否]], on样品编号, howleft ) # 中心对数比变换 def clr_transform(X): # X是二维numpy数组每行为一个样品 gm gmean(X, axis1) return np.log(X / gm[:, None]) X_clr clr_transform(raw[comp_cols].to_numpy(dtypefloat)) clr_df pd.DataFrame(X_clr, columnscomp_cols, indexraw[样品编号])这里的逻辑要点先删除全空列避免变换时连分母都是空再按列填充填充值统一取非零最小值的一半比拍脑袋填0.0001更能贴合检测下限最后把样品信息合并进来让每个样品同时拥有工艺标签和化学成分。参数上要注意三点一是raw[c].replace(0, np.nan)必须放在求最小值之前否则0会把最小值拉成0填出来的全是0二是gmean遇到0会直接报错所以填充必须在变换前完成三是合并用howleft以成分表为准信息表里多出来的样品不参与建模。2.4 清洗前后模型的差距用一个快速聚类感知数据质量拿清洗前和清洗后的数据分别做一次K2的K-Means再把聚类结果和“高钾/铅钡”这个已知标签对比你能直观看到变换的影响。from sklearn.cluster import KMeans from sklearn.metrics import adjusted_rand_score # 原始百分比直接聚类此时raw已经填过缺失 raw_filled raw[comp_cols].to_numpy(dtypefloat) labels meta.set_index(样品编号).loc[raw[样品编号], 类型].values for name, X in [(原始百分比, raw_filled), (CLR变换后, X_clr)]: km KMeans(n_clusters2, n_init25, random_state42).fit(X) ari adjusted_rand_score(labels, km.labels_) print(f{name}: ARI{ari:.3f})调整兰德指数ARI衡量聚类与真实类别的一致性1是完全吻合0是随机水平。常见结果是原始百分比空间里ARI也不低因为SiO2的绝对含量本身能区分高钾与铅钡但CLR变换后聚类不再被SiO2支配亚类之间的差异才真正浮出来。判断数据质量的要点不是看ARI谁高而是看两次聚类得到的簇中心原始空间的簇中心在SiO2上相差二十个百分点CLR空间的簇中心则在多种中低含量组分上拉开差距这才是这道题需要的分辨力。3. 玻璃亚分类用聚类找到题目里的隐藏结构3.1 亚类划分的共识做法先聚类后贴标签题目给的样品信息表里有纹饰、颜色这些工艺属性很多人会直接拿它们当亚类标签去训练分类器这是理解偏差。工艺属性反映的是制作风格化学配方才是分类的依据同一种纹饰完全可能对应两套不同的配方。这道题的常见做法是先在高钾玻璃、铅钡玻璃各自内部做无监督聚类得到化学成分上的自然分群再把分群结果和纹饰、颜色做交叉验证看工艺属性是否与配方对应。对应得好说明纹饰可以作为亚类标签对应不上就以聚类结果为准。聚类前还有一步很多人会跳过把样品按高钾、铅钡拆开。因为两类玻璃的配方完全不同PbO、K2O的含量差异本身就是最大的分类信号混在一起聚类算法只会优先区分这两个大类亚类结构彻底被掩盖。分而治之以后每类内部的数据量都不大聚类参数要相应调小后面会说。3.2 选K的两种可靠策略肘部法则与轮廓系数亚类数K不能拍脑袋定常见做法是用肘部法则先确定范围再用轮廓系数做二次确认。肘部法则看簇内误差平方和WCSSK增大WCSS必然下降但下降到某个点之后降幅变缓这个拐点就是“再增加一个簇收益不大”的位置。轮廓系数则是把每个点与自身簇的紧密度、与最近簇的分离度结合起来算取值在-1到1之间越大越好。小样本数据上两条规则都会失灵。我见过的高钾玻璃样品数通常只有十几件K取到5以上时每个簇只剩两三个点轮廓系数虚高WCSS曲线没有明显拐点。经验上K不要超过样本数的三分之一一般先试2到4铅钡玻璃样本多一些可以试到6。两条规则不一致时优先照顾业务可解释性选出来的K要让每个簇能用“高钾低钾、低铅高铅”这类配方差异描述出来而不是纯统计最优。3.3 K-Means亚分类完整代码变换后的聚类与可视化import matplotlib.pyplot as plt from sklearn.cluster import KMeans from sklearn.decomposition import PCA # 按大类拆分part_id是样品编号 for glass_type, max_k in [(高钾, 3), (铅钡, 5)]: part_id df.loc[df[类型] glass_type, 样品编号] X clr_df.loc[part_id].to_numpy(dtypefloat) # 肘部法则 wcss [] for k in range(2, max_k 1): km KMeans(n_clustersk, n_init25, random_state42).fit(X) wcss.append(km.inertia_) print(glass_type, WCSS:, [round(v, 2) for v in wcss]) # 按拐点确定K示例中高钾取3、铅钡取5以实际曲线调整 k_final 3 if glass_type 高钾 else 5 km KMeans(n_clustersk_final, n_init25, random_state42).fit(X) clr_df.loc[part_id, 亚类] km.labels_ # PCA压缩到2维做可视化 pca PCA(n_components2).fit(X) P pca.transform(X) plt.figure(figsize(6, 5)) plt.scatter(P[:, 0], P[:, 1], ckm.labels_, cmapSet2, alpha0.8) plt.title(f{glass_type} 亚类分布PCA投影) plt.xlabel(PC1) plt.ylabel(PC2) plt.show() # 把亚类标签合并回主表后续建模直接用 df df.merge(clr_df[[亚类]].reset_index(), on样品编号, howleft)这段代码把两个大类分开聚类并输出WCSS序列和PCA投影图。需要注意n_init参数K-Means的初始中心是随机的n_init25表示跑25次取最优避免局部最优random_state固定下来保证论文里写的复现结果一致。PCA只用于可视化不参与聚类如果发现投影图里某个簇完全被另一个簇吞掉说明这个K取大了回到WCSS重新选。WCSS输出的数值本身没有绝对意义看的是相邻两个K之间的降幅。3.4 把聚类结果与题目给出的纹饰、颜色分组合并聚类做完下一步是把“亚类”和工艺属性做交叉表验证分群是否对应用途。part_id df.loc[df[类型] 高钾, 样品编号] sub clr_df.loc[part_id] ct pd.crosstab( sub[亚类], meta.set_index(样品编号).loc[part_id, 纹饰] ) print(ct) print(ct.div(ct.sum(axis1), axis0).round(2)) # 按行归一化看比例交叉表按行归一化之后如果某一簇在某类纹饰上占比超过80%基本可以认定这个亚类对应一种特定的工艺配方若分布完全平均说明纹饰不影响配方论文里就不要强行把纹饰写进亚类差异分析。这里要特别提醒交叉表里的“亚类”是0、1、2这样的数字没有业务含义后续差异分析时要给每个簇起个描述性名字比如“高钾-低钙型”否则写着写着两个数字就不知道在指代什么了。4. 风化对化学成分的影响统计检验与可解释性4.1 风化不是“少了一些元素”那么简单问题2把整道题的难度拉起来了。风化玻璃在埋藏过程中长期接触土壤和地下水易溶组分流失难溶组分相对富集还可能吸附外来污染物所以“风化后成分更低”只是粗糙印象不同元素的行为差异很大。以常见的铅钡玻璃为例铅、钡在风化中容易流失SiO2在残存骨架里占比反而整体抬升高钾玻璃里K2O也容易被淋失。这套题的关键不是让模型记住“风化减少”而是找出哪些组分的流失在统计上显著、显著到足以改变后续亚类判别。分析设计上要注意性质不同的玻璃配方不同风化机理也不同绝对不能把高钾和铅钡混在一起做风化检验。常见做法是按大类分组再在每组内部比较风化与未风化样品的成分差异。另外风化状态只有“是/否”两个水平样本量又不大优先考虑两独立样本的差异检验而不是方差分析。4.2 检验方法选型正态性检查决定t检验还是Mann-Whitney两独立样本比较第一反应是t检验但成分数据大量不满足正态性百分比被压缩在0到100之间低含量组分普遍右偏。小样本下中心极限定理也救不了必须先做正态性检查。常见做法是对风化、未风化两组分别做Shapiro-Wilk检验p值小于0.05就认为不满足正态改用Mann-Whitney U检验。U检验是非参数的只比较分布位置对偏态和异常值稳健这道题里是应用频率最高的一种检验。还有一个常被忽视的点多重比较。成分列有十几个每个都做一次检验假阳性自然上升。常见的修正方案是Benjamini-HochbergBH方法控制错误发现率把p值排序后统一修正比Bonferroni宽松适合这种探索性分析。论文里写“经过BH校正后仍显著的组分”审阅印象会明显不一样。4.3 差异分析的完整代码检验表、效应量与箱线图from scipy import stats from statsmodels.stats.multitest import multipletests results [] for c in comp_cols: a df.loc[(df[类型] 高钾) (df[风化与否] 风化), c].dropna() b df.loc[(df[类型] 高钾) (df[风化与否] 未风化), c].dropna() if len(a) 3 or len(b) 3: continue # 正态性粗查 _, pa stats.shapiro(a) _, pb stats.shapiro(b) if pa 0.05 and pb 0.05: stat, p stats.ttest_ind(a, b, equal_varFalse) method Welch else: stat, p stats.mannwhitneyu(a, b, alternativetwo-sided) method Mann-Whitney results.append([c, method, stat, p, a.median(), b.median()]) res pd.DataFrame( results, columns[成分, 方法, 统计量, p值, 风化中位数, 未风化中位数] ) res[BH校正p] multipletests(res[p值], methodfdr_bh)[1] res res.sort_values(p值) print(res.to_string(indexFalse))逻辑说明先按大类切片分组做检验Shapiro p值都大于0.05才用Welch t检验否则转Mann-Whitney统计量本身不重要p值、方向、效应量才重要。参数选择上ttest_ind(equal_varFalse)是Welch校正版不假设方差齐性mannwhitneyu(alternativetwo-sided)在样本量小的时候自动走精确法样本量大时自动近似不需要手动切换multipletests的fdr_bh就是BH校正。跑完拿到结果后先去重看“风化中位数/未风化中位数”两列的方向和比例再用BH校正后的p值筛选显著组分。4.4 从“p0.05”到“结论可写进论文”的差异方向p值告诉你有差异但论文里还要说清楚差异有多大、方向是什么。我一般把结果整理成一张表组分、检验方法、风化组中位数、未风化组中位数、变化方向、BH校正后p值。中位数比均值更适合偏态数据看表格时一眼能读出“PbO中位数从18降到12”这种可解释结论。若要量化差异幅度可以补一个效应量Cliffs delta取值0.47以上被认为是中到大效应比单纯写“显著降低”有力得多。这张差异表不只为写论文还直接影响后面的判别建模风化显著影响的组分在风化样品里已经“失真”把它们作为亚类判别特征时模型学到的是“风化效应”而不是“配方差异”。所以一个值得做的中间步骤是分别用所有组分、剔除风化显著组分后的两组数据训练判别模型对比在未风化样本上的表现看特征选择是否合理。这一步做好了问题4里“为什么有些样品判别错误”的原因分析就有了数据支撑。5. 避坑判别模型里 4 个最容易翻车的地方5.1 坑1把“风化状态”当特征模型精度虚高现象训练集和测试集accuracy都超过0.96看起来完美收官。细看特征列表发现把“风化与否”编码成了0/1列喂进了模型。原因风化状态在本题里是已知的工艺背景不是待预测属性模型学会了“风化样品大概率属于某个亚类”的捷径而问题3实际要预测的未知样品根本没有风化标注模型在真实场景里等同于废掉。解决建模时只保留化学成分列工艺属性最多用来做事后分组对比不进特征矩阵。检查方法很简单——打印特征列名列表逐列过一遍凡是含义里带“是否”这种二值语义的都要问自己一句未知样品这一列能不能拿到。5.2 坑2总和归一化后直接丢进K-Means聚类跟着主成分跑现象对百分比数据做K-Means轮廓系数很高簇中心却只在SiO2这一个变量上拉开差距和铅钡、高钾已知标签对不上。原因百分比是闭合数据SiO2平均含量显著高于其他组分欧氏距离几乎被它垄断K-Means用的是欧氏距离自然把“SiO2高不高”当成唯一分组依据。解决先CLR变换再考虑是否降维最后聚类。CLR消除了闭合效应让每个组分的相对丰度都参与距离计算。检查方法输出每个簇的簇中心看差异是否分散在多个组分上而不是集中在一个变量上。5.3 坑3QDA/LDA报“协方差矩阵奇异”现象训练QDA或LDA时sklearn抛出LinAlgError或者给出警告说矩阵逆不存在个别情况下不报错但预测结果全是同一类。原因成分变量通常十几个而高钾玻璃样本只有十几件组内协方差矩阵的自由度不够求逆时矩阵奇异。解决两个思路一是先PCA降维到5到6维保留95%方差后再进判别模型二是换用带收缩的线性判别如LinearDiscriminantAnalysis(solverlsqr)或直接用正则化逻辑回归它们对奇异协方差不敏感。这个坑在论文里可以转化为“降维必要性”的论述说明你考虑了维度灾难。5.4 坑4标准化参数跨样本计算数据泄漏现象交叉验证的平均准确率0.94独立测试集上只有0.78调了半天模型没改善。原因前面做CLR、StandardScaler时用了全量数据含测试集去算几何均值和标准差测试集的信息已经渗透进训练过程这种“测试集可见”的泄漏会让交叉验证虚高而真实样本上原形毕露。解决把变换放进Pipeline让每一折训练集单独fit变换器、单独transform验证集。下面这个例子把CLR、PCA、LDA串成管道配合留一交叉验证是这个问题里最稳妥的标准写法。import numpy as np from sklearn.pipeline import Pipeline from sklearn.preprocessing import FunctionTransformer, StandardScaler from sklearn.decomposition import PCA from sklearn.discriminant_analysis import LinearDiscriminantAnalysis from sklearn.model_selection import LeaveOneOut, cross_val_score # clr_transform沿用第2章定义 pipeline Pipeline([ (clr, FunctionTransformer(clr_transform, feature_names_outone-to-one)), (scale, StandardScaler()), (pca, PCA(n_components8)), (lda, LinearDiscriminantAnalysis()), ]) X df[comp_cols].to_numpy(dtypefloat) y df[亚类].astype(str).values # 用聚类结果当标签 loo LeaveOneOut() scores cross_val_score(pipeline, X, y, cvloo) print(留一法准确率, scores.mean())这段代码把五个步骤封死在管道里交叉验证时每一折先fit管道再predictCLR的几何均值、StandardScaler的均值和方差、PCA的主成分方向全部只来自训练折泄漏通道全部堵死。参数上n_components8按经验折中若样本量紧张可降到6LinearDiscriminantAnalysis()默认solver在类数与特征数合适时不会奇异如果仍报错就换成solverlsqr。只看准确率不够建议把每一折被分错的样品编号打印出来它们很可能就是问题4要求的“误判样品分析”原料。6. 判别模型的一个进阶验证技巧判别平面上的留一法检查6.1 把留一法的误判点画到判别平面上留一法交叉验证给出一个准确率数值但数值掩盖了谁是“易错样本”。我习惯把留一法的预测结果画到判别平面上用训练数据拟合LDA并投影到两个判别轴每个点按真实亚类着色预测错误的点单独用叉号圈出。一张图能同时回答三个问题误判集中在哪两个亚类之间、误判样本是否贴着决策边界、有没有个别点被完全孤立。import numpy as np import matplotlib.pyplot as plt from sklearn.pipeline import Pipeline from sklearn.preprocessing import FunctionTransformer, StandardScaler from sklearn.decomposition import PCA from sklearn.discriminant_analysis import LinearDiscriminantAnalysis from sklearn.model_selection import LeaveOneOut pipeline Pipeline([ (clr, FunctionTransformer(clr_transform, feature_names_outone-to-one)), (scale, StandardScaler()), (pca, PCA(n_components8)), (lda, LinearDiscriminantAnalysis()), ]) X df[comp_cols].to_numpy(dtypefloat) y df[亚类].astype(str).values loo LeaveOneOut() preds, true [], [] for tr, te in loo.split(X, y): pipeline.fit(X[tr], y[tr]) preds.append(pipeline.predict(X[te])[0]) true.append(y[te][0]) # 用去掉LDA的管道做投影得到二维判别平面坐标 proj Pipeline(pipeline.steps[:-1]).fit_transform(X, y) plt.figure(figsize(7, 6)) for t in set(true): mask np.array(true) t plt.scatter(proj[mask, 0], proj[mask, 1], labelt, alpha0.8, s30) err np.array(true) ! np.array(preds) plt.scatter(proj[err, 0], proj[err, 1], facecolornone, edgecolorblack, linewidth2, s120, label误判样本) plt.legend() plt.show()这里的preds和true是留一法逐样本预测的结果proj是管道去掉最后一步LDA后得到的前置投影坐标err掩码标记哪些点预测错误scatter里s控制点大小linewidth控制误判样本圈线的粗细两个参数只影响观感不影响结论。6.2 误判样本和风化差异表配合看误判点如果贴着决策边界说明这个样本本身的成分就处在两种配方的过渡带模型判断摇摆可以解释如果误判点远离边界却被分错多半是风化样品——它们的真实配方已经被淋失改变拿未风化样本的信息去猜自然不准。这时候回去对照第4章的风化差异表看该样品里风化敏感组分是否异常问题4的“误判原因分析”就有了一段可复现的因果证据而不是笼统的“数据噪声”。我第一次做这道题时只盯准确率把误判压给“数据噪声”四个字后来把误判样本逐个投影回判别平面才发现错误大多来自风化组样品这才回头补了风化差异分析。现在这套流程已经成了我做成分类题目的固定动作清洗、CLR、聚类、差异检验、判别管道、判别平面验尸。希望帮到你。本文还有配套的精品资源点击获取
