ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

葡萄酒评价统计建模全链路:数据清洗、评分信度与AHP分级

葡萄酒评价统计建模全链路:数据清洗、评分信度与AHP分级

简介:面向备战数学建模竞赛、尤其是研究2012年A题葡萄酒质量评价的读者,这份一等奖论文PDF完整呈现从评分信度检验到葡萄分级、理化指标关联与质量预测的建模思路。压缩包仅含1个PDF文件,大小约2.41MB,便于直接查阅与打印精读。论文以K-S检验、Wilcoxon符号秩检验和肯德尔和谐系数判断两组评酒员差异与可信度,用主成分分析构建葡萄分级体系,并以典型相关分析和多元线性回归讨论理化指标对葡萄酒质量的影响,同时交代MATLAB、SPSS、SAS与EXCEL的实现路径。该文已有2908人学习,适合参赛者借鉴论文结构、模型选择、假设检验表述与结果论证方式,也可作为统计建模课程案例参考。

1. 一套被反复下载的论文,真正值钱的是那条数据处理链路

2012 年数学建模 A 题给的是同一批葡萄的两套数据:人喝出来的分,和仪器测出来的理化值。一等奖论文《葡萄酒的评价》把这两套数据接了起来——单样本 K-S 检验定分布,Wilcoxon 符号秩检验判两组评酒员有没有系统性差异,肯德尔和谐系数比谁更可信,主成分加 AHP 把二十多个葡萄指标压成一个综合分做分级,最后用典型相关和多元回归回答「理化指标能不能替掉品酒员」。它不是一篇读完就完的范文,而是一条可以照着跑的数据流水线:三份附件、四个问题、五种统计模型,每一步都有明确的输入、输出和判定阈值。适合正在做统计建模课设的人,也适合手里攥着一批评分数据、想搞清楚评分者信度到底怎么算的从业者。

2. 三份附件先过一遍:字段口径、缺失填补与异常修正

动手之前得先把数据摸清楚。这套题最容易劝退人的地方不是模型,是附件 1 到附件 3 的字段结构完全不是一个形态:一个是人打分,两个是仪器测值;一个按评酒员排,一个按样本排。不先把它们对齐成「样本 × 指标」的矩阵,后面任何一行代码都跑不通。

2.1 附件 1/2/3 的字段结构和量纲差异

三份附件的组织方式差别很大,先把结构列清楚再谈建模:

附件观测主体主要字段规模取值与量纲
附件 1评酒员 × 酒样外观(澄清度、色调)、香气(纯正度、浓度、质量)、口感(纯正度、浓度、持久性、质量)、整体评价红葡萄酒 27 个样本、白葡萄酒 28 个样本,两组各 10 名评酒员整数等级分,各分项上限不一致
附件 2酿酒葡萄样本氨基酸总量、蛋白质、VC 含量、花色苷、酒石酸、苹果酸、柠檬酸、多酚氧化酶活力、褐变度、DPPH 自由基 1/IC50、总酚、单宁、总糖、可溶性固形物、干物质、果皮质量、果穗质量、色差 a*/b*、可滴定酸、固酸比、黄酮醇等红葡萄 27 行单位混杂:g/L、百分比、色差为无量纲
附件 3葡萄酒样本花色苷、单宁、总酚、DPPH 半抑制体积、酒精度、pH 等与附件 1 一一对应部分为体积分数,部分为浓度

这张表要盯住两点。第一,附件 1 的行是评酒员、列是样本,而附件 2/3 的行就是样本,转置和对齐是第一步;第二,附件 2 的指标量纲跨度极大,总糖在 g/L 量级,色差 a* 在个位数量级,直接丢进主成分,方差全被总糖吃掉,标准化不是可选项而是必做项。

常见做法是先写一张字段映射表,把「样本编号」作为唯一主键,三份数据全部左连接到这张主键上。别小看这一步,红葡萄 27 个样本、白葡萄 28 个样本,编号只要有两位数字不一致(01 和 1),后面的配对检验就会静默出错——不报错,但结果全歪。

提示:白葡萄酒比红葡萄酒多一个样本,所有按样本循环的代码都要写成按组分别处理,别用一个写死的 27。

2.2 均值替换的两个前提:组内取均值、结果取整

原文在附件 1 里发现了两类脏数据,处理方式不一样。第一类是缺失:红葡萄酒样品 20 的色调分,4 号评酒员那一格是空的,其余 9 人打的是 6、6、4、6、6、8、6、6、8,剔掉后取均值再取整,填 6。第二类是异常:第一组白葡萄酒样品 3 的持久性分,7 号评酒员打了 77,而相邻评酒员都在 5 到 7 之间,明显是手工录入多敲了一个 7。

填补代码本身很短,关键在于「按什么范围取均值」:

import pandas as pd import numpy as np def fill_by_item_mean(df, as_int=True): """df: 行=评酒员,列=样本,值=某一分项的评分 只在该样本的其余评酒员范围内取均值,不做全表平均""" filled = df.copy() for col in filled.columns: s = filled[col] if s.isna().any(): m = s.mean(skipna=True) # 同一列 = 同一样本 filled[col] = s.fillna(round(m) if as_int else m) return filled def flag_out_of_range(df, low, high): """评分分项有明确上下界,越界直接判异常,返回待人工复核的位置""" mask = (df < low) | (df > high) return df.where(mask) # 非异常位置为 NaN,方便定位

fill_by_item_mean里最要紧的是s.mean()的作用域。这里的s是某一分项下所有样本的列,如果数据结构是「行=样本、列=评酒员」,那s.mean()就是跨评酒员取均值,方向对了;如果行列反了,就变成跨样本取均值,等于把不同酒样的水平混在一起,方差被抹平,后面 Wilcoxon 的配对差值分布会跟着变。

as_int=True不是洁癖。评酒员打的是整数等级分,如果插出一个 6.4,Wilcoxon 编秩时这一格会和别的样本产生额外的并列秩,秩和统计量变化,z 值可能和原文差零点几。原文报的红葡萄酒 z 约 -2.53、白葡萄酒 z 约 -2.23,都在 1.96 的临界线外侧,差一点点不影响结论方向,但换一批数据就未必了。

2.3 异常值走「先剔除后替换」,别用 3σ 一把梭

原文对 77 这种数据的处理是「先剔除,后替换」:先置为缺失,再走 2.2 的均值填补。这个顺序比直接改数要稳,因为它保留了「这里原来有问题」的痕迹。

我一般不会对评分数据用 3σ 准则。原因很直接:评分样本量小(一个分项就 10 个评酒员),标准差本身估不准,3σ 的边界会随样本漂;而且评分有硬上下界,用领域规则比统计规则可靠得多。分项的上限是已知的,超上限就是错,低于下限也是错,这种判断不需要统计。

真正需要小心的是「看起来不合理但其实合理」的数据。比如某位评酒员对某个样品的整体评价给了很低的分,而同组其他人给得高。这不一定是异常,可能是这位评酒员对这款酒真的不买账——评酒员之间的评价尺度差异本来就是这套题要研究的东西,把它当异常清掉,等于人为制造了一致性。

判断标准可以这样定:越界(超过分项上限、低于 0、出现两位数)的判异常;组内偏离中位数超过该分项极差一半的,标出来人工看一眼,别自动改。

2.4 把评分表压成样本级配对矩阵

问题一要比较「两组评酒员的评价结果」,比较对象是样本级的。所以中间要做一次聚合:10 名评酒员的分项分先求和成每人的总分,再取组内均值作为该样本在该组的代表值。

item_cols = ["澄清度", "色调", "纯正度_香", "浓度_香", "质量_香", "纯正度_口", "浓度_口", "持久性", "质量_口", "整体评价"] # 1) 每人总分:等权求和,和原文口径一致 score = raw.groupby(["组别", "样本", "评酒员"])[item_cols].sum().reset_index() score = score.rename(columns={0: "总分"}) if "总分" not in score else score # 2) 样本级:每组每个样本一个代表值,红白分开 red = score.query("酒种 == '红'").pivot(index="样本", columns="组别", values="总分").mean(axis=1)

第一段是分项等权求和。原文的口径就是这样,虽然香气质量、口感质量这些分项的分值上限本来就比澄清度高,等权求和实际上已经隐含了权重。如果自己再乘一套权重,问题一的 z 值会变,但方向通常不变。

第二段容易写错。pivot之后如果直接把 10 名评酒员的行铺开,得到的是 10 行 × 27 列的宽表,Wilcoxon 需要的是「同一批样本、两组各一个值」的配对序列,长度必须是样本数。所以这里先用groupby把评酒员维度压掉,再 pivot 出两组,最后两列一一对应。

注意:压掉评酒员维度会丢失组内离散程度。问题一如果还想讨论「组内分歧」,得另外保留标准差,别在这一步扔掉。

3. 评分差异与葡萄分级:Wilcoxon 编秩到 AHP 赋权的完整链路

问题一和问题二看着是两件事,实际上是一条链:先用假设检验确认两组评分不可混用,选出一组相对可信的分数作为「质量真值」;再把这个真值当作评价体系里的一级指标,和葡萄的理化指标一起做综合评价,最后落到分级。链断在哪一环,后面的结论都站不住。

3.1 单样本 K-S 检验:p=0.525 说明的是「不拒绝」

原文对四组评分数据各做了一遍正态性检验,用的是 P-P 图和单样本 K-S,第一组红葡萄酒得到 p=0.525,大于 0.05,于是认为评分服从正态分布,进而说明评分数据是科学的、合理的。

这个推理链条要打折看。K-S 检验的原假设是「样本来自指定分布」,p 大于 0.05 只能说不拒绝原假设,不等于证明了正态。而且这里的样本是 27 个样本均值,样本量小,检验功效低,本来就很难拒绝。

from scipy import stats x = red["第一组"].values # 27 个样本的组内均分 z = (x - x.mean()) / x.std(ddof=1) d, p = stats.kstest(z, "norm") # 参数由样本估计,严格说要用 Lilliefors 修正 print(f"D={d:.4f}, p={p:.4f}, n={len(x)}")

kstest这里传的是标准化后的数据和标准正态,等价于检验正态性,但临界值表是按「分布参数已知」算的。参数是从样本估出来的,p 值会偏大(偏保守),专业的做法是用 Lilliefors 修正或者直接上 Shapiro-Wilk。原文用 SPSS 的 K-S 选项,同样存在这个偏差。

不过这一步在整条链路里不是关键。问题一后面用的是 Wilcoxon 符号秩检验,本身就是非参数方法,不要求正态。K-S 在这里更多是交代「为什么选这个模型」的叙事作用。复现时对不上 p 值不用纠结,对上了反而是巧合。

3.2 Wilcoxon 符号秩检验:配对、编秩与 z 值的算法

两组评酒员评的是同一批酒样,这是天然的配对数据,所以用两配对样本的非参数检验。原文的步骤写得很规范:算成对差值 D、按绝对值编秩、恢复符号、分别求正秩和 T+ 与负秩和 T−、取较小者为统计量 T、再用正态近似算 z。

以红葡萄酒为例,n=27,E(T)=n(n+1)/4=189,D(T)=n(n+1)(2n+1)/24≈1732.5,标准差约 41.6,原文报的 z≈-2.53,反推 T 在 84 左右,落在拒绝域里。白葡萄酒 z≈-2.23,同样拒绝。结论是两组对红、白葡萄酒的评价都存在显著性差异。

Python 里一行就能跑:

from scipy.stats import wilcoxon stat, p = wilcoxon(red["第一组"], red["第二组"], zero_method="wilcox", # 差值为 0 的样本如何处理 correction=False, # 是否做连续性校正 alternative="two-sided") print(f"T={stat}, p={p:.4f}")

三个参数值得逐个说。zero_method决定差值为 0 的配对怎么算:wilcox直接丢弃、pratt保留进秩、zsplit把零差值平分到正负两侧。27 个样本里出现一两个零差值不稀奇,选哪个会改变有效 n,进而改变 z。原文没有明确交代这一点,复现时先统计一下零差值个数,写进方法说明里。

correction是连续性校正,样本量小的时候打开会让 z 往 0 缩一点,更保守。27 对数据不算大,我一般会开着跑一遍,和不开的结果对照,看结论是否翻转。

alternative用双侧,因为原文的备择假设就是「两组评价结果不同」,没预设方向。

还要提醒一句:wilcoxon返回的stat是较小的秩和,不是 z。想要和原文的 -2.53 对上,得自己按 n 换算:z = (T - n*(n+1)/4) / sqrt(n*(n+1)*(2*n+1)/24)。第一次复现时对不上数很正常,不是代码写错了。

3.3 肯德尔和谐系数 W:四组信度比较与卡方边界

差异检验只能说明两组不一样,回答不了「哪组更可信」。原文换了个角度:把评分转成秩,看组内 10 名评酒员对 27 个样本的排序有多一致。一致性强,说明这组评酒员的评价尺度稳定,可信度更高。

公式是 W = 12Q / [m²(n³−n)],其中 Q 是各样本秩和与平均秩和之差的平方和,m 是评酒员数(10),n 是样本数(27 或 28)。检验统计量 χ² = m(n−1)W。

import numpy as np from scipy.stats import rankdata def kendall_w(mat): """mat: 行=评酒员,列=样本,值=评分 返回和谐系数 W 与卡方统计量""" m, n = mat.shape R = rankdata(mat, axis=0).sum(axis=0) # 每个样本的秩和 Q = ((R - R.mean()) ** 2).sum() W = 12 * Q / (m ** 2 * (n ** 3 - n)) return W, m * (n - 1) * W

rankdata(mat, axis=0)是按列(样本)在评酒员之间编秩,方向别搞反。默认method="average",遇到并列分数取平均秩,和原文的做法一致。

原文算出的四组 W 分别是 0.044366788、0.043314652、0.040039039、0.041112055,对应 χ² 为 11.535、11.262、10.811、11.100。红葡萄酒第一组的 W 和 χ² 都略高,白葡萄酒第二组略高,于是得出「红葡萄酒第一组更可信、白葡萄酒第二组更可信」,并在后续分析中分别采用这两组分数。

这里有一个必须说清楚的边界。χ²=11.535,自由度按 n−1=26 算,0.05 水平的临界值约 38.9;按 27 算也在 40.1 左右。四组数据的 χ² 都远没有越过这道线,也就是说,没有一组达到「评分者显著一致」的统计门槛。W=0.044 这个量级本身就说明 10 名评酒员之间的排序一致性很弱。

原文的用法是横向比大小——四个值里谁高谁可信。这种相对排序在方案比较里是成立的,但不能顺着写成「评分具有显著一致性」。做工程判断时这反而是有用信息:两组评酒员的内部分歧都很大,「谁更可信」只是矮子里拔将军。复现时把这句话改准,能省掉答辩时一轮追问。

3.4 主成分降维:27 个理化指标压到 8 个的阈值

问题二要拿葡萄的理化指标做分级,但指标太多且彼此相关——总酚和单宁相关、总糖和可溶性固形物相关,直接加权会重复计权、放大某类信息。原文用主成分分析先把 27 个指标压成 8 个互不相关的主成分,累计方差贡献率 83.044%。

from sklearn.preprocessing import StandardScaler from sklearn.decomposition import PCA X = StandardScaler().fit_transform(grapes) # 行=样本,列=理化指标 pca = PCA(n_components=0.83044, svd_solver="full") Z = pca.fit_transform(X) print(pca.n_components_, pca.explained_variance_ratio_.cumsum()[-1])

n_components传浮点表示「累计方差比例达到这个值就停」,传整数表示固定个数。原文是按 83.044% 截的,正好落在第 8 个主成分上,所以浮点阈值能对上;如果传 8 也等价,但换一批数据就不一定了,浮点写法更通用。

StandardScaler必须放在 PCA 前面。附件 2 里总糖是 g/L 量级、色差是无量纲个位数,不标准化的话第一主成分基本就是总糖,其他指标全被压成噪声。

主成分解读还有个坑。原文报的是旋转后的载荷(方差最大化旋转),所以能说出「主成分 1 是总黄酮、总酚、DPPH 自由基和单宁的组合」「主成分 2 是总糖、可溶性固形物和干物质的组合」这种话。sklearn 默认不旋转,载荷矩阵和 SPSS 的输出长得不一样,同一个主成分对应的指标组合可能变。复现时要么补一步 varimax 旋转,要么只报载荷绝对值排前几位的指标,别硬套原文的主成分命名。

3.5 AHP 判断矩阵与阶梯分级:CR=0.029 之后还剩什么

指标压完了,接下来要定权重。原文用层次分析法构造了一个 8×8 的判断矩阵,最大特征根为 8,CI≈0.013,CR≈0.029,小于 0.1,通过一致性检验。归一化后的权重分布在 0.0253 到 0.0987

本文还有配套的精品资源,点击获取

RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进