1. 项目背景与核心任务拆解看到这个标题很多参加过数学建模竞赛的同学应该会心一笑。2022年高教社杯全国大学生数学建模竞赛的C题“古代玻璃制品的成分分析与鉴别”可以说是当年最具挑战性和趣味性的题目之一。它巧妙地将考古学、材料科学与数据分析结合起来要求参赛者扮演一个“科技考古”的角色通过对一批古代玻璃文物化学成分数据的分析来回答一系列科学问题。题目给出的数据是一批古代玻璃制品的化学成分检测报告包含了如二氧化硅(SiO₂)、氧化钠(Na₂O)、氧化钾(K₂O)等十几种氧化物的含量百分比。而第一问的核心任务通常是对数据进行预处理和初步的统计分析为后续的分类、风化规律研究等打下基础。具体来说第一问往往要求我们数据清洗与预处理处理原始数据中的缺失值、异常值将成分数据转换为合适的格式如和为100%的约束。描述性统计分析计算各类玻璃如高钾玻璃、铅钡玻璃在不同化学成分上的基本统计量均值、标准差等形成初步认知。差异性分析运用统计检验方法如t检验、方差分析定量分析不同类别玻璃在化学成分上是否存在显著差异。相关性分析探索不同化学成分之间的关联关系为理解玻璃的配方体系提供线索。这些步骤听起来基础但却是整个建模工作的基石。处理得当后续的模型构建事半功倍处理不当则可能带着系统性偏差一路错下去。接下来我将结合当年解题的实际经验手把手拆解第一问的完整代码实现与核心思路重点会放在为什么这么做以及实操中容易踩的坑上。2. 数据读取与初步审视避开第一个“天坑”拿到竞赛数据的第一件事绝对不是急着跑模型而是静下心来像考古学家清理文物一样仔细“清理”和审视数据。竞赛提供的通常是Excel或CSV文件我们以CSV为例。2.1 读取数据与列名处理% 假设数据文件名为 glass_data.csv data readtable(glass_data.csv, VariableNamingRule, preserve); % 查看数据前几行和基本信息 head(data) summary(data) whos data关键操作解析readtable是MATLAB读取表格数据的首选函数它会自动识别表头并将数据存储为table类型便于后续按列名操作。‘VariableNamingRule’, ‘preserve’这个参数至关重要。竞赛数据表的列名很可能包含空格、括号或中文如“二氧化硅(SiO2)”。如果不加此参数MATLAB会默认将列名修改为有效的变量名如去掉空格变成“二氧化硅_SiO2_”导致后续按列名索引时出错。这是我踩过的第一个坑务必加上。head和summary能快速查看数据概貌和每列的基本统计信息缺失值数量、范围等whos可以查看变量data的结构。2.2 识别与处理缺失值玻璃成分数据中缺失值非常常见可能因为检测限、样品污染或数据记录遗漏。不能简单地删除或填0。% 1. 统计每列的缺失值数量 missing_counts sum(ismissing(data), 1); disp(各列缺失值数量); disp(missing_counts); % 2. 可视化缺失模式可选但推荐 figure; ms missing_pattern(data); % 如果 missing_pattern 函数不可用可以用热图简单替代 % imagesc(ismissing(data)); % colorbar; % xlabel(变量); % ylabel(样本); % title(数据缺失模式热图); % 3. 针对成分数据的缺失值处理策略 % 策略A整行删除仅当缺失行极少且为随机缺失时考虑 if sum(missing_counts) / numel(data) 0.05 % 缺失比例小于5% data_clean rmmissing(data); % 删除任何列包含缺失值的行 else % 策略B基于业务逻辑的填充 - 这是更常用的方法 data_clean data; % 假设我们判断对于玻璃成分缺失可能意味着含量极低低于检测限 % 一种常见做法是用同类玻璃如“高钾”类该成分的中位数或最小值的一半填充 % 首先需要找到分类列假设列名为‘类型’ unique_types unique(data_clean.类型); for i 1:length(unique_types) type_mask strcmp(data_clean.类型, unique_types{i}); for j 1:width(data_clean) if ismember(data_clean.Properties.VariableNames{j}, {类型, 文物编号}) % 跳过非数值列 continue; end col_data data_clean.(j); missing_mask type_mask ismissing(col_data); if any(missing_mask) % 计算该类别下该列的非缺失值中位数 available_vals col_data(type_mask ~ismissing(col_data)); if ~isempty(available_vals) fill_value median(available_vals); % 或者用 min(available_vals)/2 表示“低于检测限” % fill_value min(available_vals(available_vals0)) / 2; data_clean.(j)(missing_mask) fill_value; end end end end end为什么这样处理直接删除rmmissing是最简单的方法但可能损失宝贵样本特别是小类别样本。仅在缺失率极低且随机时使用。对于成分数据缺失往往非随机。例如某种元素未检出很可能是因为其含量确实极低。用同类样本的中位数填充比用全局均值更合理因为它保留了类内分布特征。用最小值的一半填充则是一种更保守的“低于检测限”的估计在科学文献中常见。重要心得在论文中必须明确说明缺失值处理的方法及理由这是评审的得分点。3. 成分数据的归一化和为100%的约束玻璃化学成分数据通常是各种氧化物的重量百分比理论上所有成分之和应为100%。但实测数据由于误差总和往往在98%~102%之间波动。直接使用原始数据进行分析会引入误差特别是进行相关性分析或主成分分析时。% 1. 识别所有成分氧化物列 % 假设成分列名包含‘SiO2’, ‘Na2O’, ‘K2O’等且‘类型’、‘文物编号’、‘颜色’等为非成分列 component_columns {SiO2, ‘Na2O’ ‘K2O’ ‘CaO’ ‘MgO’ ‘Al2O3’ ‘Fe2O3’ ‘CuO’ ‘PbO’ ‘BaO’ ‘P2O5’ ‘SrO’ ‘SnO2’ ‘SO2’}; % 根据实际数据调整 % 更稳健的方法自动识别数值型列并排除指定的非成分列 all_vars data_clean.Properties.VariableNames; non_comp_vars {‘类型’ ‘文物编号’ ‘颜色’ ‘风化’}; % 根据实际表头调整 comp_var_mask ~ismember(all_vars, non_comp_vars); % 确保这些列是数值型 for v all_vars(comp_var_mask) if ~isnumeric(data_clean.(v{1})) comp_var_mask(strcmp(all_vars, v{1})) false; end end component_columns all_vars(comp_var_mask); % 2. 计算每个样本的原始成分总和 row_sums sum(data_clean{:, component_columns}, 2); disp([成分总和范围 ‘ num2str(min(row_sums)) ’ ~ ‘ num2str(max(row_sums))]) % 3. 执行归一化使每个样本的成分和为100% data_normalized data_clean; data_normalized{:, component_columns} data_clean{:, component_columns} ./ row_sums * 100; % 验证归一化结果 new_row_sums sum(data_normalized{:, component_columns}, 2); disp([归一化后总和范围 ‘ num2str(min(new_row_sums)) ’ ~ ‘ num2str(max(new_row_sums))])核心原理与坑点原理归一化消除了因检测总和不一带来的量纲影响使得不同样本之间的成分比例具有可比性。这对于后续任何基于距离或比例的分析如聚类、主成分分析都是必要的。大坑务必只对成分数据进行归一化千万不要把分类编号、标签等列也除一下这就是为什么需要精确识别component_columns。我见过有队伍不小心把“文物编号”也归一化了导致后续分析全部乱套。另一个坑检查归一化后的数据是否还有缺失值NaN。如果某一行所有成分都是缺失值那么row_sums为0除法会产生NaN。需要在归一化前确保没有这样的“全空”行。4. 描述性统计与可视化让数据自己“说话”在建模前用统计和图表直观感受数据特征能形成关键假设并指导后续的模型选择。4.1 按玻璃类型分组统计% 假设分类列名为‘类型’包含‘高钾’和‘铅钡’ group_stats grpstats(data_normalized, ‘类型’ {‘mean’ ‘std’ ‘min’ ‘max’ ‘median’} ‘DataVars’ component_columns); disp(group_stats); % 可以转置一下方便查看某个成分在不同类型间的对比 sio2_stats group_stats(:, {‘Group’ ‘mean_SiO2’ ‘std_SiO2’}); disp(sio2_stats);grpstats函数非常强大能一次性计算各组的多种统计量。结果group_stats是一个table行是分组高钾、铅钡列是各个成分的统计量。4.2 绘制成分对比箱线图箱线图能一眼看出分布的中心位置、离散程度和异常值。figure(‘Position’ [100, 100, 1200, 600]); % 设置大一点的图窗 for i 1:length(component_columns) subplot(3, 5, i); % 假设有15个成分排成3行5列 comp component_columns{i}; boxplot(data_normalized.(comp), data_normalized.类型); title(comp, ‘Interpreter’ ‘none’); % ‘none’防止下划线被当作下标 ylabel(‘含量 (%)’); grid on; end sgtitle(‘不同类型玻璃化学成分分布箱线图对比’) % 总标题从图中能看出什么如果某个成分如PbO、BaO在“铅钡玻璃”组的箱体明显高于“高钾玻璃”组且几乎没有重叠那它就是强区分因子。箱体的长度IQR反映了组内变异大小。变异小的成分可能配方更稳定。异常值箱须外的点需要留意可能是检测误差也可能是特殊的亚类。不要轻易删除它们可能包含重要信息。4.3 绘制成分均值对比柱状图% 提取均值 mean_highK group_stats{strcmp(group_stats.Group, ‘高钾’) startsWith(group_stats.Properties.VariableNames, ‘mean_’)} mean_PbBa group_stats{strcmp(group_stats.Group, ‘铅钡’) startsWith(group_stats.Properties.VariableNames, ‘mean_’)} mean_highK mean_highK(:)’; % 转为行向量 mean_PbBa mean_PbBa(:)’; % 绘图 figure; x 1:length(component_columns); bar(x, [mean_highK; mean_PbBa]’); set(gca, ‘XTick’ x, ‘XTickLabel’ component_columns, ‘XTickLabelRotation’ 45); legend(‘高钾玻璃’ ‘铅钡玻璃’) ylabel(‘平均含量 (%)’) title(‘不同类型玻璃化学成分平均含量对比’) grid on;这张图可以非常直观地展示两类玻璃在配方上的宏观差异比如铅钡玻璃中PbO和BaO的突出地位。5. 统计显著性检验差异是“真的”吗描述性统计显示了差异但我们需要用统计检验来确认这种差异不是由随机抽样误差造成的。对于两类玻璃高钾 vs 铅钡的每个成分比较两独立样本t检验是合适的选择。这里就涉及到热词中的一个具体问题ttest和ttest2的区别。这是MATLAB初学者常混淆的点。ttest用于单样本t检验检验一组数据的均值是否与某个假设值如0有显著差异。例如检验一批玻璃的含铅量是否显著大于0。ttest2用于两独立样本t检验检验两组独立数据的均值是否有显著差异。这正是我们当前场景需要的。% 准备两组数据 idx_highK strcmp(data_normalized.类型 ‘高钾’) idx_PbBa strcmp(data_normalized.类型 ‘铅钡’) % 初始化结果存储 p_values zeros(1, length(component_columns)); h_values zeros(1, length(component_columns)); % h1 表示拒绝原假设有显著差异 ci_cell cell(1, length(component_columns)); % 置信区间 stats_cell cell(1, length(component_columns)); % 检验统计量信息 for i 1:length(component_columns) comp component_columns{i}; data1 data_normalized.(comp)(idx_highK); data2 data_normalized.(comp)(idx_PbBa); % 进行两样本t检验默认假设两组方差不等更保守的‘Welch’s t-test’ [h, p, ci, stats] ttest2(data1, data2, ‘Vartype’ ‘unequal’) h_values(i) h; p_values(i) p; ci_cell{i} ci; stats_cell{i} stats; end % 将结果整理成表格方便查看 result_table table(component_columns’ h_values’ p_values’ ‘VariableNames’ {‘成分’ ‘显著差异’ ‘p值’}) disp(‘两独立样本t检验结果原假设两组均值无差异’) disp(result_table); % 可以标记出p值小于0.05或0.01的显著成分 sig_idx p_values 0.05; sig_components component_columns(sig_idx); disp([‘在显著性水平0.05下有以下成分存在显著差异 ‘ strjoin(sig_components, ‘ ‘)])解读与注意事项原假设两类玻璃在该成分上的均值相等。p值如果p值很小通常0.05我们就有足够证据拒绝原假设认为差异是统计显著的。‘Vartype’ ‘unequal’参数我们通常不知道两组的方差是否相等。选择‘unequal’表示使用不假设等方差的t检验Welch校正这比默认的等方差检验更稳健尤其是在样本量不等或方差明显不同时。这是实际分析中的最佳实践。结果应用找出那些p值极小的成分如PbO、BaO、K2O它们就是后续构建分类模型时最重要的特征。可以在论文中用“*”和“**”在表格中标出不同显著性水平的结果显得非常专业。6. 相关性分析与热图探索成分间的“共生关系”玻璃配方中某些元素可能一起添加或此消彼长。相关性分析可以帮助我们理解这些关系。% 计算所有成分间的相关系数矩阵 corr_matrix corrcoef(data_normalized{:, component_columns}) % 绘制相关系数热图 figure(‘Position’ [100, 100, 800, 700]) imagesc(corr_matrix) colorbar; caxis([-1, 1]); % 固定颜色轴范围 colormap(jet); % 可以使用其他配色如 parula, hot % 添加坐标轴标签 set(gca, ‘XTick’ 1:length(component_columns) ‘XTickLabel’ component_columns, ‘XTickLabelRotation’ 90) set(gca, ‘YTick’ 1:length(component_columns) ‘YTickLabel’ component_columns) title(‘古代玻璃化学成分相关系数矩阵热图’) % 为了更清晰可以只显示绝对值较大的相关性或添加数值标签 % 找出强相关例如 |r| 0.7的配对 [comp1_idx, comp2_idx] find(abs(corr_matrix) 0.7 triu(ones(size(corr_matrix)), 1)) % triu避免重复和自身相关 for k 1:length(comp1_idx) fprintf(‘%s 与 %s 的相关系数为 %.3f\n’ component_columns{comp1_idx(k)} component_columns{comp2_idx(k)} corr_matrix(comp1_idx(k), comp2_idx(k))) end从相关性中能发现什么强正相关例如Na2O和K2O若强正相关可能暗示它们来自同一种矿物原料如草木灰。强负相关例如SiO2和助熔剂Na2O, K2O之间可能存在负相关因为增加助熔剂会降低硅含量比例。铅钡玻璃中的特殊关系PbO和BaO的相关性值得特别关注它能反映这两种助熔剂的使用是固定的配方还是可变的。注意伪相关两个成分都与第三个成分如总和约束相关时它们之间也可能显示出相关性。这就是为什么先做归一化很重要它部分消除了这种“常数和”带来的伪相关。7. 第一问代码整合与进阶思考将以上步骤整合成一个脚本或函数就是第一问完整的分析代码。但作为竞赛不能只停留在跑通代码更重要的是在论文中体现你的思考。在论文中如何呈现第一问的结果数据预处理部分用一小段文字说明处理了缺失值、进行了归一化并简要说明理由。可以附上一张显示归一化前后总和分布的对比图。描述性统计制作一个清晰的表格列出两类玻璃各成分的平均值、标准差、中位数等。用箱线图或分组柱状图作为可视化支持。统计检验用另一个表格展示t检验的p值并用星号标注显著性水平。在文中指出哪些成分是显著差异的这为“高钾”和“铅钡”的分类提供了化学依据。相关性分析展示相关系数热图并挑选1-2对最具代表性的强相关或强负相关成分进行解读将其与古代玻璃制作工艺的潜在知识联系起来。进阶思考与可能的扩展风化影响第一问数据可能包含“风化”与“无风化”样本。你可以额外分析风化对成分的影响例如风化是否导致某些碱性氧化物流失这能体现你思维的深度。同样使用t检验或方差分析。子类探索即使在“高钾”或“铅钡”大类内部箱线图可能显示存在多个离群点或分布多峰。这提示可能存在亚类。可以用聚类方法如K-means对每一大类内部进行探索性分析但这通常属于第二问或第三问的范畴。特征筛选基于t检验的p值或均值差异大小可以对成分特征进行排序为后续的分类模型如第二问的SVM、随机森林提供特征选择依据。写代码只是解决了“怎么做”的问题而结合考古学背景解释数据背后的故事才是数学建模竞赛获奖的关键。你的代码和分析最终都是为了支撑一个逻辑严谨、发现新颖的科学叙事。