一.题目
确定葡萄酒质量时一般是通过聘请一批有资质的评酒员进行品评。每个评酒员在对葡萄酒进行品尝后对其分类指标打分，然后求和得到其总分，从而确定葡萄酒的质量。酿酒葡萄的好坏与所酿葡萄酒的质量有直接的关系，葡萄酒和酿酒葡萄检测的理化指标会在一定程度上反映葡萄酒和葡萄的质量。附件1给出了某一年份一些葡萄酒的评价结果，附件2和附件3分别给出了该年份这些葡萄酒的和酿酒葡萄的成分数据。请尝试建立数学模型讨论下列问题:
①分析附件1中两组评酒员的评价结果有无显著性差异，哪一组结果更可信？
②根据酿酒葡萄的理化指标和葡萄酒的质量对这些酿酒葡萄进行分级。
③分析酿酒葡萄与葡萄酒的理化指标之间的联系。
④分析酿酒葡萄和葡萄酒的理化指标对葡萄酒质量的影响，并论证能否用葡萄和葡萄酒的理化指标来评价葡萄酒的质量？
附件1：葡萄酒品尝评分表（含4个表格）
附件2：葡萄和葡萄酒的理化指标（含2个表格）
附件3：葡萄和葡萄酒的芳香物质（含4个表格）
二、问题①求解
(1)分析：问题①想确定两组品酒员的评价结果有无显著性差异，再评判哪组结果更可信。可使用统计学中的显著性检验方法。显著性检验(test of significance)的方法有t检验、F检验和卡方检验等，该问题适合用t检验。t检验分为单总体检验和双总体检验，单总体t检验：检验一个样本平均数与一个已知的总体平均数的差异是否显著。双总体检验：检验两个样本平均数与其各自所代表的总体的差异是否显著。该问题有两个样本，所以可采用双总体t检验。
T检验、F检验、卡方检验详细分析及应用场景总结
(2)方法确定后，再来确定研究对象。对每个样品，不同组酒员评价的所有指标之和应用t检验最合适。
(3)原始数据预处理：①数据质量检查与清洗，即查看数据是否有缺失，如果有缺失，则需要进行填充。②对excel表格的数据进行整理，重排，以便利用程序进行比较。
把数据都凑一块，命名为2012A_T1_processed.xls
(4)程序设计
clc
clear all
%导入数据
X1=readtable('2012A_T1_processed.xls','Sheet','T1_red_grape','Range','D3:M272');
X2=readtable('2012A_T1_processed.xls','Sheet','T2_red_grape','Range','D3:M272');
X3=readtable('2012A_T1_processed.xls','Sheet','T1_white_grape','Range','D3:M282');
X4=readtable('2012A_T1_processed.xls','Sheet','T2_white_grape','Range','D3:M282');
% 红葡萄酒T检验计算过程
[m1,n1]=size(X1); % 返回行数和列数
K1=27; % 有27个样品
% 计算每个样品的总得分
for i=1:K1
for j=1:n1
SX1(i,j)=sum(table2array(X1(10*i-9:10*i,j)));
SX2(i,j)=sum(table2array(X2(10*i-9:10*i,j)));
end
end
% 计算每组样品得分的均值
for i=1:K1
Mean1(i)=mean(SX1(i,:));
Mean2(i)=mean(SX2(i,:));
end
% 计算检验值
for i=1:K1
S1(1,i)=(sum((SX1(i,:)-Mean1(i)).^2)+sum((SX2(i,:)-Mean2(i)).^2))/(n1*(n1-1));
T1(1,i)=(Mean1(i)-Mean2(i))/(sqrt(S1(1,i)));
end
AT_R=abs(T1);
M_AT_R=mean(AT_R);
% 白葡萄酒T检验计算过程
[m2,n2]=size(X3);
K2=28;
% 计算每个样品的总得分
for i=1:K2
for j=1:n2
SX3(i,j)=sum(table2array(X3(10*i-9:10*i,j)));
SX4(i,j)=sum(table2array(X4(10*i-9:10*i,j)));
end
end
% 计算每组样品得分的均值
for i=1:K2
Mean3(i)=mean(SX3(i,:));
Mean4(i)=mean(SX4(i,:));
end
% 计算检验值
for i=1:K2
S2(1,i)=(sum((SX3(i,:)-Mean3(i)).^2)+sum((SX4(i,:)-Mean4(i)).^2))/(n2*(n2-1));
T2(1,i)=(Mean3(i)-Mean4(i))/(sqrt(S2(1,i)));
end
AT_W=abs(T2);
M_AT_W=mean(AT_W);
% 结果显示与比较
a=2.102; % T(0.05,2,18)=2.101
b=2.878; % T(0.01,2,18)=2.878
set(gca,'linewidth',2)
% 红酒结果
for i=1:K1
Ta1(i)=a;
Tb1(i)=b;
end
t1=1:K1;
subplot(2,1,1);
plot(t1,AT_R,'*k-',t1,Ta1,'r-',t1,Tb1,'-.b', 'LineWidth', 2)
title('红酒显著性检验结果','fontsize',14)
legend('T检验值', 'T(0.05)值', 'T(0.01)值')
xlabel('样品号'), ylabel('T检验值')
% 白酒结果
for i=1:K2
Ta2(i)=a;
Tb2(i)=b;
end
t2=1:K2;
subplot(2,1,2);
plot(t2,AT_W,'*k-',t2,Ta2,'r-',t2,Tb2,'-.b', 'LineWidth', 2)
title('白酒显著性检验结果','fontsize',14)
legend('T检验值', 'T(0.05)值', 'T(0.01)值')
xlabel('样品号'), ylabel('T检验值')
% 显示平均检验结果
disp(['两组品酒师对红酒的平均显著性T检验值:' num2str(M_AT_R)]);
disp(['两组品酒师对白酒的平均显著性T检验值:' num2str(M_AT_W)]);
输出结果
两组品酒师对红酒的平均显著性T检验值:1.7539
两组品酒师对白酒的平均显著性T检验值:1.1641
查表可知，t(0.05,2,18)=2.101，t(0.01,2,18)=2.878。按照t检验的步骤③，可知t<f(df)0.05(df=(10-1)+(10-1)=18)，所以可得出：
①两组品酒师对红葡萄酒和白葡萄酒的评价结果都不显著。
②对白葡萄酒评价结果的差异小于对红葡萄酒的差异。
t检验临界值表
(5)评价结果稳定性
对每一个样品而言，可以用每组品酒师的评分(对样品的各项指标评分之和)对总体样本的方差来表示各组品酒师评价结果的稳定性。
程序设计：可分别求出两组品酒师对各个红葡萄酒和白葡萄酒样品评价的方差，同时得到总方差
clc
clear all
%导入数据(白葡萄酒)
% X1=readtable('2012A_T1_processed.xls','Sheet','T1_white_grape','Range','D3:M282');
% X2=readtable('2012A_T1_processed.xls','Sheet','T2_white_grape','Range','D3:M282');
%导入数据(红葡萄酒)
X1=readtable('2012A_T1_processed.xls','Sheet','T1_red_grape','Range','D3:M272');
X2=readtable('2012A_T1_processed.xls','Sheet','T2_red_grape','Range','D3:M272');
% 计算每组品酒师对每个样品的方差
[m,n]=size(X1);
K=27; %红葡萄酒
% K=28; %白葡萄酒
% 计算每个样品的总得分
for i=1:K
for j=1:n
SX1(i,j)=sum(table2array(X1(10*i-9:10*i,j)));
SX2(i,j)=sum(table2array(X2(10*i-9:10*i,j)));
end
u0(i)=mean([SX1(i,:), SX2(i,:)]);
end
% 计算方差
for i=1:K
SD1(i,:)=(SX1(i,:)-u0(i)).*(SX1(i,:)-u0(i));
SD2(i,:)=(SX2(i,:)-u0(i)).*(SX2(i,:)-u0(i));
end
% 结果显示与比较
for i=1:K
TSD(1,i)=sum(SD1(i,:));
TSD(2,i)=sum(SD2(i,:));
end
t=1:K;
plot(t,TSD(1,:),'*k-',t,TSD(2,:),'ok--', 'LineWidth', 2)
set(gca,'linewidth',2);
legend('一组方差','二组方差')
xlabel('红葡萄酒样品编号'); ylabel('红葡萄酒评价方差');
% xlabel('白葡萄酒样品编号'); ylabel('白葡萄酒评价方差');
TSD1=sum(TSD(1,:));
TSD2=sum(TSD(2,:));
disp(['一组对红葡萄酒总方差:' num2str(TSD1)]);
disp(['二组对红葡萄酒总方差:' num2str(TSD2)]);
% disp(['一组对白葡萄酒总方差:' num2str(TSD1)]);
% disp(['二组对白葡萄酒总方差:' num2str(TSD2)]);
输出1
一组对红葡萄酒总方差:16434.675
二组对红葡萄酒总方差:10524.775
输出2
一组对白葡萄酒总方差:34961.6
二组对白葡萄酒总方差:16522.8
结论：由于方差越小，说明越稳定。综上可知，第二组对两类酒评价的稳定性都比第一组高。
三、问题②求解
(1)分析：要求对酿酒葡萄进行分级，但是并没有给出分级的标准和具体的分级数，所以该问题属于聚类问题。酿酒葡萄的编号与葡萄酒的编号是一致的，存在严格的对应关系，所以可以根据葡萄酒的理化指标进行聚类，也可以根据葡萄酒质量(评分)进行聚类。
(2)根据葡萄酒质量(评分)对酿酒葡萄进行分级
由问题①可知虽然第二组品酒师的数据更稳定，但两组品酒师的结果无显著差异，所以以两组品酒师的平均值，作为酒样的质量数据。
首先确定类别数量
clc
clear all
%导入数据
X1=readtable('2012A_T1_processed.xls','Sheet','T1_red_grape','Range','D3:M272');
X2=readtable('2012A_T1_processed.xls','Sheet','T2_red_grape','Range','D3:M272');
X3=readtable('2012A_T1_processed.xls','Sheet','T1_white_grape','Range','D3:M282');
X4=readtable('2012A_T1_processed.xls','Sheet','T2_white_grape','Range','D3:M282');
[m1,n1]=size(X1); % 返回行数和列数
K1=27; % 有27个样品
% 计算每个样品的总得分
for i=1:K1
for j=1:n1
SX1(i,j)=sum(table2array(X1(10*i-9:10*i,j)));
SX1(i,j)=(SX1(i,j)+sum(table2array(X2(10*i-9:10*i,j))))/2;
end
end
[m2,n2]=size(X3);
K2=28;
for i=1:K2
for j=1:n2
SX2(i,j)=sum(table2array(X3(10*i-9:10*i,j)));
SX2(i,j)=(SX2(i,j)+sum(table2array(X4(10*i-9:10*i,j))))/2;
end
end
% 计算红葡萄酒和白葡萄酒质量评分数据
for i=1:K1
A(i)=mean(SX1(i,:));
end
for i=1:K2
B(i)=mean(SX2(i,:));
end
% 用k-Means聚类法确定最佳的聚类数
X=A'; %换成B就是研究白葡萄酒
numC=15;
for i=1:numC
kidx = kmeans(X,i);
silh = silhouette(X,kidx); %计算轮廓值
silh_m(i) = mean(silh); %计算平均轮廓值
end
figure
plot(1:numC,silh_m,'ko-', 'linewidth',2)
set(gca,'linewidth',2);
xlabel('类别数')
ylabel('平均轮廓值')
title(' 不同类别对应的平均轮廓值')
% 绘制2至5类时的轮廓值分布图
figure
set(gca,'linewidth',2);
for i=2:5
kidx = kmeans(X,i);
subplot(2,2,i-1);
[~,h] = silhouette(X,kidx);
set(gca,'linewidth',2);
title([num2str(i), '类时的轮廓值 ' ])
snapnow
xlabel('轮廓值');
ylabel('类别数');
end
把代码X=A'换成X=B'，就是对白葡萄酒进行分析，输出效果这里省略。
由于样本数量不到30，比较适合的类别数是2~5，综上分析，类别为4比较合适。
得到类别数量，用K-Means方法、层次聚类、模糊聚类对数据进行聚类
clc
clear all
%导入数据
X1=readtable('2012A_T1_processed.xls','Sheet','T1_red_grape','Range','D3:M272');
X2=readtable('2012A_T1_processed.xls','Sheet','T2_red_grape','Range','D3:M272');
X3=readtable('2012A_T1_processed.xls','Sheet','T1_white_grape','Range','D3:M282');
X4=readtable('2012A_T1_processed.xls','Sheet','T2_white_grape','Range','D3:M282');
[m1,n1]=size(X1); % 返回行数和列数
K1=27; % 有27个样品
% 计算每个样品的总得分
for i=1:K1
for j=1:n1
SX1(i,j)=sum(table2array(X1(10*i-9:10*i,j)));
SX1(i,j)=(SX1(i,j)+sum(table2array(X2(10*i-9:10*i,j))))/2;
end
end
[m2,n2]=size(X3);
K2=28;
for i=1:K2
for j=1:n2
SX2(i,j)=sum(table2array(X3(10*i-9:10*i,j)));
SX2(i,j)=(SX2(i,j)+sum(table2array(X4(10*i-9:10*i,j))))/2;
end
end
% 计算红葡萄酒和白葡萄酒质量评分数据
for i=1:K1
A(i)=mean(SX1(i,:));
end
for i=1:K2
B(i)=mean(SX2(i,:));
end
%默认研究的是红葡萄酒，将下面的A换成B就是研究白葡萄酒
% K-means聚类过程，并将结果显示出来
[idx,ctr]=kmeans(A',4); % 用K-means法聚类
% 提取同一类别的样品号
c1=find(idx==1); c2=find(idx==2);
c3=find(idx==3); c4=find(idx==4);
figure
F1 = plot(find(idx==1), A(idx==1),'k:*', ...
find(idx==2), A(idx==2),'k:o', ...
find(idx==3), A(idx==3),'k:p', ...
find(idx==4), A(idx==4),'k:d');
set(gca,'linewidth',2);
set(F1,'linewidth',2, 'MarkerSize',8);
xlabel('编号','fontsize',12);
ylabel('得分','fontsize',12);
title('Kmeans方法聚类结果')
disp('聚类结果：');
disp(['第1类:' ,'中心点：',num2str(ctr(1)),' ','该类样品编号：', num2str(c1')]);
disp(['第2类:' ,'中心点：',num2str(ctr(2)),' ','该类样品编号：', num2str(c2')]);
disp(['第3类:' ,'中心点：',num2str(ctr(3)),' ','该类样品编号：', num2str(c3')]);
disp(['第4类:' ,'中心点：',num2str(ctr(4)),' ','该类样品编号：', num2str(c4')]);
% 层次聚类
X=A';
Y=pdist(X); %计算样品间的欧式距离
Z=linkage(Y, 'average'); % 利用类平均法创建系统聚类树
cn=size(X);
clabel=1:cn;
clabel=clabel';
figure
F2 = dendrogram(Z); % 绘制聚类树形图
set(gca,'linewidth',2);
title('层次聚类法聚类结果')
set(F2,'linewidth',2);
ylabel('标准距离');
% Fuzzy C-means聚类
X=A';
[center,U] = fcm(X,4);
Cid1 = find(U(1,:) ==max(U));
Cid2 = find(U(2,:) ==max(U));
Cid3 = find(U(3,:) ==max(U));
Cid4 = find(U(4,:) ==max(U));
figure
F3=plot(Cid1, A(Cid1),'k:*', ...
Cid2, A(Cid2),'k:o', ...
Cid3, A(Cid3),'k:p', ...
Cid4, A(Cid4),'k:d');
set(gca,'linewidth',2);
set(F3,'linewidth',2, 'MarkerSize',8);
xlabel('编号');
ylabel('得分');
title('Fuzzy C-means方法聚类结果')
3种聚类方法对该问题的聚类效果差不多，不妨以K-Means方法得到的结果作为本问题的分级依据。
聚类结果：
第1类:中心点：69.08 该类样品编号：4 6 7 8 25
第2类:中心点：78.3667 该类样品编号：2 3 9 17 20 23
第3类:中心点：63.46 该类样品编号：1 11 12 15 18
第4类:中心点：73.2409 该类样品编号：5 10 13 14 16 19 21 22 24 26 27
默认研究的是红葡萄酒，将上面代码中的A换成B就是研究白葡萄酒
(4)根据葡萄酒的理化指标对酿酒葡萄进行分级
由于葡萄酒的理化指标存在一级指标和二级指标，为了保持指标级别的一致性，统一选择一级指标作为研究对象。红葡萄酒的一级理性指标有9个，根据这9个指标进行聚类，难度较大，可用PCA(主成分分析)方法对指标进行降维，然后利用K-Means方法对降维结果进行聚类。
excel表格整一整，如图：
2012A_Table2.xls
%数据导入及处理
clc, clear all, close all
A=table2array(readtable('2012A_Table2.xls','Sheet','葡萄酒指标汇总','Range','B3:J29'));%红葡萄酒
% A=table2array(readtable('2012A_Table2.xls','Sheet','葡萄酒指标汇总','Range','C33:J60')); %白葡萄酒
%数据标准化处理
a=size(A,1);
b=size(A,2);
for i=1:b
SA(:,i)=(A(:,i)-mean(A(:,i)))/std(A(:,i));
end
%计算相关系数矩阵的特征值和特征向量
CM=corrcoef(SA); %计算相关系数矩阵(correlation matrix)
[V, D]=eig(CM); %计算特征值和特征向量
for j=1:b
DS(j,1)=D(b+1-j, b+1-j); % 对特征值按降序进行排序
end
for i=1:b
DS(i,2)=DS(i,1)/sum(DS(:,1)); %贡献率
DS(i,3)=sum(DS(1:i,1))/sum(DS(:,1)); %累积贡献率
end
%选择主成分及对应的特征向量
T=0.8; %主成分信息保留率
for K=1:b
if DS(K,3)>=T
Com_num=K;
break;
end
end
%提取主成分对应的特征向量
for j=1:Com_num
PV(:,j)=V(:,b+1-j);
end
%计算各评价对象的主成分得分
new_score=SA*PV;
for i=1:a
total_score(i,1)=sum(new_score(i,:));
total_score(i,2)=i;
end
result_report=[new_score, total_score]; %将各主成分得分与总分放在同一个矩阵中
result_report=sortrows(result_report,(K+2)); %按总分降序排序
% Kmeans聚类及结果报告
A=result_report(:,(K+1));
[idx,ctr]=kmeans(A,4);
[m,n]=size(A);
t1=ones(1,n)*30;
c1=find(idx==1); c2=find(idx==2); c3=find(idx==3); c4=find(idx==4);
h=plot(t1,A,'ko',c1,A(idx==1),'k--*', c2,A(idx==2),'k--s', c3,A(idx==3),'k--d', c4,A(idx==4),'k--p');
xlabel('红葡萄酒样品编号','fontsize',12);
ylabel('主成分得分','fontsize',12);
title('红葡萄酒理化指标聚类图','fontsize',12)
set(h,'MarkerSize',8, 'MarkerFaceColor','k');
set(gca,'linewidth',2);
disp('主成分得分(最后1列为样本编号，倒数第2列为总分，前面为各主成分得分)')
result_report
disp('分类结果：');
disp(['第1类:' ,'中心点：',num2str(ctr(1)),' ','该类样品编号：', num2str(c1')]);
disp(['第2类:' ,'中心点：',num2str(ctr(2)),' ','该类样品编号：', num2str(c2')]);
disp(['第3类:' ,'中心点：',num2str(ctr(3)),' ','该类样品编号：', num2str(c3')]);
disp(['第4类:' ,'中心点：',num2str(ctr(4)),' ','该类样品编号：', num2str(c4')]);
输出
分类结果：
第1类:中心点：-0.28266 该类样品编号：1 4 5 6 8 12 14 17 19 22
第2类:中心点：-4.8813 该类样品编号：2 3 9 21 23
第3类:中心点：3.0977 该类样品编号：10 11 15 16 20 25
第4类:中心点：1.4411 该类样品编号：7 13 18 24 26 27
在上面的程序中，默认研究的是红葡萄酒的理化指标，想研究白葡萄酒的换个数据就行。
小结：对上面两种分级方法的结果进行分析，不难发现，它们结果的一致性较差。也就是说，质量评分高的葡萄酒，理化指标得分却不一定高，两者之间没有明显关系。但这种不一致是合理的，依据酒的质量评分分级主要依据品酒师的感觉，依据理化指标进行分级单纯依据指标数值的高低，两者的评价的角度不一样，得到的结果也不一样。
四、问题③的解决
要分析酿酒葡萄与葡萄酒的理化指标之间的联系，由于两者的理化指标比较多，是多对多的关系。可从以下几个角度展开：
①只研究酿酒葡萄和葡萄酒共有的指标，比如总酚、花色苷，然后可用拟合、回归方法给出两者共有指标的关系。
②先用PCA或方差分析，分别找出酿酒葡萄和葡萄酒的几个主要指标，然后再研究几个主要指标之间的关系。
③求出指标之间的相关系数，筛选相关性比较强的指标，再拟合、回归它们的关系。
五、问题④的解决
该问题是对以上3个问题的总结。
