2012年数学建模A题葡萄酒评价：从配对检验到Lasso回归

简介2012年数学建模A题一等奖获奖论文系统研究葡萄酒质量评价问题适合数学建模参赛者、数据分析爱好者及对统计方法应用感兴趣的读者。论文完整解答四个子问题先用单样本K-S检验确定评酒员评分的概率分布再用Wilcoxon符号秩检验比较两组评酒员评分差异并通过秩相关分析与肯德尔和谐系数评价评分信度随后运用主成分分析与综合评价指标构建酿酒葡萄分级模型进一步利用典型相关分析揭示酿酒葡萄与葡萄酒理化指标间的联系最后建立多元线性回归模型论证理化指标评价葡萄酒质量的可行性。全文综合MATLAB、SPSS、SAS和EXCEL等软件实现包含问题重述、模型假设、符号说明、模型建立与求解、模型检验与推广等内容。资源为1份PDF文档约2.41MB目前已有2906人学习下载。从中可学习到秩相关分析、层次分析法、典型相关分析、多元回归等方法的综合运用是备赛数学建模与撰写高质量科技论文的实用参考。1. 2012年数学建模A题「葡萄酒的评价」一等奖论文都在解决同一个问题提到“2012年数学建模A题一等奖论文葡萄酒的评价”很多备赛的人第一反应是找一份PDF抄个模型。真按这个思路走大概率翻车。这道题表面是“品酒打分”实际考的是一整条统计决策链两组评酒员的分数能不能互相替代理化指标能不能解释感官评分解释到什么程度才算可信。当年的优秀论文里真正拉开差距的不是用了多复杂的算法而是把显著性检验、共线性处理、交叉验证这些基础工具用到了位并且每一步都留下了可核查的中间结果。这篇文章就按这个路径复现一遍从数据清洗、配对检验到回归建模和排错给出能直接运行的代码和参数取值同时把最容易出问题的几个细节单独列出来。适合正在备赛的队伍也适合以后要处理多人打分一致性、指标与主观评分关联性问题的从业者。2. 题目拆解与数据编排把两问拆成一张打分矩阵和一张理化指标表2.1 三张表其实是一条因果链先说题目给了什么。通常的版本是若干款红葡萄酒和白葡萄酒每组酒都有一张酿酒葡萄的理化指标表、一张成品葡萄酒的理化指标表另外还有两组评酒员对每一款酒的感官评分。两问也顺着这个结构来第一问比较两组评酒员的评分有没有系统性差异第二问用葡萄和酒的理化指标去解释评分差异。关键是把这三张表看成一个链条酿酒葡萄的理化指标决定了发酵后葡萄酒的理化成分而葡萄酒的理化成分又影响了品酒员的主观感受。一等奖论文通常不把三个对象混在一起直接回归而是先拆成“葡萄→酒”和“酒→评分”两个环节分别建模再串联解释。这样做的直接好处是当某个变量在“酒→评分”环节不显著时你能知道问题出在感官主观性上还是出在指标选择上而不是把所有噪声堆进一个黑匣子。我一般拿到题目第一件事是画一张对象关系图纸上画就行标注每张表的主键。这里的主键是酒样编号所有表都靠它关联。如果某张表里酒样编号有重复、缺失或格式不统一后面所有检验都会失真所以数据编排这一步值得花半小时做扎实。2.2 数据清洗长表、宽表与“评酒员编号”这个坑原始打分表往往是宽表一行是一瓶酒一列是一个评酒员的分数。这种格式给人看很直观但给统计用很别扭——pandas里的t检验、回归都偏好长表每一行是一个观测哪瓶酒、哪个评委、多少分。转换工具就是pivot / melt。import pandas as pd # 宽表wine_id 10列A组评酒员 10列B组评酒员 scores pd.read_csv(scores_wide.csv) long scores.melt( id_vars[wine_id], var_namejudge_code, value_namescore, ) # judge_code形如 A1, A2, ..., B10 long[group] long[judge_code].str.replace(r\d, , regexTrue) long[judge_id] long[judge_code].str.extract(r(\d))这里两个参数要解释一下。id_vars是保持不变的列也就是酒样编号var_name是原来列名落到的那一列起名judge_codevalue_name是分数落到的那一列。因为两组评酒员的编号都从1开始如果不加group列直接把A1和B1混在一起后面分组统计时就会出现“同一个评委id对应两个分数”的错乱——这就是评酒员id的坑。另一个高频问题是缺失值。打分表里某个评委漏评一款酒很常见。先看有没有空值再决定是删除该评委的全部数据还是用组内均值补。我倾向于删除这款酒在这一组内的数据理由很简单补一个均值进去虽然不改变组均值但会缩小组内方差让后续t检验的p值偏小等于人为制造“显著差异”。2.3 把两问翻译成统计语言数据理清之后两问就变成两个明确的统计任务。第一问同一批酒A组和B组各给了一个分数判断两组评分是否存在系统性偏差。因为是对同一款酒打的分本质是配对数据。很多人在这里就犯错——直接把两组所有分数当成两个独立样本做t检验忽略同款酒的内在关联导致标准误被低估。第二问一堆理化指标几十个里哪些能解释评分变化能解释多少预测稳定性如何。这里又分两层如果是酿酒葡萄的指标对酒评分直接建模中间隔着一道酿造转化如果是成品酒的指标则更接近直接因果。论文里一般两层都做优先看成品酒指标。下面这个表是每次给学生讲解时都会画的对应关系问题数据对象首选统计工具常见误用第一问同款酒在两组下的评分配对t检验 / Wilcoxon符号秩 / 效应量误用独立样本t检验第二问理化指标与评分相关分析、VIF、Lasso、多元回归几十个变量直接全回归然后宣称R方很高这两个问题的解法分别对应本文第3章和第4章。在做任何建模之前先把这两行对应关系确认下来能省掉后面一大半返工。3. 第一问两组评酒员的评分到底差没差配对检验与效应量3.1 为什么必须用配对设计而不是直接比平均分先讲选型逻辑。假设A组平均分比B组高0.3分能不能说“两组有差异”不能因为同一款酒本身有好坏好酒两组都给高分差酒两组都给低分如果忽略这个共同趋势把几百个评分当成独立样本组间差异的标准误会被算小p值也会被算小最后得到“显著差异”其实可能是假的。所以正确做法是先对每款酒在各自组内取平均得到每款酒的两个组平均分然后对这两个平均分的差值做配对检验。这一步把“酒与酒之间的差异”从“组与组之间的差异”里剥离开了检验的才是评酒员偏好。但配对t检验有一个前提差值近似正态。数据量小的时候通常只有几十款酒这个前提不成立就直接用t检验会出偏差所以要先做正态性检验不通过再退一步用非参的Wilcoxon符号秩检验。这也是为什么一等奖论文的检验部分通常同时报告两个方法的结果——不是炫技而是让别人能核查你的结论到底稳不稳。3.2 用scipy完成配对t检验、Wilcoxon与效应量计算直接给可运行的代码。假设long表已经按2.2整理好。import numpy as np import pandas as pd from scipy import stats # 每款酒在每个组下的平均分 avg long.pivot_table(indexwine_id, columnsgroup, valuesscore) diff avg[A] - avg[B] # 1) 差值正态性检验 # 适用场景n 50 时建议先看这个结果 shapiro_stat, shapiro_p stats.shapiro(diff) # 2) 配对t检验 t_stat, t_p stats.ttest_rel(avg[A], avg[B]) # 3) Wilcoxon符号秩检验不要求正态 w_stat, w_p stats.wilcoxon(avg[A], avg[B]) # 4) Cohens d效应量 mean_diff diff.mean() std_diff diff.std(ddof1) cohen_d mean_diff / std_diff if std_diff ! 0 else float(inf) print(fShapiro-Wilk p {shapiro_p:.3f}) print(f配对t检验: t {t_stat:.3f}, p {t_p:.3f}) print(fWilcoxon: W {w_stat:.3f}, p {w_p:.3f}) print(f差值均值 {mean_diff:.3f}, Cohens d {cohen_d:.3f})参数说明ttest_rel要求两组长度相同也就是酒样一一配对这正好是我们要的场景wilcoxon默认是检验差值的分布关于0对称如果数据里有差值为0的样本它会自动剔除并警告这不是bugddof1计算样本标准差配合小样本这个自由度更可靠。p值只能回答“有没有差异”回答不了“差异大不大”。这就是效应量的作用。Cohens d 差值均值 / 差值标准差0.2算小效应0.5中等0.8算大。比如p0.03但d0.15结论应该写成“统计上可检出差异但实际影响很小”而不是“两组评分显著不同”。3.3 效应量的置信区间与Bland-Altman一致性为了让结论更扎实我习惯再算一个效应量的95%置信区间。直接套公式在小样本下很敏感bootstrap是最稳的方案。rng np.random.default_rng(42) boot_d [] for _ in range(2000): idx rng.integers(0, len(diff), len(diff)) boot_d.append(diff[idx].mean() / diff[idx].std(ddof1)) lo, hi np.percentile(boot_d, [2.5, 97.5]) print(fCohens d 95% CI [{lo:.3f}, {hi:.3f}])default_rng(42)是固定随机种子保证任何人在任何机器上跑出来的结果一致。这是评审和复现时最看重的东西——一篇论文如果连随机种子都不写敏感数字每次跑都不一样很难让人信服。如果想看两组评分能不能互相替换再画一张Bland-Altman图横轴是两组对同一款酒平均分的均值纵轴是差值中间画均值线上下画1.96倍标准差的一致性界限。评分这一类主观数据一致性界限通常允许在±1分以内如果界限宽到±2分说明两组评酒员的偏好差异对单款酒来说相当大即使均值检验不显著也不建议直接把两组评分合并当同一个标签用。这个图用matplotlib几行就画完但关键不是图是那两条一致性界限的宽度怎么解读。4. 第二问用Lasso和逐步回归把理化指标变成有解释力的质量模型4.1 先拆因果链别把几十个指标一把梭第二问最常见的暴力做法拿几十个理化指标列加和当前评分列直接跑一个多元回归看到R方0.9就宣布“找到了评价模型”。这种模型交上去答辩环节大概率被问倒回归系数的符号反常识、变量之间高度相关导致系数不稳定、换几个样本就崩。正确打开方式是先拆链。葡萄酒的理化指标酒精度、单宁、总酚、花色苷等由葡萄原料决定感官评分又受酒的理化成分影响。所以先做“酒理化指标 → 感官评分”这一步因为它离评分最近、解释力最直接再做“葡萄理化指标 → 酒理化指标”作为补充。两个模型分开展示比混在一起更有说服力。另一个前置问题是所有回归之前指标先标准化。理化指标的数值范围差别极大pH是3点多花色苷可能是几百不标准化的话Lasso的惩罚项会对量纲大的变量不公回归系数也不可比较。这个操作要放在训练流程里不能拿到全量数据上先算好再切——否则就是数据泄漏后面第5章细说。4.2 用相关性矩阵和VIF先把共线性按下去几十个指标之间高度相关是必然的。直接上回归系数方差膨胀正负号可能都是错的。先做两件事相关系数矩阵粗筛VIF定量判断。import pandas as pd import numpy as np from sklearn.preprocessing import StandardScaler from statsmodels.stats.outliers_influence import variance_inflation_factor # chem_scores: 行为酒样列为理化指标最后一列是quality_score features chem_scores.drop(columns[quality_score]) y chem_scores[quality_score] scaler StandardScaler() X_scaled scaler.fit_transform(features) corr pd.DataFrame(X_scaled, columnsfeatures.columns).corr().round(2) print(corr.abs().unstack().sort_values(ascendingFalse).head(10)) vif pd.Series( [variance_inflation_factor(X_scaled, i) for i in range(X_scaled.shape[1])], indexfeatures.columns, ) print(vif.sort_values(ascendingFalse))variance_inflation_factor的计算逻辑是拿第i个变量对其它所有变量做回归看1/(1-R²)有多大。经验阈值是VIF10就说明这个变量几乎可以被其它变量线性表示留着只会捣乱。注意要先标准化再算VIF不然statsmodels给出的值会被量纲带偏。处理方式两种一是删除VIF高的变量保留有实际解释意义的那个二是干脆进入下一步用Lasso自动压变量。在食品、水质这类指标强相关场景里我习惯两者结合先删VIF10的再用Lasso从剩余变量里挑。4.3 Lasso选变量 OLS定系数 留一交叉验证Lasso的L1惩罚会让一部分系数变成精确的0相当于自动做变量选择特别适合“变量多、样本少”的比赛题。样本量通常就几十款酒五折交叉验证很勉强但够用变量更多时换留一法更稳。from sklearn.linear_model import LassoCV from sklearn.model_selection import LeaveOneOut, cross_val_score from sklearn.linear_model import LinearRegression from sklearn.pipeline import Pipeline import statsmodels.api as sm # alphas网格从1e-3到1对数等距取100个 lasso LassoCV(alphasnp.logspace(-3, 0, 100), cv5, max_iter20000, random_state0) lasso.fit(X_scaled, y) selected features.columns[lasso.coef_ ! 0] print(Lasso保留变量:, list(selected)) # 用保留变量建OLS模型拿p值和系数 X_sel scaler.fit_transform(features[selected]) ols sm.OLS(y, sm.add_constant(X_sel)).fit() print(ols.summary()) # 留一交叉验证看真实预测能力 # 注意标准化放进Pipeline避免数据泄漏 pipe Pipeline([ (scale, StandardScaler()), (lr, LinearRegression()), ]) cv_r2 cross_val_score(pipe, features[selected], y, cvLeaveOneOut(), scoringr2) print(fLOO R2 {cv_r2.mean():.3f} (/- {cv_r2.std():.3f}))几个参数值得说。alphas用np.logspace取对数网格是因为Lasso对惩罚系数是乘法尺度敏感等距网格会漏掉最优值。max_iter20000是防止在某些alpha下迭代不收敛报Warning实际数据里如果指标量纲特别大这个值还要加大。cv5对小数据集是折中最终结论以LOO为准因为样本太少时5折的划分方差很大。OLS是拿来定系数和显著性水平的Lasso选出来的变量再用OLS跑一遍输出的是标准误和p值这才方便写进论文的回归结果表。要注意这一步如果p值不显著也不能回头再手动塞变量进去凑数——那叫选择性报告比模型差更伤人品。4.4 残差诊断R方再高也骗不了这张图回归做完先不急着下结论。残差诊断至少看两样Q-Q图看残差正不正态残差与拟合值散点图看有没有“喇叭形”方差不齐。直接给检查代码import matplotlib.pyplot as plt from scipy import stats resid ols.resid fitted ols.fittedvalues fig, axes plt.subplots(1, 2, figsize(10, 4)) stats.probplot(resid, distnorm, plotaxes[0]) axes[1].scatter(fitted, resid, alpha0.7) axes[1].axhline(0, colorred, linestyle--) plt.show()如果Q-Q图两端明显偏离直线考虑对评分做Box-Cox变换或用稳健回归如果残差随拟合值增大而发散说明评分方差不是常数可以考虑加权最小二乘。这一步在论文里写一行“残差满足独立同分布假设”很容易但不做诊断就写这句话遇到较真的评委一眼看穿。5. 建模避坑手册显著性检验与回归里的5个翻车点5.1 翻车点一样本量不到30还硬上独立t检验现象两组评分各有百来个数据点直接ttest_ind算出p0.01结论“两组评分显著不同”结果评审指出数据不独立结论作废。原因所有评分共享同一批酒样的底层信息各评分之间不是独立的独立t检验的标准误会偏小。另外酒样数往往也就二十几款中心极限定理并不总能兜住均值差的正态近似。解决先按酒样聚合到组均值用配对t检验若差值经Shapiro-Wilk检验不是正态改用Wilcoxon符号秩。配对之后即使p值变大那也是真实情况——组间真实差异本来就不大。5.2 翻车点二只看p值不看效应量现象p0.04写进结论“有显著差异”但两组平均分只差0.2分评酒员自身对同一瓶酒重复打分都可能差出0.5分。原因p值度量的是“差异在统计上是否可检出”和“差异在业务上是否重要”不是一回事。样本量足够大时任何微小差异都能显著。解决同时报告Cohens d。d0.2时即使p显著结论也应写成“可检出但效应量很小实操上可视为同一水平”。对评酒数据通常把效应量阈值设在0.3以上再讨论合并评分。5.3 翻车点三变量几十个、样本几十行回归R方0.95交叉验证崩盘现象全变量多元回归训练集R方0.95留一交叉验证R方变成0.3甚至负数。原因指标之间强共线性模型把噪声也拟合进去了回归系数在样本外毫无泛化能力。评分本身带主观噪声能解释的真实方差本来有限。解决先看VIF删高共线变量再用LassoCV选变量最后用留一交叉验证报告泛化R方。训练R方和LOO R方一起列在论文里差距大说明模型过拟合要及时回头减变量。5.4 翻车点四标准化和特征筛选用了全量数据造成数据泄漏现象脚本里先对全量数据StandardScaler().fit_transform()再做train_test_split结果交叉验证分数虚高换成新数据一塌糊涂。原因标准化算出的均值和方差里混进了“测试集”的信息Lasso在特征筛选时也看到了全量数据的标签等于考试前看了答案。这在评酒这种小样本场景里尤其致命。解决把标准化放进Pipeline里让它和模型一起只对训练集fit。上一章4.3的pipe写法就是标准答案每一折的标准化参数都是从折内训练集学的验证分数才可信。5.5 翻车点五等级分当成连续分直接算皮尔逊相关现象评酒员打的是10、9.5、9、8.5这种半个刻度的离散分拿去做皮尔逊相关和回归相关系数总偏低回归系数也不稳定。原因皮尔逊相关假设两变量近似连续且线性相关数据被量化成少数几个档位后存在大量并列值相关性被低估这叫“离散化衰减”。解决涉及评分与其他连续指标的相关分析优先用Spearman秩相关建模时如果离散感强可以把评分秩次化再建模。至少在被问“为什么不用皮尔逊”时要能解释清楚离散分和连续分的差别。6. 把整套检验串成一条命令一份可复现的分析管线到了这一步代码功能都齐了但零散脚本没法交付。现在的习惯是开赛第一天就把骨架搭出来后面只填数据。管线分四段读表清洗、显著性检验、变量筛选回归、报告输出。# evaluate_wine.py 骨架 import pandas as pd from scipy import stats import numpy as np def load_and_clean(scores_path): df pd.read_csv(scores_path) long df.melt(id_vars[wine_id], var_namejudge_code, value_namescore) long[group] long[judge_code].str.replace(r\d, , regexTrue) return long.dropna() def pairwise_test(long): avg long.pivot_table(indexwine_id, columnsgroup, valuesscore) diff avg.iloc[:, 0] - avg.iloc[:, 1] return { shapiro_p: stats.shapiro(diff).pvalue, ttest_p: stats.ttest_rel(avg.iloc[:, 0], avg.iloc[:, 1]).pvalue, wilcoxon_p: stats.wilcoxon(diff).pvalue, mean_diff: diff.mean(), } def regression_report(chem, score_colquality_score): # 内部完成: 标准化 - LassoCV选变量 - OLS - LOO R2 ...骨架的价值不是省那几行代码而是把“检验结论”和“建模参数”固定下来形成一份每次比赛都能复用的模板。随机种子写死数据清洗版本号写进输出文件名中间结果存CSV这样复盘时能清楚知道哪一步改了什么、结论为什么变化。我踩过的教训是第一次参加这类题两天时间全耗在临时改脚本和记不清哪版结果上最后交的图自己都复现不出来。从那以后坚持先写骨架再分析流程固定后真正留给思考模型和写论文的时间反而多了。比赛和实际项目拼的都不是玄学是把每一步做到可以复现、每个结论做到可以追查。这套路径能帮你少走一段弯路也希望你最后交出去的论文里每一个p值、每一个回归系数都经得起别人照着数据重跑一遍的检验。希望帮到你。本文还有配套的精品资源点击获取