1. 从“成分”到“身份”一次数学建模实战的深度复盘去年带队参加高教社杯拿到C题《古代玻璃制品的成分分析与鉴别》时我第一反应是这题有意思但坑也不少。它不像一些纯优化或预测题那样有明确的“标准答案”更像是一次结合了化学、考古学和数据科学的跨学科探索。题目给了你一批古代玻璃文物的化学成分数据让你去分析它们的风化规律、判断它们的类别甚至推测它们的产地。听起来很酷对吧但真正做起来你会发现从一堆百分比数据里挖出有说服力的结论每一步都考验着你对数据的理解、对模型的选择以及最重要的——对问题本质的把握。很多人一上来就埋头跑PCA、聚类结果论文写出来自己都觉得牵强。今天我就以这道题为例抛开那些花哨的术语聊聊我们当时是怎么一步步拆解问题、选择方法并把代码真正跑通的。这不仅仅是一份“参考答案”更是一次完整的数据分析思维和建模流程的实战演练。2. 破题第一步理解数据与定义任务边界拿到数据通常是Excel或CSV格式千万别急着导入MATLAB或Python就开始拟合。第一步也是最重要的一步是像侦探勘察现场一样仔细审视你的数据。2.1 数据初窥与预处理陷阱题目提供的数据通常是各种氧化物如SiO₂, Na₂O, K₂O, CaO等的成分百分比。第一眼你可能会看到很多缺失值NaN或者成分总和严重偏离100%的行。这不是数据错误而是题目故意设置的现实情况——文物风化会导致部分成分流失或转化检测也有误差。我们的处理逻辑是区分有效样本首先剔除那些关键主成分如SiO₂完全缺失或总和异常比如低于80%或高于120%的极端异常样本。这些样本数据质量太差强行分析会引入巨大噪声。处理缺失值对于其他缺失值直接删除行删除是最简单但最奢侈的做法可能损失宝贵信息。更合理的策略是对于风化样品其表面成分已改变内部的“原始”成分是未知的。题目通常要求我们预测风化前的成分。这不能简单用均值填充我们当时的思路是基于未风化的同类玻璃的成分规律如高钾玻璃和铅钡玻璃的成分特征差异建立回归模型来估算风化样品缺失的原始成分。这是一个关键建模点。对于非关键成分的零星缺失可以考虑用该类玻璃高钾/铅钡该成分的中位数或均值进行填充。用均值还是中位数要看数据分布。如果数据有少数极端值中位数更稳健。我们当时会先画箱线图观察。成分数据的特殊性所有成分都是百分比是“闭合数据”。这意味着所有变量之和为100%理论上一个成分的变化会迫使其他成分相应变化。这会导致传统的相关性计算失真伪相关。因此在后续分析前通常需要进行数据转换比如采用中心对数比变换CLR来打破闭合效应这是很多新手会忽略的专业点。注意预处理没有唯一标准答案但每一步操作都必须在论文中阐明理由。比如“我们删除了总和低于85%的样本因为其可能受到严重污染或检测失误保留它们会干扰整体规律分析”这比单纯写“我们进行了数据清洗”要有说服力得多。2.2 明确四个子问题的内在联系这道题通常包含几个环环相扣的子问题玻璃类型鉴别根据成分判断玻璃是高钾玻璃还是铅钡玻璃。这是一个分类问题。风化规律分析分析风化前后成分的变化规律并预测风化前的成分。这是一个回归预测规律挖掘问题。亚类划分对每个大类高钾、铅钡进行细分。这是一个无监督聚类问题。产地推测根据成分分析不同类别玻璃的产地相关性。这涉及到差异性分析和多元统计。它们不是孤立的。例如问题1的分类结果是问题3聚类的基础必须先分开高钾和铅钡再各自聚类。问题2的风化规律分析又能为问题1中鉴别风化样品类型提供辅助依据比如发现风化后铅钡玻璃的PbO流失有特定模式。在建模之初就要有这种全局视角。3. 核心武器库模型选择与组合策略面对不同问题需要挑选合适的模型。下面是我们当时针对每个问题的思考路径和具体实现。3.1 问题一分类——如何区分高钾与铅钡这看似是一个标准的二分类问题但样本量小特征成分多且可能存在多重共线性。直接扔进神经网络容易过拟合。我们尝试并对比了三种方案逻辑回归LR 特征选择为什么用LR模型简单可解释性强能给出特征系数告诉我们哪些成分如PbO, K₂O对分类贡献大。关键步骤由于成分指标多我们先用了LASSO回归进行特征选择。LASSO在拟合逻辑回归的同时会将不重要的特征系数压缩至0从而实现自动筛选。MATLAB的lasso函数或Python的sklearn.linear_model.LogisticRegression(penaltyl1)可以轻松实现。代码要点Python示例from sklearn.linear_model import LogisticRegression from sklearn.feature_selection import SelectFromModel from sklearn.model_selection import train_test_split # 假设 X 是预处理后的成分数据y 是类型标签0/1 X_train, X_test, y_train, y_test train_test_split(X, y, test_size0.2, random_state42) # 使用L1正则化的逻辑回归进行特征选择 lr_l1 LogisticRegression(penaltyl1, solverliblinear, C0.1, random_state42) lr_l1.fit(X_train, y_train) # 查看系数不为零的特征 selected_features X.columns[lr_l1.coef_[0] ! 0] print(fSelected features by LASSO: {selected_features.tolist()}) # 用筛选后的特征重新训练模型 X_train_selected X_train[selected_features] X_test_selected X_test[selected_features] lr_final LogisticRegression(penaltyl2, solverliblinear, random_state42) lr_final.fit(X_train_selected, y_train) accuracy lr_final.score(X_test_selected, y_test) print(fTest Accuracy: {accuracy:.4f})输出结果我们得到了一个简洁的模型发现PbO和K₂O的权重绝对值最大这与化学常识铅钡玻璃含高铅高钾玻璃含高钾完美吻合增强了论文的可信度。支持向量机SVM为什么用SVM在小样本、非线性可分问题上表现稳健。我们使用了高斯核RBF。关键点需要调参惩罚系数C和核参数gamma。我们使用了网格搜索GridSearchCV寻找最优参数。心得SVM的分类准确率可能略高于LR但可解释性差。在数学建模中可解释性有时比那1-2%的准确率提升更重要。我们最终将SVM作为对比模型主模型还是用了LR因为其结论更清晰便于在论文中阐述“根据PbO和K₂O含量可有效区分两类玻璃”。决策树/随机森林优点非线性能力强能给出特征重要性排序。缺点容易过拟合小数据且树模型产生的规则如“PbO 10%”可能过于绝对忽略了成分间的协同效应。我们的用法用随机森林的特征重要性来辅助验证LASSO筛选出的特征是否合理作为双重保险。最终策略我们以LASSO-逻辑回归模型为主用SVM和随机森林的结果进行交叉验证。在论文中我们展示了特征选择的过程、模型的系数表以及交叉验证的准确率形成了一个完整的证据链。3.2 问题二风化规律与成分预测——如何“穿越”时间这是本题的难点和亮点。风化的本质是玻璃表面与环境的化学反应导致某些成分如K₂O, Na₂O流失而某些惰性成分如SiO₂, Al₂O₃相对富集。分析步骤定性规律挖掘分别对高钾和铅钡玻璃计算风化点与无风化点对应成分的均值比风化/无风化。比值明显小于1的是流失成分比值大于1的是相对富集成分。绘制风化前后成分变化的箱线图或雷达图可视化展示。例如我们发现高钾玻璃风化后K₂O大幅下降而SiO₂比例上升铅钡玻璃则是PbO和BaO显著下降。注意这里比较的是比例的变化而不是绝对含量。因为风化后总质量减少各成分百分比之和仍为100%所以某个成分百分比上升不一定代表其绝对量增加可能只是其他成分流失更严重导致的“被动富集”。在论文中需要明确指出这一点体现严谨性。定量预测建模目标根据风化后的成分数据预测其风化前的原始成分。思路将无风化样品的数据作为训练集。假设风化过程对同类玻璃的影响模式是一致的那么对于一件风化文物其原始成分应该最接近于某类无风化样品成分的“修正”版本。方法一基于风化系数的加权调整。计算训练集无风化样品中各类成分的均值向量M_original。对于每个风化样品我们假设其原始成分比例与M_original相似但经过了不同程度的风化流失。我们可以定义一个简单的线性调整模型这是最核心的假设预测原始成分 风化后成分 * 调整因子调整因子怎么来我们可以从训练集中“模拟”风化。但训练集只有原始数据。一个巧妙的思路是将训练集本身视为“未风化”状态那么整个训练集的成分均值M_original可以看作一个“虚拟标准品”的原始成分。对于风化样品我们寻找一个缩放向量k使得k * 风化后成分在成分空间中最接近M_original同时要考虑闭合数据约束各成分调整后之和为100%。这可以转化为一个带约束的优化问题。方法二多元线性回归。 将风化后的各成分作为自变量X将原始成分用同类无风化样品的均值或通过其他方法估算的“理想值”作为目标作为因变量Y建立多元线性回归模型。但这里样本量小自变量多容易过拟合必须使用岭回归Ridge或LASSO回归进行正则化。from sklearn.linear_model import Ridge from sklearn.preprocessing import StandardScaler # 假设 X_weathered 是风化样品数据Y_original 是对应的估算出的原始成分目标值 # 数据标准化非常重要 scaler_X StandardScaler() scaler_Y StandardScaler() X_scaled scaler_X.fit_transform(X_weathered) Y_scaled scaler_Y.fit_transform(Y_original) # 使用岭回归alpha是正则化强度需交叉验证选择 ridge Ridge(alpha1.0) ridge.fit(X_scaled, Y_scaled) # 预测新风化样品的原始成分 X_new_scaled scaler_X.transform(new_weathered_sample) Y_pred_scaled ridge.predict(X_new_scaled) Y_pred scaler_Y.inverse_transform(Y_pred_scaled) # 反标准化得到实际百分比重要提示无论用哪种方法预测出的原始成分各氧化物百分比之和必须为100%。如果预测结果偏离需要进行归一化处理。在论文中必须详细说明你的预测模型原理、假设和归一化后处理步骤。3.3 问题三亚类划分——聚类中的“望闻问切”在完成高钾/铅钡分类后需要对每一类内部进行细分。这属于无监督的聚类分析。方法选择K-Means聚类最常用但需要指定K类别数且对异常值敏感。层次聚类可以生成树状图谱系图直观展示样本间的层次关系有助于判断合适的类别数。DBSCAN基于密度能发现任意形状的簇并识别噪声点适用于数据分布不规则的情况。我们的实操流程数据准备使用经过预处理和转换如CLR转换后的成分数据。务必先剔除风化严重的样品因为其成分已失真会严重干扰聚类。确定最佳簇数K值肘部法则绘制不同K值对应的聚类误差平方和SSE曲线选择曲线拐点肘部对应的K。轮廓系数计算每个样本的轮廓系数并求平均。轮廓系数越接近1聚类效果越好。我们通常绘制不同K值下的平均轮廓系数图选择峰值对应的K。结合业务理解如果从考古学知识得知某类玻璃可能有2-3个主要亚型那么K的选择应在此范围内。我们将数学指标和先验知识结合来确定K。执行聚类与解读运行K-Means算法获取每个样本的簇标签。核心不是得到标签而是解读每个簇的特征。计算每个簇的成分均值剖面与整体均值对比。例如在高钾玻璃中我们可能聚类出一个“高钙高铝”亚类和一个“低钙低铝”亚类这可能对应不同的原料来源或工艺。用主成分分析PCA将高维数据降至2-3维在二维散点图上用不同颜色标记聚类结果可视化展示分类效果。给聚类结果命名根据每个簇的化学成分特征为其起一个具有物理或工艺意义的名称如“高铅钡玻璃”、“钾钙玻璃”等使分析结论更生动。踩坑提醒聚类前一定要做特征缩放标准化。因为成分含量单位都是百分比但范围差异大SiO₂可能70%某些微量元素可能不到1%。如果不标准化量级大的成分会完全主导距离计算使聚类结果失真。我们使用StandardScaler进行Z-score标准化。3.4 问题四产地推测——相关性不等于因果性最后一问通常要求分析成分与产地的关系。注意题目数据可能不直接包含产地信息或者只有部分样品有产地标签。这需要更巧妙的思路。如果数据有产地标签差异性检验对于不同产地的同类玻璃对其各项成分进行非参数检验如Mann-Whitney U检验或Kruskal-Wallis H检验判断哪些成分在产地间存在显著差异。为什么用非参数检验因为成分数据不一定服从正态分布。可视化对差异显著的成分绘制按产地分组的箱线图。判别分析使用线性判别分析LDA看看能否根据成分有效区分产地。LDA还能给出哪些成分对区分产地的贡献最大。如果数据没有产地标签要求你根据成分“推测”这实际上变成了一个无监督的探索性分析。你可以将问题三的聚类结果亚类与文献中已知的古代玻璃产地类型进行比对。例如你通过聚类发现铅钡玻璃中有两个亚类A类铅含量极高30%B类铅含量中等15-25%。结合考古文献你可能会发现A类特征与战国时期某地出土的玻璃器成分高度吻合而B类则与另一地区相近。这时你就可以提出“A亚类可能来源于X地区B亚类可能来源于Y地区”的合理推测并在论文中引用文献证据来支撑。关键这种推测必须基于充分的文献调研和严谨的成分对比不能凭空想象。在数学建模论文中即使没有确凿结论展示这个“文献比对-成分分析”的推理过程本身也是有价值的。4. 代码实现从脚本到可复现的流水线很多论文附上的代码就是一堆零散的脚本别人根本无法运行。我们当时的目标是构建一个清晰、可复现的流水线。4.1 环境与数据结构我们主要使用PythonPandas, NumPy, Scikit-learn, Matplotlib, Seaborn。MATLAB同样强大但Python在数据处理和机器学习库的生态上更丰富。import pandas as pd import numpy as np from sklearn.preprocessing import StandardScaler, LabelEncoder from sklearn.model_selection import train_test_split, GridSearchCV from sklearn.linear_model import LogisticRegression, Ridge, Lasso from sklearn.svm import SVC from sklearn.ensemble import RandomForestClassifier from sklearn.cluster import KMeans from sklearn.decomposition import PCA from sklearn.metrics import accuracy_score, silhouette_score import matplotlib.pyplot as plt import seaborn as sns # 1. 数据加载与初探 df pd.read_excel(古代玻璃数据.xlsx) print(df.info()) print(df.head()) print(df.isnull().sum())4.2 构建模块化函数将关键步骤封装成函数提高代码可读性和复用性。def preprocess_data(df, sum_threshold_low85, sum_threshold_high115): 数据预处理函数 1. 剔除成分总和异常样本 2. 分离风化/无风化样品 3. 初步处理缺失值标记 # 计算各样本成分总和 component_cols [col for col in df.columns if col not in [文物编号, 类型, 风化情况, 产地]] df[sum_components] df[component_cols].sum(axis1, skipnaFalse) # 剔除总和异常样本 df_valid df[(df[sum_components] sum_threshold_low) (df[sum_components] sum_threshold_high)].copy() # 分离数据 df_weathered df_valid[df_valid[风化情况] 风化].copy() df_unweathered df_valid[df_valid[风化情况] 无风化].copy() return df_valid, df_weathered, df_unweathered, component_cols def clr_transformation(data, component_cols): 中心对数比变换 (CLR) 处理闭合数据 # 将0值替换为一个极小值如1e-6避免对数计算错误 data_clr data[component_cols].replace(0, 1e-6) # 计算几何均值 gmean np.exp(np.mean(np.log(data_clr), axis1)) # CLR变换 data_clr np.log(data_clr.div(gmean, axis0)) return data_clr4.3 完整的分析流程示例以分类为例# 主流程 df pd.read_excel(data.xlsx) df_valid, df_weathered, df_unweathered, comp_cols preprocess_data(df) # 准备分类数据使用无风化样品训练所有样品测试预测风化样品类型 df_train df_unweathered.copy() # 假设类型列中 高钾1, 铅钡0 le LabelEncoder() df_train[type_label] le.fit_transform(df_train[类型]) X_train df_train[comp_cols].fillna(df_train[comp_cols].median()) # 用中位数填充训练集缺失值 y_train df_train[type_label] # 特征标准化 scaler StandardScaler() X_train_scaled scaler.fit_transform(X_train) # LASSO特征选择 lasso Lasso(alpha0.01, random_state42) # alpha需调优 lasso.fit(X_train_scaled, y_train) selected_idx np.where(lasso.coef_ ! 0)[0] selected_cols [comp_cols[i] for i in selected_idx] print(fSelected features: {selected_cols}) # 使用筛选后的特征训练最终逻辑回归模型 X_train_selected X_train_scaled[:, selected_idx] lr LogisticRegression(random_state42) lr.fit(X_train_selected, y_train) # 在训练集上评估或使用交叉验证 y_pred_train lr.predict(X_train_selected) acc_train accuracy_score(y_train, y_pred_train) print(fTraining Accuracy: {acc_train:.4f}) # 应用模型预测所有样品包括风化样品 # 注意预测前需要对所有样品使用相同的预处理和特征缩放 df_all df_valid.copy() X_all df_all[comp_cols].fillna(df_all[comp_cols].median()) # 填充方式需与训练集一致 X_all_scaled scaler.transform(X_all) # 使用训练集的scaler进行转换 X_all_selected X_all_scaled[:, selected_idx] df_all[predicted_type] le.inverse_transform(lr.predict(X_all_selected)) df_all[predicted_prob] np.max(lr.predict_proba(X_all_selected), axis1) # 查看预测结果 print(df_all[[文物编号, 类型, 风化情况, predicted_type, predicted_prob]].head(20))4.4 可视化与结果输出将关键结果用图表呈现并保存到文件。# 1. 特征重要性图逻辑回归系数 plt.figure(figsize(10,6)) coef_df pd.DataFrame({feature: selected_cols, coefficient: lr.coef_[0]}) coef_df coef_df.sort_values(coefficient, ascendingFalse) sns.barplot(xcoefficient, yfeature, datacoef_df) plt.title(Logistic Regression Coefficients for Glass Type Classification) plt.tight_layout() plt.savefig(classification_coefficients.png, dpi300) # 2. 聚类结果可视化PCA降维后 pca PCA(n_components2) X_pca pca.fit_transform(X_train_scaled) # 使用标准化后的数据 plt.figure(figsize(10,8)) scatter plt.scatter(X_pca[:, 0], X_pca[:, 1], cy_train, cmapviridis, alpha0.7) plt.xlabel(Principal Component 1) plt.ylabel(Principal Component 2) plt.title(PCA Visualization of Glass Samples (Colored by Type)) plt.colorbar(scatter, labelGlass Type (0:铅钡, 1:高钾)) plt.savefig(pca_visualization.png, dpi300) # 3. 将关键结果保存到Excel result_df df_all[[文物编号, 类型, 风化情况, predicted_type, predicted_prob]] result_df.to_excel(prediction_results.xlsx, indexFalse)5. 论文撰写与模型表达的技巧代码跑出结果只是成功了一半如何清晰地写在论文里是另一半。5.1 模型描述部分不要只写“我们使用了逻辑回归模型”而要写 “针对玻璃类型鉴别这一二分类问题考虑到样本量有限且特征间可能存在多重共线性我们选择了可解释性强的逻辑回归模型。为筛选关键化学成分指标我们首先引入了LASSOL1正则化回归进行特征选择其优化目标函数为min(∑(y_i - logit(βX_i))^2 λ∑|β_j|)。通过交叉验证确定正则化强度λ后得到非零系数对应的特征子集见表1。随后基于该特征子集构建逻辑回归模型其形式为log(P/(1-P)) β_0 β_1*x_1 ... β_p*x_p其中P为属于高钾玻璃的概率。”配上表格表1 LASSO特征选择结果氧化物成分回归系数是否被选中SiO₂-0.02否PbO2.15是K₂O-1.87是.........5.2 结果分析部分不要只写“准确率达到95%”而要写 “模型在测试集上准确率达到95.2%混淆矩阵见表2显示仅有一个铅钡玻璃样本被误判为高钾玻璃。进一步分析该误判样本发现其K₂O含量8.5%处于两类玻璃的临界区域且PbO含量5.1%显著低于同类铅钡玻璃的平均水平15%这可能是由于该文物经历了特殊的风化过程或原料不纯所致反映了模型在边界案例上的不确定性也与实际情况相符。”配上图表准确率曲线、混淆矩阵热力图、特征贡献度条形图、聚类散点图、风化前后成分对比雷达图等。一图胜千言。5.3 灵敏度分析与模型检验这是拿高分的关键。展示你的模型不是“黑箱”你思考过它的稳健性。改变预处理方法如果不用中位数填充缺失值而用KNN填充结果变化大吗改变模型参数LASSO的λ值、SVM的C和gamma在合理范围内变动准确率是否稳定交叉验证使用5折或10折交叉验证报告平均准确率及其标准差证明模型性能不是偶然。与简单方法对比比如如果只根据PbO是否大于10%来分类准确率有多少你的复杂模型提升有多大在论文中专门用一小节展示这些分析标题可以是“4.4 模型稳健性检验”。6. 那些我们踩过的坑与宝贵经验忽视数据闭合性最初直接对百分比数据做相关性分析和聚类结果很奇怪。后来查阅文献才知道“成分数据”需要特殊处理如CLR变换调整后规律立刻清晰了。教训看到百分比形式的成分数据第一反应就应该是“闭合数据”并考虑相应的统计方法。盲目追求复杂模型一开始用了XGBoost准确率确实比逻辑回归高一点点但模型复杂难以解释特征重要性虽然给出了排序但无法像逻辑回归系数那样明确指示成分是正相关还是负相关。评委更看重清晰的逻辑和可解释的结论。教训数学建模不是机器学习竞赛在保证性能的前提下模型简洁性与可解释性优先。聚类前未剔除风化样品第一次聚类结果一团糟各类别特征不明显。后来才意识到风化样品的成分已经严重偏离其原始组成把它们和未风化样品混在一起聚类相当于把“病人”和“健康人”混在一起分类毫无意义。教训数据分析前一定要根据问题的物理/化学背景对样本进行合理分组。代码与论文脱节论文里写的是A方法代码里实现的是B方法。答辩时被问到细节支支吾吾。教训写论文时所有的图表、数据都直接从代码运行结果中生成和导出确保绝对一致。给代码加上详细的注释并保存好每次运行的环境配置可以用requirements.txt。对“预测风化前成分”的理解偏差最初试图用一个全局模型预测所有风化样品。后来才想明白高钾玻璃和铅钡玻璃的风化机制不同必须分开建模。对于铅钡玻璃甚至还要考虑PbO和BaO流失的耦合关系。教训细分问题场景针对不同子群体建立特异性模型往往比一个通用模型效果更好。回过头看这道题之所以经典是因为它完美模拟了一个真实的数据科学项目流程从脏数据清洗、业务理解考古化学知识、到方法选择有监督、无监督、回归、模型实现与调优、结果解释与报告。它考察的不仅仅是编程和调包能力更是系统性的问题解决思维和严谨的科学表述能力。希望这份超详细的复盘能让你在下一次面对类似问题时多一份从容少踩一个坑。记住好的建模始于对数据的敬畏成于对问题的深刻理解。
数学建模实战：从数据预处理到模型选择，解析古代玻璃成分分析全流程
1. 从“成分”到“身份”一次数学建模实战的深度复盘去年带队参加高教社杯拿到C题《古代玻璃制品的成分分析与鉴别》时我第一反应是这题有意思但坑也不少。它不像一些纯优化或预测题那样有明确的“标准答案”更像是一次结合了化学、考古学和数据科学的跨学科探索。题目给了你一批古代玻璃文物的化学成分数据让你去分析它们的风化规律、判断它们的类别甚至推测它们的产地。听起来很酷对吧但真正做起来你会发现从一堆百分比数据里挖出有说服力的结论每一步都考验着你对数据的理解、对模型的选择以及最重要的——对问题本质的把握。很多人一上来就埋头跑PCA、聚类结果论文写出来自己都觉得牵强。今天我就以这道题为例抛开那些花哨的术语聊聊我们当时是怎么一步步拆解问题、选择方法并把代码真正跑通的。这不仅仅是一份“参考答案”更是一次完整的数据分析思维和建模流程的实战演练。2. 破题第一步理解数据与定义任务边界拿到数据通常是Excel或CSV格式千万别急着导入MATLAB或Python就开始拟合。第一步也是最重要的一步是像侦探勘察现场一样仔细审视你的数据。2.1 数据初窥与预处理陷阱题目提供的数据通常是各种氧化物如SiO₂, Na₂O, K₂O, CaO等的成分百分比。第一眼你可能会看到很多缺失值NaN或者成分总和严重偏离100%的行。这不是数据错误而是题目故意设置的现实情况——文物风化会导致部分成分流失或转化检测也有误差。我们的处理逻辑是区分有效样本首先剔除那些关键主成分如SiO₂完全缺失或总和异常比如低于80%或高于120%的极端异常样本。这些样本数据质量太差强行分析会引入巨大噪声。处理缺失值对于其他缺失值直接删除行删除是最简单但最奢侈的做法可能损失宝贵信息。更合理的策略是对于风化样品其表面成分已改变内部的“原始”成分是未知的。题目通常要求我们预测风化前的成分。这不能简单用均值填充我们当时的思路是基于未风化的同类玻璃的成分规律如高钾玻璃和铅钡玻璃的成分特征差异建立回归模型来估算风化样品缺失的原始成分。这是一个关键建模点。对于非关键成分的零星缺失可以考虑用该类玻璃高钾/铅钡该成分的中位数或均值进行填充。用均值还是中位数要看数据分布。如果数据有少数极端值中位数更稳健。我们当时会先画箱线图观察。成分数据的特殊性所有成分都是百分比是“闭合数据”。这意味着所有变量之和为100%理论上一个成分的变化会迫使其他成分相应变化。这会导致传统的相关性计算失真伪相关。因此在后续分析前通常需要进行数据转换比如采用中心对数比变换CLR来打破闭合效应这是很多新手会忽略的专业点。注意预处理没有唯一标准答案但每一步操作都必须在论文中阐明理由。比如“我们删除了总和低于85%的样本因为其可能受到严重污染或检测失误保留它们会干扰整体规律分析”这比单纯写“我们进行了数据清洗”要有说服力得多。2.2 明确四个子问题的内在联系这道题通常包含几个环环相扣的子问题玻璃类型鉴别根据成分判断玻璃是高钾玻璃还是铅钡玻璃。这是一个分类问题。风化规律分析分析风化前后成分的变化规律并预测风化前的成分。这是一个回归预测规律挖掘问题。亚类划分对每个大类高钾、铅钡进行细分。这是一个无监督聚类问题。产地推测根据成分分析不同类别玻璃的产地相关性。这涉及到差异性分析和多元统计。它们不是孤立的。例如问题1的分类结果是问题3聚类的基础必须先分开高钾和铅钡再各自聚类。问题2的风化规律分析又能为问题1中鉴别风化样品类型提供辅助依据比如发现风化后铅钡玻璃的PbO流失有特定模式。在建模之初就要有这种全局视角。3. 核心武器库模型选择与组合策略面对不同问题需要挑选合适的模型。下面是我们当时针对每个问题的思考路径和具体实现。3.1 问题一分类——如何区分高钾与铅钡这看似是一个标准的二分类问题但样本量小特征成分多且可能存在多重共线性。直接扔进神经网络容易过拟合。我们尝试并对比了三种方案逻辑回归LR 特征选择为什么用LR模型简单可解释性强能给出特征系数告诉我们哪些成分如PbO, K₂O对分类贡献大。关键步骤由于成分指标多我们先用了LASSO回归进行特征选择。LASSO在拟合逻辑回归的同时会将不重要的特征系数压缩至0从而实现自动筛选。MATLAB的lasso函数或Python的sklearn.linear_model.LogisticRegression(penaltyl1)可以轻松实现。代码要点Python示例from sklearn.linear_model import LogisticRegression from sklearn.feature_selection import SelectFromModel from sklearn.model_selection import train_test_split # 假设 X 是预处理后的成分数据y 是类型标签0/1 X_train, X_test, y_train, y_test train_test_split(X, y, test_size0.2, random_state42) # 使用L1正则化的逻辑回归进行特征选择 lr_l1 LogisticRegression(penaltyl1, solverliblinear, C0.1, random_state42) lr_l1.fit(X_train, y_train) # 查看系数不为零的特征 selected_features X.columns[lr_l1.coef_[0] ! 0] print(fSelected features by LASSO: {selected_features.tolist()}) # 用筛选后的特征重新训练模型 X_train_selected X_train[selected_features] X_test_selected X_test[selected_features] lr_final LogisticRegression(penaltyl2, solverliblinear, random_state42) lr_final.fit(X_train_selected, y_train) accuracy lr_final.score(X_test_selected, y_test) print(fTest Accuracy: {accuracy:.4f})输出结果我们得到了一个简洁的模型发现PbO和K₂O的权重绝对值最大这与化学常识铅钡玻璃含高铅高钾玻璃含高钾完美吻合增强了论文的可信度。支持向量机SVM为什么用SVM在小样本、非线性可分问题上表现稳健。我们使用了高斯核RBF。关键点需要调参惩罚系数C和核参数gamma。我们使用了网格搜索GridSearchCV寻找最优参数。心得SVM的分类准确率可能略高于LR但可解释性差。在数学建模中可解释性有时比那1-2%的准确率提升更重要。我们最终将SVM作为对比模型主模型还是用了LR因为其结论更清晰便于在论文中阐述“根据PbO和K₂O含量可有效区分两类玻璃”。决策树/随机森林优点非线性能力强能给出特征重要性排序。缺点容易过拟合小数据且树模型产生的规则如“PbO 10%”可能过于绝对忽略了成分间的协同效应。我们的用法用随机森林的特征重要性来辅助验证LASSO筛选出的特征是否合理作为双重保险。最终策略我们以LASSO-逻辑回归模型为主用SVM和随机森林的结果进行交叉验证。在论文中我们展示了特征选择的过程、模型的系数表以及交叉验证的准确率形成了一个完整的证据链。3.2 问题二风化规律与成分预测——如何“穿越”时间这是本题的难点和亮点。风化的本质是玻璃表面与环境的化学反应导致某些成分如K₂O, Na₂O流失而某些惰性成分如SiO₂, Al₂O₃相对富集。分析步骤定性规律挖掘分别对高钾和铅钡玻璃计算风化点与无风化点对应成分的均值比风化/无风化。比值明显小于1的是流失成分比值大于1的是相对富集成分。绘制风化前后成分变化的箱线图或雷达图可视化展示。例如我们发现高钾玻璃风化后K₂O大幅下降而SiO₂比例上升铅钡玻璃则是PbO和BaO显著下降。注意这里比较的是比例的变化而不是绝对含量。因为风化后总质量减少各成分百分比之和仍为100%所以某个成分百分比上升不一定代表其绝对量增加可能只是其他成分流失更严重导致的“被动富集”。在论文中需要明确指出这一点体现严谨性。定量预测建模目标根据风化后的成分数据预测其风化前的原始成分。思路将无风化样品的数据作为训练集。假设风化过程对同类玻璃的影响模式是一致的那么对于一件风化文物其原始成分应该最接近于某类无风化样品成分的“修正”版本。方法一基于风化系数的加权调整。计算训练集无风化样品中各类成分的均值向量M_original。对于每个风化样品我们假设其原始成分比例与M_original相似但经过了不同程度的风化流失。我们可以定义一个简单的线性调整模型这是最核心的假设预测原始成分 风化后成分 * 调整因子调整因子怎么来我们可以从训练集中“模拟”风化。但训练集只有原始数据。一个巧妙的思路是将训练集本身视为“未风化”状态那么整个训练集的成分均值M_original可以看作一个“虚拟标准品”的原始成分。对于风化样品我们寻找一个缩放向量k使得k * 风化后成分在成分空间中最接近M_original同时要考虑闭合数据约束各成分调整后之和为100%。这可以转化为一个带约束的优化问题。方法二多元线性回归。 将风化后的各成分作为自变量X将原始成分用同类无风化样品的均值或通过其他方法估算的“理想值”作为目标作为因变量Y建立多元线性回归模型。但这里样本量小自变量多容易过拟合必须使用岭回归Ridge或LASSO回归进行正则化。from sklearn.linear_model import Ridge from sklearn.preprocessing import StandardScaler # 假设 X_weathered 是风化样品数据Y_original 是对应的估算出的原始成分目标值 # 数据标准化非常重要 scaler_X StandardScaler() scaler_Y StandardScaler() X_scaled scaler_X.fit_transform(X_weathered) Y_scaled scaler_Y.fit_transform(Y_original) # 使用岭回归alpha是正则化强度需交叉验证选择 ridge Ridge(alpha1.0) ridge.fit(X_scaled, Y_scaled) # 预测新风化样品的原始成分 X_new_scaled scaler_X.transform(new_weathered_sample) Y_pred_scaled ridge.predict(X_new_scaled) Y_pred scaler_Y.inverse_transform(Y_pred_scaled) # 反标准化得到实际百分比重要提示无论用哪种方法预测出的原始成分各氧化物百分比之和必须为100%。如果预测结果偏离需要进行归一化处理。在论文中必须详细说明你的预测模型原理、假设和归一化后处理步骤。3.3 问题三亚类划分——聚类中的“望闻问切”在完成高钾/铅钡分类后需要对每一类内部进行细分。这属于无监督的聚类分析。方法选择K-Means聚类最常用但需要指定K类别数且对异常值敏感。层次聚类可以生成树状图谱系图直观展示样本间的层次关系有助于判断合适的类别数。DBSCAN基于密度能发现任意形状的簇并识别噪声点适用于数据分布不规则的情况。我们的实操流程数据准备使用经过预处理和转换如CLR转换后的成分数据。务必先剔除风化严重的样品因为其成分已失真会严重干扰聚类。确定最佳簇数K值肘部法则绘制不同K值对应的聚类误差平方和SSE曲线选择曲线拐点肘部对应的K。轮廓系数计算每个样本的轮廓系数并求平均。轮廓系数越接近1聚类效果越好。我们通常绘制不同K值下的平均轮廓系数图选择峰值对应的K。结合业务理解如果从考古学知识得知某类玻璃可能有2-3个主要亚型那么K的选择应在此范围内。我们将数学指标和先验知识结合来确定K。执行聚类与解读运行K-Means算法获取每个样本的簇标签。核心不是得到标签而是解读每个簇的特征。计算每个簇的成分均值剖面与整体均值对比。例如在高钾玻璃中我们可能聚类出一个“高钙高铝”亚类和一个“低钙低铝”亚类这可能对应不同的原料来源或工艺。用主成分分析PCA将高维数据降至2-3维在二维散点图上用不同颜色标记聚类结果可视化展示分类效果。给聚类结果命名根据每个簇的化学成分特征为其起一个具有物理或工艺意义的名称如“高铅钡玻璃”、“钾钙玻璃”等使分析结论更生动。踩坑提醒聚类前一定要做特征缩放标准化。因为成分含量单位都是百分比但范围差异大SiO₂可能70%某些微量元素可能不到1%。如果不标准化量级大的成分会完全主导距离计算使聚类结果失真。我们使用StandardScaler进行Z-score标准化。3.4 问题四产地推测——相关性不等于因果性最后一问通常要求分析成分与产地的关系。注意题目数据可能不直接包含产地信息或者只有部分样品有产地标签。这需要更巧妙的思路。如果数据有产地标签差异性检验对于不同产地的同类玻璃对其各项成分进行非参数检验如Mann-Whitney U检验或Kruskal-Wallis H检验判断哪些成分在产地间存在显著差异。为什么用非参数检验因为成分数据不一定服从正态分布。可视化对差异显著的成分绘制按产地分组的箱线图。判别分析使用线性判别分析LDA看看能否根据成分有效区分产地。LDA还能给出哪些成分对区分产地的贡献最大。如果数据没有产地标签要求你根据成分“推测”这实际上变成了一个无监督的探索性分析。你可以将问题三的聚类结果亚类与文献中已知的古代玻璃产地类型进行比对。例如你通过聚类发现铅钡玻璃中有两个亚类A类铅含量极高30%B类铅含量中等15-25%。结合考古文献你可能会发现A类特征与战国时期某地出土的玻璃器成分高度吻合而B类则与另一地区相近。这时你就可以提出“A亚类可能来源于X地区B亚类可能来源于Y地区”的合理推测并在论文中引用文献证据来支撑。关键这种推测必须基于充分的文献调研和严谨的成分对比不能凭空想象。在数学建模论文中即使没有确凿结论展示这个“文献比对-成分分析”的推理过程本身也是有价值的。4. 代码实现从脚本到可复现的流水线很多论文附上的代码就是一堆零散的脚本别人根本无法运行。我们当时的目标是构建一个清晰、可复现的流水线。4.1 环境与数据结构我们主要使用PythonPandas, NumPy, Scikit-learn, Matplotlib, Seaborn。MATLAB同样强大但Python在数据处理和机器学习库的生态上更丰富。import pandas as pd import numpy as np from sklearn.preprocessing import StandardScaler, LabelEncoder from sklearn.model_selection import train_test_split, GridSearchCV from sklearn.linear_model import LogisticRegression, Ridge, Lasso from sklearn.svm import SVC from sklearn.ensemble import RandomForestClassifier from sklearn.cluster import KMeans from sklearn.decomposition import PCA from sklearn.metrics import accuracy_score, silhouette_score import matplotlib.pyplot as plt import seaborn as sns # 1. 数据加载与初探 df pd.read_excel(古代玻璃数据.xlsx) print(df.info()) print(df.head()) print(df.isnull().sum())4.2 构建模块化函数将关键步骤封装成函数提高代码可读性和复用性。def preprocess_data(df, sum_threshold_low85, sum_threshold_high115): 数据预处理函数 1. 剔除成分总和异常样本 2. 分离风化/无风化样品 3. 初步处理缺失值标记 # 计算各样本成分总和 component_cols [col for col in df.columns if col not in [文物编号, 类型, 风化情况, 产地]] df[sum_components] df[component_cols].sum(axis1, skipnaFalse) # 剔除总和异常样本 df_valid df[(df[sum_components] sum_threshold_low) (df[sum_components] sum_threshold_high)].copy() # 分离数据 df_weathered df_valid[df_valid[风化情况] 风化].copy() df_unweathered df_valid[df_valid[风化情况] 无风化].copy() return df_valid, df_weathered, df_unweathered, component_cols def clr_transformation(data, component_cols): 中心对数比变换 (CLR) 处理闭合数据 # 将0值替换为一个极小值如1e-6避免对数计算错误 data_clr data[component_cols].replace(0, 1e-6) # 计算几何均值 gmean np.exp(np.mean(np.log(data_clr), axis1)) # CLR变换 data_clr np.log(data_clr.div(gmean, axis0)) return data_clr4.3 完整的分析流程示例以分类为例# 主流程 df pd.read_excel(data.xlsx) df_valid, df_weathered, df_unweathered, comp_cols preprocess_data(df) # 准备分类数据使用无风化样品训练所有样品测试预测风化样品类型 df_train df_unweathered.copy() # 假设类型列中 高钾1, 铅钡0 le LabelEncoder() df_train[type_label] le.fit_transform(df_train[类型]) X_train df_train[comp_cols].fillna(df_train[comp_cols].median()) # 用中位数填充训练集缺失值 y_train df_train[type_label] # 特征标准化 scaler StandardScaler() X_train_scaled scaler.fit_transform(X_train) # LASSO特征选择 lasso Lasso(alpha0.01, random_state42) # alpha需调优 lasso.fit(X_train_scaled, y_train) selected_idx np.where(lasso.coef_ ! 0)[0] selected_cols [comp_cols[i] for i in selected_idx] print(fSelected features: {selected_cols}) # 使用筛选后的特征训练最终逻辑回归模型 X_train_selected X_train_scaled[:, selected_idx] lr LogisticRegression(random_state42) lr.fit(X_train_selected, y_train) # 在训练集上评估或使用交叉验证 y_pred_train lr.predict(X_train_selected) acc_train accuracy_score(y_train, y_pred_train) print(fTraining Accuracy: {acc_train:.4f}) # 应用模型预测所有样品包括风化样品 # 注意预测前需要对所有样品使用相同的预处理和特征缩放 df_all df_valid.copy() X_all df_all[comp_cols].fillna(df_all[comp_cols].median()) # 填充方式需与训练集一致 X_all_scaled scaler.transform(X_all) # 使用训练集的scaler进行转换 X_all_selected X_all_scaled[:, selected_idx] df_all[predicted_type] le.inverse_transform(lr.predict(X_all_selected)) df_all[predicted_prob] np.max(lr.predict_proba(X_all_selected), axis1) # 查看预测结果 print(df_all[[文物编号, 类型, 风化情况, predicted_type, predicted_prob]].head(20))4.4 可视化与结果输出将关键结果用图表呈现并保存到文件。# 1. 特征重要性图逻辑回归系数 plt.figure(figsize(10,6)) coef_df pd.DataFrame({feature: selected_cols, coefficient: lr.coef_[0]}) coef_df coef_df.sort_values(coefficient, ascendingFalse) sns.barplot(xcoefficient, yfeature, datacoef_df) plt.title(Logistic Regression Coefficients for Glass Type Classification) plt.tight_layout() plt.savefig(classification_coefficients.png, dpi300) # 2. 聚类结果可视化PCA降维后 pca PCA(n_components2) X_pca pca.fit_transform(X_train_scaled) # 使用标准化后的数据 plt.figure(figsize(10,8)) scatter plt.scatter(X_pca[:, 0], X_pca[:, 1], cy_train, cmapviridis, alpha0.7) plt.xlabel(Principal Component 1) plt.ylabel(Principal Component 2) plt.title(PCA Visualization of Glass Samples (Colored by Type)) plt.colorbar(scatter, labelGlass Type (0:铅钡, 1:高钾)) plt.savefig(pca_visualization.png, dpi300) # 3. 将关键结果保存到Excel result_df df_all[[文物编号, 类型, 风化情况, predicted_type, predicted_prob]] result_df.to_excel(prediction_results.xlsx, indexFalse)5. 论文撰写与模型表达的技巧代码跑出结果只是成功了一半如何清晰地写在论文里是另一半。5.1 模型描述部分不要只写“我们使用了逻辑回归模型”而要写 “针对玻璃类型鉴别这一二分类问题考虑到样本量有限且特征间可能存在多重共线性我们选择了可解释性强的逻辑回归模型。为筛选关键化学成分指标我们首先引入了LASSOL1正则化回归进行特征选择其优化目标函数为min(∑(y_i - logit(βX_i))^2 λ∑|β_j|)。通过交叉验证确定正则化强度λ后得到非零系数对应的特征子集见表1。随后基于该特征子集构建逻辑回归模型其形式为log(P/(1-P)) β_0 β_1*x_1 ... β_p*x_p其中P为属于高钾玻璃的概率。”配上表格表1 LASSO特征选择结果氧化物成分回归系数是否被选中SiO₂-0.02否PbO2.15是K₂O-1.87是.........5.2 结果分析部分不要只写“准确率达到95%”而要写 “模型在测试集上准确率达到95.2%混淆矩阵见表2显示仅有一个铅钡玻璃样本被误判为高钾玻璃。进一步分析该误判样本发现其K₂O含量8.5%处于两类玻璃的临界区域且PbO含量5.1%显著低于同类铅钡玻璃的平均水平15%这可能是由于该文物经历了特殊的风化过程或原料不纯所致反映了模型在边界案例上的不确定性也与实际情况相符。”配上图表准确率曲线、混淆矩阵热力图、特征贡献度条形图、聚类散点图、风化前后成分对比雷达图等。一图胜千言。5.3 灵敏度分析与模型检验这是拿高分的关键。展示你的模型不是“黑箱”你思考过它的稳健性。改变预处理方法如果不用中位数填充缺失值而用KNN填充结果变化大吗改变模型参数LASSO的λ值、SVM的C和gamma在合理范围内变动准确率是否稳定交叉验证使用5折或10折交叉验证报告平均准确率及其标准差证明模型性能不是偶然。与简单方法对比比如如果只根据PbO是否大于10%来分类准确率有多少你的复杂模型提升有多大在论文中专门用一小节展示这些分析标题可以是“4.4 模型稳健性检验”。6. 那些我们踩过的坑与宝贵经验忽视数据闭合性最初直接对百分比数据做相关性分析和聚类结果很奇怪。后来查阅文献才知道“成分数据”需要特殊处理如CLR变换调整后规律立刻清晰了。教训看到百分比形式的成分数据第一反应就应该是“闭合数据”并考虑相应的统计方法。盲目追求复杂模型一开始用了XGBoost准确率确实比逻辑回归高一点点但模型复杂难以解释特征重要性虽然给出了排序但无法像逻辑回归系数那样明确指示成分是正相关还是负相关。评委更看重清晰的逻辑和可解释的结论。教训数学建模不是机器学习竞赛在保证性能的前提下模型简洁性与可解释性优先。聚类前未剔除风化样品第一次聚类结果一团糟各类别特征不明显。后来才意识到风化样品的成分已经严重偏离其原始组成把它们和未风化样品混在一起聚类相当于把“病人”和“健康人”混在一起分类毫无意义。教训数据分析前一定要根据问题的物理/化学背景对样本进行合理分组。代码与论文脱节论文里写的是A方法代码里实现的是B方法。答辩时被问到细节支支吾吾。教训写论文时所有的图表、数据都直接从代码运行结果中生成和导出确保绝对一致。给代码加上详细的注释并保存好每次运行的环境配置可以用requirements.txt。对“预测风化前成分”的理解偏差最初试图用一个全局模型预测所有风化样品。后来才想明白高钾玻璃和铅钡玻璃的风化机制不同必须分开建模。对于铅钡玻璃甚至还要考虑PbO和BaO流失的耦合关系。教训细分问题场景针对不同子群体建立特异性模型往往比一个通用模型效果更好。回过头看这道题之所以经典是因为它完美模拟了一个真实的数据科学项目流程从脏数据清洗、业务理解考古化学知识、到方法选择有监督、无监督、回归、模型实现与调优、结果解释与报告。它考察的不仅仅是编程和调包能力更是系统性的问题解决思维和严谨的科学表述能力。希望这份超详细的复盘能让你在下一次面对类似问题时多一份从容少踩一个坑。记住好的建模始于对数据的敬畏成于对问题的深刻理解。