葡萄酒质量评价建模:主成分分析与聚类回归实战解析 简介2012年数学建模A题一等奖论文完整呈现“葡萄酒的评价”赛题建模与求解全过程。资源包为单个PDF文档共1个文件大小约2.41MB可清晰阅读论文图表、公式与完整正文。论文围绕葡萄酒质量评价综合运用单样本K-S检验、符号秩检验、秩相关分析、主成分分析、典型相关分析与多元线性回归等方法依次解决两组评酒员评分差异检验、酿酒葡萄分级、葡萄与葡萄酒理化指标关联分析以及用理化指标评价葡萄酒质量的可行性论证等问题建模中借助多种统计软件完成计算并涵盖问题重述、模型假设、符号说明、模型建立与求解等完整环节。已有2906人学习下载适合数学建模参赛学生、备战国赛美赛的选手以及统计建模与综合评价方法学习者参考可从中学习问题拆解、模型选择、软件实现、论文撰写等完整思路。1. 葡萄酒评价这道题为什么值得琢磨从竞赛A题到实用分析模型做数据分析的人第一次看到2012年数学建模A题通常会被两样东西吓住一是酿酒葡萄和成品葡萄酒各自那几十项理化指标二是评酒员给出的感官评分表。这道题表面上问“葡萄酒的质量怎么评价”实际考的是两件事把理化成分和专家打分联系起来并让评价结果能解释、能复现。当年能拿到一等奖的论文共性往往不是用了多高深的算法而是把样品对应关系、量纲处理、降维和聚类这几步梳理得足够干净。这篇文章不还原任何一篇获奖原文只讲这个题目从数据清洗到模型落地的可靠路径以及参数怎么调、坑在哪里。适合参加类似建模比赛的学生也适合在质检、食品饮料行业做感官数据分析的从业者。2. 先看懂题数据里藏着哪几类信息以及一等奖论文的建模主线2.1 数据结构与评价目标理化指标、感官评分、葡萄与酒的对应关系这道题提供给参赛者的通常是三张表酿酒葡萄的理化指标表、成品葡萄酒的理化指标表、评酒员对葡萄酒的感官评分表。葡萄理化指标往往包含几十列比如多种氨基酸、糖酸比、矿物质含量、花色苷等葡萄酒理化指标则涉及酒精度、总酚、单宁、色泽等。感官评分表更麻烦通常会有两组评酒员每个人对每个样品的外观、香气、口感和整体评价打分有时还是多次重复测量。所以在动笔写模型之前第一步不是选算法而是把这三张表的关系理清。每行是一个样品编号葡萄样品和葡萄酒样品一一对应感官评分表里的样品编号也必须能和前两张表对上。只有把“哪种葡萄酿出哪种酒、酒得到多少分”这三者绑在一起后面所有分析才有意义。很多翻车现场都是因为复制粘贴时丢了一个编号或者合并表时顺序错位后面跑出的模型再漂亮也是错的。目标可以拆成两层一是给葡萄酒质量打一个可量化的分二是解释哪些理化指标对这个分影响最大。如果只看感官评分那就用聚类或判别分析把样品分成优、良、差几级如果想把葡萄品质和酒品质打通就需要回归建模用葡萄理化指标预测酒的口感得分。大多数一等奖论文不会只用一种方法而是先用主成分分析降维再聚类分级最后用回归或判别定位关键指标形成一个闭环。2.2 建模主线的选择理由为什么降维、聚类和回归要搭配着用直接拿几十列理化指标去预测感官评分很容易过拟合因为样本量通常只有几十个而指标可能有四十甚至更多。指标之间还存在多重共线性氨基酸总量和多种氨基酸组分高度相关花色苷和色泽指标也相关。这种情况下最小二乘回归的系数极不稳定训练集表现很好交叉验证一测就垮。所以常见做法是先做主成分分析用几个综合变量代替原来的几十个变量把相关性和冗余去掉。聚类在这里的作用是解决“等级怎么分”的问题。感官评分虽然给出了分数但分数本身受评酒员尺度影响相同质量的酒在不同评酒员手里可能差出三五分。与其直接用原始分当回归目标不如先对理化因子做聚类再用每个聚类的平均感官分去判断哪个等级更高。这样既能规避评分噪声又能让模型输出“优、良、差”这类可解读的结果。回归的作用则是把等级映射回连续分数或者找到影响评分的关键因子。可以用主成分得分作为自变量感官平均分作为因变量建立线性回归也可以用Lasso回归直接在原始指标上做变量选择。两种路径各有各的用途PCA回归适合解释主成分贡献Lasso适合指出具体是哪几项理化指标在起作用。下表概括了常用方法在本题中的定位方便按自己的数据情况选。方法解决什么问题在本题中的输出边界主成分分析降维、消除共线性主成分得分、累计方差贡献率载荷解释可能不直观K/层次聚类确定质量等级样品的等级标签聚类数量需要结合评分校准多元线性回归建立连续分数预测模型预测得分、关键变量系数样本量小容易过拟合Lasso回归变量选择稀疏系数筛出重要指标需要调节正则强度判别分析分类及变量重要性分类准确率、判别函数需要已知等级标签这套组合在历年竞赛里都很常见因为每一步都有明确的现实意义降维解决“指标太多”聚类解决“等级模糊”回归解决“哪项指标说了算”。后面几章我就按这条主线展开先做数据清洗再跑模型最后聊踩过的坑。3. 把数据洗干净三表对齐、缺失值和量纲统一的具体做法3.1 第一步把葡萄理化、葡萄酒理化和感官评分三张表对齐很多人在这一步就翻车。三张表来自不同文件行顺序未必一致甚至样品编号格式都不同有的叫“红葡萄1”有的叫“酒样1”。如果直接按行号合并等于给自己埋雷。我通常先把每张表的样品编号列统一成同一格式然后做内连接确保只保留三张表都有的样品。import pandas as pd # 读取三张表 df_grape pd.read_excel(grape_physicochemical.xlsx) df_wine pd.read_excel(wine_physicochemical.xlsx) df_score pd.read_excel(sensory_score.xlsx) # 查看编号列的数据类型和唯一性 print(df_grape[样品号].dtype, df_grape[样品号].nunique()) print(df_wine[样品号].dtype, df_wine[样品号].nunique()) print(df_score[样品号].dtype, df_score[样品号].nunique())如果发现格式不一致比如葡萄表编号是“样1”酒表编号是“Sample1”要先把它们改成同一个规范。常见做法是写一个映射函数统一成数字编号或者干脆用“样1”这种短字符串。合并不需要用全部变量只要保留后面的建模列即可。先做一次内连接看看还剩多少行如果比最小表行数还少说明编号有错位。# 统一编号为字符串并去掉多余空格 df_grape[样品号] df_grape[样品号].astype(str).str.strip() df_wine[样品号] df_wine[样品号].astype(str).str.strip() df_score[样品号] df_score[样品号].astype(str).str.strip() # 三表做内连接 df_all df_grape.merge(df_wine, on样品号, suffixes(_葡萄, _葡萄酒)) df_all df_all.merge(df_score, on样品号) print(df_all.shape)这一步的逻辑很直白先确保样品维度对齐再谈数据清洗。合并后要检查是否有重复列比如两边都叫“总酚”后面会自动加上后缀。我一般会检查一遍df_all.columns把不需要的编号列和重复列删掉。需要特别注意的是感官评分表如果是长格式也就是每行是一个评酒员对某个样品的评分要先聚合或转成宽格式否则合并后会变成每个样品对应多行把模型训练集搞乱。3.2 缺失值和异常值能中位数填充的就不要随手删行理化指标表里经常有少量缺失值原因包括检测设备故障、某批次样品量不足等。如果一个指标的缺失率超过30%我倾向于直接去掉这一列因为后续不管用什么插补方法这一列都不太可能是稳定信号。如果缺失率不高用中位数填充比均值填充更稳尤其当指标分布偏斜时均值会把异常值带进数据。# 统计缺失率剔除缺失率高的列 miss_rate df_all.isnull().mean() cols_high_miss miss_rate[miss_rate 0.3].index.tolist() df_all.drop(columnscols_high_miss, inplaceTrue) # 数值列用中位数填充 numeric_cols df_all.select_dtypes(includenumber).columns df_all[numeric_cols] df_all[numeric_cols].fillna(df_all[numeric_cols].median())填充之后再检查异常值。食材类理化数据常常有极端值比如某一批葡萄的花色苷含量是其他样品的五六倍。用IQR四分位距筛一遍会很有用但注意不要直接把整个样本删掉因为竞赛数据只有几十个样品删一个可能就少一个关键类别。我一般先把超过上下限3倍IQR的值标记出来再回到原始记录里看是不是录入错误。如果是真实数值就保留因为异常值本身可能代表一个特殊等级。# 基于IQR标记异常值 def mark_outliers(s): q1, q3 s.quantile(0.25), s.quantile(0.75) iqr q3 - q1 lower_bound q1 - 3.0 * iqr upper_bound q3 3.0 * iqr return (s lower_bound) | (s upper_bound) outlier_flags df_all[numeric_cols].apply(mark_outliers) print(outlier_flags.sum())这段代码的作用是生成一个与数据框同形状的布尔矩阵方便定位哪一行哪一列超界。参数3.0是经验值竞赛数据通常可以放宽到2.5或3.0因为样本量太小IQR本身就受异常值影响。不建议把超过阈值的行全部删除更好的做法是看具体数值是否合理或者用Winsorization把极端值缩到合理范围内。3.3 标准化与归一化做PCA之前必须统一量纲主成分分析对变量尺度极其敏感。如果直接用原始数据酒精度的单位是“%vol”氨基酸的单位是“mg/L”两者数值差异巨大PCA会不自觉把数值大的变量当成主成分的主导完全失去解释意义。因此标准化是必做的一步。from sklearn.preprocessing import StandardScaler # 去掉样品号之后选取数值列 feature_cols numeric_cols.drop(样品号, errorsignore) X df_all[feature_cols] scaler StandardScaler() X_scaled scaler.fit_transform(X)这里我选择StandardScaler而不是MinMaxScaler理由是想保留原始分布的形状让每个变量的均值归零、方差为1。MinMax归一化会把所有值压缩到0到1之间如果分布有极端值大多数点会挤在很窄的区间里反而不利于后续PCA。需要注意的是StandardScaler要在分成训练集和测试集之前对全部数据拟合还是先拟合训练集再变换测试集这里有个常见误区如果后面要做交叉验证最好在每次CV fold内部单独做标准化而不是用全局均值和标准差。不过在竞赛场景里样本量小通常先全量标准化再建模但要在心里清楚这个操作会轻微泄漏测试集信息。我一般会用管道Pipeline来处理避免自己手动做错。4. 建立葡萄酒质量评价模型主成分降维、聚类分级与回归预测4.1 主成分分析累计方差贡献率选到85%还是90%主成分个数的选择没有唯一答案但在这道题里常见做法是保留累计方差贡献率85%左右。因为理化指标之间相关性强通常前四五个主成分就够解释大部分信息。取85%而不是90%目的是让模型更简洁减少噪声。当然如果前几个主成分在业务上很难解释可以适度增加主成分个数但不要为了凑累计贡献率而留到十几个。from sklearn.decomposition import PCA # 先跑一次完整PCA查看累计方差贡献率 pca_full PCA() pca_full.fit(X_scaled) cum_ratio pca_full.explained_variance_ratio_.cumsum() print(cum_ratio) # 找到累计贡献率首次超过0.85的主成分个数 n_components (cum_ratio 0.85).argmax() 1 print(n_components , n_components)解释一下argmax的用法布尔数组里首个True的索引加1就是所需个数。如果累计贡献率曲线一直很平缓说明数据本身没有特别强的结构这时候可以考虑用碎石图拐点来选或者直接固定3到5个主成分。主成分数量不要超过样本量的五分之一否则后面回归容易过拟合。比如只有30个样品主成分选到10个以上就不太合适了。选定个数后重新拟合PCA并获取主成分得分pca PCA(n_componentsn_components) X_pca pca.fit_transform(X_scaled) # 查看每个主成分的方差占比 print(pca.explained_variance_ratio_) # 打印前两个主成分的载荷矩阵行名 print(pca.components_)主成分得分矩阵X_pca每一行对应一个样品每一列是一个主成分。后面聚类和回归都用这个得分矩阵而不是原始指标。载荷矩阵pca.components_的维度是主成分数×原始特征数每一行对应一个主成分值的大小和正负表示原始指标在该主成分上的贡献方向和强度。4.2 聚类分级KMeans还是层次聚类轮廓系数选K拿到主成分得分后可以聚类来给样品分级。KMeans简单稳定适合样本量不大的情况层次聚类能画树状图但在竞赛论文里反而不如KMeans直观。这里用KMeans重点是如何决定分成几类。感官评分表一般会暗示质量等级可能是优、良、差三类但理化数据不一定能分成三类所以先用轮廓系数挑K再结合感官评分校准。from sklearn.cluster import KMeans from sklearn.metrics import silhouette_sample_score, silhouette_score best_k 2 best_sil -1 sil_scores [] for k in range(2, 7): km KMeans(n_clustersk, random_state42, n_init10).fit(X_pca) sil silhouette_score(X_pca, km.labels_) sil_scores.append(sil) if sil best_sil: best_sil sil best_k k print(best_k , best_k, best_sil , best_sil)range(2, 7)意味着从2类试到6类竞赛数据一般不会超过6类。random_state42是为了固定初始质心让结果可复现后面分析不会因为随机数不同而改变。n_init10是让KMeans跑10次不同的初始质心选最小SSE的一次能减少陷入局部最优的概率。确定K后给每个样品打上等级标签km KMeans(n_clustersbest_k, random_state42, n_init10) cluster_labels km.fit_predict(X_pca) df_all[等级标签] cluster_labels注意这里聚类得到的类别编号是无序的编号0不一定对应质量最差的等级。要把聚类结果和感官评分对应起来常见的做法是算每个类别的感官平均分然后按平均分从低到高映射到“差、中、优”等级。# 假设感官评分列是整体评分 group_score df_all.groupby(等级标签)[整体评分].mean().sort_values() rank_map {old_label: new_label for new_label, old_label in enumerate(group_score.index)} df_all[质量等级] df_all[等级标签].map(rank_map)后面所有分析都建议用质量等级而不是等级标签因为前者是有顺序的后者是任意编号。4.3 建立回归预测模型从主成分得分到感官评分等级是离散的如果题目要求给每个样品一个具体分数还需要回归模型。最简单可靠的做法是将主成分得分作为自变量将评酒员打分的均值作为因变量做多元线性回归。因为样本量有限我建议用留一交叉验证而不是普通的5折交叉验证。from sklearn.linear_model import LinearRegression from sklearn.model_selection import LeaveOneOut, cross_val_score # 因变量感官评分的平均值假设已经算好 y df_all[整体评分] # 主成分得分作为自变量 reg LinearRegression() loo LeaveOneOut() cv_scores cross_val_score(reg, X_pca, y, cvloo, scoringr2) print(LOO-CV R2 , cv_scores.mean(), /-, cv_scores.std())留一交叉验证适合几十个小样本场景每轮只留一个样品做测试最大程度利用训练数据。但它的方差大所以要看均值和方差不要单看某一个折的结果。如果R2均值不理想常见改进是改用Lasso回归在原始指标上做变量选择或者用主成分得分配合随机森林不过随机森林在小样本上很容易过拟合需要谨慎。from sklearn.linear_model import Lasso from sklearn.pipeline import make_pipeline from sklearn.preprocessing import StandardScaler # 如果直接用原始指标pipeline自动做标准化和Lasso lasso_pipe make_pipeline(StandardScaler(), Lasso(alpha0.05, max_iter10000)) lasso_cv_scores cross_val_score(lasso_pipe, X, y, cvloo, scoringr2) print(Lasso LOO-CV R2 , lasso_cv_scores.mean())alpha0.05是正则强度需要调参。如果用交叉验证来选alpha比如用LassoCV会得到更稳妥的结果。但要注意LassoCV在内部会再做一次交叉验证最外层再套LOO计算量不大样本量小也能接受。最后回归模型训练完成后要看系数标准化后的系数绝对值越大说明该指标对感官评分的影响越大。这就是论文里“找关键指标”的那张表。5. 避坑与常见问题数据对应、评分尺度、载荷解释和过拟合5.1 样品编号错位导致模型结论全错现象模型跑出来很漂亮主成分贡献率、聚类轮廓系数都不错但回归系数的符号和常识相反比如酒精度越高评分越低。回头检查时发现葡萄表、酒表和评分表行顺序不一致合并后样品编号没有对齐。原因合并时用了concat或者只按行索引合并没有用样品号做key。复制粘贴时也可能把某一行删掉导致后续所有行错位。解决在建任何模型之前先检查三张表的样品号列表是否一致用set差集找丢失或多余的编号。合并时只用merge不用concat保证以样品号为主键。每一步清洗后保存一次带编号的中间文件方便回溯。5.2 感官评分尺度差异大直接平均不可靠现象同一款葡萄酒评酒员A打了85分评酒员B只打了72分两组的整体均值全被这两人的尺度左右。聚类结果和感官评分对不上某类样品平均分反而比另一类低。原因不同评酒员对好坏的判断标准不同有人天生手松有人手紧。直接对原始分取平均相当于把个人偏差当成了质量差异。解决先对评分做标准化。常用做法是对每个评酒员的打分单独做z-score再按样品取平均。还有一种做法是只比较同一评酒员对不同样品的相对排序不依赖绝对值。处理后的评分再作为回归目标模型会更稳定。如果评分表有多个评酒员建议先画一张箱线图看看是否存在明显异常分。5.3 主成分载荷方向解释不清楚现象第一个主成分的线性组合里一半指标系数为正一半为负说不清PC1代表什么东西。论文里硬写“综合品质因子”但读者看了载荷表也不知道在说什么。原因主成分是最大方差方向的线性组合没有考虑业务含义。当指标数量多且相关性复杂时载荷会分散在许多变量上解释力下降。解决可以用Varimax旋转把载荷向0或1两端靠更容易解释也可以不强行解释主成分而是直接用它做聚类和回归只把主成分得分作为中间变量使用。论文里把主成分解释放在次要位置重点展示聚类结果和回归预测能力反而更实在。5.4 训练集R2接近1交叉验证R2却是负数现象线性回归在训练集上R2高达0.99但LOO交叉验证R2只有0.2甚至为负预测值和真实值几乎没关系。原因指标数远多于样本数模型学到了噪声和个别样本的独特模式。主成分虽然降到了几个维度但若主成分数依然太多比如样本量30选了15个主成分仍会过拟合。解决降低主成分个数比如从85%改为70%或者固定为3。还可以改用Lasso回归或岭回归用正则化惩罚大系数。最直接的办法是把样本量与特征数的比例控制在10比1以上即30个样本最多用3个主成分。这样子虽然训练集R2会下降但交叉验证分数更真实。5.5 聚类结果不稳定每次跑出来等级标签都变现象同一份数据换台电脑跑KMeans聚类中心和标签编号对不上画图时发现类别顺序完全变了论文里没法复现。原因KMeans初始质心是随机选择的没有固定随机种子。标签编号本身也没有业务含义每次运行编号会变。解决固定random_state比如random_state42。同时将聚类结果映射到感官评分后使用有序的质量等级列不要直接用原始标签列。写代码时把random_state和聚类数量写进脚本注释方便自己或合作者复现。6. 让模型结果更可信稳定性验证与可视化技巧6.1 用灵敏度分析检查主成分个数对结果的影响主成分个数从2调到6聚类结果和回归R2会跟着变。如果某个结论只在某个固定个数下成立那它就是脆弱的如果结论在3到5个主成分范围内都稳定说服力就强得多。我在建模时会把主成分个数作为一个循环变量跑一遍记录每个个数下交叉验证R2的变化。results [] for nc in range(2, 7): pca PCA(n_componentsnc).fit(X_scaled) X_pc pca.transform(X_scaled) reg LinearRegression() loo LeaveOneOut() r2 cross_val_score(reg, X_pc, y, cvloo, scoringr2).mean() results.append((nc, r2)) print(results)这段代码的作用不是找最优而是看趋势。如果R2从0.52跳到0.55又跳回0.51说明模型对主成分数不敏感结论可信如果R2从0.2一路爬到0.8那说明你选的主成分个数对结果影响太大要谨慎解释。6.2 用碎石图和载荷热力图快速定位关键因子的可视化技巧碎石图能直观展示每个主成分的方差贡献适合放在论文里支撑主成分个数选择。载荷热力图则能展示哪些指标在哪些主成分上权重更大。注意不要把载荷矩阵直接画成热力图要先对载荷取绝对值或者做归一化否则颜色全被最大载荷主导。import matplotlib.pyplot as plt import seaborn as sns # 碎石图 plt.figure(figsize(8, 4)) plt.plot(range(1, len(cum_ratio) 1), cum_ratio, o-) plt.xlabel(主成分个数) plt.ylabel(累计方差贡献率) plt.axhline(0.85, colorred, linestyle--) plt.show() # 载荷热力图 loadings pd.DataFrame( pca.components_.T, indexfeature_cols, columns[fPC{i1} for i in range(pca.n_components_)] ) plt.figure(figsize(10, 8)) sns.heatmap(loadings, cmapcoolwarm, center0, annotTrue) plt.show()热力图里annotTrue会显示每个格子的数值如果指标太多建议关掉annot只保留颜色。画图时还要注意如果主成分数为0会报错所以要确保pca已经选定有效的n_components。这类图在报告里比一大段文字更有说服力。6.3 一个值得养成的工作习惯保存中间过程给后悔药留个出口我做这类小样本建模时最怕的不是模型效果差而是隔一天回来想改某个参数发现原始处理过后的数据已经覆盖无法回头。后来我养成了一个习惯每一步清洗和转换后都保存一版中间数据文件名带处理阶段比如grape_scaled.csv、X_pca.csv。代码里固定所有随机种子并且把关键的参数写在脚本开头的配置变量里。这样即使实验结果不理想也不用全部重跑。有一次我因为没有固定KMeans的随机种子第二次运行时聚类标签全部变了图也乱了整个早上都在查为什么结果对不上。后来发现只是忘了写random_state。从那时起凡是涉及随机数的算法我都会在一开始就固定种子并将这个习惯写进模板。做这类评价建模复现比精度更重要一个能稳定复现的模型哪怕R2低一些也比一个运气好的黑匣子更有实用价值。希望这些细节能帮你在自己的数据和题目上少走几段弯路。本文还有配套的精品资源点击获取