2022数学建模C题复盘:PCA、因子分析与灰色关联度实战
发布时间:2026/10/2 14:42:16来源:尧图网络
简介这份资源是2022年全国大学生数学建模竞赛C题的完整解题资料面向备战数模竞赛的高校学生及指导教师聚焦古代玻璃文物的成分分析与鉴别这一典型赛题。压缩包共1个PDF文件约3.55MB内容涵盖赛题文档与配套代码便于对照阅读与复现。资料围绕数据预处理、数据探索、特征工程展开系统运用主成分分析降维、因子分析构建风化到未风化成分的转换矩阵并借助敏感性分析验证模型鲁棒性同时以灰色关联度分析揭示不同类别玻璃化学成分的关联差异。读者可据此掌握从风化系数定义、亚类划分到未知样品分类预测的完整建模链路理解高钾玻璃与铅钡玻璃在风化效应上的规律并获得可直接参考的代码实现与结果分析思路。目前已有926人学习适合需要系统复盘赛题、查漏补缺的参赛者。1. 2022全国大学生数学建模C题一份文档加代码的复盘为什么值得你花时间2022年全国大学生数学建模竞赛C题题目围绕古代玻璃制品的成分分析与鉴别展开。这道题在赛后讨论度一直很高原因不复杂它把数据预处理、降维、分类、关联分析这几件事串成了一条完整的链路而且数据里埋了不少“看起来能用、直接跑就翻车”的坑。你手里如果只有一份“文档代码”的压缩包真正的问题不是有没有资料而是这份资料能不能让你复现出合理结果、能不能讲清楚每一步为什么这么做。我见过太多人拿到代码直接run结果第一问的分类准确率看着还行第二问的关联规则一塌糊涂第三问的敏感性分析干脆没做。这篇笔记就按“文档代码”这个组合把C题从读题到出图的全流程拆开重点放在主成分分析、因子分析、灰色关联度这几个热搜词对应的实操环节上顺带把Python环境配置和常见报错也带一笔。适合正在准备华为杯数学建模大赛、研究生数学建模或者想拿往年国赛题练手的人。你不需要先看完一堆理论跟着走一遍至少能判断手里的代码包值不值得改、往哪个方向改。2. 先搞清楚C题的数据长什么样玻璃成分表里的三个陷阱2.1 表单结构、缺失值和“未检测”的真实含义2022年C题给的是古代玻璃制品的化学成分数据分两个表单一个记录玻璃类型、纹饰、颜色、表面风化等基本信息另一个是各化学成分的含量百分比。表面上看是标准的结构化数据但实际打开会发现几个问题。第一成分列有大量空值不是随机缺失而是“未检测”或“低于检出限”。直接dropna会丢掉大量样本直接fillna(0)又会把“未检测”和“含量为零”混为一谈。我一般会先统计每列的缺失比例超过30%的列考虑整列舍弃低于10%的用中位数填充中间区间的单独标记。第二风化与否这个字段在后续分类里权重很高但原始数据里可能有拼写不一致比如“风化”和“已风化”混用需要先做值统一。第三成分百分比加起来不一定等于100%有的样本因为四舍五入或检测方法差异总和在95%到105%之间波动做归一化之前要先检查。import pandas as pd import numpy as np # 读取两个表单注意sheet名按实际文件调整 df_info pd.read_excel(2022C.xlsx, sheet_name表单1) df_chem pd.read_excel(2022C.xlsx, sheet_name表单2) # 查看缺失比例 missing_ratio df_chem.isnull().sum() / len(df_chem) print(missing_ratio[missing_ratio 0].sort_values(ascendingFalse)) # 成分列求和检查是否接近100 chem_cols [c for c in df_chem.columns if c not in [文物编号, 表面风化]] df_chem[sum_check] df_chem[chem_cols].sum(axis1) print(df_chem[sum_check].describe()) # 缺失填充策略低缺失列用中位数高缺失列标记 low_missing missing_ratio[missing_ratio 0.1].index.tolist() for col in low_missing: df_chem[col] df_chem[col].fillna(df_chem[col].median())这段代码先做缺失统计和求和校验目的是让你在动手建模之前对数据质量有个底。missing_ratio帮你决定哪些列能留、哪些列该扔sum_check的分布如果集中在100附近说明数据整体可信如果出现大量偏离就要考虑是不是单位不统一。填充策略里中位数比均值稳因为成分数据常有偏态。注意这里没有直接对高缺失列做填充而是留到后面用模型自带缺失处理或单独分析避免引入太大偏差。2.2 类型变量编码风化、纹饰、颜色怎么进模型玻璃类型、纹饰、颜色这些字段是字符串直接扔给sklearn会报错。常见做法是独热编码但独热之后维度会膨胀尤其是颜色这种取值较多的列。我一般会先看每个类别变量的取值分布取值超过10个的考虑做频次编码或者合并稀有类别。表面风化只有两个值用map转成0和1最省事。纹饰和颜色可以用pd.get_dummies但记得加drop_firstTrue避免共线性。这里有个容易忽略的点训练集和测试集要一起做编码或者用同一个编码器否则类别对齐会出问题。如果你手里的代码是分开编码的检查一下有没有出现训练集有、测试集没有的类别那种情况预测时会直接报错。# 风化字段二值化 df_info[表面风化] df_info[表面风化].map({无风化: 0, 风化: 1}) # 纹饰和颜色做独热编码合并到主表 df_info pd.get_dummies(df_info, columns[纹饰, 颜色], drop_firstTrue) # 合并信息表和成分表按文物编号对齐 df pd.merge(df_info, df_chem, on文物编号, howinner) print(df.shape)合并之后要检查行数有没有异常减少如果inner join丢了很多样本说明两个表单的编号对不上需要手动核对。drop_firstTrue是为了避免虚拟变量陷阱线性模型里尤其要注意。如果你后面用的是树模型其实可以保留全部哑变量但为了统一处理我习惯都加上。2.3 高维成分数据的可视化初探PCA之前先画个相关矩阵在跑主成分分析之前先看成分之间的相关性。玻璃成分里SiO2、Na2O、K2O、PbO这些往往此消彼长相关矩阵能帮你判断哪些变量适合降维。如果大部分相关系数都在0.3以下PCA的效果可能不明显这时候要考虑是不是数据本身就不适合降维或者需要先做标准化。我一般会画一个热力图快速扫一眼有没有强相关块。这一步不需要写太多代码但能省掉后面很多无效调参。import seaborn as sns import matplotlib.pyplot as plt corr df[chem_cols].corr() plt.figure(figsize(12, 10)) sns.heatmap(corr, cmapRdBu_r, center0, annotFalse) plt.title(成分相关矩阵) plt.show()热力图里如果看到几个变量抱团比如PbO和BaO高度相关那PCA的第一主成分大概率由它们主导。这时候你要想清楚降维是为了可视化还是为了去共线性如果是为了后续分类保留太多主成分反而会丢信息。我一般会保留累计方差贡献率到85%左右但不会死守这个数具体看碎石图的拐点。3. 主成分分析在C题里的正确打开方式别把PCA当万能降维3.1 标准化、KMO检验和碎石图PCA前的三个必做动作主成分分析PCA是热搜词里出现频率最高的但很多人直接对原始数据跑PCA结果第一主成分方差贡献率90%以上看着很漂亮实际上是因为量纲没统一。成分数据虽然都是百分比但有的列方差大、有的列方差小不标准化的话方差大的列会主导主成分方向。所以第一步一定是StandardScaler。第二步做KMO检验和Bartlett球形检验判断数据是否适合因子分析或PCA。KMO低于0.6一般就不建议硬做了说明变量间相关性太弱。第三步画碎石图看特征值下降的拐点决定保留几个主成分。这三个动作做完你才能说“我做了PCA”而不是“我调了个包”。from sklearn.preprocessing import StandardScaler from factor_analyzer import FactorAnalyzer from factor_analyzer.factor_analyzer import calculate_kmo, calculate_bartlett_sphericity # 标准化 scaler StandardScaler() X_scaled scaler.fit_transform(df[chem_cols].fillna(0)) # KMO和Bartlett检验 kmo_all, kmo_model calculate_kmo(X_scaled) chi_square, p_value calculate_bartlett_sphericity(X_scaled) print(fKMO: {kmo_model:.3f}, Bartlett p-value: {p_value:.3e}) # 碎石图 fa FactorAnalyzer(n_factorslen(chem_cols), rotationNone) fa.fit(X_scaled) ev, v fa.get_eigenvalues() plt.plot(range(1, len(ev)1), ev, o-) plt.axhline(y1, colorr, linestyle--) plt.xlabel(主成分个数) plt.ylabel(特征值) plt.show()KMO值如果低于0.6说明变量间共同方差偏低PCA降维后主成分解释力有限这时候要么换方法要么重新筛选变量。Bartlett检验的p值小于0.05才适合做因子分析。碎石图里特征值大于1的主成分个数是常用标准但更稳妥的是看累计方差贡献率。factor_analyzer这个库不是sklearn自带的需要pip install如果环境里没有用sklearn的PCA也能做只是没有KMO检验。3.2 主成分个数的选择累计方差85%还是特征值大于1这两个标准经常打架。特征值大于1是Kaiser准则适合变量数不多的情况累计方差85%更直观但主成分个数可能偏多。我的做法是两个都算取交集。如果特征值大于1的主成分有4个累计方差85%需要5个那就选4个然后看第5个的方差贡献率是不是接近10%如果接近就加到5个。C题的数据里成分变量大概十来个一般保留3到4个主成分就能解释80%以上的方差。保留太多后续分类模型容易过拟合保留太少信息损失大分类边界模糊。这里没有绝对正确的答案但你可以通过对比不同主成分个数下的分类准确率来反推。from sklearn.decomposition import PCA pca PCA(n_components0.85) # 保留累计方差85% X_pca pca.fit_transform(X_scaled) print(f保留主成分数: {pca.n_components_}) print(f累计方差贡献率: {pca.explained_variance_ratio_.sum():.3f}) # 如果要用固定个数 pca_fixed PCA(n_components4) X_pca_fixed pca_fixed.fit_transform(X_scaled) print(pca_fixed.explained_variance_ratio_)n_components0.85会让PCA自动选择达到85%方差所需的主成分数省去手动判断。但要注意这个85%是经验值不是铁律。如果后续分类模型对降维敏感可以试0.8、0.85、0.9三个档看哪个效果最好。explained_variance_ratio_输出每个主成分的方差贡献率第一个通常最大后面递减。3.3 主成分得分的业务解释第一主成分到底代表什么PCA跑完不是终点你得能解释每个主成分的含义。看载荷矩阵components_每个主成分在原始变量上的权重。比如第一主成分在PbO和BaO上载荷很高在SiO2上载荷为负那它可能代表“铅钡玻璃体系”的强度。第二主成分在K2O上载荷高可能代表“钾玻璃体系”。这种解释在写论文时很重要评委看的就是你能不能把数学结果翻译成材料学语言。如果载荷矩阵很乱每个主成分都均匀分布说明数据本身没有明显的因子结构这时候硬解释就是玄学不如老老实实说“降维后主成分物理意义不明确仅用于后续分类”。loadings pd.DataFrame( pca_fixed.components_.T, columns[fPC{i1} for i in range(4)], indexchem_cols ) print(loadings.round(3).sort_values(PC1, ascendingFalse))载荷矩阵的绝对值越大说明该变量对主成分的贡献越大。正负号代表方向不影响重要性。我一般会挑每个主成分里载荷绝对值前3的变量结合玻璃考古的背景知识给个名字。如果实在解释不通就在论文里注明“主成分命名仅供参考”不要强行编。4. 因子分析与灰色关联度C题后两问的实操路径4.1 因子旋转让载荷矩阵更好看但别过度解读因子分析和PCA经常被混用区别在于因子分析有旋转能让载荷矩阵更稀疏每个变量只在一个因子上载荷高。C题里如果要做变量归类因子分析比PCA更合适。旋转方法常用方差最大化Varimax旋转后因子之间仍然正交。但旋转有个副作用因子方差贡献率会重新分配原来第一因子贡献30%旋转后可能变成25%这是正常的。不要因为旋转后贡献率下降就认为模型变差了。我一般会对比旋转前后的载荷矩阵如果旋转后每个因子都有明确的变量簇那就用旋转后的结果如果旋转后更乱说明数据不适合因子分析退回PCA。fa_varimax FactorAnalyzer(n_factors4, rotationvarimax) fa_varimax.fit(X_scaled) loadings_varimax pd.DataFrame( fa_varimax.loadings_, columns[fFactor{i1} for i in range(4)], indexchem_cols ) print(loadings_varimax.round(3))rotationvarimax是最常用的正交旋转。如果因子之间允许相关可以用promax斜交旋转但解释起来更复杂。C题这种成分数据我一般先用varimax不行再换promax。载荷矩阵里绝对值大于0.5的可以认为是显著载荷小于0.3的可以忽略。4.2 灰色关联度小样本分类的备选方案灰色关联度分析适合样本量小、信息不完全的情况。C题里如果某些玻璃类型样本很少用分类模型容易过拟合这时候灰色关联度可以作为一个补充。它的核心是计算每个样本与参考序列的关联度关联度越大说明样本与参考越相似。参考序列一般选每个类别的均值或者理想样本。计算步骤先无量纲化再求差序列然后算关联系数最后求关联度。听起来绕代码其实不长。注意分辨系数ρ通常取0.5但也可以调ρ越小关联度区分度越大。def grey_relation(data, ref, rho0.5): # data: 样本矩阵ref: 参考序列 data_norm (data - data.min()) / (data.max() - data.min()) ref_norm (ref - data.min()) / (data.max() - data.min()) diff np.abs(data_norm - ref_norm) min_diff diff.min() max_diff diff.max() xi (min_diff rho * max_diff) / (diff rho * max_diff) return xi.mean(axis1) # 假设取第一类样本的均值作为参考 ref df[df[类型] 铅钡][chem_cols].mean() grey_scores grey_relation(df[chem_cols].fillna(0), ref) df[灰色关联度] grey_scores print(df.groupby(类型)[灰色关联度].mean())这段代码先做min-max归一化然后计算每个样本与参考序列的差最后用公式算关联系数并取均值。rho0.5是默认值如果关联度区分不明显可以试0.3或0.7。灰色关联度的结果可以跟分类模型的概率输出做对比如果两者一致性高说明分类结果可信如果差异大就要检查是不是某类样本太少导致模型偏了。4.3 分类模型选型逻辑回归、随机森林还是SVMC题第一问通常要求分类第二问可能要求关联分析。分类模型里逻辑回归可解释性强但假设线性边界随机森林抗过拟合能输出特征重要性SVM在小样本高维数据上表现稳。我一般会三个都跑用交叉验证对比准确率。如果样本量小于100优先SVM或逻辑回归加正则如果样本量几百随机森林更省心。注意分类之前一定要做类别平衡检查如果某类样本只有几个要么用SMOTE过采样要么在交叉验证里用分层抽样。C题的数据里铅钡玻璃和钾玻璃的比例可能不均衡直接跑准确率会虚高。from sklearn.model_selection import cross_val_score, StratifiedKFold from sklearn.ensemble import RandomForestClassifier from sklearn.linear_model import LogisticRegression from sklearn.svm import SVC X df[chem_cols].fillna(0) y df[类型] cv StratifiedKFold(n_splits5, shuffleTrue, random_state42) models { LR: LogisticRegression(max_iter1000), RF: RandomForestClassifier(n_estimators100, random_state42), SVM: SVC(kernelrbf, probabilityTrue) } for name, model in models.items(): scores cross_val_score(model, X, y, cvcv, scoringf1_macro) print(f{name}: {scores.mean():.3f} ± {scores.std():.3f})StratifiedKFold保证每折里类别比例一致避免某折缺类。f1_macro比准确率更适合不平衡数据。如果三个模型差距不大选可解释性强的逻辑回归如果随机森林明显高就用随机森林但记得输出特征重要性方便写论文。SVM的probabilityTrue会慢一些但能输出概率方便做后续分析。5. 避坑与排查C题代码跑不通时先看这五条5.1 现象PCA方差贡献率异常高第一主成分超过95%原因没有做标准化量纲大的变量主导了主成分方向。解决在PCA之前加StandardScaler确保每个变量均值为0、方差为1。如果标准化之后还是很高检查是不是有常量列或者近似常量列这种列方差接近0标准化后会变成NaN或极大值需要提前删除。5.2 现象分类模型准确率很高但混淆矩阵里某一类全错原因类别不平衡模型倾向于预测多数类。解决用class_weightbalanced或者用SMOTE过采样少数类。另外检查交叉验证有没有用分层抽样普通KFold在不平衡数据上会不稳定。5.3 现象灰色关联度所有样本得分几乎一样原因分辨系数ρ取值不合适或者数据没有做无量纲化。解决先确认归一化步骤有没有漏然后调整ρ从0.5试到0.3和0.7看区分度有没有改善。如果还是不行说明参考序列选得不好换一个参考样本或者用类别均值重新算。5.4 现象因子分析报错“矩阵不是正定的”原因变量之间存在完全共线性或者样本数少于变量数。解决先检查相关矩阵把相关系数大于0.95的变量删掉一个。如果样本数确实少减少因子个数或者改用PCA。factor_analyzer对数据要求比sklearn的PCA严格遇到报错不要硬跑先做共线性诊断。5.5 现象Python环境里factor_analyzer装不上pip报错原因这个库依赖scipy和numpy的特定版本版本冲突常见。解决先升级pip然后pip install factor_analyzer --no-deps再手动装依赖。如果还不行用sklearn的PCA替代KMO检验可以用pingouin库的calculate_kmo。环境问题没有后悔药建议用conda建独立环境别在base里折腾。6. 从代码包到论文三个让评委多看两眼的技巧第一个技巧把主成分得分和原始成分放在同一张图里做双标图biplot。双标图能同时展示样本分布和变量方向评委一眼就能看出你的降维结果有没有业务含义。代码不长用matplotlib的quiver画箭头就行。第二个技巧分类模型的混淆矩阵不要只画热力图加上准确率、召回率、F1的数值标注尤其是你重点关心的那一类。第三个技巧敏感性分析不要只调一个参数把主成分个数、分类模型、训练集比例三个因素做交叉对比用表格呈现。哪怕结果波动不大也能说明你考虑过稳健性。# 双标图示例 fig, ax plt.subplots(figsize(10, 8)) ax.scatter(X_pca_fixed[:, 0], X_pca_fixed[:, 1], cy.map({铅钡: 0, 钾: 1}), cmapcoolwarm, alpha0.6) for i, col in enumerate(chem_cols): ax.arrow(0, 0, pca_fixed.components_[0, i]*3, pca_fixed.components_[1, i]*3, colorgreen, alpha0.5, head_width0.1) ax.text(pca_fixed.components_[0, i]*3.2, pca_fixed.components_[1, i]*3.2, col, fontsize8) ax.set_xlabel(PC1) ax.set_ylabel(PC2) plt.show()箭头长度代表变量在主成分上的载荷方向代表正负相关。样本点的颜色区分玻璃类型如果两类样本在PC1或PC2上分得开说明降维有效。这张图放在论文里比单纯放碎石图更有说服力。我自己的习惯是拿到任何一份数学建模的文档加代码先不看代码先把数据读进来跑一遍缺失值和分布心里有数了再去对照代码逻辑。2022年C题的数据不算大但坑不少尤其是缺失值和类别不平衡这两块翻车的人最多。你要是能把上面这些步骤走通再根据自己队伍的特长调整模型拿个不错的成绩是有希望的。希望帮到你。本文还有配套的精品资源点击获取
网站建设高端定制企业官网