2012年全国大学生数学建模竞赛A题国家二等奖论文

葡萄酒的评价

摘要

本文主要采用了方差分析、聚类分析、BP神经网络、偏最小二乘回归分析以及主成分分析方法，对葡萄酒的质量评价问题进行了深入探讨。

对于问题1，本文选用方差分析模型,首先分别求出了两组评酒员的组间评价结果，再通过计算得出第一、二组评酒员对红葡萄酒的质量评价无显著性差异，对白葡萄酒的评价有显著性差异。最后通过计算两组评酒员分别对于红酒和白酒评分的方差，得出第二组的方差更小，从而结果更可靠。

对于问题2，在问题1的求解基础上，采用结果更可靠的第二组专家的评分结果，作为评价葡萄酒质量的依据。先用R型聚类分析法对酿酒葡萄的理化指标进行分类，选出其中较具有代表的理化指标，在此基础上，再根据选出的理化指标用Q型聚类分析模型对各样品葡萄进行分类。然后利用葡萄酒的得分建立 BP神经网络模型对酿酒葡萄进行分级，最后综合前面两种模型的结论得到酿酒葡萄的分级表，将两种酿酒葡萄分为I、II、III、IV四个等级，即红葡萄4、13、25、26、27 为I级；3、5、7、9、17、20、21、22、23、24为II级；1、2、8、10、14、16为III级；6、11、12、15、18、19为IV级，白葡萄1、3、4、6、9、10、15、18、19、21、23、27、28 为I级；2、5、14、17、20、22、24、25、26为II级；7、8、11、12、13为III级；16为IV级。

对于问题3，考虑酿酒葡萄与葡萄酒的各项理化指标个数较多，我们建立偏最小二乘回归模型，求得酿酒葡萄与葡萄酒的主要理化指标两者之间的联系。经计算可得出结论：在红色的酿酒葡萄中，花色苷对红色葡萄酒的影响最大，其他的影响都不大；而在白色的酿酒葡萄中，果皮颜色对白色葡萄酒各项指标都存在着较大的影响。

对于问题4，建立主成分分析模型分别分析酿酒葡萄和葡萄酒的理化指标对葡萄酒质量的影响，然后再用基于主成分分析的综合评价模型论证能否用葡萄酒和葡萄酒的理化指标对葡萄酒的质量进行评价。得出结论为：酿酒葡萄中的单宁、蛋白质和果皮质量与澄清度有关；酿酒葡萄中的酯类和酚类与香气有关；葡萄酒的单宁、酒石酸、苹果酸和柠檬酸与口感有密切联系。 葡萄酒的质量与葡萄及葡萄酒的理化指标有密切联系，可以用它们来部分评价葡萄酒的质量。

本文中用到了多种统计学方法，将统计学方法应用于葡萄酒质量分析与评价中，可以更加清楚地了解葡萄酒成分与感官质量之间的相互关系，为葡萄酒的质量控制、预测、预报、等级区分提供了一种有效的途径。最后，对模型的优缺点进行了分析，并提出了合理的模型改进方向。

关键词：BP神经网络模型 偏最小二乘回归分析 聚类分析 葡萄酒

1

一、问题重述

1.问题的背景

葡萄是我们日常生活中常见的水果之一，随着人民生活水平的提高，用葡萄酿造的葡萄酒正被越来越多的人所享用。葡萄、葡萄酒源自伊朗，葡萄酒之历史，是人类孜孜不倦努力创造的文化史的一部分。在远古时代的黑暗中，她步伐缓慢，而且显得暧昧模糊。这第一步大概就开始于某一天，人类偶然创造出葡萄酒；从此以后，葡萄酒随着时代的变迁而变化。人类创造发明了葡萄酒，葡萄酒相应地给人以影响，这种影响又进一步驱使人们奔向新的工作目标的原因。葡萄酒就是在这种漫长的因果关系中，不断地与每个时代的人们的感官知觉相适应，并一直成长到今天。

发展到今天，随着酿造工艺和技术的提高，葡萄酒的种类越来越多，质量越来越好，但如何对葡萄酒的质量进行确定呢？确定葡萄酒质量时一般是通过聘请一批有资质的评酒员进行品评。每个评酒员在对葡萄酒进行品尝后对其分类指标打分，然后求和得到其总分，从而确定葡萄酒的质量。酿酒葡萄的好坏与所酿葡萄酒的质量有直接的关系，葡萄酒和酿酒葡萄检测的理化指标会在一定程度上反映葡萄酒和葡萄的质量。附件1给出了某一年份一些葡萄酒的评价结果，附件2和附件3分别给出了该年份这些葡萄酒的和酿酒葡萄的成分数据。 2.问题的提出

请根据所给资料尝试建立数学模型讨论下列问题：

1. 分析附件1中两组评酒员的评价结果有无显著性差异，哪一组结果更可信？

2. 根据酿酒葡萄的理化指标和葡萄酒的质量对这些酿酒葡萄进行分级。 3. 分析酿酒葡萄与葡萄酒的理化指标之间的联系。

4．分析酿酒葡萄和葡萄酒的理化指标对葡萄酒质量的影响，并论证能否用葡萄和葡萄酒的理化指标来评价葡萄酒的质量？

二、问题分析

对问题1的分析，由附件1所给数据可知，附件中有两个数据有明显错误，故先对其进行修正。此外，另有一空数据缺失，补全该空的空缺数据。对于解决显著性差异问题，常采用方差分析模型和t检验法求解，在本题的第一问中，我们采用了方差分析法，发现该方法可以很好地解决本题所提出的问题。除此之外，第一问还要求确定哪组结果更可信，通过计算各组评酒员分别对红酒和白酒打分的方差可以判断两组评酒员打分的波动性大小，取波动性较小的一组，说明该组品酒员对各样品酒打分比较合理，即评价结果更可信。

2

对问题2的分析，问题要求根据葡萄酒的理化指标和葡萄酒的质量对葡萄酒进行分级，可先用第二组专家对各样品葡萄酒所打的分数排出各葡萄酒的得分顺序，即各葡萄酒的质量高低。对于数据量大的多指标的分级问题，已有很多成型的算法，如：主成分分析、层次分析、聚类分析、BP神经网络。由于在本题中采用主成分分析法和层次分析法步骤太过于繁琐，不易得出结论。因此，我们选用聚类分析和神经网络模型求解本题并得到了比较精确的结论。

对问题3的分析，本题要求分析酿酒葡萄与葡萄酒的理化指标之间的关系，解决这类问题常用的方法有灰色关联分析、典型相关分析、线性回归分析、偏最小二乘回归分析。由于偏最小二乘回归分析集中了主成分分析、典型相关分析和线性回归分析等各种方法的特点，在分析结果中，除了可以提供一个更为合理的回归模型外，还可以同时完成一些类似于主成分分析和典型相关分析的研究内容，提供一些更丰富、深入的信息，故在本题中采用偏最小二乘回归分析法。

对问题4的分析，由于酿酒葡萄和葡萄酒的理化指标对葡萄酒的影响程度难以直接量化，影响葡萄酒质量的因素种类不确定。本文采用主成分发分析法分别确定酿酒葡萄的理化指标对葡萄酒质量的影响以及葡萄酒的理化指标对葡萄酒质量的影响，同时用基于主成分分析法的综合评价法对能否用葡萄和葡萄酒的理化指标来评价葡萄酒的质量进行论证。

对于本文中的各个问题，分别采用了不同的模型和方法，并对各模型进行了评价和推广。

三、模型假设

1、假设第一组的评酒员为A组评酒员，第二组的评酒员为B组评酒员； 2、假设每种样品葡萄酒的最后得分为各位品酒员打分的均值；

3、除个别数据外，文中所给的其他数据均是真实可信的，能反映各品酒员对各葡萄酒样品的客观评价；

4、从附件2中可以得知，酿酒葡萄的理化指标比较多，且有些一级指标已经包含了二级指标，故不考虑二级指标对解题结果的影响； 5、假设评酒员对样品酒所打的分数主要服从正态分布。

四、符号说明

Aij：表示第一组的评酒员对i种酒样品j所打的总分；

Bij：表示第二组的评酒员对i种酒样品j所打的总分（i?1,2;j?1,2?,27）如A12表示第一组评酒员对白酒样品2打的总分；

3

xij：第i组对第j种酿酒葡萄样品的总体观测值； yij：第i组对第j种葡萄酒样品的总体观测值；

ST:偏差得分方程；

pi：评酒员i对酒的评价得分；

S

2

：评分方差值；

：各个样本之间的欧式距离；

个神经元的连接权值；

dijwij：输入层第i个神经元与隐含层第j?j：隐含层第j个神经元的阀值；

Lt：输入第t个神经元； xj：自变量组酿酒葡萄指标值向量，j?1,2,?6； ：因变量组葡萄酒的指标向量，j?1,2,?4; :指标观测值；

yjaijui；自变量组第i个主成分； vi：因变量组第i个主成分； ~xj：数据标准化后的样本自变量； ~yj：数据标准化后的样本因变量。

五、模型的建立与求解

问题1

模型I 方差分析模型 1.1 数据处理

4

先对附件1中的错误数据进行修正，由于对持久性指标打的分值不能超过8分，故附录一中第一组白葡萄酒品尝评分的第233行第10列以及第298行第12列的数据均是错误的，分别取该错误数字的末位，即分别修改为7和6。而对于第76行第6列中空缺的数据，常用的方法是插值或拟合，由于各评酒员对同一种样品酒的同一指标所打的分数没有明显的线性组合关系或其他显著规律，故不宜采用插值或拟合的方法。通过观察数据可以发现10个评酒员中已有6个评酒员对该指标打了6分，故可考虑将该空缺值补为6。另外，求得已知9个评酒员对该指标的平均分为6.22，由于所打分数均为整数，可取为6，由此可以说明将该空缺值补为6是比较合理的。数据完整后，先求出各评酒员对各样品酒所打的总分，然后再求得各组评酒员对各样品酒打分的均值，如表1所示。 表1 各样品酒的得分表 评 分 样 A组红酒评分 A组白酒评分 Ｂ组红酒评分 B组白酒评分 品 酒样品1 酒样品2 酒样品3 酒样品4 酒样品5 酒样品6 酒样品7 酒样品8 酒样品9 酒样品10 酒样品11 酒样品12 酒样品13 酒样品14 酒样品15 酒样品16 酒样品17 酒样品18 酒样品19 酒样品20 酒样品21 酒样品22 酒样品23 酒样品24 酒样品25 酒样品26 酒样品27 酒样品28

62.7 80.3 80.4 68.6 73.3 72.2 71.5 72.3 81.5 74.2 70.1 53.9 74.6 73 58.7 74.9 79.3 59.9 78.6 78.6 77.1 77.2 85.6 78 69.2 73.8 73 82 74.2 78.3 79.4 71 68.4 77.5 70.4 72.9 74.3 72.3 63.3 65.9 72 72.4 74 78.8 73.1 72.2 77.8 76.4 71 75.9 73.3 77.1 81.3 64.8 81.3 68.1 74 74.6 71.2 72.1 66.3 65.3 66 78.2 68.8 61.6 68.3 68.8 72.6 65.7 69.9 74.5 65.4 72.6 75.8 72.2 71.6 77.1 71.5 68.2 72 71.5 77.9 75.8 75.6 76.9 81.5 75.5 74.2 72.3 80.4 79.8 71.4 72.4 73.9 77.1 78.4 67.3 80.3 76.7 76.4 76.6 79.2 79.4 77.4 76.1 79.5 74.3 77 79.6 5

表4 红色酿酒葡萄分级表 等级 葡萄样品 I 4、13、25、26、27 3、5、7、9、17、 II 20、21、22、23、24 III 1、2、8、10、14、16 IV 6、11、12、15、18、19 表5 白色酿酒葡萄分级表 等级 葡萄样品 1、3、4、6、9、10、15、 I 18、19、21、23、27、28 2、5、14、17、20、22、II 24、25、26 III 7、8、11、12、13 IV 16 从表中可以看出，对于红色酿酒葡萄样品，各等级分布比较均匀，且 主要分布在第二等级，而对于白色酿酒样品，分布不怎么均匀，只有少数样品酒分布在第三、第四等级，且集中在第一等级，这与葡萄酒的质量排名吻合度较高，进一步说明了该分类方法的正确性。 问题3

模型III：偏最小二乘回归分析

考虑p个变量y1,y2,?yp与m个自变量x1,x2,?xm的建模问题。偏最小二乘回归分析的基本做法是首先在自变量集中提出第一成分t1（ t1是x1,x2,?xm的线性组合，且尽可能多地提取原自变量集中的变异信息）；同时在因变量集中也提取第一成分u1，并要求t1与u1相关程度达到最大。然后建立因变量y1,y2,?yp与

t1 的回归，如果回归方程已达到满意的精度，则算法中止。否则继续第二对自

变量与因变量成分的提取，直到能达到满意的精度为止。最终对自变量集提取r个成分u1,u2,?ur，然后再表示为yp与原自变量xm的回归方程式，即偏最小二乘回归方程式。

16

根据问题2中葡萄酒的理化指标，我们通过R型聚类分析与综合选定了附件2中酿酒葡萄的基本理化指标为：VC含量、花色苷、白藜芦醇、可滴定酸、固酸比、果皮颜色。我们分别令它们作为自变量x1,x2,?x6；附件2中葡萄酒的各项指标显然存在着关联，为了精确简练的凸显它们之间的联系，依照问题2的方法，我们也对其进行了一次R型聚类（程序见附录）。通过分析，我们选定了葡萄酒的基本理化指标：花色苷，白藜芦醇 、DPPH半抑制体积、色泽L*。分别令其为因变量y1,y2,?y4。 模型的基本步骤如下： 步骤1：数据的处理 根据附表2的数据，我们应用excel基本知识对其多组测量的数据进行求平均，然后再提取作为应用。 自变量的观测数据矩阵记为A??aij?27?6，因变量的观测数据矩阵记为B??bij?~，有 aij转化成标准化指标a。将各指标值ij27?4~?aijaij??jsj(1)j(1)(1)，i?1,2,?27,j?1,2,?6， 127ij其中：?(1)j??a，s27iji?1127??(a27?1i?1??j)，j?1,2,?6,即?(1)2(1)j、s(j1)为的j个自变量xj的样本均值和样本方差。 类似的，将bij转化成标准化指标值bij，有 bij??j~bij?(2)sj(2)~，i?1,2,?27,j?1,2,3,4， 127ij其中：?(2)j?12727?i?1bij，s(2)j?(a?27?1i?1??j)，j?1,2,3,4,即?(2)2(2)j、s(j2) 为的j个因变量xj的样本均值和样本方差。 步骤2：求相关系数矩阵

17

表6 相关系数矩阵表

VC含量

VC含量

花色苷 白藜芦醇 可滴定酸 固酸比 果皮颜色 花色苷 白藜芦醇 DPPH 色泽L*

如表所示，给出了这10个变量的简单相关矩阵。从相关系数矩阵可以看出，酿酒葡萄的VC含量、果皮颜色与红葡萄酒中的花色苷，白藜芦醇 、DPPH半抑制体积成负相关，与色泽L*成正相关；酿酒葡萄的花色苷与红葡萄酒中的花色苷，白藜芦醇 、DPPH半抑制体积成正相关，与色泽L*成负相关；酿酒葡萄的白藜芦醇与红葡萄酒中的花色苷成负相关，与白藜芦醇 、DPPH半抑制体积、色泽L*成正相关；酿酒葡萄的可滴定酸与红葡萄酒中的花色苷、色泽L*成负相关，与白藜芦醇 、DPPH半抑制体积成正相关；酿酒葡萄的固酸比与红葡萄酒中的花色苷、DPPH半抑制体积成负相关，与白藜芦醇、色泽L*成正相关。

步骤3：分别提出自变量和因变量组的成分 使用matlab软件编程求得各对成分分别为：

对于自变量xj，有主成分如下：

?u1?0.26x1?0.46x2?0.3x3?0.43x4?0.09x5?0.65x6??u2?0.03x1?0.10x2?0.39x3?0.69x4?0.57x5?0.14x6 ???1 -0.11 0.27 0.1313 -0.019 0.1313 -0.089 -0.028 -0.122 0.1224

花色苷 -0.11

1 -0.06 0.129 -0.25 -0.52 0.923 0.2 0.671 -0.83 白藜芦醇 0.27 -0.06

1 0.182 -0.41 0.309 -0.04 0.014 0.073 0.162 可滴定酸 0.131 0.129 0.182

1 -0.51 0.039 -0.02 0.179 0.232 -0.14 固酸比 -0.02 -0.25 -0.41 -0.51

1 -0.25 -0.22 0.109 -0.06 0.198 果皮颜色 0.13 -0.5 0.31 0.04 -0.2

1 -0.4 -0.01 -0.4 0.49 花色苷 -0.09 0.923 -0.04 -0.02 -0.22 -0.37

1 0.124 0.675 -0.82 白藜芦醇 -0.03 0.2 0.014 0.179 0.109 -0.03 0.124

1 0.529 -0.35

DPPH -0.1 0.67 0.07 0.23 -0.1 -0.4 0.68 0.53

1 -0.8

色泽L* 0.12 -0.8 0.16 -0.1 0.2 0.49 -0.8 -0.4 -0.8

1

对于因变量yj，有主成分如下：

?v1?0.002y1?0.0163y2?0.999y3?0.0019y4??v2??0.0068y1?0.9939y2?0.0160y3?0..109y4 ???由于前两个主成分的累计贡献率都超过了95%，只要取两个成分就可以了 步骤4：求两个成分对时标准化指标变量与成分变量之间的回归方程 求得自变量组与因变量组与u1,u2之间的回归方程分别为；

18

~x1??0.1167u1?0.2398u2 ~x2?0.8657u1?0.6370u2 ~x3??0.0508u1?0.2907u2 ~x4?0.1247u1?0.4102u2 ~x5??0.1639u1?0.0604u2 ~x6??0.4381u1?0.5659u2 ~y1?0.6626u1?0.5673u2 ~y2?0.1286u1?0.0653u2 ~y3?0.5294u1?0.1622u2 ~y4??0.6665u1?0.1785u2 步骤5：确定因变量组与自变量组之间的回归方程 yi的回归方程，得到标准化指标之间的把步骤（3）中的ui代入步骤（4）中~回归方程 ~y1?0.06~x1?0.93~x2?0.13~x3?0.15~x4?0.07~x5?0.03~x6~y2?0.007~x1?0.15~x2?0.01~x3?0.01~x4?0.02~x5?0.02~x6 ~y??0.0229~x?0.56~x?0.02~x?0.0005~x?0.077~x?0.14~x3123456~y4?0.035~x1?0.69~x2?0.018~x3?0.01~x4?0.0985~x5?0.191~x6 将标准化变量~xj,~yj分别还原成原10.80.60.4始变量xj,yj，得到回归方程： y1?397.507?6.916x1??2.4x2?5.5x3?1400.2.74x4?12.51x5?6.1174x600.1271x?0.0363x?0.0484xy2?5.0373?0.001x1?0.0049x2?0.0066x3?456-0.2y3?0.5926?0.0015x1?0.0008x2?0.0005x?0.0003x4?0.0072x5?0.0155x63y4??41.3713?0.3826x1?0.1647x2?0.0703x3?0.8619x4?1.5401x5?3.5302x6-0.4-0.6-0.8 123 图9 回归系数图 4 19

步骤6：模型的解释与检验 为了更直观、迅速地观察各个自变 量在解释yj时的边际作用，不妨绘制回归系数图，见右图。这个图是针对标准化数据的回归方程的。 从回归系数图中可以立刻观察到，酿酒葡萄中的花色苷对红葡萄酒影响最大，其他的都影响不大，换言之，酿酒葡萄中花色苷一类对葡萄酒的指标联系最为紧密。 为了考察这四个回归方程的模型精度，我们以yj为参考，对所有样本值进行绘制预测图。在这个预测图上，如果所有点能在图的对角线附近均匀分布，则方程的拟合度较高，比较令人满意。如图10所示： 10008006004002000050010000051015101550.80.60.40.2000.20.40.60.8100806040200-50050100 图10 影响程度预测图 从图10中，我们可以发现第二个图像的拟合效果最差，即葡萄酒中的白藜芦醇一类酿酒葡萄指标之间的联系不是很紧密。 同理，我们可以分析白葡萄酒的酿酒葡萄与0.3酒之间的联系

0.2我们通过matlab编程可以得出回归方程

?y1??y2??y3?y?4??0.9527?0.0732x1?0.0710x2?0.0509x3?0.0962x4?0.0217x5?0.0604x6?1.5804?0.0699x1?0.0166x2?0.0448x3?0.13031x4?0.0218x5?0.0221x6??0.0285?0.0091x1?0.0034x2?0.0004x3?0.-0.10022x4?0.0024x6?103.4179?0.2220x1?0.0009x2?0.0849x3?0-0.2.2772x4?0.0430x5?0.021x6-0.300.1

20

-0.41234图11 回归系数图

tm=reshape(tm,1,length(tm)); %变成行向量 fprintf('第%d类的有%s\\n',i,int2str(tm)); %显示分类结果 end

程序13 bar(xishu') load mydata num

ch0=repmat(ch0,num,1);

yhat=ch0+x0*xish; %计算y 的预测值 y1max=max(yhat); y2max=max(y0);

ymax=max([y1max;y2max]) cancha=yhat-y0; %计算残差 subplot(2,2,1)

plot(0:ymax(1),0:ymax(1),yhat(:,1),y0(:,1),'*') subplot(2,2,2)

plot(0:ymax(2),0:ymax(2),yhat(:,2),y0(:,2),'O') subplot(2,2,3)

plot(0:ymax(3),0:ymax(3),yhat(:,3),y0(:,3),'H') subplot(2,2,4)

plot(0:ymax(4),0:ymax(4),yhat(:,4),y0(:,4),'+')

程序14

function result=cwfac(vector); fprintf('协方差矩阵:\\n')

std=cov(vector) %计算相关系数矩阵CORRCOEF fprintf('特征向量(vec)及特征值(val)：\\n')

[vec,val]=eig(std) %求特征值(val)及特征向量(vec)

newval=diag(val) ;

[y,i]=sort(newval) ; %对特征根进行排序，y为排序结果，i为索引 fprintf('特征根排序：\\n') for z=1:length(y)

newy(z)=y(length(y)+1-z); end

fprintf('%g\\n',newy) rate=y/sum(y);

fprintf('\\n贡献率：\\n') newrate=newy/sum(newy) sumrate=0; newi=[];

for k=length(y):-1:1

sumrate=sumrate+rate(k); newi(length(y)+1-k)=i(k);

31

if sumrate>0.85 break;

end

end %记下累积贡献率大85%的特征值的序号放入newi中 fprintf('主成分数：%g\\n\\n',length(newi)); fprintf('主成分载荷：\\n') for p=1:length(newi) for q=1:length(y)

result(q,p)=sqrt(newval(newi(p)))*vec(q,newi(p)); end

end disp(result)

%计算载荷 32

break end

fprintf('**********************************\\n'); end

程序5

zscore(C); %标准化处理数据

d=pdist(ans','correlation'); %计算相关系数导出的距离 z=linkage(d,'average'); %按类平均法聚类 h=dendrogram(z); %画聚类图

set(h,'Color','k','LineWidth',1.3) %把聚类图线的颜色改成黑色，线宽加粗 T=cluster(z,'maxclust',6) ; %把变量划分成12类 for i=1:12

tm=find(T==i); %求第i类的对象

tm=reshape(tm,1,length(tm)); %变成行向量

fprintf('第%d类的有%s\\n',i,int2str(tm)); %显示分类结果 end

程序6

gj=C(:,1:3);gj(:,4:5)=C(:,7:7:14);gj(:,6)=C(:,23); gj=zscore(gj); %数据标准化

y=pdist(gj); %求对象间的欧氏距离,每行是一个对象 z=linkage(y,'average'); %按类平均法聚类 h=dendrogram(z); %画聚类图 for k=6:12

fprintf('划分成%d类的结果如下：\\n',k)

T=cluster(z,'maxclust',k); %把样本点划分成k类 for i=1:k

tm=find(T==i); %求第i类的对象

tm=reshape(tm,1,length(tm)); %变成行向量

fprintf('第%d类的有%s\\n',i,int2str(tm)); %显示分类结果 end if k==5 break end

fprintf('**********************************\\n'); end

程序7

p_test=p_test'; input=size(p,1); output=size(t,1);

net=newff(minmax(p),[61,4],{'tansig','logsig'},'trainlm'); net.trainparam.epochs=1000;

26

net.trainparam.goal=0.01; lp.lr=0.1;

net=train(net,p,t);

disp(['1.使用权值和阀值']) disp('测试样本预测结果:' ) Y1=sim(net,p_test) err=norm(Y1-t_test); err1=norm(sim(net,p)-t);

disp(['测试样本的仿真误差: 'num2str(err)]) disp(['训练样本的仿真误差: 'num2str(err1)])

程序8

p_test=p_test'; input=size(p,1); output=size(t,1);

net=newff(minmax(p),[61,4],{'tansig','logsig'},'trainlm'); net.trainparam.epochs=1000; net.trainparam.goal=0.01; lp.lr=0.1;

net=train(net,p,t);

disp(['1.使用权值和阀值']) disp('测试样本预测结果:' ) Y1=sim(net,p_test) err=norm(Y1-t_test); err1=norm(sim(net,p)-t);

disp(['测试样本的仿真误差: 'num2str(err)]) disp(['训练样本的仿真误差: 'num2str(err1)])

程序9 clc,clear

load pz.txt %原始数据存放在纯文本文件pz.txt 中 mu=mean(pz);sig=std(pz); %求均值和标准差 rr=corrcoef(pz); %求相关系数矩阵 data=zscore(pz); %数据标准化

n=3;m=3; %n 是自变量的个数,m 是因变量的个数 x0=pz(:,1:n);y0=pz(:,n+1:end); e0=data(:,1:n);f0=data(:,n+1:end); num=size(e0,1);%求样本点的个数

chg=eye(n); %w 到w*变换矩阵的初始化 for i=1:n

%以下计算w，w*和t 的得分向量， matrix=e0'*f0*f0'*e0;

[vec,val]=eig(matrix); %求特征值和特征向量

27

val=diag(val); %提出对角线元素 [val,ind]=sort(val,'descend');

w(:,i)=vec(:,ind(1)); %提出最大特征值对应的特征向量 w_star(:,i)=chg*w(:,i); %计算w*的取值 t(:,i)=e0*w(:,i); %计算成分ti 的得分

alpha=e0'*t(:,i)/(t(:,i)'*t(:,i)); %计算alpha_i

chg=chg*(eye(n)-w(:,i)*alpha'); %计算w 到w*的变换矩阵 e=e0-t(:,i)*alpha'; %计算残差矩阵 e0=e;

%以下计算ss(i)的值

beta=[t(:,1:i),ones(num,1)]\\f0; %求回归方程的系数 beta(end,:)=[]; %删除回归分析的常数项 cancha=f0-t(:,1:i)*beta; %求残差矩阵

ss(i)=sum(sum(cancha.^2)); %求误差平方和 %以下计算press(i) for j=1:num t1=t(:,1:i);f1=f0;

she_t=t1(j,:);she_f=f1(j,:); %把舍去的第j 个样本点保存起来 t1(j,:)=[];f1(j,:)=[]; %删除第j 个观测值

beta1=[t1,ones(num-1,1)]\\f1; %求回归分析的系数 beta1(end,:)=[]; %删除回归分析的常数项 cancha=she_f-she_t*beta1; %求残差向量 press_i(j)=sum(cancha.^2); end

press(i)=sum(press_i); if i>1

Q_h2(i)=1-press(i)/ss(i-1); else

Q_h2(1)=1; end

if Q_h2(i)<0.0975

fprintf('提出的成分个数r=%d',i); r=i; break end end

beta_z=[t(:,1:r),ones(num,1)]\\f0; %求Y 关于t 的回归系数 beta_z(end,:)=[]; %删除常数项

xishu=w_star(:,1:r)*beta_z; %求Y关于X的回归系数，且是针对标准数据的回归系数，

每一列是一个回归方程

mu_x=mu(1:n);mu_y=mu(n+1:end); sig_x=sig(1:n);sig_y=sig(n+1:end); for i=1:m

28

ch0(i)=mu_y(i)-mu_x./sig_x*sig_y(i)*xishu(:,i); %计算原始数据的回归方程的常数 项 end

for i=1:m

xish(:,i)=xishu(:,i)./sig_x'*sig_y(i); %计算原始数据的回归方程的系数，每一列是一个回归方程 end

sol=[ch0;xish] %显示回归方程的系数，每一列是一个方程，每一列的第一个数是 常数项

save mydata x0 y0 num xishu ch0 xish

程序10

Pz=[ 191 36 50 5 162 60 189 37 52 2 110 60 193 38 58 12 101 101 162 35 62 12 105 37 189 35 46 13 155 58 182 36 56 4 101 42 211 38 56 8 101 38 167 34 60 6 125 40 176 31 74 15 200 40 154 33 56 17 251 250 169 34 50 17 120 38 166 33 52 13 210 115 154 34 64 14 215 105 247 46 50 1 50 50 193 36 46 6 70 31 202 37 62 12 210 120 176 37 54 4 60 25 157 32 52 11 230 80 156 33 54 15 225 73 138 33 68 2 110 43];

save mydata x0 y0 num xishu ch0 xish

程序11 bar(xishu') load mydata num

ch0=repmat(ch0,num,1);

yhat=ch0+x0*xish; %计算y 的预测值y1max=max(yhat); y2max=max(y0);

ymax=max([y1max;y2max]) cancha=yhat-y0; %计算残差

29

subplot(2,2,1)

plot(0:ymax(1),0:ymax(1),yhat(:,1),y0(:,1),'*') subplot(2,2,2)

plot(0:ymax(2),0:ymax(2),yhat(:,2),y0(:,2),'O') subplot(2,2,3)

plot(0:ymax(3),0:ymax(3),yhat(:,3),y0(:,3),'H')

程序12

A=[973.878 11.030 9.983 8.020 2.4382 0.358 517.581 11.078 9.560 13.300 3.6484 0.460 398.770 13.259 8.549 7.368 5.2456 0.396 183.519 6.477 5.982 4.306 2.9337 0.177 280.190 5.849 6.034 3.644 4.9969 0.207 117.026 7.354 5.858 4.445 4.4311 0.211 90.825 4.014 3.858 2.765 1.8205 0.112 918.688 12.028 10.137 7.748 1.0158 0.346 387.765 12.933 11.313 9.905 3.8599 0.386 138.714 5.567 4.343 3.145 3.2459 0.136 11.838 4.588 4.023 2.103 0.3816 0.105 84.079 6.458 4.817 2.986 2.1628 0.141 200.080 6.385 4.930 3.957 1.3388 0.166 251.570 6.073 5.013 3.068 2.1659 0.163 122.592 3.985 4.064 1.836 0.8886 0.068 171.502 4.832 4.044 2.668 1.1620 0.117 234.420 9.170 6.168 4.912 1.6504 0.310 71.902 4.447 4.353 3.531 1.7396 0.138 198.614 5.981 5.157 3.875 9.0269 0.167 74.377 5.864 4.858 4.044 0.9641 0.158 313.784 10.090 8.941 4.440 8.7937 0.358 251.017 7.105 6.199 5.827 4.4666 0.231 413.940 10.888 12.529 12.144 12.6821 0.566 270.108 5.747 5.394 3.731 6.8689 0.165 158.569 5.406 4.425 3.022 2.5789 0.165 151.481 3.615 3.889 2.154 2.7369 0.076 138.455 5.961 4.734 3.284 4.7758 0.151 ]; zscore(A); %标准化处理数据 d=pdist(ans','correlation'); %计算相关系数导出的距离 z=linkage(d,'average'); %按类平均法聚类 h=dendrogram(z); %画聚类图

set(h,'Color','k','LineWidth',1.3) %把聚类图线的颜色改成黑色，线宽加粗 T=cluster(z,'maxclust',3) %把变量划分成3类 for i=1:3

tm=find(T==i); %求第i类的对象

30

本文来源：https://www.bwwdw.com/article/8i9f.html