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_vals>0)) / 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)); % h=1 表示拒绝原假设（有显著差异） 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、随机森林）提供特征选择依据。
写代码只是解决了“怎么做”的问题，而结合考古学背景解释数据背后的故事，才是数学建模竞赛获奖的关键。你的代码和分析，最终都是为了支撑一个逻辑严谨、发现新颖的科学叙事。
