1. 项目概述从一道赛题到一套完整的数据分析实战方案拿到“古代玻璃制品的成分分析与鉴别”这个题目很多同学的第一反应可能是懵的。这看起来像是一个考古学或者材料科学的问题跟数学建模有什么关系这正是这道赛题的巧妙之处也是它成为经典的原因。它本质上是一个多源、高维、不完整数据的分类与溯源问题完美地将数学统计、机器学习与一个具体的交叉学科应用场景结合了起来。简单来说题目提供了一批古代玻璃文物样本的化学成分检测数据如二氧化硅、氧化钠、氧化钾等氧化物的含量百分比以及这些样本的部分背景信息如出土环境、风化情况等。你的核心任务是利用这些数据通过数学建模的方法回答一系列科学问题如何根据成分区分玻璃的类型如高钾玻璃、铅钡玻璃如何判断一个严重风化的玻璃片原本属于哪种类型不同化学成分的玻璃在风化规律上有什么不同能否根据成分推测文物的产地或制作工艺这绝不仅仅是套用一个现成的分类算法那么简单。数据里充满了陷阱成分数据是“定和约束”数据所有氧化物含量加起来是100%一个成分的变化必然引起其他成分的变化存在大量的缺失值且变量间存在严重的多重共线性。直接把这组数据扔进SVM或者随机森林效果往往很差甚至得到违背化学常识的结果。因此这道题考察的是你数据预处理、特征工程、模型选择与结果解释的全链条能力。而Python凭借其强大的科学计算库如Pandas, NumPy和机器学习库如Scikit-learn, Statsmodels成为了实现这一系列想法最得心应手的工具。接下来我将以这道C题为例拆解从数据清洗到模型构建再到结果可视化的完整流程并提供可直接运行、修改的Python代码源码。无论你是备战数模的新手还是希望提升数据分析实战能力的朋友这篇详解都能提供一条清晰的路径。2. 核心思路拆解如何将考古问题转化为数学模型面对一个跨学科问题建立数学模型的第一步是完成“问题翻译”。我们需要把模糊的自然语言描述转化为清晰的、可量化的数学任务。2.1 问题一玻璃类型的成分分析与判别任务本质这是一个有监督的分类问题。我们拥有已知类型的玻璃样本训练集需要建立一个模型能够根据未知样本的化学成分预测其类型高钾或铅钡。关键难点与思路数据非独立性成分数据是定和约束Compositional Data。直接使用欧氏距离或相关系数会得出错误结论。例如二氧化硅含量从70%上升到80%看似变化很大但这可能仅仅是因为其他成分如氧化钠减少了而非二氧化硅的绝对量增加。处理这类数据通常需要进行对数比变换将数据从单纯形空间映射到欧氏空间。常用的有加性对数比变换ALR或等距对数比变换ILR。特征选择与降维氧化物种类多达十几种且很多存在强相关性如氧化钾和氧化钠可能此消彼长。直接全部放入模型会导致维度灾难和过拟合。我们需要通过主成分分析PCA或因子分析在去除相关性的同时提取主要信息维度。注意PCA应在对数比变换后的数据上进行。模型选择由于样本量通常不大国赛数据一般几十到上百个应选择适合小样本、可解释性强的模型。线性判别分析LDA和逻辑回归Logistic Regression是很好的起点它们能给出每个特征成分对分类的贡献度便于物理解释。如果线性模型效果不佳再考虑支持向量机SVM或简单的决策树。注意很多论文直接使用原始百分比数据进行K-Means聚类来“探索”分类这从统计学上是错误的。对于成分数据聚类应在变换后的空间进行或者使用专门针对成分数据的聚类方法如基于Aitchison距离。2.2 问题二风化玻璃的类型鉴别任务本质这是一个数据缺失下的分类问题或者说是“迁移学习”问题。风化玻璃的某些易流失成分如氧化钾、氧化钠含量严重缺失或异常但其他稳定成分如二氧化硅、氧化铝可能保留较好。关键难点与思路识别风化影响首先需要定量分析风化对每种化学成分的影响。可以比较同类型玻璃中风化与未风化样本各成分含量的分布差异如使用箱线图、T检验。找出风化后系统性升高的成分如二氧化硅因其他成分流失而相对富集和系统性降低的成分如碱金属氧化物。构建稳健特征放弃受风化影响剧烈的成分作为主要判别特征。转而依赖那些在风化过程中相对稳定的成分或者构建新的“抗风化”特征比率例如SiO2/(Al2O3Fe2O3)硅铝铁比这个比值可能更能反映玻璃的基础基质受风化影响较小。模型调整与预测使用问题一中构建的模型时需要剔除或填补风化敏感特征。可以采用两种策略策略一删除在预测风化样本时直接从特征向量中移除已识别出的风化敏感变量仅使用稳定特征进行预测。这要求问题一的模型在构建时就有对应的“稳定特征子集”版本。策略二填补基于未风化数据建立风化敏感成分与稳定成分的回归模型。对于风化样本先用稳定成分预测出其原始敏感成分的可能范围再进行填补最后用完整特征的模型预测。这种方法不确定性更大。2.3 问题三风化程度与成分关联分析任务本质这是一个相关性分析与预测问题。需要探究化学成分如何影响玻璃的抗风化能力。关键难点与思路定义“风化程度”指标这是建模的关键。不能简单用“是/否”风化。可以构建一个连续的风化指数例如风化指数 (风化后SiO2含量 - 同类未风化SiO2平均含量) / 同类未风化SiO2平均含量或者使用主成分分析将多种成分在风化前后的变化综合成一个或两个“风化主成分”得分作为指标。建立关联模型全局分析将所有玻璃样本含类型信息的风化指数作为因变量其原始化学成分作为自变量进行多元线性回归或岭回归解决共线性。通过回归系数的大小和正负判断哪种成分有助于抗风化系数为负即该成分高则风化指数低哪种成分加剧风化。分类型分析对高钾玻璃和铅钡玻璃分别建立上述回归模型。比较两个模型中相同成分的系数是否显著不同。这可以回答“不同类型玻璃的风化机理是否不同”。可视化验证绘制关键成分与风化指数的散点图并按玻璃类型着色可以直观地验证模型结论。2.4 问题四玻璃亚类划分与规律分析任务本质这是一个无监督的聚类分析问题并结合聚类结果进行模式总结。关键难点与思路聚类方法选择在完成对数比变换和PCA降维后可以在主成分得分上进行聚类。K-Means和层次聚类Hierarchical Clustering是常用方法。关键在于确定最佳聚类数K。可以使用肘部法则看误差平方和SSE的拐点或轮廓系数来辅助选择。聚类特征选择不一定使用全部成分进行聚类。可以结合问题三的发现选取与风化行为、或与类型判别密切相关的关键成分进行聚类这样得到的亚类可能具有更明确的化学或考古学意义。规律分析与解释聚类完成后需要描述每个亚类的特征成分剖面计算每个亚类在各化学成分上的均值绘制雷达图或条形图形成“成分指纹”。关联背景将亚类标签与文物的其他信息如出土地区、颜色、纹饰进行交叉分析使用卡方检验或可视化探索亚类与这些背景信息之间是否存在显著关联从而为考古学研究提供线索。3. 数据预处理与特征工程实战理论思路清晰后我们进入实战环节。数据预处理的质量直接决定了模型的上限。这是最繁琐但也最能体现建模者功底的一步。3.1 原始数据清洗与探索假设我们有一个glass_data.csv文件包含字段ID,类型,风化,SiO2,Na2O,K2O,CaO,MgO,Al2O3,Fe2O3,CuO,PbO,BaO... 等以及出土地区等信息。import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns from scipy import stats # 1. 加载数据 df pd.read_csv(glass_data.csv) print(数据形状:, df.shape) print(\n前5行数据:) print(df.head()) print(\n数据基本信息:) print(df.info()) print(\n描述性统计:) print(df.describe()) # 2. 处理缺失值 # 成分数据中的缺失可能以“ND”未检出或空值形式存在 df.replace(ND, np.nan, inplaceTrue) # 将成分列转换为数值型 comp_cols [SiO2, Na2O, K2O, CaO, MgO, Al2O3, Fe2O3, CuO, PbO, BaO] # 根据实际数据调整 for col in comp_cols: df[col] pd.to_numeric(df[col], errorscoerce) print(f\n缺失值统计:\n{df[comp_cols].isnull().sum()}) # 3. 缺失值填补策略需谨慎 # 策略1对于同一类型、同一风化状态的样本用中位数填补 df_filled df.copy() for col in comp_cols: for glass_type in df[类型].unique(): for weather_state in [0, 1]: # 假设0未风化1风化 mask (df[类型] glass_type) (df[风化] weather_state) median_val df.loc[mask, col].median() df_filled.loc[mask df_filled[col].isnull(), col] median_val # 策略2对于风化样本的易流失成分如K2O, Na2O如果缺失过多考虑不填补而是将其作为风化的标志在特征工程中处理。 # 检查填补后是否还有缺失 print(f\n填补后缺失值统计:\n{df_filled[comp_cols].isnull().sum().sum()}) # 4. 数据探索可视化 # 4.1 成分分布箱线图按类型和风化 fig, axes plt.subplots(2, 5, figsize(20, 8)) # 假设有10种成分 axes axes.ravel() for idx, col in enumerate(comp_cols[:10]): # 画前10种 sns.boxplot(x类型, ycol, hue风化, datadf_filled, axaxes[idx]) axes[idx].set_title(f{col} Distribution) axes[idx].tick_params(axisx, rotation45) plt.tight_layout() plt.show() # 4.2 成分相关性热图仅使用未风化数据避免风化干扰 df_unweathered df_filled[df_filled[风化] 0] corr_matrix df_unweathered[comp_cols].corr() plt.figure(figsize(10,8)) sns.heatmap(corr_matrix, annotTrue, fmt.2f, cmapcoolwarm, center0) plt.title(成分相关性热图 (未风化样本)) plt.show()实操心得不要盲目填补对于风化样本中完全流失的成分如K2O含量为0填补为0可能比用中位数更合理因为0具有明确的物理意义完全流失。关注异常值箱线图能快速发现异常点。对于极端异常值需要结合化学知识判断是检测误差还是特殊样品决定是否剔除。相关性分析很重要高热图能直观展示成分间的共生或拮抗关系如Na2O和K2O常呈负相关这为后续的特征工程和模型选择需处理多重共线性提供了依据。3.2 成分数据的对数比变换这是处理本题数据最核心、最专业的一步也是很多入门者容易忽略的一步。from sklearn.preprocessing import StandardScaler # 1. 准备数据选择需要变换的成分列确保数据为正且无零值可用一个小值替代零值 comp_data df_filled[comp_cols].copy() # 为防止取对数时出现无穷大将0值替换为一个极小值如1e-6 comp_data comp_data.replace(0, 1e-6) # 2. 加性对数比变换ALR- 选择一种成分为分母通常是含量最高的SiO2 # ALR: log(X_i / X_D) 其中X_D是分母成分 denominator SiO2 # 选择SiO2作为参考成分 alr_data pd.DataFrame() for col in comp_cols: if col ! denominator: alr_data[flog({col}/{denominator})] np.log(comp_data[col] / comp_data[denominator]) print(ALR变换后的数据前5行:) print(alr_data.head()) # 3. 等距对数比变换ILR- 更复杂但性质更优可以使用第三方库如compositions或scikit-bio # 这里演示一个简单的手动实现基于一组正交对比 # 假设我们有三种成分 A, B, C, D (D为参考) # ILR1 sqrt(3/4) * ln( (A*B*C)^(1/3) / D ) # ILR2 sqrt(2/3) * ln( A^(1/2) / B^(1/2) ) # ILR3 sqrt(1/2) * ln( A / B ) # 简化 # 实际中建议使用现成库此处仅示意原理。 # 4. 变换后数据的标准化对于许多模型很重要 scaler StandardScaler() alr_data_scaled pd.DataFrame(scaler.fit_transform(alr_data), columnsalr_data.columns)注意事项分母的选择ALR变换的结果依赖于分母的选择。选择含量高、变化相对稳定的成分如SiO2作为分母通常更稳健。可以尝试不同分母观察对后续分类效果的影响。ILR的优势ILR变换生成的是正交坐标消除了成分间的冗余信息更适用于PCA和基于欧氏距离的模型。如果条件允许优先使用ILR。解释性经过对数比变换后特征的含义变成了“某种成分相对于参考成分的比例的对数”。在解释模型系数时需要回溯到这个含义。3.3 特征选择与降维处理好多重共线性后我们进行降维提取核心信息。from sklearn.decomposition import PCA # 使用ALR变换并标准化后的数据 X alr_data_scaled # 1. 执行PCA pca PCA() X_pca pca.fit_transform(X) # 2. 查看方差解释比例 explained_variance_ratio pca.explained_variance_ratio_ cumulative_variance np.cumsum(explained_variance_ratio) plt.figure(figsize(10, 4)) plt.subplot(1, 2, 1) plt.bar(range(1, len(explained_variance_ratio)1), explained_variance_ratio, alpha0.8) plt.xlabel(Principal Component) plt.ylabel(Explained Variance Ratio) plt.title(Scree Plot) plt.subplot(1, 2, 2) plt.plot(range(1, len(cumulative_variance)1), cumulative_variance, markero) plt.axhline(y0.85, colorr, linestyle--, label85% Variance) plt.xlabel(Number of Principal Components) plt.ylabel(Cumulative Explained Variance) plt.title(Cumulative Variance Plot) plt.legend() plt.tight_layout() plt.show() print(f前3个主成分解释的方差比例: {cumulative_variance[2]:.2%}) # 3. 查看主成分的载荷Loading理解其物理意义 loadings pd.DataFrame(pca.components_.T, columns[fPC{i1} for i in range(X.shape[1])], indexX.columns) print(\n主成分载荷矩阵前3个主成分:) print(loadings.iloc[:, :3]) # 4. 选择保留的主成分数例如累计贡献率85% n_components np.argmax(cumulative_variance 0.85) 1 print(f\n建议保留的主成分数量: {n_components} (解释{cumulative_variance[n_components-1]:.2%}方差)) X_pca_selected X_pca[:, :n_components] # 5. 可视化主成分得分按类型 pca_df pd.DataFrame(X_pca_selected, columns[fPC{i1} for i in range(n_components)]) pca_df[类型] df_filled[类型].values pca_df[风化] df_filled[风化].values plt.figure(figsize(8,6)) sns.scatterplot(xPC1, yPC2, hue类型, style风化, datapca_df, s100) plt.title(PCA Score Plot (PC1 vs PC2)) plt.show()实操心得主成分的意义一定要查看载荷矩阵。例如PC1可能在log(K2O/SiO2)上有很高的正载荷在log(PbO/SiO2)上有很高的负载荷那么PC1本质上可能就代表了“高钾-低铅”与“低钾-高铅”的对立轴这与玻璃类型的划分可能高度相关。降维不是必须的如果主成分的前2-3个已经能解释绝大部分方差并且散点图能很好地区分类别那么后续的分类模型可以基于主成分得分进行。如果主成分区分度不明显可能需要回到原始特征ALR变换后或尝试其他特征构造方法。4. 模型构建、训练与评估数据准备就绪后我们开始针对不同问题构建模型。4.1 问题一模型玻璃类型分类我们使用未风化的样本作为训练集。from sklearn.model_selection import train_test_split, cross_val_score, GridSearchCV from sklearn.discriminant_analysis import LinearDiscriminantAnalysis from sklearn.linear_model import LogisticRegression from sklearn.svm import SVC from sklearn.tree import DecisionTreeClassifier from sklearn.metrics import classification_report, confusion_matrix, accuracy_score from sklearn.pipeline import Pipeline # 准备数据使用未风化样本特征为PCA降维后的结果或ALR特征 df_train df_filled[df_filled[风化] 0].copy() X_raw df_train[comp_cols] y df_train[类型].map({高钾: 0, 铅钡: 1}) # 编码为0/1 # 特征工程管道在实际应用中应将变换应用于整个数据集再拆分此处为演示简化 # 假设我们有一个预先拟合好的PCA模型 pca_fitted (用训练集拟合的) # X_train_pca pca_fitted.transform(apply_alr_and_scale(X_raw)) # 为了流程完整这里演示从原始数据开始在训练集上拟合变换然后评估 from sklearn.preprocessing import StandardScaler def preprocess_for_clf(X): 对成分数据执行ALR变换和标准化 X_no_zero X.replace(0, 1e-6) alr_features [] for col in comp_cols: if col ! SiO2: alr_features.append(np.log(X_no_zero[col] / X_no_zero[SiO2])) alr_df pd.DataFrame(np.column_stack(alr_features), columns[flog({c}/SiO2) for c in comp_cols if c ! SiO2]) scaler StandardScaler() return scaler.fit_transform(alr_df) X_processed preprocess_for_clf(X_raw) # 划分训练集和测试集 X_train, X_test, y_train, y_test train_test_split(X_processed, y, test_size0.2, random_state42, stratifyy) # 1. 线性判别分析 (LDA) lda LinearDiscriminantAnalysis() lda.fit(X_train, y_train) y_pred_lda lda.predict(X_test) print( LDA 分类报告 ) print(classification_report(y_test, y_pred_lda, target_names[高钾, 铅钡])) print(f准确率: {accuracy_score(y_test, y_pred_lda):.4f}) # 查看判别函数的系数用于解释 print(f\nLDA 系数 (对应特征重要性):) coef_df pd.DataFrame(lda.coef_.T, index[flog({c}/SiO2) for c in comp_cols if c ! SiO2], columns[判别函数系数]) print(coef_df) # 2. 逻辑回归 logreg LogisticRegression(max_iter1000, random_state42) logreg.fit(X_train, y_train) y_pred_log logreg.predict(X_test) print(\n 逻辑回归 分类报告 ) print(classification_report(y_test, y_pred_log, target_names[高钾, 铅钡])) # 3. 支持向量机 (SVM) - 使用网格搜索优化参数 svm_pipeline Pipeline([ (clf, SVC(random_state42)) ]) param_grid { clf__C: [0.1, 1, 10], clf__kernel: [linear, rbf], clf__gamma: [scale, auto, 0.01, 0.1] } grid_search GridSearchCV(svm_pipeline, param_grid, cv5, scoringaccuracy, n_jobs-1) grid_search.fit(X_train, y_train) print(f\n SVM 最佳参数: {grid_search.best_params_}) print(f最佳交叉验证准确率: {grid_search.best_score_:.4f}) y_pred_svm grid_search.predict(X_test) print(classification_report(y_test, y_pred_svm, target_names[高钾, 铅钡])) # 4. 模型比较与选择 models {LDA: lda, Logistic Regression: logreg, SVM: grid_search.best_estimator_} for name, model in models.items(): if hasattr(model, predict): y_pred model.predict(X_test) acc accuracy_score(y_test, y_pred) print(f{name} 测试集准确率: {acc:.4f})注意事项数据泄露PCA、标准化等步骤必须在训练集上拟合fit然后同时应用于训练集和测试集transform。绝对不能在整个数据集上先fit再拆分。使用Pipeline可以很好地避免这个问题。模型解释性LDA和逻辑回归的系数可以直接解释。例如log(K2O/SiO2)的系数为正且很大说明该比值越高模型越倾向于预测为高钾玻璃。这符合化学常识。类别不平衡检查两类样本数量是否均衡。如果不均衡需要在评估时使用F1-score、精确率、召回率等指标而非只看准确率并考虑在模型中设置class_weight参数。4.2 问题二模型风化样本鉴别这里演示基于稳定特征的分类策略。# 假设通过之前的分析我们确定以下成分在风化中相对稳定 stable_features [SiO2, Al2O3, Fe2O3, CaO] # 示例需根据实际数据分析确定 # 1. 为分类器构建仅包含稳定特征的训练数据 df_train_stable df_filled[df_filled[风化] 0].copy() X_train_stable_raw df_train_stable[stable_features] y_train_stable df_train_stable[类型].map({高钾: 0, 铅钡: 1}) # 对稳定特征进行适当的预处理如标准化由于不是成分数据可能不需要对数比变换 scaler_stable StandardScaler() X_train_stable scaler_stable.fit_transform(X_train_stable_raw) # 2. 训练一个基于稳定特征的分类器例如逻辑回归 clf_stable LogisticRegression(max_iter1000, random_state42) clf_stable.fit(X_train_stable, y_train_stable) # 3. 预测风化样本 df_weathered df_filled[df_filled[风化] 1].copy() X_pred_stable_raw df_weathered[stable_features] X_pred_stable scaler_stable.transform(X_pred_stable_raw) # 使用训练集的scaler predictions clf_stable.predict(X_pred_stable) df_weathered[预测类型] [高钾 if p 0 else 铅钡 for p in predictions] print(风化样本的预测结果:) print(df_weathered[[ID, 类型, 预测类型]].head()) # 4. 评估如果有部分风化样本已知类型可用于验证 # 假设 df_weathered 中有一部分‘类型’是已知的 if df_weathered[类型].notnull().any(): known_mask df_weathered[类型].notnull() y_true df_weathered.loc[known_mask, 类型].map({高钾: 0, 铅钡: 1}) y_pred df_weathered.loc[known_mask, 预测类型].map({高钾: 0, 铅钡: 1}) print(f\n风化样本已知类型的预测准确率: {accuracy_score(y_true, y_pred):.4f})4.3 问题三模型风化程度回归分析from sklearn.linear_model import Ridge from sklearn.metrics import mean_squared_error, r2_score # 1. 构建风化指数示例使用SiO2的相对变化 df_all df_filled.copy() # 计算每类玻璃未风化样本的SiO2平均含量 mean_sio2_by_type df_all[df_all[风化]0].groupby(类型)[SiO2].mean() df_all[未风化_SiO2_均值] df_all[类型].map(mean_sio2_by_type) df_all[风化指数_SiO2] (df_all[SiO2] - df_all[未风化_SiO2_均值]) / df_all[未风化_SiO2_均值] # 只对风化样本进行回归分析 df_reg df_all[df_all[风化]1].copy() # 特征原始化学成分或ALR变换后的特征 X_reg_raw df_reg[comp_cols] # 目标风化指数 y_reg df_reg[风化指数_SiO2] # 预处理特征例如使用ALR变换 X_reg_processed preprocess_for_clf(X_reg_raw) # 复用之前的函数 # 2. 划分数据集 X_train_reg, X_test_reg, y_train_reg, y_test_reg train_test_split(X_reg_processed, y_reg, test_size0.2, random_state42) # 3. 使用岭回归处理共线性 ridge Ridge(alpha1.0) # alpha是正则化强度 ridge.fit(X_train_reg, y_train_reg) y_pred_reg ridge.predict(X_test_reg) print( 风化指数回归模型岭回归) print(f测试集 R^2 分数: {r2_score(y_test_reg, y_pred_reg):.4f}) print(f测试集均方根误差 (RMSE): {np.sqrt(mean_squared_error(y_test_reg, y_pred_reg)):.4f}) # 4. 查看模型系数分析成分影响 coef_series pd.Series(ridge.coef_, index[flog({c}/SiO2) for c in comp_cols if c ! SiO2]) coef_series_sorted coef_series.sort_values(ascendingFalse) print(\n回归系数绝对值越大影响越大:) print(coef_series_sorted) # 5. 分类型建模比较 print(\n 分类型回归比较 ) for glass_type in df_reg[类型].unique(): df_type df_reg[df_reg[类型]glass_type] if len(df_type) 5: # 样本量太小则跳过 continue X_type preprocess_for_clf(df_type[comp_cols]) y_type df_type[风化指数_SiO2] ridge_type Ridge(alpha1.0).fit(X_type, y_type) # 比较关键成分的系数... # 这里可以提取特定成分的系数进行对比4.4 问题四模型亚类聚类分析from sklearn.cluster import KMeans from sklearn.metrics import silhouette_score # 使用所有未风化样本进行聚类特征使用PCA降维后的主成分得分或ALR特征 df_cluster df_filled[df_filled[风化]0].copy() X_cluster preprocess_for_clf(df_cluster[comp_cols]) # ALR标准化 # 1. 确定最佳聚类数K inertia [] silhouette_scores [] K_range range(2, 8) for k in K_range: kmeans KMeans(n_clustersk, random_state42, n_init10) kmeans.fit(X_cluster) inertia.append(kmeans.inertia_) silhouette_scores.append(silhouette_score(X_cluster, kmeans.labels_)) plt.figure(figsize(12,4)) plt.subplot(1,2,1) plt.plot(K_range, inertia, markero) plt.xlabel(Number of clusters (K)) plt.ylabel(Inertia (SSE)) plt.title(Elbow Method for Optimal K) plt.subplot(1,2,2) plt.plot(K_range, silhouette_scores, markero) plt.xlabel(Number of clusters (K)) plt.ylabel(Silhouette Score) plt.title(Silhouette Score for Optimal K) plt.tight_layout() plt.show() # 2. 选择K例如根据轮廓系数最大或肘部拐点 optimal_k 3 # 假设我们选择3 kmeans_final KMeans(n_clustersoptimal_k, random_state42, n_init10) cluster_labels kmeans_final.fit_predict(X_cluster) df_cluster[亚类] cluster_labels print(f聚类完成共分为 {optimal_k} 个亚类。) print(df_cluster[[ID, 类型, 亚类]].head()) # 3. 分析每个亚类的成分特征 cluster_profiles df_cluster.groupby(亚类)[comp_cols].mean() print(\n各亚类平均成分剖面:) print(cluster_profiles) # 可视化成分剖面雷达图 from math import pi categories comp_cols N len(categories) angles [n / float(N) * 2 * pi for n in range(N)] angles angles[:1] # 闭合 fig, axes plt.subplots(1, optimal_k, subplot_kwdict(projectionpolar), figsize(5*optimal_k, 4)) if optimal_k 1: axes [axes] for idx, (cluster, profile) in enumerate(cluster_profiles.iterrows()): values profile.values.tolist() values values[:1] ax axes[idx] ax.plot(angles, values, linewidth2, linestylesolid) ax.fill(angles, values, alpha0.25) ax.set_xticks(angles[:-1]) ax.set_xticklabels(categories, size8) ax.set_title(f亚类 {cluster} 成分剖面, size12, y1.1) plt.tight_layout() plt.show() # 4. 分析亚类与类型、出土背景的关联 contingency_table_type pd.crosstab(df_cluster[亚类], df_cluster[类型]) print(\n亚类与玻璃类型的列联表:) print(contingency_table_type) # 可以进行卡方检验 chi2, p, dof, ex stats.chi2_contingency(contingency_table_type) print(f卡方检验 p-value: {p:.4f})5. 常见问题、避坑指南与技巧实录在实际操作和论文写作中你会遇到各种各样的问题。以下是我总结的一些高频“坑点”和应对技巧。5.1 数据预处理中的陷阱问题成分数据中存在零值导致对数变换无穷大。解决方法用一个小正数如1e-6, 1e-10或检测限的一半来替换零值。在论文中需要说明这种处理方式及其合理性。对于成分数据更严谨的方法是使用零值替换方法如zCompositionsR包中的cmultRepl函数但Python中需手动实现或寻找类似库。问题缺失值过多尤其是风化样本的某些成分。解决方法不要对所有列使用同一种填补方法。对于明确因风化流失的成分如K2O, Na2O在风化样本中填补为0或一个极小的值可能比用均值/中位数更合理。可以分别对“高钾-未风化”、“高钾-风化”、“铅钡-未风化”、“铅钡-风化”这四组数据分别进行缺失值处理。问题PCA后主成分的含义难以解释。技巧在运行PCA时设置pca PCA(n_componentsNone)然后查看pca.components_载荷矩阵。将载荷绝对值大的几个原始特征找出来结合其正负号给主成分赋予一个物理意义例如“富钾贫铅因子”、“钙镁相关因子”等。这能极大提升论文的可读性。5.2 模型选择与评估的误区问题直接使用原始百分比数据做K-Means聚类或计算欧氏距离。避坑这是成分数据分析中最常见的错误。必须强调对于成分数据任何基于欧氏距离的方法包括K-Means、层次聚类、KNN都必须在对数比变换后的空间进行。在论文中需要引用统计学中关于成分数据的理论如Aitchison几何并说明你进行了对数比变换。问题模型在训练集上准确率很高但一预测新样本就出错。排查数据泄露检查是否在划分训练测试集之前就做了全局的标准化或PCA。务必使用Pipeline或在训练集上fit再应用到测试集。过拟合小样本数据集上使用了过于复杂的模型如深度神经网络、没有正则化的复杂SVM。优先使用简单模型LDA、逻辑回归并采用交叉验证来评估泛化能力。特征工程不一致确保预测新样本时使用的特征工程流程如对数比变换的参考成分、PCA的载荷矩阵、标准化的均值和方差与训练时完全一致。问题分类结果很好但无法从化学或考古学角度解释。技巧模型的可解释性与预测性能同等重要。多使用LDA、逻辑回归、决策树这类可解释模型。在论文中不仅要给出准确率更要分析“是哪些化学成分起到了关键判别作用”并将这个结论与考古学文献中关于古代玻璃配方的知识相互印证。5.3 论文写作与结果呈现技巧技巧一图胜千言。PCA散点图用不同颜色和形状表示类型和风化状态是展示数据整体结构和模型判别效果的利器。成分剖面雷达图/堆叠柱状图用于直观对比不同类别或亚类玻璃的化学成分差异。模型系数条形图横向条形图展示LDA或逻辑回归的系数一目了然地看出哪些成分是正相关哪些是负相关。聚类结果热图将样本行和成分列以热图形式呈现并用聚类标签标注行可以清晰展示亚类内的成分一致性。技巧量化分析避免主观描述。不要说“A成分在两类玻璃中差异很大”要说“独立样本t检验表明A成分在高钾玻璃和铅钡玻璃中的平均含量存在显著差异p 0.01”。不要说“我们的模型效果很好”要说“在测试集上该模型的准确率达到92.5%F1-score为0.93显著优于基线模型准确率70%”。技巧敏感性分析。在论文中增加一个“敏感性分析”小节讨论你的关键选择对结果的影响。例如对数比变换选择不同参考成分如SiO2 vs Al2O3对分类结果的影响有多大PCA保留的主成分数从3个变成4个聚类结果是否稳定对于风化样本使用不同的稳定特征组合预测准确率的变化范围是多少这能体现你思考的全面性和模型的稳健性是加分项。5.4 代码实现与复现性问题代码杂乱无法复现。规范设置随机种子在代码开头使用np.random.seed(42)和random.seed(42)确保每次运行结果一致。封装函数将数据预处理、特征工程等步骤封装成函数使主流程清晰。使用Pipeline将标准化、降维、模型训练等步骤组合成sklearn Pipeline避免数据泄露也便于网格搜索。注释与文档关键步骤和参数选择理由要写清楚注释。提供数据与代码在提交的论文附录或附件中提供清洗后的数据文件和核心代码方便评委复现。最后我想分享的一点个人体会是数学建模竞赛的本质是用数学工具讲一个逻辑自洽、证据充分的故事。对于这道C题你的故事线应该是面对成分数据的内在约束定和约束我采用了严谨的数据变换方法对数比变换来正确处理它为了从高维、相关的特征中提取信息我使用了降维技术PCA并解释了其物理意义针对不同的子问题分类、回归、聚类我选择了最合适且可解释的模型并给出了符合化学和考古学常识的解读我对模型的关键假设和参数进行了测试敏感性分析证明了结论的稳健性。当你把代码、图表和这些文字叙述有机结合起来时一篇优秀的数模论文就诞生了。 成分数据分析实战：从对数比变换到机器学习建模的完整流程
1. 项目概述从一道赛题到一套完整的数据分析实战方案拿到“古代玻璃制品的成分分析与鉴别”这个题目很多同学的第一反应可能是懵的。这看起来像是一个考古学或者材料科学的问题跟数学建模有什么关系这正是这道赛题的巧妙之处也是它成为经典的原因。它本质上是一个多源、高维、不完整数据的分类与溯源问题完美地将数学统计、机器学习与一个具体的交叉学科应用场景结合了起来。简单来说题目提供了一批古代玻璃文物样本的化学成分检测数据如二氧化硅、氧化钠、氧化钾等氧化物的含量百分比以及这些样本的部分背景信息如出土环境、风化情况等。你的核心任务是利用这些数据通过数学建模的方法回答一系列科学问题如何根据成分区分玻璃的类型如高钾玻璃、铅钡玻璃如何判断一个严重风化的玻璃片原本属于哪种类型不同化学成分的玻璃在风化规律上有什么不同能否根据成分推测文物的产地或制作工艺这绝不仅仅是套用一个现成的分类算法那么简单。数据里充满了陷阱成分数据是“定和约束”数据所有氧化物含量加起来是100%一个成分的变化必然引起其他成分的变化存在大量的缺失值且变量间存在严重的多重共线性。直接把这组数据扔进SVM或者随机森林效果往往很差甚至得到违背化学常识的结果。因此这道题考察的是你数据预处理、特征工程、模型选择与结果解释的全链条能力。而Python凭借其强大的科学计算库如Pandas, NumPy和机器学习库如Scikit-learn, Statsmodels成为了实现这一系列想法最得心应手的工具。接下来我将以这道C题为例拆解从数据清洗到模型构建再到结果可视化的完整流程并提供可直接运行、修改的Python代码源码。无论你是备战数模的新手还是希望提升数据分析实战能力的朋友这篇详解都能提供一条清晰的路径。2. 核心思路拆解如何将考古问题转化为数学模型面对一个跨学科问题建立数学模型的第一步是完成“问题翻译”。我们需要把模糊的自然语言描述转化为清晰的、可量化的数学任务。2.1 问题一玻璃类型的成分分析与判别任务本质这是一个有监督的分类问题。我们拥有已知类型的玻璃样本训练集需要建立一个模型能够根据未知样本的化学成分预测其类型高钾或铅钡。关键难点与思路数据非独立性成分数据是定和约束Compositional Data。直接使用欧氏距离或相关系数会得出错误结论。例如二氧化硅含量从70%上升到80%看似变化很大但这可能仅仅是因为其他成分如氧化钠减少了而非二氧化硅的绝对量增加。处理这类数据通常需要进行对数比变换将数据从单纯形空间映射到欧氏空间。常用的有加性对数比变换ALR或等距对数比变换ILR。特征选择与降维氧化物种类多达十几种且很多存在强相关性如氧化钾和氧化钠可能此消彼长。直接全部放入模型会导致维度灾难和过拟合。我们需要通过主成分分析PCA或因子分析在去除相关性的同时提取主要信息维度。注意PCA应在对数比变换后的数据上进行。模型选择由于样本量通常不大国赛数据一般几十到上百个应选择适合小样本、可解释性强的模型。线性判别分析LDA和逻辑回归Logistic Regression是很好的起点它们能给出每个特征成分对分类的贡献度便于物理解释。如果线性模型效果不佳再考虑支持向量机SVM或简单的决策树。注意很多论文直接使用原始百分比数据进行K-Means聚类来“探索”分类这从统计学上是错误的。对于成分数据聚类应在变换后的空间进行或者使用专门针对成分数据的聚类方法如基于Aitchison距离。2.2 问题二风化玻璃的类型鉴别任务本质这是一个数据缺失下的分类问题或者说是“迁移学习”问题。风化玻璃的某些易流失成分如氧化钾、氧化钠含量严重缺失或异常但其他稳定成分如二氧化硅、氧化铝可能保留较好。关键难点与思路识别风化影响首先需要定量分析风化对每种化学成分的影响。可以比较同类型玻璃中风化与未风化样本各成分含量的分布差异如使用箱线图、T检验。找出风化后系统性升高的成分如二氧化硅因其他成分流失而相对富集和系统性降低的成分如碱金属氧化物。构建稳健特征放弃受风化影响剧烈的成分作为主要判别特征。转而依赖那些在风化过程中相对稳定的成分或者构建新的“抗风化”特征比率例如SiO2/(Al2O3Fe2O3)硅铝铁比这个比值可能更能反映玻璃的基础基质受风化影响较小。模型调整与预测使用问题一中构建的模型时需要剔除或填补风化敏感特征。可以采用两种策略策略一删除在预测风化样本时直接从特征向量中移除已识别出的风化敏感变量仅使用稳定特征进行预测。这要求问题一的模型在构建时就有对应的“稳定特征子集”版本。策略二填补基于未风化数据建立风化敏感成分与稳定成分的回归模型。对于风化样本先用稳定成分预测出其原始敏感成分的可能范围再进行填补最后用完整特征的模型预测。这种方法不确定性更大。2.3 问题三风化程度与成分关联分析任务本质这是一个相关性分析与预测问题。需要探究化学成分如何影响玻璃的抗风化能力。关键难点与思路定义“风化程度”指标这是建模的关键。不能简单用“是/否”风化。可以构建一个连续的风化指数例如风化指数 (风化后SiO2含量 - 同类未风化SiO2平均含量) / 同类未风化SiO2平均含量或者使用主成分分析将多种成分在风化前后的变化综合成一个或两个“风化主成分”得分作为指标。建立关联模型全局分析将所有玻璃样本含类型信息的风化指数作为因变量其原始化学成分作为自变量进行多元线性回归或岭回归解决共线性。通过回归系数的大小和正负判断哪种成分有助于抗风化系数为负即该成分高则风化指数低哪种成分加剧风化。分类型分析对高钾玻璃和铅钡玻璃分别建立上述回归模型。比较两个模型中相同成分的系数是否显著不同。这可以回答“不同类型玻璃的风化机理是否不同”。可视化验证绘制关键成分与风化指数的散点图并按玻璃类型着色可以直观地验证模型结论。2.4 问题四玻璃亚类划分与规律分析任务本质这是一个无监督的聚类分析问题并结合聚类结果进行模式总结。关键难点与思路聚类方法选择在完成对数比变换和PCA降维后可以在主成分得分上进行聚类。K-Means和层次聚类Hierarchical Clustering是常用方法。关键在于确定最佳聚类数K。可以使用肘部法则看误差平方和SSE的拐点或轮廓系数来辅助选择。聚类特征选择不一定使用全部成分进行聚类。可以结合问题三的发现选取与风化行为、或与类型判别密切相关的关键成分进行聚类这样得到的亚类可能具有更明确的化学或考古学意义。规律分析与解释聚类完成后需要描述每个亚类的特征成分剖面计算每个亚类在各化学成分上的均值绘制雷达图或条形图形成“成分指纹”。关联背景将亚类标签与文物的其他信息如出土地区、颜色、纹饰进行交叉分析使用卡方检验或可视化探索亚类与这些背景信息之间是否存在显著关联从而为考古学研究提供线索。3. 数据预处理与特征工程实战理论思路清晰后我们进入实战环节。数据预处理的质量直接决定了模型的上限。这是最繁琐但也最能体现建模者功底的一步。3.1 原始数据清洗与探索假设我们有一个glass_data.csv文件包含字段ID,类型,风化,SiO2,Na2O,K2O,CaO,MgO,Al2O3,Fe2O3,CuO,PbO,BaO... 等以及出土地区等信息。import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns from scipy import stats # 1. 加载数据 df pd.read_csv(glass_data.csv) print(数据形状:, df.shape) print(\n前5行数据:) print(df.head()) print(\n数据基本信息:) print(df.info()) print(\n描述性统计:) print(df.describe()) # 2. 处理缺失值 # 成分数据中的缺失可能以“ND”未检出或空值形式存在 df.replace(ND, np.nan, inplaceTrue) # 将成分列转换为数值型 comp_cols [SiO2, Na2O, K2O, CaO, MgO, Al2O3, Fe2O3, CuO, PbO, BaO] # 根据实际数据调整 for col in comp_cols: df[col] pd.to_numeric(df[col], errorscoerce) print(f\n缺失值统计:\n{df[comp_cols].isnull().sum()}) # 3. 缺失值填补策略需谨慎 # 策略1对于同一类型、同一风化状态的样本用中位数填补 df_filled df.copy() for col in comp_cols: for glass_type in df[类型].unique(): for weather_state in [0, 1]: # 假设0未风化1风化 mask (df[类型] glass_type) (df[风化] weather_state) median_val df.loc[mask, col].median() df_filled.loc[mask df_filled[col].isnull(), col] median_val # 策略2对于风化样本的易流失成分如K2O, Na2O如果缺失过多考虑不填补而是将其作为风化的标志在特征工程中处理。 # 检查填补后是否还有缺失 print(f\n填补后缺失值统计:\n{df_filled[comp_cols].isnull().sum().sum()}) # 4. 数据探索可视化 # 4.1 成分分布箱线图按类型和风化 fig, axes plt.subplots(2, 5, figsize(20, 8)) # 假设有10种成分 axes axes.ravel() for idx, col in enumerate(comp_cols[:10]): # 画前10种 sns.boxplot(x类型, ycol, hue风化, datadf_filled, axaxes[idx]) axes[idx].set_title(f{col} Distribution) axes[idx].tick_params(axisx, rotation45) plt.tight_layout() plt.show() # 4.2 成分相关性热图仅使用未风化数据避免风化干扰 df_unweathered df_filled[df_filled[风化] 0] corr_matrix df_unweathered[comp_cols].corr() plt.figure(figsize(10,8)) sns.heatmap(corr_matrix, annotTrue, fmt.2f, cmapcoolwarm, center0) plt.title(成分相关性热图 (未风化样本)) plt.show()实操心得不要盲目填补对于风化样本中完全流失的成分如K2O含量为0填补为0可能比用中位数更合理因为0具有明确的物理意义完全流失。关注异常值箱线图能快速发现异常点。对于极端异常值需要结合化学知识判断是检测误差还是特殊样品决定是否剔除。相关性分析很重要高热图能直观展示成分间的共生或拮抗关系如Na2O和K2O常呈负相关这为后续的特征工程和模型选择需处理多重共线性提供了依据。3.2 成分数据的对数比变换这是处理本题数据最核心、最专业的一步也是很多入门者容易忽略的一步。from sklearn.preprocessing import StandardScaler # 1. 准备数据选择需要变换的成分列确保数据为正且无零值可用一个小值替代零值 comp_data df_filled[comp_cols].copy() # 为防止取对数时出现无穷大将0值替换为一个极小值如1e-6 comp_data comp_data.replace(0, 1e-6) # 2. 加性对数比变换ALR- 选择一种成分为分母通常是含量最高的SiO2 # ALR: log(X_i / X_D) 其中X_D是分母成分 denominator SiO2 # 选择SiO2作为参考成分 alr_data pd.DataFrame() for col in comp_cols: if col ! denominator: alr_data[flog({col}/{denominator})] np.log(comp_data[col] / comp_data[denominator]) print(ALR变换后的数据前5行:) print(alr_data.head()) # 3. 等距对数比变换ILR- 更复杂但性质更优可以使用第三方库如compositions或scikit-bio # 这里演示一个简单的手动实现基于一组正交对比 # 假设我们有三种成分 A, B, C, D (D为参考) # ILR1 sqrt(3/4) * ln( (A*B*C)^(1/3) / D ) # ILR2 sqrt(2/3) * ln( A^(1/2) / B^(1/2) ) # ILR3 sqrt(1/2) * ln( A / B ) # 简化 # 实际中建议使用现成库此处仅示意原理。 # 4. 变换后数据的标准化对于许多模型很重要 scaler StandardScaler() alr_data_scaled pd.DataFrame(scaler.fit_transform(alr_data), columnsalr_data.columns)注意事项分母的选择ALR变换的结果依赖于分母的选择。选择含量高、变化相对稳定的成分如SiO2作为分母通常更稳健。可以尝试不同分母观察对后续分类效果的影响。ILR的优势ILR变换生成的是正交坐标消除了成分间的冗余信息更适用于PCA和基于欧氏距离的模型。如果条件允许优先使用ILR。解释性经过对数比变换后特征的含义变成了“某种成分相对于参考成分的比例的对数”。在解释模型系数时需要回溯到这个含义。3.3 特征选择与降维处理好多重共线性后我们进行降维提取核心信息。from sklearn.decomposition import PCA # 使用ALR变换并标准化后的数据 X alr_data_scaled # 1. 执行PCA pca PCA() X_pca pca.fit_transform(X) # 2. 查看方差解释比例 explained_variance_ratio pca.explained_variance_ratio_ cumulative_variance np.cumsum(explained_variance_ratio) plt.figure(figsize(10, 4)) plt.subplot(1, 2, 1) plt.bar(range(1, len(explained_variance_ratio)1), explained_variance_ratio, alpha0.8) plt.xlabel(Principal Component) plt.ylabel(Explained Variance Ratio) plt.title(Scree Plot) plt.subplot(1, 2, 2) plt.plot(range(1, len(cumulative_variance)1), cumulative_variance, markero) plt.axhline(y0.85, colorr, linestyle--, label85% Variance) plt.xlabel(Number of Principal Components) plt.ylabel(Cumulative Explained Variance) plt.title(Cumulative Variance Plot) plt.legend() plt.tight_layout() plt.show() print(f前3个主成分解释的方差比例: {cumulative_variance[2]:.2%}) # 3. 查看主成分的载荷Loading理解其物理意义 loadings pd.DataFrame(pca.components_.T, columns[fPC{i1} for i in range(X.shape[1])], indexX.columns) print(\n主成分载荷矩阵前3个主成分:) print(loadings.iloc[:, :3]) # 4. 选择保留的主成分数例如累计贡献率85% n_components np.argmax(cumulative_variance 0.85) 1 print(f\n建议保留的主成分数量: {n_components} (解释{cumulative_variance[n_components-1]:.2%}方差)) X_pca_selected X_pca[:, :n_components] # 5. 可视化主成分得分按类型 pca_df pd.DataFrame(X_pca_selected, columns[fPC{i1} for i in range(n_components)]) pca_df[类型] df_filled[类型].values pca_df[风化] df_filled[风化].values plt.figure(figsize(8,6)) sns.scatterplot(xPC1, yPC2, hue类型, style风化, datapca_df, s100) plt.title(PCA Score Plot (PC1 vs PC2)) plt.show()实操心得主成分的意义一定要查看载荷矩阵。例如PC1可能在log(K2O/SiO2)上有很高的正载荷在log(PbO/SiO2)上有很高的负载荷那么PC1本质上可能就代表了“高钾-低铅”与“低钾-高铅”的对立轴这与玻璃类型的划分可能高度相关。降维不是必须的如果主成分的前2-3个已经能解释绝大部分方差并且散点图能很好地区分类别那么后续的分类模型可以基于主成分得分进行。如果主成分区分度不明显可能需要回到原始特征ALR变换后或尝试其他特征构造方法。4. 模型构建、训练与评估数据准备就绪后我们开始针对不同问题构建模型。4.1 问题一模型玻璃类型分类我们使用未风化的样本作为训练集。from sklearn.model_selection import train_test_split, cross_val_score, GridSearchCV from sklearn.discriminant_analysis import LinearDiscriminantAnalysis from sklearn.linear_model import LogisticRegression from sklearn.svm import SVC from sklearn.tree import DecisionTreeClassifier from sklearn.metrics import classification_report, confusion_matrix, accuracy_score from sklearn.pipeline import Pipeline # 准备数据使用未风化样本特征为PCA降维后的结果或ALR特征 df_train df_filled[df_filled[风化] 0].copy() X_raw df_train[comp_cols] y df_train[类型].map({高钾: 0, 铅钡: 1}) # 编码为0/1 # 特征工程管道在实际应用中应将变换应用于整个数据集再拆分此处为演示简化 # 假设我们有一个预先拟合好的PCA模型 pca_fitted (用训练集拟合的) # X_train_pca pca_fitted.transform(apply_alr_and_scale(X_raw)) # 为了流程完整这里演示从原始数据开始在训练集上拟合变换然后评估 from sklearn.preprocessing import StandardScaler def preprocess_for_clf(X): 对成分数据执行ALR变换和标准化 X_no_zero X.replace(0, 1e-6) alr_features [] for col in comp_cols: if col ! SiO2: alr_features.append(np.log(X_no_zero[col] / X_no_zero[SiO2])) alr_df pd.DataFrame(np.column_stack(alr_features), columns[flog({c}/SiO2) for c in comp_cols if c ! SiO2]) scaler StandardScaler() return scaler.fit_transform(alr_df) X_processed preprocess_for_clf(X_raw) # 划分训练集和测试集 X_train, X_test, y_train, y_test train_test_split(X_processed, y, test_size0.2, random_state42, stratifyy) # 1. 线性判别分析 (LDA) lda LinearDiscriminantAnalysis() lda.fit(X_train, y_train) y_pred_lda lda.predict(X_test) print( LDA 分类报告 ) print(classification_report(y_test, y_pred_lda, target_names[高钾, 铅钡])) print(f准确率: {accuracy_score(y_test, y_pred_lda):.4f}) # 查看判别函数的系数用于解释 print(f\nLDA 系数 (对应特征重要性):) coef_df pd.DataFrame(lda.coef_.T, index[flog({c}/SiO2) for c in comp_cols if c ! SiO2], columns[判别函数系数]) print(coef_df) # 2. 逻辑回归 logreg LogisticRegression(max_iter1000, random_state42) logreg.fit(X_train, y_train) y_pred_log logreg.predict(X_test) print(\n 逻辑回归 分类报告 ) print(classification_report(y_test, y_pred_log, target_names[高钾, 铅钡])) # 3. 支持向量机 (SVM) - 使用网格搜索优化参数 svm_pipeline Pipeline([ (clf, SVC(random_state42)) ]) param_grid { clf__C: [0.1, 1, 10], clf__kernel: [linear, rbf], clf__gamma: [scale, auto, 0.01, 0.1] } grid_search GridSearchCV(svm_pipeline, param_grid, cv5, scoringaccuracy, n_jobs-1) grid_search.fit(X_train, y_train) print(f\n SVM 最佳参数: {grid_search.best_params_}) print(f最佳交叉验证准确率: {grid_search.best_score_:.4f}) y_pred_svm grid_search.predict(X_test) print(classification_report(y_test, y_pred_svm, target_names[高钾, 铅钡])) # 4. 模型比较与选择 models {LDA: lda, Logistic Regression: logreg, SVM: grid_search.best_estimator_} for name, model in models.items(): if hasattr(model, predict): y_pred model.predict(X_test) acc accuracy_score(y_test, y_pred) print(f{name} 测试集准确率: {acc:.4f})注意事项数据泄露PCA、标准化等步骤必须在训练集上拟合fit然后同时应用于训练集和测试集transform。绝对不能在整个数据集上先fit再拆分。使用Pipeline可以很好地避免这个问题。模型解释性LDA和逻辑回归的系数可以直接解释。例如log(K2O/SiO2)的系数为正且很大说明该比值越高模型越倾向于预测为高钾玻璃。这符合化学常识。类别不平衡检查两类样本数量是否均衡。如果不均衡需要在评估时使用F1-score、精确率、召回率等指标而非只看准确率并考虑在模型中设置class_weight参数。4.2 问题二模型风化样本鉴别这里演示基于稳定特征的分类策略。# 假设通过之前的分析我们确定以下成分在风化中相对稳定 stable_features [SiO2, Al2O3, Fe2O3, CaO] # 示例需根据实际数据分析确定 # 1. 为分类器构建仅包含稳定特征的训练数据 df_train_stable df_filled[df_filled[风化] 0].copy() X_train_stable_raw df_train_stable[stable_features] y_train_stable df_train_stable[类型].map({高钾: 0, 铅钡: 1}) # 对稳定特征进行适当的预处理如标准化由于不是成分数据可能不需要对数比变换 scaler_stable StandardScaler() X_train_stable scaler_stable.fit_transform(X_train_stable_raw) # 2. 训练一个基于稳定特征的分类器例如逻辑回归 clf_stable LogisticRegression(max_iter1000, random_state42) clf_stable.fit(X_train_stable, y_train_stable) # 3. 预测风化样本 df_weathered df_filled[df_filled[风化] 1].copy() X_pred_stable_raw df_weathered[stable_features] X_pred_stable scaler_stable.transform(X_pred_stable_raw) # 使用训练集的scaler predictions clf_stable.predict(X_pred_stable) df_weathered[预测类型] [高钾 if p 0 else 铅钡 for p in predictions] print(风化样本的预测结果:) print(df_weathered[[ID, 类型, 预测类型]].head()) # 4. 评估如果有部分风化样本已知类型可用于验证 # 假设 df_weathered 中有一部分‘类型’是已知的 if df_weathered[类型].notnull().any(): known_mask df_weathered[类型].notnull() y_true df_weathered.loc[known_mask, 类型].map({高钾: 0, 铅钡: 1}) y_pred df_weathered.loc[known_mask, 预测类型].map({高钾: 0, 铅钡: 1}) print(f\n风化样本已知类型的预测准确率: {accuracy_score(y_true, y_pred):.4f})4.3 问题三模型风化程度回归分析from sklearn.linear_model import Ridge from sklearn.metrics import mean_squared_error, r2_score # 1. 构建风化指数示例使用SiO2的相对变化 df_all df_filled.copy() # 计算每类玻璃未风化样本的SiO2平均含量 mean_sio2_by_type df_all[df_all[风化]0].groupby(类型)[SiO2].mean() df_all[未风化_SiO2_均值] df_all[类型].map(mean_sio2_by_type) df_all[风化指数_SiO2] (df_all[SiO2] - df_all[未风化_SiO2_均值]) / df_all[未风化_SiO2_均值] # 只对风化样本进行回归分析 df_reg df_all[df_all[风化]1].copy() # 特征原始化学成分或ALR变换后的特征 X_reg_raw df_reg[comp_cols] # 目标风化指数 y_reg df_reg[风化指数_SiO2] # 预处理特征例如使用ALR变换 X_reg_processed preprocess_for_clf(X_reg_raw) # 复用之前的函数 # 2. 划分数据集 X_train_reg, X_test_reg, y_train_reg, y_test_reg train_test_split(X_reg_processed, y_reg, test_size0.2, random_state42) # 3. 使用岭回归处理共线性 ridge Ridge(alpha1.0) # alpha是正则化强度 ridge.fit(X_train_reg, y_train_reg) y_pred_reg ridge.predict(X_test_reg) print( 风化指数回归模型岭回归) print(f测试集 R^2 分数: {r2_score(y_test_reg, y_pred_reg):.4f}) print(f测试集均方根误差 (RMSE): {np.sqrt(mean_squared_error(y_test_reg, y_pred_reg)):.4f}) # 4. 查看模型系数分析成分影响 coef_series pd.Series(ridge.coef_, index[flog({c}/SiO2) for c in comp_cols if c ! SiO2]) coef_series_sorted coef_series.sort_values(ascendingFalse) print(\n回归系数绝对值越大影响越大:) print(coef_series_sorted) # 5. 分类型建模比较 print(\n 分类型回归比较 ) for glass_type in df_reg[类型].unique(): df_type df_reg[df_reg[类型]glass_type] if len(df_type) 5: # 样本量太小则跳过 continue X_type preprocess_for_clf(df_type[comp_cols]) y_type df_type[风化指数_SiO2] ridge_type Ridge(alpha1.0).fit(X_type, y_type) # 比较关键成分的系数... # 这里可以提取特定成分的系数进行对比4.4 问题四模型亚类聚类分析from sklearn.cluster import KMeans from sklearn.metrics import silhouette_score # 使用所有未风化样本进行聚类特征使用PCA降维后的主成分得分或ALR特征 df_cluster df_filled[df_filled[风化]0].copy() X_cluster preprocess_for_clf(df_cluster[comp_cols]) # ALR标准化 # 1. 确定最佳聚类数K inertia [] silhouette_scores [] K_range range(2, 8) for k in K_range: kmeans KMeans(n_clustersk, random_state42, n_init10) kmeans.fit(X_cluster) inertia.append(kmeans.inertia_) silhouette_scores.append(silhouette_score(X_cluster, kmeans.labels_)) plt.figure(figsize(12,4)) plt.subplot(1,2,1) plt.plot(K_range, inertia, markero) plt.xlabel(Number of clusters (K)) plt.ylabel(Inertia (SSE)) plt.title(Elbow Method for Optimal K) plt.subplot(1,2,2) plt.plot(K_range, silhouette_scores, markero) plt.xlabel(Number of clusters (K)) plt.ylabel(Silhouette Score) plt.title(Silhouette Score for Optimal K) plt.tight_layout() plt.show() # 2. 选择K例如根据轮廓系数最大或肘部拐点 optimal_k 3 # 假设我们选择3 kmeans_final KMeans(n_clustersoptimal_k, random_state42, n_init10) cluster_labels kmeans_final.fit_predict(X_cluster) df_cluster[亚类] cluster_labels print(f聚类完成共分为 {optimal_k} 个亚类。) print(df_cluster[[ID, 类型, 亚类]].head()) # 3. 分析每个亚类的成分特征 cluster_profiles df_cluster.groupby(亚类)[comp_cols].mean() print(\n各亚类平均成分剖面:) print(cluster_profiles) # 可视化成分剖面雷达图 from math import pi categories comp_cols N len(categories) angles [n / float(N) * 2 * pi for n in range(N)] angles angles[:1] # 闭合 fig, axes plt.subplots(1, optimal_k, subplot_kwdict(projectionpolar), figsize(5*optimal_k, 4)) if optimal_k 1: axes [axes] for idx, (cluster, profile) in enumerate(cluster_profiles.iterrows()): values profile.values.tolist() values values[:1] ax axes[idx] ax.plot(angles, values, linewidth2, linestylesolid) ax.fill(angles, values, alpha0.25) ax.set_xticks(angles[:-1]) ax.set_xticklabels(categories, size8) ax.set_title(f亚类 {cluster} 成分剖面, size12, y1.1) plt.tight_layout() plt.show() # 4. 分析亚类与类型、出土背景的关联 contingency_table_type pd.crosstab(df_cluster[亚类], df_cluster[类型]) print(\n亚类与玻璃类型的列联表:) print(contingency_table_type) # 可以进行卡方检验 chi2, p, dof, ex stats.chi2_contingency(contingency_table_type) print(f卡方检验 p-value: {p:.4f})5. 常见问题、避坑指南与技巧实录在实际操作和论文写作中你会遇到各种各样的问题。以下是我总结的一些高频“坑点”和应对技巧。5.1 数据预处理中的陷阱问题成分数据中存在零值导致对数变换无穷大。解决方法用一个小正数如1e-6, 1e-10或检测限的一半来替换零值。在论文中需要说明这种处理方式及其合理性。对于成分数据更严谨的方法是使用零值替换方法如zCompositionsR包中的cmultRepl函数但Python中需手动实现或寻找类似库。问题缺失值过多尤其是风化样本的某些成分。解决方法不要对所有列使用同一种填补方法。对于明确因风化流失的成分如K2O, Na2O在风化样本中填补为0或一个极小的值可能比用均值/中位数更合理。可以分别对“高钾-未风化”、“高钾-风化”、“铅钡-未风化”、“铅钡-风化”这四组数据分别进行缺失值处理。问题PCA后主成分的含义难以解释。技巧在运行PCA时设置pca PCA(n_componentsNone)然后查看pca.components_载荷矩阵。将载荷绝对值大的几个原始特征找出来结合其正负号给主成分赋予一个物理意义例如“富钾贫铅因子”、“钙镁相关因子”等。这能极大提升论文的可读性。5.2 模型选择与评估的误区问题直接使用原始百分比数据做K-Means聚类或计算欧氏距离。避坑这是成分数据分析中最常见的错误。必须强调对于成分数据任何基于欧氏距离的方法包括K-Means、层次聚类、KNN都必须在对数比变换后的空间进行。在论文中需要引用统计学中关于成分数据的理论如Aitchison几何并说明你进行了对数比变换。问题模型在训练集上准确率很高但一预测新样本就出错。排查数据泄露检查是否在划分训练测试集之前就做了全局的标准化或PCA。务必使用Pipeline或在训练集上fit再应用到测试集。过拟合小样本数据集上使用了过于复杂的模型如深度神经网络、没有正则化的复杂SVM。优先使用简单模型LDA、逻辑回归并采用交叉验证来评估泛化能力。特征工程不一致确保预测新样本时使用的特征工程流程如对数比变换的参考成分、PCA的载荷矩阵、标准化的均值和方差与训练时完全一致。问题分类结果很好但无法从化学或考古学角度解释。技巧模型的可解释性与预测性能同等重要。多使用LDA、逻辑回归、决策树这类可解释模型。在论文中不仅要给出准确率更要分析“是哪些化学成分起到了关键判别作用”并将这个结论与考古学文献中关于古代玻璃配方的知识相互印证。5.3 论文写作与结果呈现技巧技巧一图胜千言。PCA散点图用不同颜色和形状表示类型和风化状态是展示数据整体结构和模型判别效果的利器。成分剖面雷达图/堆叠柱状图用于直观对比不同类别或亚类玻璃的化学成分差异。模型系数条形图横向条形图展示LDA或逻辑回归的系数一目了然地看出哪些成分是正相关哪些是负相关。聚类结果热图将样本行和成分列以热图形式呈现并用聚类标签标注行可以清晰展示亚类内的成分一致性。技巧量化分析避免主观描述。不要说“A成分在两类玻璃中差异很大”要说“独立样本t检验表明A成分在高钾玻璃和铅钡玻璃中的平均含量存在显著差异p 0.01”。不要说“我们的模型效果很好”要说“在测试集上该模型的准确率达到92.5%F1-score为0.93显著优于基线模型准确率70%”。技巧敏感性分析。在论文中增加一个“敏感性分析”小节讨论你的关键选择对结果的影响。例如对数比变换选择不同参考成分如SiO2 vs Al2O3对分类结果的影响有多大PCA保留的主成分数从3个变成4个聚类结果是否稳定对于风化样本使用不同的稳定特征组合预测准确率的变化范围是多少这能体现你思考的全面性和模型的稳健性是加分项。5.4 代码实现与复现性问题代码杂乱无法复现。规范设置随机种子在代码开头使用np.random.seed(42)和random.seed(42)确保每次运行结果一致。封装函数将数据预处理、特征工程等步骤封装成函数使主流程清晰。使用Pipeline将标准化、降维、模型训练等步骤组合成sklearn Pipeline避免数据泄露也便于网格搜索。注释与文档关键步骤和参数选择理由要写清楚注释。提供数据与代码在提交的论文附录或附件中提供清洗后的数据文件和核心代码方便评委复现。最后我想分享的一点个人体会是数学建模竞赛的本质是用数学工具讲一个逻辑自洽、证据充分的故事。对于这道C题你的故事线应该是面对成分数据的内在约束定和约束我采用了严谨的数据变换方法对数比变换来正确处理它为了从高维、相关的特征中提取信息我使用了降维技术PCA并解释了其物理意义针对不同的子问题分类、回归、聚类我选择了最合适且可解释的模型并给出了符合化学和考古学常识的解读我对模型的关键假设和参数进行了测试敏感性分析证明了结论的稳健性。当你把代码、图表和这些文字叙述有机结合起来时一篇优秀的数模论文就诞生了。