简介这份2022全国大学生数学建模C题资料面向备赛国赛、需要复盘古代玻璃文物成分分析与鉴别案例的高校学生及建模爱好者核心是用数学建模与数据分析方法解决玻璃风化效应、成分预测和类别鉴别问题。压缩包共1个文件为3.55MB的PDF文档集中呈现赛题完整求解方案可对照学习问题重述、数据探索、预处理与四问建模思路。内容覆盖主成分分析降维、因子分析转换未风化成分、敏感性分析、灰色关联度分析及分类预测等关键环节并结合高钾玻璃与铅钡玻璃的风化规律、亚类划分和关联度比较给出步骤与结果便于读者梳理建模框架、理解方法选择并借鉴论文表达。已有928人学习适合希望系统掌握赛题分析、模型构建与结果解释的参赛者查阅。1. 风化层下的成分账本这道题到底在考什么2022 全国大学生数学建模 C 题给的是一张“玻璃文物身份证”和一张“化学成分账本”前者记录纹饰、类型、颜色和表面是否风化后者记录十几个氧化物的质量百分比附件末尾还附了一批类型未知的样品等你判断归属。它表面上是一道统计题实质是一次完整的成分数据处理训练——缺失怎么定义、百分比怎么变换、小样本怎么验证、判别阈值怎么定四件事一件都躲不掉。做完这套流程你拿到的不只是四个问题的答案而是一套能迁移到材料检测、文物断源、工业配方归类的分析方法。适合正在准备数学建模、或者手头就有一堆成分百分比表格要分析的人。真正让人犯难的在于“风化”这两个字把数据分成了两种性质不同的样本。风化层的成分和本体不是一回事铅钡玻璃风化后氧化铅会明显流失高钾玻璃的钾也会跟着环境走。如果一上来就把风化样品和未风化样品混在一起算均值后面所有结论都是虚的。所以这道题的第一步是把“数据可信边界”画清楚再谈模型。2. 把成分表洗成可建模矩阵缺失值、合计量与对数比拿到成分表后不要急着画图先把表结构和缺失机制搞清楚。这一步做扎实后面四问的结论才有落脚点。我一般按“读表 → 对齐编号 → 判断缺失类型 → 变换”四步走每一步都留一个可复查的中间结果。2.1 先分清“未检测”和“含量为零”坑就在这一列。赛题附件里部分氧化物在某些样品上是空值这代表当时没有测这一项不代表含量真的是零。如果你用fillna(0)一把梭等于把大量“没测”当成“没有”算均值时分子不变、分母变大结果系统性偏低。import pandas as pd import numpy as np # 三张表文物基本信息、化学成分、未知类别样品 info pd.read_excel(附件1.xlsx) chem pd.read_excel(附件2.xlsx) unknown pd.read_excel(附件3.xlsx) # 列名去空格避免合并时对不上 for df in (info, chem, unknown): df.columns [str(c).strip().replace( , ) for c in df.columns] # 采样点编码里带部位信息抽出前导数字得到文物编号 chem[文物编号] chem[文物采样点].astype(str).str.extract(r^(\d))[0] base chem.merge(info, on文物编号, howleft, validatemany_to_one) # 识别氧化物列以氧化/二氧化/三氧化/五氧化开头 oxide_cols [c for c in chem.columns if c.startswith((氧化, 二氧化, 三氧化, 五氧化))] # 强制转数值非数字一律变 NaN保留“未检测”语义 base[oxide_cols] base[oxide_cols].apply(pd.to_numeric, errorscoerce) print(base[oxide_cols].isna().sum().sort_values(ascendingFalse).head(8))这段的关键动作有三个str.extract从采样点里抽文物编号解决“一张表按样品、一张表按采样点”的粒度错位validatemany_to_one让合并时如果出现一对多会直接报错而不是悄悄产生笛卡尔积errorscoerce把“0.1”这类文本或异常符号统一变成 NaN防止后面求均值时崩掉。参数上唯一需要斟酌的是howleft。用左连接保留化学成分表全部行是因为基本信息表可能缺某些样品但成分数据不能丢。如果反过来用右连接你会莫名其妙少掉几条记录而且不一定马上发现。接下来看合计量。每种玻璃检测的氧化物项目不完全一样把所有列直接加起来通常会偏离 100偏离多少反映的是检测覆盖度# min_count1 保证整行全空时结果为 NaN而不是 0 base[检测合计] base[oxide_cols].sum(axis1, min_count1) print(base[检测合计].describe()) print((base[检测合计] 95).sum(), 条记录检测项不完整)如果相当一部分样品合计明显低于 100说明这些样品有好几项没测。这时候做跨样品比较就必须限定在同一组“测了哪些项”的样品之间或者先把缺失项多的样品单独标出来。2.2 成分数据的闭合效应与 CLR 变换成分数据有一个绕不开的性质所有分量加起来是常数100% 或某个检测总量这叫闭合效应。一个分量升高其他分量被动降低于是直接对百分比做相关分析会得到一堆伪相关——氧化铅和氧化钡看起来强相关其实只是“此消彼长”的算术结果。常见的处理是做中心对数比变换CLR把单纯形空间里的数据映射到欧氏空间再做聚类或回归def clr_transform(df, cols, eps1e-6): 对成分列做 CLR先加极小值避免 log(0)再减去行几何均值 x df[cols].fillna(0) eps logx np.log(x) # 每行减去该行所有分量对数的均值得到零和约束下的坐标 return logx.sub(logx.mean(axis1), axis0) clr_df clr_transform(base, oxide_cols)eps的取值是这里唯一的自由参数。取 1e-6 是因为原始百分比量级在 0.01 到 90 之间太大会扭曲低含量成分的相对关系太小则在浮点精度上容易出问题。需要注意的是CLR 里fillna(0)只是为了算几何均值、让变换可执行真正做统计检验时原始表里的 NaN 不能替换两套数据要分开存放。另一个常见做法是只用“主要成分子集”做 CLR比如只保留 SiO₂、K₂O、PbO、BaO、Al₂O₃ 这几个区分度高的分量。分量选得越多微量元素的噪声占比越大小样本下反而不稳。提示CLR 之后再接标准化要谨慎。CLR 结果本身已经是零和约束再套一层 StandardScaler 等于把约束信息又打散一次两种处理二选一即可别叠着用。到这里你手上应该有三份东西一份保留原始缺失的明细表、一份用于聚类的 CLR 矩阵、一份用于描述统计的按类型分组的汇总表。后面所有分析都从这三份出发不要临时再从 Excel 里重读一遍否则口径迟早对不上。3. 风化与类型的关联检验卡方、效应量和配对差异第一问要回答“表面风化跟类型、纹饰、颜色有没有关系”第二问要回答“风化前后成分怎么变”。两个问题的性质完全不同前者是分类变量之间的关联后者是连续变量的组间差异。用错检验方法是这一章最常见的翻车点。3.1 用卡方和 Cramérs V 判断类型、纹饰、颜色与风化的关系类型、纹饰、颜色、风化状态都是分类变量常规做法是列联表加卡方独立性检验。但卡方只给 p 值p 值小只说明“不太可能是巧合”不说明关联有多强。样本量几十条的时候p 值对样本量极其敏感必须补一个效应量否则你没法判断这个关联是“强到能用作判别依据”还是“统计显著但实际没用”。from scipy.stats import chi2_contingency import numpy as np def cramers_v(x, y): tab pd.crosstab(x, y) chi2, p, dof, expected chi2_contingency(tab) n tab.values.sum() # 自由度取行列数较小者减一 v np.sqrt(chi2 / (n * (min(tab.shape) - 1))) # 小样本偏差修正 v_corrected max(0, v - (tab.shape[0]-1)*(tab.shape[1]-1)/(2*n-2)) small (expected 5).sum() return {chi2: chi2, p: p, V: v, V_adj: v_corrected, 期望频数5的格子数: small} for col in [类型, 纹饰, 颜色]: print(col, cramers_v(base[表面风化], base[col]))代码里有两个细节值得说。一是expected 5的格子数卡方检验的前提是期望频数不能太小如果超过两成的格子期望频数低于 5p 值就不可信了。这时候要么合并类别比如颜色里的深浅相近项合并要么改用 Fisher 精确检验的扩展版本。二是V_adj的修正公式样本量小于一百时未修正的 Cramérs V 会系统性高估关联强度修正后才好做横向比较。解读时的经验门槛V 在 0.1 以下基本可以认为无关0.3 左右算中等0.5 以上才谈得上强关联。如果某个纹饰类别只有两三个样品哪怕 V 很大也不要写进结论那是单个样本撑起来的假象。另一个容易忽略的点是风化状态本身可能是“类型”的下游变量。铅钡玻璃和高钾玻璃的耐候性不同所以“类型与风化相关”这件事本质上是材料性质差异的体现不是风化导致了类型。写结论时要区分相关和因果别把顺序说反。3.2 风化前后成分怎么比非参数检验加效应量第二问要求分析风化对成分的影响并预测风化前含量。这里第一个要澄清的概念是同一件文物的风化部位和未风化部位不是配对样本因为采样点不同、部位不同甚至可能来自不同器物。所以用配对 t 检验是错的应该按“类型 × 风化状态”分组做独立样本比较。成分百分比分布通常右偏、样本量小、还带缺失正态性靠不住我一般直接用 Mann-Whitney U 检验并配上秩双列相关系数作为效应量from scipy.stats import mannwhitneyu rows [] for t in base[类型].dropna().unique(): sub base[base[类型] t] for c in oxide_cols: a sub.loc[sub[表面风化] 风化, c].dropna() b sub.loc[sub[表面风化] 未风化, c].dropna() # 每组至少 3 个有效值才做检验否则结果没有意义 if len(a) 3 and len(b) 3: u, p mannwhitneyu(a, b, alternativetwo-sided) # 秩双列相关范围 -1 到 1绝对值越大差异越明显 r 1 - 2*u/(len(a)*len(b)) rows.append({类型: t, 成分: c, n_风化: len(a), n_未风化: len(b), 均值_风化: a.mean(), 均值_未风化: b.mean(), 差值: a.mean() - b.mean(), p: p, 效应量r: r}) res pd.DataFrame(rows).sort_values([类型, p]) print(res.head(15))参数上alternativetwo-sided是默认选择因为风化既可能让某个成分富集也可能让它流失没有先验方向。len(a) 3这个门槛是我自己加的两组各只有一两个值时U 检验的 p 值只能是 0.1 以上的几个离散值报出来纯属自欺欺人。这类成分在论文里应该写成“样本不足不做统计推断”只给描述性均值。真正决定结论质量的是效应量那一列。p 值受样本量影响效应量反映差异的实际幅度。氧化铅在高钾玻璃里可能只有 1% 到 2% 的绝对差p 值却很显著这种差异对判别毫无帮助而在铅钡玻璃里风化前后氧化铅差十几个百分点才是真正的规律。关于“预测风化前含量”稳妥路径不是去拟合一个从风化到未风化的回归而是用同类未风化样品的稳健统计量中位数或截尾均值作为该类型的参考值再对风化样品给出区间估计。原因很直接风化样品数量有限回归模型的外推能力没有保障而中位数在小样本、偏态分布下更稳。如果一定要建模用分位数回归比普通最小二乘更合适因为它对离群值不敏感。到这里第一、二问的骨架就有了一张关联检验表带效应量、一张分组差异表带效应量和样本量、一组按类型给出的风化前参考值。表格里的数字必须带样本量这是小样本分析的底线。4. 从聚类到判别两类玻璃的分类边界与未知样品归属第三问是全题的核心找出高钾玻璃和铅钡玻璃的分类规律再把未知样品判进去。很多人一上来就丢个 KMeans 或者随机森林跑出准确率 0.95 就收工这在答辩环节基本会被问穿——样本量只有几十条交叉验证的方差比模型间的差异还大。4.1 先用 KMeans 看自然分群再决定是否相信标签聚类的作用不是给答案而是检查“已知的类型标签”和“数据自身的结构”是否一致。如果 KMeans 在 k2 时的分群和类型标签高度吻合说明成分确实能分开这两类后面的判别模型有物理基础如果吻合度很差那要么是特征选错了要么是风化样品把结构搅乱了。from sklearn.preprocessing import StandardScaler from sklearn.cluster import KMeans from sklearn.metrics import silhouette_score, adjusted_rand_score # 只用未风化样品探索结构避免风化层干扰 clean base[base[表面风化] 未风化].dropna(subset[类型]).copy() X clr_transform(clean, oxide_cols).values for k in range(2, 5): km KMeans(n_clustersk, n_init20, random_state42).fit(X) sil silhouette_score(X, km.labels_) ari adjusted_rand_score(clean[类型], km.labels_) if k 2 else float(nan) print(fk{k} 轮廓系数{sil:.3f} 与标签的ARI{ari:.3f})n_init20是因为 KMeans 对初始质心敏感小样本下多跑几次取最优能显著降低结果随机性。random_state42保证结果可复现这在提交论文时很重要。adjusted_rand_score用来衡量聚类结果和真实标签的一致程度取值在 0 附近表示跟随机分组差不多接近 1 表示几乎完全一致。如果 k2 的轮廓系数只有 0.2 出头说明结构不明显这时候不要硬套聚类结论转而检查是不是该先做特征筛选——把氧化铅、氧化钡、氧化钾这几个判别力最强的分量单独拿出来看二维散点往往比十几维空间里的轮廓系数直观得多。这里有个反直觉的点风化样品在聚类里经常自成一簇。因为它所有成分都被稀释或富集了整体偏离本体。如果你的聚类把风化状态当成了主要分群依据那这个模型对“类型判别”毫无价值。所以探索阶段必须把风化样品排除或者把风化状态作为分层变量。4.2 决策树和 SVM 的判别规则以及未知样品的预测判别模型要解决两个需求一是准确二是可解释。赛题要求“给出分类规律”所以模型必须能翻译成人话比如“氧化铅大于某阈值判为铅钡”。决策树天然满足这一点SVM 作为对照验证非线性边界是否存在。from sklearn.tree import DecisionTreeClassifier, export_text from sklearn.svm import SVC from sklearn.pipeline import make_pipeline from sklearn.model_selection import LeaveOneOut, cross_val_score X clr_transform(base, oxide_cols).values y base[类型].values tree make_pipeline( DecisionTreeClassifier(max_depth3, min_samples_leaf3, class_weightbalanced, random_state42) ) svm make_pipeline(StandardScaler(), SVC(kernelrbf, C1.0, gammascale, class_weightbalanced)) loo LeaveOneOut() print(决策树 LOO 准确率:, cross_val_score(tree, X, y, cvloo, scoringaccuracy).mean()) print(SVM LOO 准确率:, cross_val_score(svm, X, y, cvloo, scoringaccuracy).mean()) tree.fit(X, y) print(export_text(tree[-1], feature_namesoxide_cols, decimals2))max_depth3和min_samples_leaf3是刻意压低的。样本几十条的情况下不限制深度的树会把每个样品都记住训练集准确率 100%换个数据集就崩。class_weightbalanced应对两类样品数量不平衡否则模型会偏向多数类。交叉验证用留一法而不是 5 折是因为总量太小5 折时每折验证集只有十几条方差大到没法比较模型。export_text打印出来的规则就是第三问的核心答案。你会看到类似“氧化铅 ≤ 某值 且 氧化钾 某值 → 高钾”的路径。把这些路径整理成阈值表比报一个准确率有说服力得多。对未知样品的预测要在同一套变换下做千万别对未知样品单独做 CLRXu clr_transform(unknown, oxide_cols).values pred tree.predict(Xu) prob tree.predict_proba(Xu) unknown_out unknown[[文物采样点]].copy() unknown_out[预测类型] pred unknown_out[置信度] prob.max(axis1) print(unknown_out.sort_values(置信度).head(10))注意prob.max(axis1)低的那几条它们大概率落在两类交界处是论文里必须单独讨论的对象。把它们简单归到某一类而不加说明是答辩时最容易被追问的地方。合理做法是给出“倾向铅钡但置信度 0.6建议补充检测氧化钡”这样的表述。5. 五处最容易翻车的细节填零、泄漏与阈值拍脑袋这一章是本篇的血泪经验。下面五条几乎每一届都有人踩而且踩了之后结果表面上还挺好看属于典型的黑匣子。5.1 数据层的三个坑第一条缺失值填零。现象某类玻璃的氧化铅均值只有 3%比文献里的典型值低一大截。 原因大量“未检测”被填成了 0均值被严重稀释。 解决缺失保持 NaN统计时用dropna()分组计算只有在必须构造完整矩阵比如聚类时才用极小值或中位数替代并且替代后的副本单独命名别覆盖原始表。第二条把多次采样当成多个独立样品。现象某件文物因为采了三个部位在统计里权重变成别人的三倍卡方检验结果全歪。 原因附件按采样点记录直接value_counts会让采样多的文物占主导。 解决先按文物编号聚合成分取均值类别取唯一值或者在做关联检验时用文物级数据在做成分分析时用采样点级数据两套口径分开维护。第三条百分比直接做相关分析。现象氧化铅和氧化钡的相关系数飙到 0.9 以上写成“两者共生”。 原因闭合效应一个涨另一个必然被动跌或者反过来。 解决先做 CLR 或对数比变换再算相关如果坚持用原始百分比就改用 Spearman 秩相关并明确说明局限。5.2 建模层的两个坑第四条训练集和验证集不分离。现象模型准确率 0.98换一批数据就掉到 0.7。 原因在全体样本上做了特征筛选或标准化然后再交叉验证信息从验证折泄漏进了训练过程。 解决把特征选择、变换、标准化全部塞进 Pipeline交叉验证只作用于 Pipeline 整体。或者更稳妥用嵌套交叉验证。第五条判别阈值拍脑袋。现象氧化铅阈值取 5%未知样品里有两个被判成高钾但它们的氧化铅是 4.8%。 原因阈值没有经过敏感性检验而且风化样品的氧化铅可能已被淋失同一类型风化前后的值差异很大。 解决阈值按风化状态分层设定或者改用未风化样品建立阈值同时做阈值扫描看结论在多大范围内稳定下一章展开。这五条的共同特征是出错时不会报错只会让你得到一个“看起来很合理”的错误结论。所以每一处都要留一个中间检查点把处理前后的样本量、均值、分布画出来对比别等最后写论文时才发现数字对不上。6. 敏感性分析把“判为铅钡”变成一句有底气的话第四问要求比较两类玻璃的成分关联性并对新样品提出鉴别方法真正拉开差距的是敏感性分析。评委想知道的不是“你的方法准确率多少”而是“换个假设你的结论还站得住吗”。我一般做三组扫描。第一组是阈值扫描把氧化铅的判别阈值从 1% 扫到 10%步长 0.5%记录每个阈值下未知样品的归属变化。第二组是重采样对训练集做 200 次 bootstrap每次重训决策树统计“氧化铅阈值”这条规则的取值分布。第三组是缺失处理方案对比分别用丢弃、中位数填补、极小值填补三种策略跑完整流程看最终归类是否一致。判别阈值PbO判为铅钡的未知样品数与基准方案一致数结论稳定性2.0%97偏低易把高钾误判5.0%66基准8.0%44偏保守可能漏判10.0%33明显漏判表格里的规律是阈值在 4% 到 6% 之间时结论几乎不变超过 8% 开始漏判低于 2% 开始误判。把这段区间写进结论比单给一个数字强得多——“建议阈值取 5%在 4% 至 6% 范围内判别结论稳定”。bootstrap 那组还有一个额外好处你能看到哪些样品在重采样中反复变类。这些样品就是判别边界上的模糊样本应该在论文里单独列出说明“其成分同时具备两类特征建议结合纹饰与颜色辅助判断”。至于两类玻璃的成分关联性比较我的做法是分别对高钾组和铅钡组算 CLR 后的相关矩阵再做矩阵差异检验而不是把两类混在一起算一个总相关。混算的结果会被类型这个混杂变量主导看起来所有成分都相关实际上只是两类玻璃的均值不同。最后说一句自己的习惯这道题我做完第一版之后会把所有中间表导成 CSV 存到同一个目录命名带上日期和处理策略比如chem_未风化_中位数填补.csv。因为做到第四问时你一定会想回头改前面的某个假设没有这些中间产物就得从头再跑一遍那种感觉比返工还难受。建模这件事后悔药就是中间文件。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/31 7:13:23 · 阅读 1.2k 玻璃文物成分数据建模：缺失值处理、CLR变换与分类判别
简介这份2022全国大学生数学建模C题资料面向备赛国赛、需要复盘古代玻璃文物成分分析与鉴别案例的高校学生及建模爱好者核心是用数学建模与数据分析方法解决玻璃风化效应、成分预测和类别鉴别问题。压缩包共1个文件为3.55MB的PDF文档集中呈现赛题完整求解方案可对照学习问题重述、数据探索、预处理与四问建模思路。内容覆盖主成分分析降维、因子分析转换未风化成分、敏感性分析、灰色关联度分析及分类预测等关键环节并结合高钾玻璃与铅钡玻璃的风化规律、亚类划分和关联度比较给出步骤与结果便于读者梳理建模框架、理解方法选择并借鉴论文表达。已有928人学习适合希望系统掌握赛题分析、模型构建与结果解释的参赛者查阅。1. 风化层下的成分账本这道题到底在考什么2022 全国大学生数学建模 C 题给的是一张“玻璃文物身份证”和一张“化学成分账本”前者记录纹饰、类型、颜色和表面是否风化后者记录十几个氧化物的质量百分比附件末尾还附了一批类型未知的样品等你判断归属。它表面上是一道统计题实质是一次完整的成分数据处理训练——缺失怎么定义、百分比怎么变换、小样本怎么验证、判别阈值怎么定四件事一件都躲不掉。做完这套流程你拿到的不只是四个问题的答案而是一套能迁移到材料检测、文物断源、工业配方归类的分析方法。适合正在准备数学建模、或者手头就有一堆成分百分比表格要分析的人。真正让人犯难的在于“风化”这两个字把数据分成了两种性质不同的样本。风化层的成分和本体不是一回事铅钡玻璃风化后氧化铅会明显流失高钾玻璃的钾也会跟着环境走。如果一上来就把风化样品和未风化样品混在一起算均值后面所有结论都是虚的。所以这道题的第一步是把“数据可信边界”画清楚再谈模型。2. 把成分表洗成可建模矩阵缺失值、合计量与对数比拿到成分表后不要急着画图先把表结构和缺失机制搞清楚。这一步做扎实后面四问的结论才有落脚点。我一般按“读表 → 对齐编号 → 判断缺失类型 → 变换”四步走每一步都留一个可复查的中间结果。2.1 先分清“未检测”和“含量为零”坑就在这一列。赛题附件里部分氧化物在某些样品上是空值这代表当时没有测这一项不代表含量真的是零。如果你用fillna(0)一把梭等于把大量“没测”当成“没有”算均值时分子不变、分母变大结果系统性偏低。import pandas as pd import numpy as np # 三张表文物基本信息、化学成分、未知类别样品 info pd.read_excel(附件1.xlsx) chem pd.read_excel(附件2.xlsx) unknown pd.read_excel(附件3.xlsx) # 列名去空格避免合并时对不上 for df in (info, chem, unknown): df.columns [str(c).strip().replace( , ) for c in df.columns] # 采样点编码里带部位信息抽出前导数字得到文物编号 chem[文物编号] chem[文物采样点].astype(str).str.extract(r^(\d))[0] base chem.merge(info, on文物编号, howleft, validatemany_to_one) # 识别氧化物列以氧化/二氧化/三氧化/五氧化开头 oxide_cols [c for c in chem.columns if c.startswith((氧化, 二氧化, 三氧化, 五氧化))] # 强制转数值非数字一律变 NaN保留“未检测”语义 base[oxide_cols] base[oxide_cols].apply(pd.to_numeric, errorscoerce) print(base[oxide_cols].isna().sum().sort_values(ascendingFalse).head(8))这段的关键动作有三个str.extract从采样点里抽文物编号解决“一张表按样品、一张表按采样点”的粒度错位validatemany_to_one让合并时如果出现一对多会直接报错而不是悄悄产生笛卡尔积errorscoerce把“0.1”这类文本或异常符号统一变成 NaN防止后面求均值时崩掉。参数上唯一需要斟酌的是howleft。用左连接保留化学成分表全部行是因为基本信息表可能缺某些样品但成分数据不能丢。如果反过来用右连接你会莫名其妙少掉几条记录而且不一定马上发现。接下来看合计量。每种玻璃检测的氧化物项目不完全一样把所有列直接加起来通常会偏离 100偏离多少反映的是检测覆盖度# min_count1 保证整行全空时结果为 NaN而不是 0 base[检测合计] base[oxide_cols].sum(axis1, min_count1) print(base[检测合计].describe()) print((base[检测合计] 95).sum(), 条记录检测项不完整)如果相当一部分样品合计明显低于 100说明这些样品有好几项没测。这时候做跨样品比较就必须限定在同一组“测了哪些项”的样品之间或者先把缺失项多的样品单独标出来。2.2 成分数据的闭合效应与 CLR 变换成分数据有一个绕不开的性质所有分量加起来是常数100% 或某个检测总量这叫闭合效应。一个分量升高其他分量被动降低于是直接对百分比做相关分析会得到一堆伪相关——氧化铅和氧化钡看起来强相关其实只是“此消彼长”的算术结果。常见的处理是做中心对数比变换CLR把单纯形空间里的数据映射到欧氏空间再做聚类或回归def clr_transform(df, cols, eps1e-6): 对成分列做 CLR先加极小值避免 log(0)再减去行几何均值 x df[cols].fillna(0) eps logx np.log(x) # 每行减去该行所有分量对数的均值得到零和约束下的坐标 return logx.sub(logx.mean(axis1), axis0) clr_df clr_transform(base, oxide_cols)eps的取值是这里唯一的自由参数。取 1e-6 是因为原始百分比量级在 0.01 到 90 之间太大会扭曲低含量成分的相对关系太小则在浮点精度上容易出问题。需要注意的是CLR 里fillna(0)只是为了算几何均值、让变换可执行真正做统计检验时原始表里的 NaN 不能替换两套数据要分开存放。另一个常见做法是只用“主要成分子集”做 CLR比如只保留 SiO₂、K₂O、PbO、BaO、Al₂O₃ 这几个区分度高的分量。分量选得越多微量元素的噪声占比越大小样本下反而不稳。提示CLR 之后再接标准化要谨慎。CLR 结果本身已经是零和约束再套一层 StandardScaler 等于把约束信息又打散一次两种处理二选一即可别叠着用。到这里你手上应该有三份东西一份保留原始缺失的明细表、一份用于聚类的 CLR 矩阵、一份用于描述统计的按类型分组的汇总表。后面所有分析都从这三份出发不要临时再从 Excel 里重读一遍否则口径迟早对不上。3. 风化与类型的关联检验卡方、效应量和配对差异第一问要回答“表面风化跟类型、纹饰、颜色有没有关系”第二问要回答“风化前后成分怎么变”。两个问题的性质完全不同前者是分类变量之间的关联后者是连续变量的组间差异。用错检验方法是这一章最常见的翻车点。3.1 用卡方和 Cramérs V 判断类型、纹饰、颜色与风化的关系类型、纹饰、颜色、风化状态都是分类变量常规做法是列联表加卡方独立性检验。但卡方只给 p 值p 值小只说明“不太可能是巧合”不说明关联有多强。样本量几十条的时候p 值对样本量极其敏感必须补一个效应量否则你没法判断这个关联是“强到能用作判别依据”还是“统计显著但实际没用”。from scipy.stats import chi2_contingency import numpy as np def cramers_v(x, y): tab pd.crosstab(x, y) chi2, p, dof, expected chi2_contingency(tab) n tab.values.sum() # 自由度取行列数较小者减一 v np.sqrt(chi2 / (n * (min(tab.shape) - 1))) # 小样本偏差修正 v_corrected max(0, v - (tab.shape[0]-1)*(tab.shape[1]-1)/(2*n-2)) small (expected 5).sum() return {chi2: chi2, p: p, V: v, V_adj: v_corrected, 期望频数5的格子数: small} for col in [类型, 纹饰, 颜色]: print(col, cramers_v(base[表面风化], base[col]))代码里有两个细节值得说。一是expected 5的格子数卡方检验的前提是期望频数不能太小如果超过两成的格子期望频数低于 5p 值就不可信了。这时候要么合并类别比如颜色里的深浅相近项合并要么改用 Fisher 精确检验的扩展版本。二是V_adj的修正公式样本量小于一百时未修正的 Cramérs V 会系统性高估关联强度修正后才好做横向比较。解读时的经验门槛V 在 0.1 以下基本可以认为无关0.3 左右算中等0.5 以上才谈得上强关联。如果某个纹饰类别只有两三个样品哪怕 V 很大也不要写进结论那是单个样本撑起来的假象。另一个容易忽略的点是风化状态本身可能是“类型”的下游变量。铅钡玻璃和高钾玻璃的耐候性不同所以“类型与风化相关”这件事本质上是材料性质差异的体现不是风化导致了类型。写结论时要区分相关和因果别把顺序说反。3.2 风化前后成分怎么比非参数检验加效应量第二问要求分析风化对成分的影响并预测风化前含量。这里第一个要澄清的概念是同一件文物的风化部位和未风化部位不是配对样本因为采样点不同、部位不同甚至可能来自不同器物。所以用配对 t 检验是错的应该按“类型 × 风化状态”分组做独立样本比较。成分百分比分布通常右偏、样本量小、还带缺失正态性靠不住我一般直接用 Mann-Whitney U 检验并配上秩双列相关系数作为效应量from scipy.stats import mannwhitneyu rows [] for t in base[类型].dropna().unique(): sub base[base[类型] t] for c in oxide_cols: a sub.loc[sub[表面风化] 风化, c].dropna() b sub.loc[sub[表面风化] 未风化, c].dropna() # 每组至少 3 个有效值才做检验否则结果没有意义 if len(a) 3 and len(b) 3: u, p mannwhitneyu(a, b, alternativetwo-sided) # 秩双列相关范围 -1 到 1绝对值越大差异越明显 r 1 - 2*u/(len(a)*len(b)) rows.append({类型: t, 成分: c, n_风化: len(a), n_未风化: len(b), 均值_风化: a.mean(), 均值_未风化: b.mean(), 差值: a.mean() - b.mean(), p: p, 效应量r: r}) res pd.DataFrame(rows).sort_values([类型, p]) print(res.head(15))参数上alternativetwo-sided是默认选择因为风化既可能让某个成分富集也可能让它流失没有先验方向。len(a) 3这个门槛是我自己加的两组各只有一两个值时U 检验的 p 值只能是 0.1 以上的几个离散值报出来纯属自欺欺人。这类成分在论文里应该写成“样本不足不做统计推断”只给描述性均值。真正决定结论质量的是效应量那一列。p 值受样本量影响效应量反映差异的实际幅度。氧化铅在高钾玻璃里可能只有 1% 到 2% 的绝对差p 值却很显著这种差异对判别毫无帮助而在铅钡玻璃里风化前后氧化铅差十几个百分点才是真正的规律。关于“预测风化前含量”稳妥路径不是去拟合一个从风化到未风化的回归而是用同类未风化样品的稳健统计量中位数或截尾均值作为该类型的参考值再对风化样品给出区间估计。原因很直接风化样品数量有限回归模型的外推能力没有保障而中位数在小样本、偏态分布下更稳。如果一定要建模用分位数回归比普通最小二乘更合适因为它对离群值不敏感。到这里第一、二问的骨架就有了一张关联检验表带效应量、一张分组差异表带效应量和样本量、一组按类型给出的风化前参考值。表格里的数字必须带样本量这是小样本分析的底线。4. 从聚类到判别两类玻璃的分类边界与未知样品归属第三问是全题的核心找出高钾玻璃和铅钡玻璃的分类规律再把未知样品判进去。很多人一上来就丢个 KMeans 或者随机森林跑出准确率 0.95 就收工这在答辩环节基本会被问穿——样本量只有几十条交叉验证的方差比模型间的差异还大。4.1 先用 KMeans 看自然分群再决定是否相信标签聚类的作用不是给答案而是检查“已知的类型标签”和“数据自身的结构”是否一致。如果 KMeans 在 k2 时的分群和类型标签高度吻合说明成分确实能分开这两类后面的判别模型有物理基础如果吻合度很差那要么是特征选错了要么是风化样品把结构搅乱了。from sklearn.preprocessing import StandardScaler from sklearn.cluster import KMeans from sklearn.metrics import silhouette_score, adjusted_rand_score # 只用未风化样品探索结构避免风化层干扰 clean base[base[表面风化] 未风化].dropna(subset[类型]).copy() X clr_transform(clean, oxide_cols).values for k in range(2, 5): km KMeans(n_clustersk, n_init20, random_state42).fit(X) sil silhouette_score(X, km.labels_) ari adjusted_rand_score(clean[类型], km.labels_) if k 2 else float(nan) print(fk{k} 轮廓系数{sil:.3f} 与标签的ARI{ari:.3f})n_init20是因为 KMeans 对初始质心敏感小样本下多跑几次取最优能显著降低结果随机性。random_state42保证结果可复现这在提交论文时很重要。adjusted_rand_score用来衡量聚类结果和真实标签的一致程度取值在 0 附近表示跟随机分组差不多接近 1 表示几乎完全一致。如果 k2 的轮廓系数只有 0.2 出头说明结构不明显这时候不要硬套聚类结论转而检查是不是该先做特征筛选——把氧化铅、氧化钡、氧化钾这几个判别力最强的分量单独拿出来看二维散点往往比十几维空间里的轮廓系数直观得多。这里有个反直觉的点风化样品在聚类里经常自成一簇。因为它所有成分都被稀释或富集了整体偏离本体。如果你的聚类把风化状态当成了主要分群依据那这个模型对“类型判别”毫无价值。所以探索阶段必须把风化样品排除或者把风化状态作为分层变量。4.2 决策树和 SVM 的判别规则以及未知样品的预测判别模型要解决两个需求一是准确二是可解释。赛题要求“给出分类规律”所以模型必须能翻译成人话比如“氧化铅大于某阈值判为铅钡”。决策树天然满足这一点SVM 作为对照验证非线性边界是否存在。from sklearn.tree import DecisionTreeClassifier, export_text from sklearn.svm import SVC from sklearn.pipeline import make_pipeline from sklearn.model_selection import LeaveOneOut, cross_val_score X clr_transform(base, oxide_cols).values y base[类型].values tree make_pipeline( DecisionTreeClassifier(max_depth3, min_samples_leaf3, class_weightbalanced, random_state42) ) svm make_pipeline(StandardScaler(), SVC(kernelrbf, C1.0, gammascale, class_weightbalanced)) loo LeaveOneOut() print(决策树 LOO 准确率:, cross_val_score(tree, X, y, cvloo, scoringaccuracy).mean()) print(SVM LOO 准确率:, cross_val_score(svm, X, y, cvloo, scoringaccuracy).mean()) tree.fit(X, y) print(export_text(tree[-1], feature_namesoxide_cols, decimals2))max_depth3和min_samples_leaf3是刻意压低的。样本几十条的情况下不限制深度的树会把每个样品都记住训练集准确率 100%换个数据集就崩。class_weightbalanced应对两类样品数量不平衡否则模型会偏向多数类。交叉验证用留一法而不是 5 折是因为总量太小5 折时每折验证集只有十几条方差大到没法比较模型。export_text打印出来的规则就是第三问的核心答案。你会看到类似“氧化铅 ≤ 某值 且 氧化钾 某值 → 高钾”的路径。把这些路径整理成阈值表比报一个准确率有说服力得多。对未知样品的预测要在同一套变换下做千万别对未知样品单独做 CLRXu clr_transform(unknown, oxide_cols).values pred tree.predict(Xu) prob tree.predict_proba(Xu) unknown_out unknown[[文物采样点]].copy() unknown_out[预测类型] pred unknown_out[置信度] prob.max(axis1) print(unknown_out.sort_values(置信度).head(10))注意prob.max(axis1)低的那几条它们大概率落在两类交界处是论文里必须单独讨论的对象。把它们简单归到某一类而不加说明是答辩时最容易被追问的地方。合理做法是给出“倾向铅钡但置信度 0.6建议补充检测氧化钡”这样的表述。5. 五处最容易翻车的细节填零、泄漏与阈值拍脑袋这一章是本篇的血泪经验。下面五条几乎每一届都有人踩而且踩了之后结果表面上还挺好看属于典型的黑匣子。5.1 数据层的三个坑第一条缺失值填零。现象某类玻璃的氧化铅均值只有 3%比文献里的典型值低一大截。 原因大量“未检测”被填成了 0均值被严重稀释。 解决缺失保持 NaN统计时用dropna()分组计算只有在必须构造完整矩阵比如聚类时才用极小值或中位数替代并且替代后的副本单独命名别覆盖原始表。第二条把多次采样当成多个独立样品。现象某件文物因为采了三个部位在统计里权重变成别人的三倍卡方检验结果全歪。 原因附件按采样点记录直接value_counts会让采样多的文物占主导。 解决先按文物编号聚合成分取均值类别取唯一值或者在做关联检验时用文物级数据在做成分分析时用采样点级数据两套口径分开维护。第三条百分比直接做相关分析。现象氧化铅和氧化钡的相关系数飙到 0.9 以上写成“两者共生”。 原因闭合效应一个涨另一个必然被动跌或者反过来。 解决先做 CLR 或对数比变换再算相关如果坚持用原始百分比就改用 Spearman 秩相关并明确说明局限。5.2 建模层的两个坑第四条训练集和验证集不分离。现象模型准确率 0.98换一批数据就掉到 0.7。 原因在全体样本上做了特征筛选或标准化然后再交叉验证信息从验证折泄漏进了训练过程。 解决把特征选择、变换、标准化全部塞进 Pipeline交叉验证只作用于 Pipeline 整体。或者更稳妥用嵌套交叉验证。第五条判别阈值拍脑袋。现象氧化铅阈值取 5%未知样品里有两个被判成高钾但它们的氧化铅是 4.8%。 原因阈值没有经过敏感性检验而且风化样品的氧化铅可能已被淋失同一类型风化前后的值差异很大。 解决阈值按风化状态分层设定或者改用未风化样品建立阈值同时做阈值扫描看结论在多大范围内稳定下一章展开。这五条的共同特征是出错时不会报错只会让你得到一个“看起来很合理”的错误结论。所以每一处都要留一个中间检查点把处理前后的样本量、均值、分布画出来对比别等最后写论文时才发现数字对不上。6. 敏感性分析把“判为铅钡”变成一句有底气的话第四问要求比较两类玻璃的成分关联性并对新样品提出鉴别方法真正拉开差距的是敏感性分析。评委想知道的不是“你的方法准确率多少”而是“换个假设你的结论还站得住吗”。我一般做三组扫描。第一组是阈值扫描把氧化铅的判别阈值从 1% 扫到 10%步长 0.5%记录每个阈值下未知样品的归属变化。第二组是重采样对训练集做 200 次 bootstrap每次重训决策树统计“氧化铅阈值”这条规则的取值分布。第三组是缺失处理方案对比分别用丢弃、中位数填补、极小值填补三种策略跑完整流程看最终归类是否一致。判别阈值PbO判为铅钡的未知样品数与基准方案一致数结论稳定性2.0%97偏低易把高钾误判5.0%66基准8.0%44偏保守可能漏判10.0%33明显漏判表格里的规律是阈值在 4% 到 6% 之间时结论几乎不变超过 8% 开始漏判低于 2% 开始误判。把这段区间写进结论比单给一个数字强得多——“建议阈值取 5%在 4% 至 6% 范围内判别结论稳定”。bootstrap 那组还有一个额外好处你能看到哪些样品在重采样中反复变类。这些样品就是判别边界上的模糊样本应该在论文里单独列出说明“其成分同时具备两类特征建议结合纹饰与颜色辅助判断”。至于两类玻璃的成分关联性比较我的做法是分别对高钾组和铅钡组算 CLR 后的相关矩阵再做矩阵差异检验而不是把两类混在一起算一个总相关。混算的结果会被类型这个混杂变量主导看起来所有成分都相关实际上只是两类玻璃的均值不同。最后说一句自己的习惯这道题我做完第一版之后会把所有中间表导成 CSV 存到同一个目录命名带上日期和处理策略比如chem_未风化_中位数填补.csv。因为做到第四问时你一定会想回头改前面的某个假设没有这些中间产物就得从头再跑一遍那种感觉比返工还难受。建模这件事后悔药就是中间文件。希望帮到你。本文还有配套的精品资源点击获取