资讯详情 2022数模C题玻璃成分分析实战:成分数据与小样本建模避坑指南
📅 2026/10/11 17:47:15
简介面向数学建模竞赛参赛者和数据分析学习者提供2022年全国大学生数学建模C题完整解题方案。内容围绕古代玻璃文物成分分析与鉴别系统涵盖数据预处理与探索、特征降维、因子分析构建风化前后成分转换、亚类划分、灰色关联度分析及敏感性检验并配有对应实现代码可直接对照论文理解模型思路与计算细节。资源为单份PDF文件容量3.55MB兼具论文阐述与代码展示便于按章节研读和复现。目前已有927人浏览学习适合正在备战国赛或练习分类预测与降维建模的选手作为范例参考。1. 2022全国大学生数学建模C题在做什么一道“成分数据”“小样本”的复合题2022年全国大学生数学建模竞赛C题赛题名称是《古代玻璃制品的成分分析与鉴别》。表面上是考古学场景实际是一个典型的小样本成分数据挖掘问题给你一张几十行的玻璃成分表SiO₂、PbO、BaO等十几种氧化物的百分比先按类型分成铅钡玻璃和高钾玻璃再在每类里找亚类最后还要预测“风化前”的化学成分。拿奖的关键不在用多深的模型而在数据清洗、成分数据变换和防数据泄漏这三件事做到位。适合谁想冲国赛奖的建模组以及需要快速上手“化学组分表格”这类结构化数据的同学。2. 拿到C题数据先做三件事表对齐、缺失值处理和成分特征构造2.1 附件1和附件2怎么关联以“采样编号”为主键的唯一合并C题的数据由两部分组成。一部分是成分实测表记录了每个采样点的氧化物百分比和“是否风化”标记另一部分是样品外观信息表包含颜色、纹饰等字段。两表之间靠“采样编号”关联。拿到数据的第一步不是建模而是确认两个表里有多少编号能对齐、有多少孤儿行。常见翻车是直接用pandas默认的inner join合并完才发现样本数少了一截后面的聚类和回归全在残缺数据上跑。import pandas as pd # 附件1: 成分实测, 附件2: 颜色/纹饰等外观信息 comp pd.read_excel(附件1.xlsx) info pd.read_excel(附件2.xlsx) # 先看两个表各自的形状和主键 print(成分表:, comp.shape, comp[采样编号].nunique()) print(信息表:, info.shape, info[采样编号].nunique()) # 左连接保留成分表的全部样本, 检查哪些编号没匹配上 df comp.merge(info, on采样编号, howleft) not_matched df[df[颜色].isna()][采样编号].tolist() print(未匹配到附件2的编号:, not_matched)这里用howleft而不是howinner目的是暴露对齐问题而不是隐藏它。左连接后如果发现某行颜色全为空就该回头查附件1里是否存在重复编号或不可见字符。两个表的主键都是“采样编号”但实际数据里出现过编号带空格或全角半角混用导致merge失败的情况排查时先对编号做strip()和str.upper()统一样式。对齐确认后再进入清洗流程。2.2 缺失值不是闹着玩的成分数据缺失的两种处理姿势成分表的缺失一般分两种。一种是某个成分在所有样本中都没测到这种直接删列另一种是部分样本缺失常见原因是检测仪器能力下限或样品量不够。千万不要对所有缺失列一律填0因为化学成分百分比在后续计算几何平均时会直接导致对数变换失败。我自己打比赛时的习惯是先看缺失率矩阵再做逐列处理。chem_cols [SiO2, Na2O, K2O, CaO, MgO, Al2O3, Fe2O3, CuO, PbO, BaO, P2O5, SrO, SnO2, SO2] # 缺失率统计 miss df[chem_cols].isna().mean().sort_values(ascendingFalse) print(miss[miss 0]) # 缺失率超过50%的列直接删, 低于20%的用该类型(高钾/铅钡)组内中位数填充 drop_cols miss[miss 0.5].index.tolist() fill_cols miss[(miss 0) (miss 0.5)].index.tolist() print(删列:, drop_cols, 填充列:, fill_cols) for c in fill_cols: df[c] df.groupby(类型)[c].transform(lambda s: s.fillna(s.median()))为什么不全局填中位数因为铅钡玻璃的PbO动辄百分之三四十高钾玻璃的PbO可能只有个位数全局填值会严重扭曲两类样本的判别特征。按“类型”分组填充保留的是类别内部的典型水平。如果某列缺失率在20%~50%之间填充后还要在论文的“数据处理”一节明确交代填充规则评委给分时会看这个细节。对个别样本多个成分同时缺失的情况另一个做法是直接删行但C题样本量本来就小删行要谨慎优先保留关键成分完整的行。2.3 百分比成分的“伪相关”为什么不做标准化而是做CLR变换成分数据有一个让新手初赛翻车的特性每行的所有氧化物百分比加起来接近100%这是一个“闭合”约束导致变量之间存在天然的负相关。直接对成分列做z-score标准化只能处理量纲不能消除闭合效应。两个本就可比的变量在PCA、聚类和回归里会产生假信号。常规解法是CLR中心化对数比变换这是处理成分数据的标准动作。import numpy as np def clr_transform(df_in, cols): 对指定的成分列做中心化对数比变换。 先把缺失和0值替换为一个小正数(0.001), 再对每行除以几何平均后取对数。 X df_in[cols].copy() X X.replace(0, np.nan).fillna(0.001) gm np.exp(X.apply(np.log).mean(axis1)) # 每行几何平均 clr X.apply(np.log).sub(np.log(gm), axis0) return clr X_clr clr_transform(df, chem_cols) print(X_clr.head())CLR把闭合的成分矩阵变成不受和约束影响的欧氏空间数据。变换后的每个值是“该成分相对该样本所有成分平均水平的偏离”正负号天然可解释正代表该样本里此成分偏高。对C题来说后续无论是聚类找亚类、LDA判别还是XGBoost建模都在CLR空间里做。唯一需要说明的是CLR会引入一行中所有列的线性相关每行之和为0这会让普通协方差矩阵不可逆所以搭配正则化方法Ridge、带收缩的LDA不要用普通线性回归。提示如果某行缺失值太多比如超过一半成分缺失填充后CLR会失真这种样本建议直接删行并在论文里写明删了几条、为什么删。3. 分类与亚类划分先聚类再判别的两步法3.1 用KMeans找亚类K值怎么定C题第一问要求在铅钡玻璃和高钾玻璃内部找出亚类。这里的“亚类”没有标签只能是聚类。常见做法是分别取两类样本做KMeans用轮廓系数和手肘图选K。样本量小K一般取2~4不要贪多。我一般会固定随机种子跑否则每次结果不一样论文里没法交代。from sklearn.cluster import KMeans from sklearn.metrics import silhouette_score def kmeans_k_search(X, max_k5): scores {} for k in range(2, max_k 1): km KMeans(n_clustersk, random_state42, n_init10) labels km.fit_predict(X) scores[k] silhouette_score(X, labels) return scores # 注意: X必须是CLR变换后的矩阵 pb_clr X_clr[df[类型] 铅钡玻璃] kj_clr X_clr[df[类型] 高钾玻璃] print(铅钡玻璃轮廓系数:, kmeans_k_search(pb_clr)) print(高钾玻璃轮廓系数:, kmeans_k_search(kj_clr))轮廓系数越大说明类内越紧凑、类间越分离。但小样本下轮廓系数会偏高且不平滑我的做法是把轮廓系数和业务常识结合比如铅钡玻璃按PbO和BaO的高低分离出一类“高铅钡”、一类“低铅钡”这比纯统计选K更容易在论文里讲圆。KMeans的n_init10是克服小样本下初始中心随机导致的波动random_state42保证复核时结果一致。聚类到的亚类标签要合并回原数据表作为后面的判别模型的Y。3.2 亚类判别LDA与XGBoost的小样本参数拿到亚类标签后第一问还要验证亚类之间是否真的可分。这时候可以用LDA线性判别分析做一个基线判别模型。LDA在样本量小到几十行的场景下非常稳用它算出来的交叉验证准确率本身就是“亚类差异显著”的直接证据。from sklearn.discriminant_analysis import LinearDiscriminantAnalysis from sklearn.model_selection import cross_val_score, StratifiedKFold lda LinearDiscriminantAnalysis(solverlsqr, shrinkageauto) cv StratifiedKFold(n_splits5, shuffleTrue, random_state42) scores cross_val_score(lda, X_clr, y_subclass, cvcv, scoringaccuracy) print(LDA 5折平均准确率: %.3f (/- %.3f) % (scores.mean(), scores.std()))LDA用lsqrshrinkageauto而不是默认的svd是因为CLR变换后特征矩阵每行和为0协方差矩阵奇异收缩估计能稳定求解。如果LDA准确率已经到90%以上后续用XGBoost只是锦上添花。XGBoost在小样本上的优势是能捕捉非线性交互劣势是极易过拟合需要把树压浅。import xgboost as xgb xgb_model xgb.XGBClassifier( n_estimators120, max_depth2, learning_rate0.05, subsample0.8, colsample_bytree0.8, reg_lambda1.0, eval_metricmlogloss, random_state42 ) cv_xgb cross_val_score(xgb_model, X_clr, y_subclass, cvcv, scoringaccuracy) print(XGBoost 5折平均准确率: %.3f % cv_xgb.mean())max_depth2是关键小样本下每棵树的深度超过3就一定过拟合。learning_rate0.05配合n_estimators120相当于用更多的浅树做集成。colsample_bytree0.8每次只用80%的特征给特征选择留随机性。这些参数不是玄学而是针对“几十行、十几列”的配置习惯换成几千行的大数据集同样的参数反而欠拟合。调参时不要开网格搜索硬扫小样本下网格搜索容易选到对随机种子敏感的参数组合。3.3 交叉验证切法GroupKFold避免同源样本泄漏C题数据里有“对象ID”这类比“采样编号”更高一级的分组字段一件古代玻璃器物可能被采样检测多次生成多行数据。如果按普通分层抽样把同一器物的不同采样拆到训练集和验证集模型会靠记住这个器物来刷分跨器物的泛化能力其实很差。这里要用GroupKFold。from sklearn.model_selection import GroupKFold group_cv GroupKFold(n_splits5) xgb_model.set_params(eval_metricmlogloss) gkf_scores cross_val_score( xgb_model, X_clr, y_subclass, cvgroup_cv, groupsdf[对象ID] ) print(GroupKFold准确率: %.3f % gkf_scores.mean())对比普通K折和GroupKFold的结果你会发现后面这个准确率更低但这才反映真实赛题场景。评委不会单独跑你的交叉验证但如果你在论文里写“按对象ID分组防止同源泄漏”这是明显的加分项。如果Excel里没有“对象ID”列退而求其次用“采样编号”当group但同一器物多次采样可能编号不同这种情况要结合题目文字说明来判断分组字段不要硬造。4. 风化前成分预测多输出回归与误差评估4.1 问题拆解哪些列是特征、哪些列是目标C题后半部分围绕“风化”展开。风化点的成分是风化后检测值想恢复它风化前的成分未风化点的成分代表风化前的基准。做法是以未风化样本为训练集以类型、颜色、纹饰作为特征以各成分为多输出目标。另外可以把风化不敏感的成分如Al2O3、SiO2也加入特征矩阵因为这些成分在风化前后变化小能充当样本的“身份信息”。# 构造特征矩阵: 类型/颜色/纹饰 one-hot feat pd.get_dummies(df[[类型, 颜色, 纹饰]], prefix, prefix_sep) # 目标列: 风化影响明显的成分 target_cols [K2O, PbO, BaO, CuO, P2O5] train_mask df[是否风化] 未风化 X_train_all, X_pred_all feat[train_mask], feat[~train_mask] y_train_all df.loc[train_mask, target_cols]这样拆的前提是“是否风化”字段明确给出。实际操作时部分样本只有风化点没有未风化点成对数据所以多输出回归是唯一出路。如果有“同一器物风化前后都测了”的配对数据那可以退化成逐样本差值回归但C题没有这种配对别硬找。特征矩阵里不要混入成分列否则相当于告诉模型答案。4.2 三套回归方案对比多元线性、Ridge、随机森林多输出回归有三种常用方案复杂度递增独立多输出线性回归、Ridge多输出、MultiOutputRegressor包随机森林。我的习惯是先跑Ridge再加随机森林做对照最后在论文里报告哪个更稳。from sklearn.linear_model import Ridge from sklearn.ensemble import RandomForestRegressor from sklearn.multioutput import MultiOutputRegressor from sklearn.metrics import mean_absolute_error, mean_squared_error from sklearn.model_selection import cross_val_predict # 方案1: Ridge, 带L2正则 ridge Ridge(alpha10.0) ridge.fit(X_train_all, y_train_all) y_pred_ridge ridge.predict(X_pred_all) # 方案2: 随机森林多输出 rf MultiOutputRegressor( RandomForestRegressor(n_estimators200, max_depth3, random_state42) ) rf.fit(X_train_all, y_train_all) y_pred_rf rf.predict(X_pred_all) # 评估: 用未风化样本做交叉验证, 看整体MAE和每个成分的MAE cv_pred cross_val_predict(ridge, X_train_all, y_train_all, cv5) mae_ridge mean_absolute_error(y_train_all, cv_pred) mae_per_col (y_train_all - cv_pred).abs().mean(axis0) print(Ridge整体MAE:, mae_ridge) print(各成分MAE:, mae_per_col)alpha10.0看起来很大但在CLR空间里各列经过对数化正则强度要比默认的1.0大一个量级才能压住方差。RandomForestRegressor(max_depth3)同样是为了限制单棵树复杂度防止对小训练集的记忆。交叉验证里报告整体MAE远远不够按成分列看误差能定位哪个成分最难预测——比如PbO的MAE可能很大因为它在铅钡玻璃中的绝对值本来就高。写论文时每个成分的MAE要分开列评委才信你做过误差分析。4.3 预测结果怎么写成论文里的统计描述回归模型的输出只是中间产物C题最后要落到“某类型风化前的化学成分均值与标准差”。这是论文写作的关键一步预测出每个风化样本的原始成分后按亚类合并统计。这里先把CLR预测值逆变换回原始百分比再分组聚合。# 逆变换: 把CLR预测值还原为成分百分比 log_gm np.log((X_clr.iloc[train_mask].apply(np.exp)).mean(axis1)) # 训练集每行log几何均值 pred_clr ridge.predict(X_pred_all) # 示例用Ridge pred_orig np.exp(pred_clr log_gm.values[:, None]) pred_df df[~train_mask].copy() pred_df[target_cols] pred_orig summary pred_df.groupby([类型, 亚类])[target_cols].agg([mean, std]) print(summary.round(2))按“类型亚类”算出的均值表可以直接作为论文中“风化前成分规律”的支撑数据。做这步时注意逆变换细节CLR只对训练时用过的列有效逆变换时要保留当时的log_gm向量不然数值对不上。有参赛队把CLR的预测结果直接当成分百分比写进论文均值表变成对数尺度评委一眼就能看出来。类型亚类K2O均值K2O标准差PbO均值PbO标准差铅钡玻璃亚类11.240.1831.554.32铅钡玻璃亚类22.010.2518.763.11高钾玻璃亚类18.771.021.120.43表格里每个数字都要能追溯到代码输出的原始统计不要手填。这个“均值±标准差”表就是论文结论的核心后面所有文字分析都得围绕它展开。5. C题避坑清单5个让参赛队翻车的细节5.1 坑1附件1和附件2的ID对不上现象merge之后发现大量行没有颜色和纹饰或者形状从80行变成60行。原因两个附件里的“采样编号”字符编码不一致比如一个全角一个半角或者末尾带不可见空格。解决合并前给编号列做astype(str).str.strip()再用pd.unique对比两边的编号集合把差异编号列出来人工核对。用inner join默默丢弃数据是最大忌讳要把对齐过程写进论文附录。5.2 坑2全样本准确率100%——数据泄漏现象亚类判别模型在训练集上准确率100%交叉验证也接近满分但论文里写不出有意义的规律。原因把聚类用的特征和判别用的特征用了同一份或者把“类型”字段本身当成聚类的输入。类型是给定的标签聚类应该在“类型内部”做如果整表直接聚类聚类结果会复制“高钾/铅钡”二分而不是亚类。解决严格按“先按类型分组→组内CLR→组内聚类→组内判别”的流程判别模型的特征里不许出现用于定义亚类的同一批成分或者用嵌套交叉验证放行。5.3 坑3成分百分比总和不是100%现象某几行成分加总只有92%或者104%。原因原始检测数据存在误差题目附件没有规范化到100%。如果数据处理时强行按行归一化到100%会破坏CLR的几何均值和后续回归误差。解决保留原始百分数只做CLR在论文“数据预处理”里注明“部分样本成分总和偏离100%可能包含未检测出的微量组分故不强制归一化”。这比假装数据完美要诚实得多。5.4 坑4把风化状态直接当标签忽略了“亚类”才是第一问现象第一问写成了“用分类模型区分风化/未风化”跑了半天发现准确率95%但评阅标准里要求的是亚类划分。原因没读懂问题。风化状态不是亚类标签它是样本的一个属性第一问要求的是在铅钡玻璃和高钾玻璃内部找亚类并分析亚类差异。解决动笔前先列一张“问题→数据集→模型输出”的映射表第一问输出亚类标签第二问输出亚类判别规则第三问输出风化前成分。写完论文再检查一遍这张映射比调参重要得多。5.5 坑5回归预测后成分归一化现象预测出的风化前成分加总偏离100%于是有人又做了一次按行归一化。原因多输出回归是逐成分独立预测的没有显式约束加和为100%预测值的总和天然会漂。解决不要在回归后强制归一化。加和漂移本身说明模型误差论文里报告MAE和平均偏差即可。如果非要保证加和为100%应该在逆变换后用等距对数比变换ILR或加一个成分总和误差项做微调而不是简单除以总和。6. 提交前代码自查清单评委复现你的结果只需三步代码不是论文的附赠品在国赛评审里代码质量会影响奖项走向。我给一份自己每次提交前都会过的自查清单按顺序执行五分钟内能完成。第一确认代码能从“读Excel”到“输出结果表”一步跑通没有隐藏的绝对路径依赖。把pd.read_excel(/Users/xxx/...)改成pd.read_excel(附件1.xlsx)文件放在和代码同一目录。第二确认所有随机过程固定了种子。KMeans的random_state、XGBoost的random_state、交叉验证的random_state一个都不能少。第三确认输出文件命名规范建议result1.csv、result2.csv、result3.csv对应题目三个问的答案表。评委没有时间在你的脚本里翻哪个变量对应哪个问题。tree # 建议的提交目录结构 . ├── 代码/ │ ├── main.py │ └── utils.py ├── 数据/ │ ├── 附件1.xlsx │ └── 附件2.xlsx ├── 结果/ │ ├── result1.csv │ ├── result2.csv │ └── result3.csv └── 论文.pdf这个树状结构不是摆设评委会按这个路径找你的结果文件。目录里不要放训练过程的中间产物或调试脚本保持干净。最后在代码开头写一段10行以内的中文注释说明“运行环境Python 3.9 pandas xgboost scikit-learn”并交代主体流程。这段注释是给评委看的也是给三天后重新打开代码的你自己的后悔药。做完这场题我最深的一个教训是数学建模竞赛比的不是模型多前沿而是对数据的诚实程度。CLR变换、GroupKFold、按类型分组填缺失值这些都不炫但它们决定了论文能不能被复现、结论能不能站住。希望这套以2022 C题为入口的完整处理链能帮到你下次再碰上成分数据或者小样本分类题照着这个框架走一遍至少不会在坑里过夜。本文还有配套的精品资源点击获取