简介这份资源是2022年全国大学生数学建模竞赛C题的完整解题资料面向备战数模竞赛的高校学生及指导教师聚焦古代玻璃文物的成分分析与鉴别这一典型赛题。压缩包共1个PDF文件约3.55MB内容涵盖赛题文档与配套代码便于对照阅读与复现。资料围绕数据预处理、数据探索、特征工程展开重点讲解主成分分析降维、因子分析构建风化到未风化成分的转换矩阵、敏感性分析验证模型鲁棒性以及灰色关联度分析比较不同类别玻璃化学成分的关联差异并给出高钾玻璃与铅钡玻璃的分类预测思路。目前已有926人学习下载适合需要系统复盘赛题建模流程、掌握降维与关联分析方法的读者参考借鉴。1. 2022全国大学生数学建模C题一份文档加代码到底能帮你复现什么2022年全国大学生数学建模竞赛C题题目背景是古代玻璃制品的成分分析与鉴别。这道题在赛后流传的资料里最常见的组合就是一份论文文档加一套代码。很多人拿到手第一反应是收藏第二反应是打开看一眼摘要然后就再也没动过。真正的问题是这份文档加代码能不能让你在本地把整条分析链路跑通能不能让你理解每一步为什么这么做能不能让你在遇到同类数据时自己改参数。这道题的核心工作可以拆成三块第一块是数据清洗和缺失值处理玻璃成分数据里有大量零值和缺失处理方式直接决定后续结论第二块是降维和分类主成分分析、因子分析、灰色关联度这些方法在题目里都有出场空间第三块是判别与预测用化学成分判断玻璃类型和风化状态。适合的读者是正在准备数学建模竞赛的学生、需要复现经典案例的数据分析新手以及想找一套完整流程练手的Python使用者。如果你只想要一个能跑的脚本网上很多如果你想搞清楚每一步的参数为什么这么设、换一组数据该怎么调那这份文档加代码值得你花时间拆开看。2. 从原始成分表到可分析数据清洗与缺失值处理的完整链路2.1 玻璃成分数据的三个特殊结构2022年C题的数据表里每一行是一个玻璃样本每一列是一种化学成分的含量比如二氧化硅、氧化钠、氧化钾、氧化铅、氧化钡等。表面上看是一张普通的数值表但实际打开会发现三个麻烦。第一很多成分列存在大量零值这些零值不是真的含量为零而是检测手段没有测到或者低于检出限。第二不同玻璃类型高钾、铅钡、普通的缺失模式不一样铅钡玻璃的钡含量列缺失少高钾玻璃的钾含量列缺失少如果直接按列填充会把类型信息抹掉。第三成分之间存在约束所有主要成分加起来应该接近百分之百但原始数据因为缺失和检测误差加总经常偏离。这三个结构决定了你不能上来就dropna()或者fillna(0)。常见做法是先把零值和缺失值分开标记再按玻璃类型分组看缺失比例最后决定是删除样本还是填充。我一般会先做一张缺失热力图把行和列都排序肉眼确认缺失是不是随机分布。2.2 用Python做分组缺失统计和条件填充下面这段代码做三件事读数据、按类型统计缺失、对成分列做基于类型中位数的填充。注意这里用的是中位数而不是均值因为成分数据里存在极端值均值会被拉偏。import pandas as pd import numpy as np # 读取原始数据假设文件名为 glass.csv df pd.read_csv(glass.csv, encodingutf-8) # 把零值替换为NaN因为零值在成分数据里通常代表未检出 component_cols [c for c in df.columns if c not in [样本编号, 玻璃类型, 风化状态]] df[component_cols] df[component_cols].replace(0, np.nan) # 按玻璃类型统计每列的缺失比例 missing_by_type df.groupby(玻璃类型)[component_cols].apply(lambda x: x.isna().mean()) print(missing_by_type.round(3)) # 按玻璃类型分组用该类型的中位数填充 df[component_cols] df.groupby(玻璃类型)[component_cols].transform( lambda x: x.fillna(x.median()) ) # 如果整列在某个类型下全缺失中位数填充会留下NaN用全局中位数兜底 df[component_cols] df[component_cols].fillna(df[component_cols].median()) # 检查填充后是否还有缺失 print(剩余缺失值数量, df[component_cols].isna().sum().sum())逻辑说明第一步把零值转成NaN是为了让后续的缺失统计和填充统一处理。第二步按玻璃类型分组统计是为了看清缺失是不是和类型相关。第三步用transform配合fillna(x.median())保证填充值来自同类型样本不会把铅钡玻璃的钡含量中位数填到高钾玻璃上。最后一步兜底是防止某个类型下某列全空导致中位数也是NaN。参数说明replace(0, np.nan)里的0可以根据实际数据调整如果某些成分确实可能为零比如某些微量元素就不要替换。groupby(玻璃类型)的列名要和你的数据表一致如果表头是英文改成对应的英文列名。中位数填充适合成分数据因为成分数据分布偏斜中位数比均值稳健。2.3 成分加和校验与异常样本标记填充完之后要做一次加和校验。把所有主要成分列的数值加起来看每个样本的总和是不是在合理范围内比如95%到105%之间。超出这个范围的样本要么是填充引入了偏差要么是原始数据有录入错误。我一般会把加和异常的样本单独标出来在后续分析里做敏感性测试看删掉它们结论会不会变。# 计算每个样本的主要成分加和 df[成分总和] df[component_cols].sum(axis1) # 标记加和异常样本 df[加和异常] (df[成分总和] 95) | (df[成分总和] 105) # 输出异常样本数量和占比 print(加和异常样本数, df[加和异常].sum()) print(占比, df[加和异常].mean().round(3)) # 查看异常样本的类型分布 print(df[df[加和异常]][玻璃类型].value_counts())这段代码的输出会告诉你异常样本集中在哪个玻璃类型。如果某个类型异常比例特别高说明该类型的缺失模式更复杂可能需要单独处理而不是统一用中位数填充。到这一步你手里就有一张干净且带质量标记的数据表可以进入降维和分类环节。3. 主成分分析与因子分析降维之前先想清楚你要解释什么3.1 PCA和因子分析在玻璃数据上的分工主成分分析PCA和因子分析因子分析经常被放在一起提但在2022年C题里它们的作用不一样。PCA是把原始成分变量线性组合成几个互不相关的主成分目的是压缩维度、去掉共线性适合在分类之前做特征提取。因子分析假设原始变量背后有几个潜在的公共因子比如“助熔剂因子”“稳定剂因子”目的是解释成分之间的相关性结构适合做机理解释。很多参赛论文把两个都做了但没写清楚为什么两个都做。我的建议是如果你后续要用判别分析或者聚类先用PCA降维把主成分得分作为输入如果你要写“铅钡玻璃的钡含量和铅含量共同反映了一个什么工艺特征”用因子分析做旋转后的载荷矩阵来解释。两个方法的输入都是标准化后的成分矩阵因为成分量纲不同不标准化的话二氧化硅的大数值会主导主成分方向。3.2 用sklearn跑PCA并确定保留几个主成分下面这段代码对标准化后的成分矩阵做PCA并输出方差解释率和载荷矩阵。关键参数是n_components我一般先设成和变量数一样看方差解释率再决定保留几个。from sklearn.preprocessing import StandardScaler from sklearn.decomposition import PCA # 提取成分列排除标记列 X df[component_cols].values # 标准化PCA对量纲敏感必须做 scaler StandardScaler() X_scaled scaler.fit_transform(X) # 先保留所有主成分看方差解释率 pca_full PCA(n_componentslen(component_cols)) pca_full.fit(X_scaled) # 输出每个主成分的方差解释率和累计解释率 explained pca_full.explained_variance_ratio_ cum_explained np.cumsum(explained) for i, (e, c) in enumerate(zip(explained, cum_explained), 1): print(fPC{i}: 解释率{e:.3f}, 累计{c:.3f}) # 一般保留累计解释率达到85%以上的主成分 n_keep np.argmax(cum_explained 0.85) 1 print(保留主成分数, n_keep) # 用保留的主成分做降维 pca PCA(n_componentsn_keep) X_pca pca.fit_transform(X_scaled) # 输出载荷矩阵看每个主成分主要受哪些成分影响 loadings pd.DataFrame( pca.components_.T, indexcomponent_cols, columns[fPC{i1} for i in range(n_keep)] ) print(loadings.round(3))逻辑说明StandardScaler把每个成分列变成均值0方差1消除量纲影响。PCA(n_componentslen(component_cols))先拟合全部主成分拿到方差解释率曲线。np.argmax(cum_explained 0.85) 1找到累计解释率首次达到85%的位置这个阈值可以根据题目要求调整有的论文用90%。载荷矩阵告诉你每个主成分的物理含义比如PC1在二氧化硅上载荷高、在氧化铅上载荷低可能代表“硅质-铅质”对比轴。参数说明n_components如果设成小数比如0.85sklearn会自动保留达到85%解释率的主成分数但为了输出载荷矩阵方便我习惯先手动算。StandardScaler默认按列标准化如果你的数据里某些成分列方差极小标准化后会放大噪声这时候要考虑先删掉近似常数列。3.3 因子分析旋转后的载荷怎么读因子分析和PCA的代码结构类似但多了旋转步骤。旋转的目的是让载荷矩阵更稀疏每个变量只在一个因子上有高载荷方便命名。下面用factor_analyzer库做因子分析如果你没装这个库可以用pip install factor_analyzer。from factor_analyzer import FactorAnalyzer # 先做KMO检验和Bartlett球形检验判断数据是否适合因子分析 from factor_analyzer.factor_analyzer import calculate_kmo, calculate_bartlett_sphericity kmo_all, kmo_model calculate_kmo(X_scaled) chi_square, p_value calculate_bartlett_sphericity(X_scaled) print(fKMO值{kmo_model:.3f}) print(fBartlett检验p值{p_value:.3f}) # 设因子数为3用最大方差法旋转 fa FactorAnalyzer(n_factors3, rotationvarimax) fa.fit(X_scaled) # 输出旋转后的载荷矩阵 loadings_fa pd.DataFrame( fa.loadings_, indexcomponent_cols, columns[fFactor{i1} for i in range(3)] ) print(loadings_fa.round(3)) # 输出每个变量的共同度 communalities pd.DataFrame( fa.get_communalities(), indexcomponent_cols, columns[共同度] ) print(communalities.round(3))逻辑说明KMO值衡量变量间的偏相关性一般大于0.6才适合做因子分析小于0.5说明变量间相关性太弱不适合。Bartlett检验的p值小于0.05说明相关矩阵不是单位矩阵有因子结构。rotationvarimax是最大方差正交旋转让每个因子上的载荷两极分化。共同度表示每个变量被公共因子解释的比例共同度太低比如小于0.4说明这个变量不适合放进因子模型。参数说明n_factors3是我根据玻璃成分的工艺背景预设的实际做的时候可以先看特征值大于1的因子个数或者看碎石图拐点。rotation还可以选promax做斜交旋转如果因子之间允许相关斜交旋转更合适但解释起来比正交旋转复杂。4. 灰色关联度与判别分析把分类结果落到具体样本上4.1 灰色关联度在玻璃类型判别中的用法灰色关联度分析适合样本量小、信息不完全的场景。在2022年C题里你可以把已知类型的玻璃样本作为参考序列把待判样本作为比较序列计算每个待判样本和各类参考序列的关联度关联度最大的类就是判别结果。和判别分析相比灰色关联度不要求数据服从特定分布对样本量要求低但它的结果受分辨系数影响。具体做法是先按玻璃类型把已知样本的成分均值算出来作为该类型的参考序列。然后把待判样本的成分向量和每个参考序列做关联系数计算最后加权平均得到关联度。分辨系数一般取0.5这个值影响关联系数的区分度取太小会导致关联度都接近1取太大会导致关联度都接近0。4.2 灰色关联度的Python实现与分辨系数调参下面这段代码实现灰色关联度判别输入是训练集和测试集输出是测试集样本的预测类型和关联度矩阵。def grey_relational_degree(reference, compare, rho0.5): 计算比较序列与参考序列的灰色关联度 reference: 参考序列形状 (n_features,) compare: 比较序列矩阵形状 (n_samples, n_features) rho: 分辨系数通常取0.5 # 计算绝对差矩阵 diff np.abs(compare - reference) # 两级最小差和最大差 min_diff diff.min() max_diff diff.max() # 关联系数 xi (min_diff rho * max_diff) / (diff rho * max_diff) # 关联度取均值 degree xi.mean(axis1) return degree # 按类型计算参考序列成分均值 reference_seqs df[df[加和异常] False].groupby(玻璃类型)[component_cols].mean() # 对待判样本计算与各类型的关联度 test_samples df[df[加和异常] True][component_cols].values results {} for glass_type, ref in reference_seqs.iterrows(): degree grey_relational_degree(ref.values, test_samples, rho0.5) results[glass_type] degree # 整理成DataFrame每行是一个样本每列是一个类型的关联度 result_df pd.DataFrame(results) result_df[预测类型] result_df.idxmax(axis1) print(result_df.round(3))逻辑说明diff矩阵是每个待判样本的每个成分与参考序列对应成分的绝对差。min_diff和max_diff是所有差里的最小值和最大值对应灰色关联理论里的两级最小差和最大差。关联系数公式(min_diff rho * max_diff) / (diff rho * max_diff)保证关联系数在0到1之间。最后对每个样本的所有成分关联系数取均值得到该样本与参考序列的关联度。参数说明rho是分辨系数取值区间(0,1)我一般先用0.5跑一遍然后试0.3和0.7看预测结果稳不稳定。如果三个值下预测类型一致说明结果稳健如果变化大说明样本在类型边界上需要结合其他方法判断。reference_seqs用的是均值你也可以用中位数或者去掉异常样本后的均值取决于你对参考序列稳健性的要求。4.3 判别分析做交叉验证的完整流程判别分析是更标准的分类方法线性判别分析LDA假设各类协方差矩阵相同二次判别分析QDA不要求这个假设。玻璃成分数据里不同类型玻璃的成分协方差可能不同所以QDA有时比LDA更合适但QDA参数多小样本下容易过拟合。我一般两个都跑用交叉验证比较准确率。from sklearn.discriminant_analysis import LinearDiscriminantAnalysis, QuadraticDiscriminantAnalysis from sklearn.model_selection import cross_val_score, StratifiedKFold # 准备特征和标签只用加和正常的样本做训练 train_df df[df[加和异常] False] X_train train_df[component_cols].values y_train train_df[玻璃类型].values # 标准化 scaler_cv StandardScaler() X_train_scaled scaler_cv.fit_transform(X_train) # 分层交叉验证保证每折里各类比例一致 cv StratifiedKFold(n_splits5, shuffleTrue, random_state42) # LDA交叉验证 lda LinearDiscriminantAnalysis() lda_scores cross_val_score(lda, X_train_scaled, y_train, cvcv, scoringaccuracy) print(fLDA准确率{lda_scores.mean():.3f} ± {lda_scores.std():.3f}) # QDA交叉验证 qda QuadraticDiscriminantAnalysis() qda_scores cross_val_score(qda, X_train_scaled, y_train, cvcv, scoringaccuracy) print(fQDA准确率{qda_scores.mean():.3f} ± {qda_scores.std():.3f}) # 用全部训练数据拟合LDA输出混淆矩阵 from sklearn.metrics import confusion_matrix lda.fit(X_train_scaled, y_train) y_pred lda.predict(X_train_scaled) print(LDA混淆矩阵) print(confusion_matrix(y_train, y_pred))逻辑说明StratifiedKFold保证每一折里各玻璃类型的比例和整体一致避免某一折里某个类型样本太少导致评估偏差。cross_val_score返回5折的准确率取均值和标准差标准差大说明模型对数据划分敏感。混淆矩阵告诉你哪些类型容易被混淆比如高钾玻璃和普通玻璃可能在某个成分上重叠。参数说明n_splits5是常用选择样本量少可以改成3样本量多可以改成10。random_state42固定随机种子保证结果可复现。scoringaccuracy是准确率如果类别不平衡可以改成f1_macro。LDA和QDA的准确率如果差距不大优先选LDA因为参数少、解释性强。5. 避坑与排查复现2022年C题时最容易翻车的五个地方5.1 零值当缺失处理导致成分加和异常现象填充完缺失值后成分总和普遍超过105%PCA载荷矩阵第一主成分几乎全是正载荷看不出对比关系。原因原始数据里的零值有两种含义一种是未检出一种是真为零。如果把所有零值都当缺失填充原本真为零的成分被填成了中位数加和自然偏高。解决先看成分的检测背景对于玻璃主要成分二氧化硅、氧化钠、氧化铅等零值大概率是未检出可以填充对于微量元素零值可能是真为零保留。填充后做加和校验超过105%的样本标记出来在后续分析里做敏感性测试。5.2 PCA之前忘记标准化导致主成分被大数值成分主导现象PCA方差解释率显示PC1解释了90%以上的方差载荷矩阵里PC1在二氧化硅上有极高载荷其他成分载荷接近零。原因二氧化硅的含量通常在60%到80%之间而其他成分可能在个位数甚至小数点后不标准化的话PCA会优先沿着方差最大的方向投影也就是二氧化硅的方向。解决PCA之前必须用StandardScaler做列标准化。标准化之后每个成分的方差都是1PCA找的是相关性结构而不是量纲差异。如果你用的是自己写的PCA记得先减均值再除以标准差。5.3 因子分析因子数选太多导致解释困难现象因子分析保留5个因子旋转后每个因子上都有两三个高载荷变量没法给因子起一个合理的工艺名称。原因因子数太多每个因子解释的方差少旋转后载荷分散。或者因子数太少变量共同度低信息丢失多。解决先用特征值大于1准则初选因子数再看碎石图拐点最后结合工艺背景。玻璃成分数据一般3到4个因子就够了比如“硅质因子”“助熔剂因子”“着色剂因子”。如果旋转后因子含义不清试试斜交旋转promax或者删掉共同度低于0.4的变量重新做。5.4 灰色关联度分辨系数取极端值导致判别失效现象分辨系数取0.1时所有待判样本与各类型的关联度都在0.95以上无法区分取0.9时关联度都在0.3以下也难区分。原因分辨系数rho控制关联系数的区分度。rho太小rho * max_diff项太小关联系数趋近于1rho太大关联系数趋近于0。解决rho取0.5是默认值先用0.5跑一遍然后试0.3、0.4、0.6、0.7看预测类型是否稳定。如果不同rho下预测类型变化大说明样本在类型边界上需要结合判别分析或者增加特征。5.5 交叉验证时没做分层导致某折缺少某个类型现象5折交叉验证里某一折的准确率特别低或者报错提示某个类型在测试集里没有出现。原因玻璃类型分布不均匀比如高钾玻璃样本少随机划分时可能全部被分到训练集测试集里没有高钾玻璃准确率计算就不完整。解决用StratifiedKFold代替KFold保证每一折里各类型的比例和整体一致。如果某个类型样本数少于折数比如只有3个样本但做5折那就把折数降到3或者用留一法交叉验证。6. 把2022年C题的代码改造成你自己的建模模板这套文档加代码最大的价值不是让你复现一遍2022年的结果而是让你把它拆成可复用的模块。我自己的习惯是建三个文件夹data_cleaning放清洗和缺失值处理脚本feature_engineering放PCA、因子分析、灰色关联度的函数modeling放判别分析和交叉验证的流程。每个模块的输入输出用DataFrame约定好换一道题只需要改列名和参数。具体来说清洗模块的入口是一个原始CSV和一个配置字典配置字典里写清楚哪些列是成分列、哪些列是标签列、零值是否替换。降维模块的入口是标准化后的特征矩阵和保留主成分数的阈值输出是降维后的矩阵和载荷矩阵。判别模块的入口是特征矩阵和标签输出是交叉验证准确率和混淆矩阵。这样拆完之后你拿到2023年或者2024年的C题只需要改配置字典和特征列流程代码基本不用动。还有一个技巧把每次跑的方差解释率、交叉验证准确率、混淆矩阵写到一个日志文件里带上时间戳和参数。我吃过亏调了半天的参数结果忘了之前哪组参数效果最好只能重新跑。后来养成习惯每次跑完自动追加一行记录参数和指标一目了然。这个习惯在比赛期间尤其重要因为时间紧你不可能记住每一组实验。最后说一个我自己的教训不要等到论文快写完才去整理代码。比赛的时候代码和论文是同步迭代的你改一个参数论文里的数字就要跟着改。如果代码没有模块化改一处要动好几个文件很容易漏改。我一般会在代码里把关键结果直接输出成Markdown表格复制到论文里就行减少手工抄错的机会。希望帮到你。本文还有配套的精品资源点击获取