简介:这份2012年全国大学生数学建模竞赛A题一等奖论文,围绕葡萄酒质量评价问题展开,适合数模参赛者、统计建模学习者及需要主观评分与多指标分析案例的研究者参考。全文以一份PDF文档呈现,压缩包内共1个PDF文件,约2.41MB,篇幅完整,便于直接查阅与打印。论文针对评酒员评分差异、酿酒葡萄分级、葡萄与葡萄酒理化指标关联及质量预测四个问题,综合运用K-S检验、Wilcoxon符号秩检验、肯德尔和谐系数、主成分分析、典型相关分析和多元线性回归等方法,并借助MATLAB、SPSS、SAS与EXCEL完成计算。读者可从中获取完整的赛题建模框架、指标筛选与模型检验思路,以及一等奖论文的摘要撰写、问题重述与模型假设组织方式。目前已有2907人学习下载,适合作为赛前研读、方法复现和统计建模思路迁移的参考范例。
1. 评酒员打分的黑匣子里,藏着一条完整的统计链路
十个人给同一杯酒打分,最高 9 分最低 4 分,直接取平均排名基本是玄学。2012 年数学建模 A 题要处理的正是这堆主观评分:红葡萄酒 27 个样本、白葡萄酒 28 个样本,每组 10 名评酒员,附件里还有二十多项酿酒葡萄理化指标和近十项葡萄酒理化指标。这份一等奖论文《葡萄酒的评价》真正值钱的地方不是最后那句"样本 23 最优",而是它把"人打的分"和"仪器测的指标"用一条统计链路串了起来:K-S 检验定分布、Wilcoxon 符号秩检验判差异、肯德尔和谐系数评信度、主成分加层次分析做分级、典型相关和多元回归打通葡萄与酒。如果你手头正好有主观评分表、又需要给出可复现的一致性结论,这套流程可以直接拆开用。
2. 评分数据先洗后检:K-S 定分布、Wilcoxon 判差异、和谐系数评信度
2.1 原始评分表的两处脏点:缺失值与突变值
附件 1 是典型的"多人打分矩阵",行是酒样、列是评酒员,单元格里是各分类指标得分。这类表在录入阶段一定会出两种问题:一是漏填,二是手误多敲一位。论文里举的两个例子很典型——红葡萄酒样品 20 的色调数据里,4 号评酒员的格子是空的;第一组白葡萄酒样品 3 的持久性评分里,7 号评酒员给出 77 分,而同一行其他人都是 5 到 7 分。
处理缺失用的是列内均值替换,而不是整表均值。原因不难想:不同评酒员的打分尺度本来就不一样,有人偏爱给高分,有人手紧,跨评酒员取均值会把尺度差异混进来。异常值走的是"先剔除、后替换"——把那个格子置空,再用同一列剩余数据的均值回填,样本本身不删。删样本是下策,27 个红葡萄样本删掉一个,后面主成分的样本量就少一个,累计贡献率都会跟着变。
| 品酒员 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|---|---|---|---|---|---|---|---|---|---|---|
| 修补前 | 6 | 6 | 4 | --- | 6 | 6 | 8 | 6 | 6 | 8 |
| 修补后 | 6 | 6 | 4 | 6 | 6 | 6 | 8 | 6 | 6 | 8 |
替换值取整是有讲究的,评分是整数尺度,回填 6.375 这种小数会让后续的秩次计算出现并列,Wilcoxon 里处理并列秩要额外取平均秩,平白增加复杂度。
2.2 单样本 K-S 检验:P-P 图和渐近显著性怎么读
洗完之后第一件事不是直接做差异检验,而是确认这批数据能不能用参数检验。论文的做法是先算每个酒样 10 名评酒员的评分均值,得到 27 个(红)或 28 个(白)样本均值,再对这批均值做正态性检验。
SPSS 里走"分析 → 非参数检验 → 旧对话框 → 1-样本 K-S",勾选正态分布,输出的 P-P 图看两点:散点是否沿着对角线分布、偏离对角线最大处的垂直距离就是 K-S 统计量 D。论文给出的红葡萄酒第一组结果是 p = 0.525,远大于 0.05,不拒绝正态假设。
这里有个容易被忽略的细节:K-S 检验的标准临界值表假设分布参数已知,而实际是拿样本估计的均值和标准差,直接查表会偏保守。SPSS 内部做了 Lilliefors 修正,Python 里对应的是statsmodels的lilliefors,不是scipy.stats.kstest。
import numpy as np import pandas as pd from scipy import stats from statsmodels.stats.diagnostic import lilliefors # df:行 = 酒样,列 = 10 名评酒员,值是品尝总分 df = pd.read_excel("red_group1.xlsx", index_col=0) # 第一步:按列(评酒员)回填缺失,避免跨评酒员的尺度污染 df = df.apply(lambda c: c.fillna(round(c.mean())), axis=0) # 第二步:用中位数绝对偏差识别突变值,比 3σ 更抗异常点 def repair(col, k=3.5): med = col.median() mad = (col - med).abs().median() if mad == 0: return col out = (col - med).abs() > k * 1.4826 * mad col = col.mask(out) # 先剔除 return col.fillna(round(col.mean())) # 再按列均值替换 df = df.apply(repair, axis=0) # 第三步:每个酒样取组内均值,作为该样本的代表性评分 sample_score = df.mean(axis=1) # 第四步:Lilliefors 修正版 K-S 检验,等价于 SPSS 的单样本 K-S d_stat, p_val = lilliefors(sample_score.values, dist="norm") print(f"D={d_stat:.4f}, p={p_val:.4f}")参数说明:k=3.5是 MAD 的经验阈值,对应正态下约 3σ 的等价强度;1.4826是把 MAD 换算成正态标准差的一致性系数。如果某个评酒员整列打分都很平(mad 为 0),函数直接跳过,这种情况在只打整数分的表里很常见。
为什么不用 3σ?因为那个 77 分本身就把该列的标准差撑大了,27 分的偏离在 10 个观测里未必超过 3 倍标准差,异常值会把自己"藏"进正常区间。MAD 基于中位数,不受极端值影响,这是血泪经验。
2.3 Wilcoxon 符号秩检验:配对数据的非参数解法
分布确认之后进差异检验。两组评酒员评的是同一批酒样,这是标准的配对结构,理论上该用配对 t 检验。但论文放弃 t 检验的理由写得很清楚:配对 t 检验要求差值服从正态分布,而评酒员打分是感官评价,主观性强,差值分布没理由假设正态。
Wilcoxon 符号秩检验只要求差值关于中位数对称,门槛低得多。计算流程是:先算每对观测的差值 D,取绝对值排序编秩(并列取平均秩),恢复正负号后分别求正秩和 T+ 与负秩和 T-,取较小者作为统计量 T。
红葡萄酒 n = 27 时,E(T) = n(n+1)/4 = 189,D(T) = n(n+1)(2n+1)/24 = 1732.5,标准差约 41.62。论文算出的 z ≈ -2.53,反推 T ≈ 84。白葡萄酒 n = 28,E(T) = 203,标准差约 43.9,z ≈ -2.23。两组都超过 1.96 的临界值,结论是红、白葡萄酒的两组评分均存在显著性差异。
from scipy import stats g1 = df_red_group1["总分"] # 第一组评酒员对 27 个红葡萄样本的总分 g2 = df_red_group2["总分"] # 第二组,样本顺序必须与 g1 一一对应 stat, p = stats.wilcoxon( g1, g2, zero_method="wilcox", # 差值为 0 的样本直接丢弃,不参与编秩 alternative="two-sided", # 双侧检验,对应 α=0.05 下 |z|>1.96 method="approx" # 正态近似;样本量小于 20 时改成 "exact" ) print(f"T={stat}, p={p:.4f}")参数说明:zero_method有三个取值,wilcox是标准做法;pratt会把零差值也编秩,样本里零差值多时结论会变,不要随手改。method在新版 scipy 里替代了旧的mode,老代码跑不通多半是这个原因。
2.4 肯德尔和谐系数:判断哪一组评分更可信
差异检验只能回答"两组不一样",回答不了"哪组更可信"。论文的做法是把评分的具体数值丢掉,只保留秩次信息,用肯德尔和谐系数 W 度量组内 10 名评酒员的一致性。
W 的计算是 W = 12S / [m²(n³ − n)],其中 m 是评酒员数,n 是样本数,S 是各样本秩和与平均秩和之差的平方和。检验统计量 χ² = m(n − 1)W,自由度 n − 1。
| 组别 | W 值 | χ² | 自由度 |
|---|---|---|---|
| 第一组红葡萄酒 | 0.04437 | 11.535 | 26 |
| 第二组红葡萄酒 | 0.04331 | 11.262 | 26 |
| 第一组白葡萄酒 | 0.04004 | 10.811 | 27 |
| 第二组白葡萄酒 | 0.04111 | 11.100 | 27 |
数值可以对上:10 × 26 × 0.04437 ≈ 11.54,和表里的 11.535 一致。红葡萄酒第一组的 W 和 χ² 都大于第二组,白葡萄酒则反过来,所以论文选了"红用第一组、白用第二组"作为后续建模的评分输入。
def kendall_w(mat): # mat:行 = 评酒员,列 = 酒样,值是总分 m, n = mat.shape ranks = np.apply_along_axis(stats.rankdata, 1, mat) # 每个评酒员内部单独编秩 R = ranks.sum(axis=0) # 各酒样的秩和 S = ((R - R.mean()) ** 2).sum() W = 12 * S / (m ** 2 * (n ** 3 - n)) chi2 = m * (n - 1) * W # 自由度 n-1 p = stats.chi2.sf(chi2, n - 1) return W, chi2, pnp.apply_along_axis(stats.rankdata, 1, mat)这一行是关键,沿行方向编秩意味着每个评酒员的 27 个打分自己排一遍,消除个人尺度差异。如果图省事对整个矩阵展平编秩,算出来的 W 会完全失真。
还有一点必须说清楚:W 的绝对值只有 0.04 量级,这意味着评酒员组内的一致性其实很弱,χ² 也没超过 0.05 显著性水平的临界值(26 自由度约 38.9,27 自由度约 40.1)。论文用 W 的相对大小来选组是合理的,但把它读成"评分高度一致"就过头了。这个坑我后面单独拎出来讲。
3. 27 个指标压成 8 个主成分:PCA 遴选加 AHP 定权的分级模型
3.1 先标准化,再算相关系数矩阵
问题二要给酿酒葡萄分级,输入是 27 项理化指标:氨基酸总量、蛋白质、VC、花色苷、酒石酸、苹果酸、柠檬酸、多酚氧化酶活力、褐变度、DPPH 自由基、总酚、单宁、总糖、可溶性固形物、干物质含量、果皮质量、果穗质量、色差 a*、b*、可滴定酸、固酸比、黄酮醇等。这些指标量纲五花八门,有的是百分比,有的是 mg/L,有的是色差单位,不标准化直接做主成分,方差大的指标会独吞所有载荷。
标准化用 z-score:F~ij = (Fij − Fi) / si,减均值除标准差。标准化不改变变量之间的相关系数,只消除量纲。
然后算 27×27 的相关系数矩阵 R。论文在相关性系数表里发现多个指标高度相关,比如总酚和 DPPH 自由基、总糖和可溶性固形物,直接用 27 个指标做评价会造成信息重叠,同一件事被重复计权两次。
3.2 主成分提取与正交旋转后的指标归并
论文以累计方差贡献率 83.044% 为门槛,红葡萄提取出 8 个主成分。原始载荷矩阵里各主成分含义很模糊,做最大方差正交旋转(varimax)之后才清晰:
| 主成分 | 归并的主要指标 | 解释方向 |
|---|---|---|
| 1 | 葡萄总黄酮、总酚、DPPH 自由基、单宁 | 抗氧化与酚类物质 |
| 2 | 总糖、可溶性固形物、干物质含量 | 成熟度与糖分 |
| 3 | 苹果酸、褐变度 | 酸度与氧化程度 |
| 4 | 果皮质量、果穗质量 | 果穗结构 |
| 5 | 色差 a*、b* | 色泽 |
| 6 | 可滴定酸、固酸比 | 总酸水平 |
| 7 | 黄酮醇 | 次生代谢物 |
| 8 | 酒石酸 | 特定有机酸 |
这张表的用途是给后续 AHP 判断矩阵提供依据——同一主成分里的指标可以视为同向,跨主成分的指标才需要认真比较重要度。
import pandas as pd from sklearn.preprocessing import StandardScaler from factor_analyzer import FactorAnalyzer # grape_df:行 = 27 个红葡萄样本,列 = 27 项理化指标 X = grape_df.values Z = StandardScaler().fit_transform(X) # z-score,等价于论文里的 F~ fa = FactorAnalyzer( n_factors=8, # 由累计贡献率 ≥ 80% 反推得到 rotation="varimax", # 最大方差正交旋转,让载荷向 0/1 两极分化 method="principal" # 主成分法,对应 SPSS 的"主成分"提取 ) fa.fit(Z) loadings = pd.DataFrame(fa.loadings_, index=grape_df.columns) variance = pd.DataFrame(fa.get_factor_variance(), index=["方差", "方差贡献率", "累计贡献率"]).T print(variance) print(loadings.round(3))参数说明:n_factors=8不是拍脑袋定的,标准做法是从相关矩阵特征值大于 1 的成分数起步,再结合累计贡献率 80% 以上的要求微调。rotation="varimax"会改变载荷矩阵但不改变累计贡献率,所以贡献率表在旋转前后可以对照看。如果环境里装不上factor_analyzer,用sklearn.decomposition.PCA也能拿到主成分得分,但它不做旋转,载荷没法直接命名,这是复现时最常见的翻车点。
3.3 指标一致化与极差无量纲化
进入综合评价之前还有两道预处理。第一道是一致化,27 项指标里混着极大型、极小型和中间型。极小型用x' = M − x平移翻转;中间型更麻烦,要按指标最优值分两段做线性变换,低于最优值和高于最优值分别映射。第二道是极差化,x* = (x − m) / (M − m),把值压到 [0,1]。
极差化有个副作用要提前知道:它对极端值极其敏感,一个离群样本会把整个区间拉开,其余样本全挤在中间一小段。所以第 2 章那套异常值处理在这里是有回报的,前面不洗,后面所有样本的综合评价值都会失真。
3.4 AHP 判断矩阵、一致性检验与综合评价值
权重的确定用层次分析法。论文构造了一个 8 阶判断矩阵,行列对应 8 个主成分,元素aij表示第 i 个因子相对第 j 个因子的重要程度,对角线为 1,aij = 1/aji。矩阵里因子 1 对因子 8 的比值是 2.7346,因子 2 对因子 1 是 0.8580。
求最大特征根对应的特征向量,归一化之后就是权重。论文给出的 27 项指标权重里,口感"质量"是 0.0987,"整体评价"是 0.0915,是所有三级指标里最高的两个;理化指标里"苹果酸"0.0316、"总糖"0.0302 相对靠前,权重最低的是"单宁"0.0247。
一致性检验是这一步不能省的动作。CI = (λmax − n)/(n − 1),CR = CI / RI,8 阶矩阵的 RI 取 1.41。论文算出 CI ≈ 0.013,单层 CR 远小于 0.1,三层组合一致性比率 CR = 0.029 < 0.1,判断矩阵可用。
import numpy as np RI = {1: 0.0, 2: 0.0, 3: 0.58, 4: 0.90, 5: 1.12, 6: 1.24, 7: 1.32, 8: 1.41, 9: 1.45} A = np.array([ [1.0000, 1.1655, 1.6571, 1.9159, 1.9311, 2.0250, 2.3020, 2.7346], [0.8580, 1.0000, 1.4218, 1.6439, 1.6569, 1.7375, 1.9751, 2.3463], [0.6035, 0.7033, 1.0000, 1.1562, 1.1654, 1.2221, 1.3892, 1.6502], [0.5219, 0.6083, 0.8649, 1.0000, 1.0079, 1.0569, 1.2015, 1.4273], [0.5178, 0.6035, 0.8581, 0.9921, 1.0000, 1.0486, 1.1921, 1.4161], [0.4938, 0.5755, 0.8183, 0.9461, 0.9536, 1.0000, 1.1368, 1.3504], [0.4344, 0.5063, 0.7198, 0.8323, 0.8389, 0.8797, 1.0000, 1.1879], [0.3657, 0.4262, 0.6060, 0.7006, 0.7062, 0.7405, 0.8418, 1.0000], ]) n = A.shape[0] eigvals, eigvecs = np.linalg.eig(A) k = np.argmax(eigvals.real) # 最大特征根的位置 lam = eigvals.real[k] w = np.abs(eigvecs[:, k].real) w = w / w.sum() # 归一化得到权重向量 CI = (lam - n) / (n - 1) CR = CI / RI[n] print(f"lambda_max={lam:.4f}, CI={CI:.4f}, CR={CR:.4f}")np.abs()这步不能漏,特征向量方向不唯一,numpy 求出来可能是负的。CR < 0.1才继续,不满足就得回去改矩阵元素重新比较,硬算下去权重没有意义。
综合评价值用线性加权求和,三层权重逐级合成:三级指标加权成二级,二级加权成一级,一级加权成总分 y。评价结果和等级正相关,y 越大等级越高。分级阶梯按 0.1 一档切:[0.8,1]特级,[0.7,0.8)一级,依此类推到[0,0.2)七级。
论文的求解结果里,27 个红葡萄样本大部分落在二、三级,样本 23 的评价值 0.80 正好卡在特级下界,样本 12 只有 0.23,属于六级。样本 23 之所以能上特级,看主成分 1 和 2 的载荷就明白了——它的酚类物质和糖分指标都处在高位,两个权重最大的方向同时占优。
import pandas as pd bins = [0, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0] labels = ["7th", "6th", "5th", "4th", "3rd", "2nd", "1st", "super"] grade = pd.cut(comprehensive_score, bins=bins, labels=labels, right=False, include_lowest=True)right=False对应论文区间表里[0.8,1]这种左闭右开的写法,如果写成默认的右闭,0.80 会落进一级而不是特级,样本 23 的结论直接反过来。
4. 两组理化指标的桥:典型相关分析与多元回归的接续建模
4.1 简单相关为什么不够用
问题三要分析酿酒葡萄理化指标和葡萄酒理化指标之间的联系。两组变量各自都有十几到二十几项,如果一项一项算 Pearson 相关系数,会得到几百个数字,其中大量是伪相关——葡萄的总糖和酒的总糖当然相关,因为发酵本身就在转化它们,这种相关没有解释价值。
更要命的是组内相关。葡萄侧的指标彼此高度相关,酒侧的指标也彼此高度相关,简单相关分析孤立看待单个 X 和单个 Y,完全忽略组内结构。典型相关分析的做法是分别对两组变量做线性组合,找出使两组线性组合之间相关系数最大的那对权重向量,这样一次就能把组内冗余全部吸收掉。
4.2 典型相关分析:矩阵构造与显著性检验
设葡萄侧指标为 X(p 维),酒侧指标为 Y(q 维),样本 n 个。先对两组数据各自标准化,算相关系数矩阵 R,并按块划分:
R = [[R11, R12], [R21, R22]]R11 是葡萄侧指标内部相关,R22 是酒侧内部相关,R12 和 R21 是跨组相关。构造M1 = R11⁻¹ R12 R22⁻¹ R21,求它的特征值 λ²,特征值开方就是典型相关系数 ρ;对应的特征向量归一化后就是典型变量的权重。
典型相关系数不是算完就完事,还要做 Bartlett 渐近 χ² 检验,原假设是第 k 个及以后的所有典型相关系数为 0。只有通过检验的那几对典型变量才有解释意义。
import numpy as np from scipy import stats from sklearn.preprocessing import StandardScaler def canonical_corr(X, Y): Xs = StandardScaler().fit_transform(X) Ys = StandardScaler().fit_transform(Y) p, q = X.shape[1], Y.shape[1] Rxx = np.corrcoef(Xs, rowvar=False) # R11 Ryy = np.corrcoef(Ys, rowvar=False) # R22 Rxy = np.corrcoef(Xs, Ys, rowvar=False)[:p, p:] # R12 Ryx = Rxy.T # R21 M = np.linalg.inv(Rxx) @ Rxy @ np.linalg.inv(Ryy) @ Ryx vals = np.sort(np.linalg.eigvals(M).real)[::-1] return np.sqrt(np.clip(vals, 0, 1)) def bartlett_test(rho, n, p, q): # 依次检验第 k 对及之后的典型相关是否为零 out = [] for k in range(len(rho)): lam = np.prod(1 - rho[k:] ** 2) # Wilks' Lambda chi2 = -(n - 1 - (p + q + 1) / 2) * np.log(lam) df = (p - k) * (q - k) out.append((k + 1, chi2, df, stats.chi2.sf(chi2, df))) return out rho = canonical_corr(grape_df.values, wine_df.values) for row in bartlett_test(rho, n=27, p=27, q=9): print(f"第{row[0]}对: chi2={row[1]:.3f}, df={row[2]}, p={row[3]:.4f}")np.clip(vals, 0, 1)是防御性写法,数值误差会让特征值出现 -1e-16 这种负数,开方直接报错。Bartlett 检验里(p + q + 1) / 2是修正项,样本量小时不能省。
论文通过典型相关得到三条具体结论:一是增大酿酒葡萄果皮含量,对葡萄酒中 DPPH 半抑制体积含量的增加有重要影响;二是苹果酸既促进发酵,又能保护红葡萄酒中起主要呈色作用的花色苷和起辅助辅色作用的单宁;三是葡萄总黄酮清除自由基的抗氧化作用与总酚的保护清除作用共同作用,把 DPPH 自由基转化为葡萄酒中的 DPPH 半抑制体积。这三条结论如果只做简单相关,是挖不出来的。
4.3 指标剔除与多元线性回归
问题四要论证"能不能用理化指标评价葡萄酒质量"。论文的处理顺序值得学:先在问题三的典型相关结果基础上,保留葡萄酒指标,剔除葡萄指标中那些被认为"已经被对应酒指标完全代表"的部分,降低共线性;筛选后的指标进多元线性回归,因变量是葡萄酒质量评分。
这里的关键判断是剔除哪些。典型相关里如果某一对典型变量的载荷显示葡萄的某个指标和酒的某个指标几乎一一对应(比如总糖对总糖),那这对变量贡献的是"同一件事的两种测法",保留两个只会让回归系数不稳定,标准误膨胀。
import statsmodels.api as sm X_sel = sm.add_constant(selected_df) # 筛选后的葡萄 + 酒理化指标 model = sm.OLS(y_score, X_sel).fit() print(model.summary()) # 看 R²、F、各系数 t 值与 VIF print(sm.stats.variance_inflation_factor(X_sel.values, 1)) # 逐个查 VIFsm.add_constant必须加,statsmodels 默认不包含截距项,漏了会导致 R² 虚高、系数偏移。VIF 超过 10 说明还有严重共线性,得回去继续剔指标。
论文的检验方式是对比样本组的模型评价值与真实评分值的显著性差异,结论是两者差异不显著,用理化指标评价葡萄酒质量可行。这个"可行"的边界要说清楚:它成立的前提是样本来自同一产区、同一年份、同一品种体系,换一个年份,回归系数就得重估。
5. 复现时最容易翻车的五处:数据、旋转、一致性与验证
5.1 数据层:缺失值、异常值和红白混检
现象:回填缺失值之后,第一组评酒员的整体均分比第二组高出零点几分,后面 Wilcoxon 检验的 z 值怎么调都对不上论文的 -2.53。
原因:用了整表均值填充。评酒员的打分尺度不同,有人习惯性给高分,整表均值会把这种个人偏移抹平,人为缩小或放大组间差异。
解决:一律按列(评酒员维度)算均值回填,缺失值所在的列是谁的就用谁自己的均值。异常值
本文还有配套的精品资源,点击获取