简介：这份资源是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', encoding='utf-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(axis=1) # 标记加和异常样本 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_components=len(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(f'PC{i}: 解释率={e:.3f}, 累计={c:.3f}') # 一般保留累计解释率达到85%以上的主成分 n_keep = np.argmax(cum_explained >= 0.85) + 1 print('保留主成分数：', n_keep) # 用保留的主成分做降维 pca = PCA(n_components=n_keep) X_pca = pca.fit_transform(X_scaled) # 输出载荷矩阵，看每个主成分主要受哪些成分影响 loadings = pd.DataFrame( pca.components_.T, index=component_cols, columns=[f'PC{i+1}' for i in range(n_keep)] ) print(loadings.round(3))逻辑说明：StandardScaler把每个成分列变成均值0方差1，消除量纲影响。PCA(n_components=len(component_cols))先拟合全部主成分，拿到方差解释率曲线。np.argmax(cum_explained >= 0.85) + 1找到累计解释率首次达到85%的位置，这个阈值可以根据题目要求调整，有的论文用90%。载荷矩阵告诉你每个主成分的物理含义，比如PC1在二氧化硅上载荷高、在氧化铅上载荷低，可能代表“硅质-铅质”对比轴。
参数说明：n_components如果设成小数比如0.85，sklearn会自动保留达到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(f'KMO值：{kmo_model:.3f}') print(f'Bartlett检验p值：{p_value:.3f}') # 设因子数为3，用最大方差法旋转 fa = FactorAnalyzer(n_factors=3, rotation='varimax') fa.fit(X_scaled) # 输出旋转后的载荷矩阵 loadings_fa = pd.DataFrame( fa.loadings_, index=component_cols, columns=[f'Factor{i+1}' for i in range(3)] ) print(loadings_fa.round(3)) # 输出每个变量的共同度 communalities = pd.DataFrame( fa.get_communalities(), index=component_cols, columns=['共同度'] ) print(communalities.round(3))逻辑说明：KMO值衡量变量间的偏相关性，一般大于0.6才适合做因子分析，小于0.5说明变量间相关性太弱，不适合。Bartlett检验的p值小于0.05说明相关矩阵不是单位矩阵，有因子结构。rotation='varimax'是最大方差正交旋转，让每个因子上的载荷两极分化。共同度表示每个变量被公共因子解释的比例，共同度太低比如小于0.4，说明这个变量不适合放进因子模型。
参数说明：n_factors=3是我根据玻璃成分的工艺背景预设的，实际做的时候可以先看特征值大于1的因子个数，或者看碎石图拐点。rotation还可以选promax做斜交旋转，如果因子之间允许相关，斜交旋转更合适，但解释起来比正交旋转复杂。
4. 灰色关联度与判别分析：把分类结果落到具体样本上
4.1 灰色关联度在玻璃类型判别中的用法
灰色关联度分析适合样本量小、信息不完全的场景。在2022年C题里，你可以把已知类型的玻璃样本作为参考序列，把待判样本作为比较序列，计算每个待判样本和各类参考序列的关联度，关联度最大的类就是判别结果。和判别分析相比，灰色关联度不要求数据服从特定分布，对样本量要求低，但它的结果受分辨系数影响。
具体做法是：先按玻璃类型把已知样本的成分均值算出来，作为该类型的参考序列。然后把待判样本的成分向量和每个参考序列做关联系数计算，最后加权平均得到关联度。分辨系数一般取0.5，这个值影响关联系数的区分度，取太小会导致关联度都接近1，取太大会导致关联度都接近0。
4.2 灰色关联度的Python实现与分辨系数调参
下面这段代码实现灰色关联度判别，输入是训练集和测试集，输出是测试集样本的预测类型和关联度矩阵。
def grey_relational_degree(reference, compare, rho=0.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(axis=1) 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, rho=0.5) results[glass_type] = degree # 整理成DataFrame，每行是一个样本，每列是一个类型的关联度 result_df = pd.DataFrame(results) result_df['预测类型'] = result_df.idxmax(axis=1) 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_splits=5, shuffle=True, random_state=42) # LDA交叉验证 lda = LinearDiscriminantAnalysis() lda_scores = cross_val_score(lda, X_train_scaled, y_train, cv=cv, scoring='accuracy') print(f'LDA准确率：{lda_scores.mean():.3f} ± {lda_scores.std():.3f}') # QDA交叉验证 qda = QuadraticDiscriminantAnalysis() qda_scores = cross_val_score(qda, X_train_scaled, y_train, cv=cv, scoring='accuracy') print(f'QDA准确率：{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_splits=5是常用选择，样本量少可以改成3，样本量多可以改成10。random_state=42固定随机种子，保证结果可复现。scoring='accuracy'是准确率，如果类别不平衡，可以改成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做列标准化。标准化之后，每个成分的方差都是1，PCA找的是相关性结构而不是量纲差异。如果你用的是自己写的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项太小，关联系数趋近于1；rho太大，关联系数趋近于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表格，复制到论文里就行，减少手工抄错的机会。希望帮到你。
本文还有配套的精品资源，点击获取
