Matlab数据挖掘实战：从古代玻璃成分分析到分类预测模型构建

1. 项目概述从一堆玻璃碎片到一段尘封历史去年带队参加国赛拿到C题“古代玻璃制品的成分分析与鉴别”时团队里几个搞算法的同学第一反应是兴奋——数据、分类、预测这不正是我们熟悉的领域吗但当我们真正开始处理那些来自考古现场的、成分复杂且残缺不全的数据时才发现这远不是一个简单的数学建模问题。它更像是一次穿越时空的侦探工作我们手里的数据不是冷冰冰的数字而是古代工匠无意中留下的“化学指纹”我们的任务是通过Matlab这把“放大镜”去解读这些指纹背后的文化迁徙、技术传播与贸易往来。这道题的核心是要求我们根据提供的古代玻璃制品化学成分数据对其所属类型高钾玻璃或铅钡玻璃进行鉴别并进一步分析其风化情况、关联规律甚至推测其亚类划分。这听起来像是一个标准的模式识别问题但难点在于数据本身的特性成分百分比数据存在大量缺失、总和不为100%、且存在显著的理化特性如风化导致某些成分流失。直接套用现成的分类算法效果往往惨不忍睹。我们需要做的是先用化学和考古学的逻辑去理解数据再用数学和编程的工具去处理数据最后构建出稳健的鉴别模型。整个过程就是一次典型的“数据驱动”与“领域知识”深度融合的实战。2. 核心思路与解题框架设计面对这道赛题一个清晰的解题框架是成功的一半。盲目地扎进代码里调参很容易迷失方向。我们的整体思路遵循了“数据理解 - 数据预处理 - 特征工程 - 模型构建 - 结果分析”的经典数据科学流程但每一步都深深烙上了本题特有的印记。2.1 问题拆解与任务定义首先我们必须把赛题冗长的描述转化为几个明确、可执行、可评估的数学任务。题目要求可以分解为以下四个核心子问题分类判别根据玻璃文物的化学成分数据判断其属于“高钾玻璃”还是“铅钡玻璃”。这是典型的二分类问题但样本量小数据质量差是核心挑战。风化分析针对风化前后的玻璃成分数据研究风化前后化学成分含量的变化规律并预测风化前的化学成分含量。这涉及到变化趋势建模、缺失数据预测和敏感性分析。关联分析探究玻璃类型、纹饰、颜色等表面属性与化学成分之间的关联关系。这属于多变量关联分析或统计推断的范畴。亚类划分对高钾玻璃和铅钡玻璃分别进行合理的亚类划分并分析亚类间的化学成分差异。这可以看作是无监督的聚类问题。这四个任务环环相扣。例如准确的分类是后续所有分析的基础而对风化规律的理解又能反过来指导我们如何更合理地填补缺失数据提升分类模型的稳健性。2.2 技术路线选型与工具栈基于以上任务我们选择了Matlab作为核心工具。原因有三其一Matlab在矩阵运算、统计分析和可视化方面具有天然优势处理这种“表格型”化学数据非常高效其二其统计与机器学习工具箱功能完整从传统的判别分析到现代的SVM、决策树都能快速实现原型其三团队对Matlab熟悉能在有限时间内将精力集中于建模本身而非环境配置。我们的技术路线图如下数据预处理层以Matlab的表格table数据类型为核心载体进行缺失值处理基于化学规则填补、归一化应对成分和约束、异常值检测。特征工程层利用Matlab的矩阵运算构造比值特征如K2O/PbO、统计特征如各成分的变异系数、以及基于主成分分析PCA的降维特征。模型构建层分类任务主要尝试支持向量机SVM和Fisher线性判别分析LDA风化预测采用多元线性回归或偏最小二乘回归PLSR亚类划分使用系统聚类法或K-Means聚类。评估与可视化层利用Matlab强大的绘图功能plot,scatter,heatmap进行结果展示和规律挖掘用混淆矩阵、均方误差等指标进行模型评估。注意很多队伍一开始就想上复杂的深度学习模型这在样本量仅百来条的数据集上是极其危险的极易过拟合。我们的原则是“先用简单的模型抓住主要矛盾再用复杂模型精细调整”。3. 数据预处理清洗、填补与重构的艺术原始数据通常是一团乱麻而预处理就是将其梳理成模型能“读懂”的语言的过程。对于古代玻璃成分数据这一步至关重要直接决定了后续所有模型的天花板。3.1 缺失值处理不仅仅是填个平均数数据中存在大量“ND”未检测标记传统做法是直接删除或填均值。但对于化学成分数据这两种方法都可能引入严重偏差。我们的策略是基于化学知识进行条件填补识别缺失模式使用Matlab的ismissing()函数快速定位缺失值。我们发现风化后的玻璃中K2O、Na2O等易溶成分常缺失这与它们易被风化淋滤的化学性质相符。规则填补对于高钾玻璃若K2O缺失但PbO和BaO含量极低则用高钾玻璃样本中K2O的中位数或分位数填补。对于铅钡玻璃若PbO或BaO缺失则用铅钡玻璃样本的相应统计量填补。对于其他成分若同一样本的其他主要成分含量很高则该缺失成分可能原本含量就极低考虑填补为0或一个极小值如0.01。迭代优化初步填补后计算所有成分百分比之和。由于检测误差和未测成分存在总和通常不等于100%。我们采用归一化到100%的方法修正后成分 (原始成分 / 所有成分之和) * 100。这一步在Matlab中一行代码即可完成data_normalized data ./ sum(data, 2) * 100;。% 示例基于类型的条件缺失值填补 % 假设 data 是table包含‘Type’列‘高钾’或‘铅钡’和成分列 high_k_idx strcmp(data.Type, ‘高钾’); lead_barium_idx strcmp(data.Type, ‘铅钡’); % 填补高钾玻璃的K2O缺失 k2o_median_highK median(data.K2O(high_k_idx ~ismissing(data.K2O)), ‘omitnan’); data.K2O(high_k_idx ismissing(data.K2O)) k2o_median_highK; % 填补铅钡玻璃的PbO缺失 pbo_median_lb median(data.PbO(lead_barium_idx ~ismissing(data.PbO)), ‘omitnan’); data.PbO(lead_barium_idx ismissing(data.PbO)) pbo_median_lb; % 成分百分比归一化 component_cols {‘SiO2’, ‘Na2O’, ‘K2O’, ‘CaO’, ‘MgO’, ‘Al2O3’, ‘Fe2O3’, ‘CuO’, ‘PbO’, ‘BaO’}; % 所有成分列名 comp_matrix table2array(data(:, component_cols)); row_sums sum(comp_matrix, 2, ‘omitnan’); comp_matrix_normalized comp_matrix ./ row_sums * 100; data(:, component_cols) array2table(comp_matrix_normalized);3.2 特征构造从原始数据中挖掘“金矿”直接使用原始成分百分比作为特征模型可能难以学习到深层规律。我们需要根据化学知识构造更有判别力的特征。比值特征这是最强的一类特征。例如K2O/(PbOBaO)直接区分高钾和铅钡玻璃的核心指标。高钾玻璃此值远大于1铅钡玻璃则远小于1。SiO2/(Na2OK2O)反映玻璃网络形成体与修饰体的比例与玻璃稳定性相关。(CaOMgO)/Al2O3与玻璃的耐风化性可能有关。风化敏感指数针对风化分析任务构造如(风化前成分 - 风化后成分) / 风化前成分的相对变化率特征可以量化不同成分的抗风化能力。统计特征对于每个样本可以计算其所有成分的熵值成分越均匀熵越大或主要成分集中度前三大成分占比这些特征可能隐含了制作工艺的标准化程度。% 示例构造比值特征 data.K_PbBa_Ratio data.K2O ./ (data.PbO data.BaO eps); % 加eps防止除零 data.Si_Alkali_Ratio data.SiO2 ./ (data.Na2O data.K2O eps); data.CaMg_Al_Ratio (data.CaO data.MgO) ./ (data.Al2O3 eps); % 将无穷大Inf或无效值替换为一个大数或NaN data.K_PbBa_Ratio(isinf(data.K_PbBa_Ratio)) 1000; % 假设高钾玻璃该比值极大3.3 异常样本检测通过箱线图boxplot或马氏距离我们可以找出在化学成分上明显偏离群体的样本。例如一个同时含有较高K2O和PbO的样本可能在分类上会造成混淆。我们需要结合考古背景判断这是一个罕见的过渡类型还是数据记录有误在竞赛中通常将其作为待研究的特殊点或在初步建模时暂时剔除以观察模型性能变化。4. 核心模型构建与实现数据准备就绪后就进入了核心的建模环节。我们针对不同子问题选择了最合适且可解释性强的模型。4.1 任务一玻璃类型分类模型我们对比了三种模型Fisher线性判别分析LDA原理是寻找一个投影方向使得两类样本在该方向上投影的类间散度最大类内散度最小。它的优势在于有严格的数学推导且结果稳定特别适合线性可分的小样本数据。Matlab中通过fitcdiscr函数实现。支持向量机SVM寻找一个最优超平面来分割两类样本对于非线性情况可以使用核函数如高斯核。我们使用fitcsvm函数并通过交叉验证选择惩罚参数C和核参数。决策树易于理解和解释可以直观看到哪些化学成分是主要决策依据如K_PbBa_Ratio 5则判为高钾玻璃。使用fitctree函数。实操步骤数据分割将已知类型的样本按7:3分为训练集和测试集。特征选择使用训练集通过LDA的系数权重或决策树的重要性排序筛选出最重要的3-5个特征如K_PbBa_Ratio,PbO,SiO2。模型训练在训练集上分别训练LDA、SVM和决策树模型。模型评估在测试集上计算准确率、召回率、F1分数并绘制混淆矩阵。特别注意对“铅钡玻璃”的召回率因为其样本数可能较少。% 示例使用LDA进行分类 predictors data{:, {‘K_PbBa_Ratio’, ‘PbO’, ‘SiO2’, ‘BaO’}}; % 选择特征 response data.Type; % 响应变量 % 划分训练集和测试集80%训练20%测试 cv cvpartition(response, ‘HoldOut’, 0.2); trainIdx cv.training; testIdx cv.test; % 训练LDA模型 ldaModel fitcdiscr(predictors(trainIdx, :), response(trainIdx), ‘DiscrimType’, ‘linear’); % 预测并评估 predictedType predict(ldaModel, predictors(testIdx, :)); trueType response(testIdx); % 计算混淆矩阵和准确率 C confusionmat(trueType, predictedType); accuracy sum(diag(C)) / sum(C, ‘all’); fprintf(‘LDA模型测试集准确率%.2f%%\n’, accuracy*100); % 可视化混淆矩阵 confusionchart(C, categories(trueType)); title(‘LDA分类结果混淆矩阵’);实操心得在实际测试中LDA和线性SVM的表现通常不相上下且都优于决策树。因为化学成分与类型间存在较强的线性关系。决策树容易过拟合小样本。我们的策略是以LDA作为基准模型用SVM进行微调最后用决策树的结果来辅助解释。最终提交时可以选择在测试集上表现最稳定的那个模型。4.2 任务二风化规律分析与成分预测这是本题的难点和创新点。风化本质上是化学成分的流失但流失比例并非固定。规律分析定性分析绘制风化前后成分含量的簇状柱状图或折线图直观对比。可以明显看到K2O、Na2O等碱金属氧化物含量在风化后大幅下降而SiO2、Al2O3等骨架成分相对稳定PbO在铅钡玻璃风化后也可能流失。定量分析计算每种成分的平均流失率(风化前-风化后)/风化前。进一步可以按玻璃类型分别计算会发现高钾玻璃的碱金属流失更显著。相关性分析计算风化前后成分含量的相关系数矩阵corrcoef寻找哪些成分的流失是同步发生的。预测模型 目标根据风化后的成分X_weathered预测风化前的成分X_original。多元线性回归假设X_original A * X_weathered b。对于每种成分分别建立回归方程。但问题在于风化后某些成分可能已检测不出缺失且成分间存在共线性。偏最小二乘回归PLSR这是更优的选择。PLSR能同时处理自变量风化后成分和因变量风化前成分的建模特别擅长处理多重共线性数据和样本量少于变量数的情况。Matlab中使用plsregress函数。% 示例使用PLSR预测风化前PbO含量 % 假设 weathered_data 是风化后成分表 original_data 是风化前成分表 X table2array(weathered_data(:, {‘SiO2’, ‘Na2O’, ‘K2O’, ‘CaO’, ‘PbO’, ‘BaO’})); % 风化后特征 Y original_data.PbO; % 要预测的风化前PbO含量 % 选择PLSR成分数潜在变量数 - 通过交叉验证确定 ncomp 3; % 例如选择3个主成分 [Xloadings, Yloadings, Xscores, Yscores, beta, PLSPctVar] plsregress(X, Y, ncomp); % 进行预测 Y_pred [ones(size(X,1),1), X] * beta; % 评估预测效果计算均方根误差(RMSE)和决定系数(R^2) rmse sqrt(mean((Y - Y_pred).^2)); ss_tot sum((Y - mean(Y)).^2); ss_res sum((Y - Y_pred).^2); r2 1 - (ss_res / ss_tot); fprintf(‘PLSR模型预测风化前PbO: RMSE %.2f, R^2 %.4f\n’, rmse, r2); % 绘制预测值与真实值散点图 figure; scatter(Y, Y_pred, ‘filled’); hold on; plot([min(Y), max(Y)], [min(Y), max(Y)], ‘r--‘, ‘LineWidth’, 2); % 绘制yx参考线 xlabel(‘风化前PbO真实值’); ylabel(‘风化前PbO预测值’); title(sprintf(‘PLSR预测效果 (R^2%.3f)’, r2)); grid on;4.3 任务三与四关联分析与亚类划分关联分析探究类型、纹饰、颜色与成分的关系。这更像一个探索性数据分析EDA过程。方法对于类别变量如纹饰、颜色可以使用方差分析ANOVA检验不同纹饰下某化学成分的平均含量是否有显著差异。在Matlab中用anovan函数。也可以绘制分组箱线图进行直观比较。发现我们当时分析出某种特定纹饰可能与较高的PbO含量相关这或许反映了当时工匠的配方偏好。亚类划分这是一个无监督的聚类问题。方法对“高钾玻璃”和“铅钡玻璃”的数据分别进行聚类。常用的有K-Means聚类和系统聚类层次聚类。关键点数据标准化聚类前必须对成分数据进行Z-score标准化zscore函数消除量纲影响。确定聚类数K使用“肘部法则”观察不同K值下簇内误差平方和的变化拐点或“轮廓系数”来帮助确定。解释聚类结果计算每个簇的成分均值剖面图给每个簇一个化学意义上的解释例如“高硅高钾型”、“高铅低钡型”等。% 示例对高钾玻璃进行K-Means聚类 high_k_data data(strcmp(data.Type, ‘高钾’), :); features_for_cluster table2array(high_k_data(:, {‘SiO2’, ‘K2O’, ‘CaO’, ‘Al2O3’})); features_scaled zscore(features_for_cluster); % 标准化 % 使用肘部法则寻找最佳K值 distortions []; for k 1:6 [idx, C, sumd] kmeans(features_scaled, k); distortions(k) sum(sumd); end figure; plot(1:6, distortions, ‘bo-‘); xlabel(‘聚类数量 K’); ylabel(‘簇内误差平方和’); title(‘肘部法则’); grid on; % 假设根据肘部法则选择 K3 bestK 3; [idx, C] kmeans(features_scaled, bestK, ‘Replicates’, 10); % 重复10次以避免局部最优 % 将聚类标签添加回原数据表 high_k_data.Cluster idx; % 分析每个簇的特征 cluster_profile grpstats(high_k_data(:, {‘SiO2’, ‘K2O’, ‘CaO’, ‘Al2O3’, ‘Cluster’}), ‘Cluster’, ‘mean’); disp(cluster_profile);5. 模型优化、验证与结果整合单一的模型往往有局限性我们需要通过集成和交叉验证来提升稳健性并将各个子问题的结果有机整合形成一个完整的分析报告。5.1 分类模型的优化与集成在确定了LDA/SVM作为基模型后我们可以进一步优化特征再精选使用递归特征消除RFE方法配合交叉验证自动找出最优特征子集。Matlab中可以通过自定义循环调用sequentialfs顺序特征选择函数来实现。模型集成采用投票法集成。例如让LDA、线性SVM和高斯核SVM三个模型同时对未知样本进行分类采用“少数服从多数”的原则决定最终类型。这能有效降低单一模型误判的风险。处理类别不平衡如果铅钡玻璃样本较少在训练SVM时可以通过‘Prior’, ‘empirical’参数设置先验概率或使用fitcsvm的‘Weights’参数给少数类样本更高权重。5.2 交叉验证与稳健性评估所有模型尤其是预测模型如风化成分预测必须经过严格的交叉验证以确保其泛化能力。方法使用K折交叉验证K5或10。将数据分成K份轮流将其中一份作为测试集其余作为训练集重复K次最终取性能指标的平均值。Matlab实现对于分类可以使用cvpartition和crossval函数对于回归可以使用crossval结合‘mse’损失函数。% 示例对LDA分类器进行5折交叉验证 cv cvpartition(response, ‘KFold’, 5); % 5折划分 ldaCVModel fitcdiscr(predictors, response, ‘DiscrimType’, ‘linear’, ‘CVPartition’, cv); % 计算交叉验证准确率 cvAccuracy 1 - kfoldLoss(ldaCVModel, ‘LossFun’, ‘classiferror’); fprintf(‘LDA模型5折交叉验证平均准确率%.2f%%\n’, cvAccuracy*100);5.3 结果整合与考古学解释数学建模的终点不是得到一个高精度的模型而是产出有考古学意义的结论。在论文中我们需要将数字结果“翻译”成历史语言分类结论不仅给出鉴别正确率还要分析哪些化学成分是鉴别的关键。例如“我们建立的模型表明K2O/(PbOBaO)比值是区分两类玻璃最有效的指标其阈值约为5这与两类玻璃的基础配方理论完全吻合。”风化规律结论“定量分析显示高钾玻璃中K2O的平均流失率高达60%而SiO2的流失率不足5%。这证实了碱金属氧化物在埋藏环境中极不稳定而玻璃骨架网络相对耐久。基于PLSR建立的预测模型可以为严重风化样品原始成分的复原提供定量参考。”亚类划分结论“通过聚类分析我们将高钾玻璃进一步分为三个亚类I类高硅钾型、II类高钙铝型、III类过渡型。结合出土地点信息发现I类多集中于A地区可能代表了当地的主流工艺II类则与B地区关联暗示了不同的原料来源或技术传统。”6. 常见问题、避坑指南与参赛心得回顾整个备战和解题过程我们踩过不少坑也积累了一些宝贵的经验。6.1 数据处理中的“天坑”坑1忽视成分和约束。直接对原始百分比数据做聚类或回归会因为“定和约束”导致虚假的相关性。务必先进行归一化或使用对数比变换。坑2粗暴处理缺失值。对所有缺失值用同一均值填充会模糊两类玻璃的化学界限。必须分类型、分成分进行有条件填补。坑3过度依赖复杂模型。在不足百条的数据上尝试神经网络几乎必然过拟合。坚持“简单模型优先可解释性优先”的原则。6.2 模型选择与调参陷阱陷阱1不验证就相信。在训练集上表现完美的模型可能在测试集上一塌糊涂。严格区分训练集、验证集和测试集并始终用测试集做最终评估。陷阱2盲目追求高精度。对于分类问题如果铅钡玻璃样本很少即使模型把所有样本都预测为高钾玻璃整体准确率也可能很高。一定要关注召回率、精确率和F1分数特别是对少数类的识别能力。陷阱3忽略模型假设。LDA假设数据服从正态分布且各类协方差相等。虽然对于本题数据有一定鲁棒性但最好先做一下正态性检验如lillietest或直接使用更稳健的线性SVM。6.3 论文写作与可视化要点图表胜千言多用高质量的图表展示你的发现。用散点图矩阵展示特征间关系。用热图展示相关系数矩阵或成分含量矩阵。用聚类树状图展示亚类划分过程。所有图表务必清晰标注坐标轴、图例使用易于区分的颜色和标记。逻辑清晰论文结构应严格对应题目要求。每个部分先简述方法再展示结果图表数据最后给出分析和结论。避免把代码和中间过程堆砌上去。体现思考过程在论文中适当描述你尝试过但最终放弃的方法如某种复杂的填补方法效果不佳并简要说明原因这能体现你思维的深度和严谨性。6.4 团队协作与时间管理分工明确一人主攻数据预处理和特征工程一人主攻模型构建与调优一人负责论文撰写与图表制作。但每天必须集中讨论同步进展。版本控制使用Git或至少用带日期的文件夹来管理代码和论文版本避免混乱。设置里程碑第一天结束前完成数据清洗和探索性分析第二天中午前完成核心模型的初步构建第三天专注于优化、整合和论文写作。留出最后半天进行全文检查和格式调整。这道赛题的魅力在于它完美地诠释了如何用现代计算工具解决古老的科学问题。它考验的不仅仅是编程和数学能力更是跨学科的理解力、对数据的敬畏心以及将复杂问题层层拆解的思维能力。当你看到自己构建的模型能够清晰地将一堆看似混乱的化学数据还原为两条不同的古代技术路线时那种成就感是无与伦比的。最终我们提交的不仅仅是一份答案更是一份用数据和逻辑写就的、关于古代科技与文明的微型研究报告。