格式参考：应用统计方法解决典型数学模型.doc

东北石油大学本科生毕业设计（论文） 东北石油大学东北石油大学毕业设计（论文）任务书毕业设计（论文）任务书题目 应用统计方法求解典型数学模型专业 信息与计算科学 学号 111001140109 姓名 武杰主要内容、基本要求、主要参考资料等：主要内容：本文拟首先简单介绍应用统计方法基本理论，再通过实例来掌握如何应用统计方法求解典型数学模型。基本要求：1、查阅不少于 10 篇参考文献，其中至少 2 篇外文文献；2、完成不少于 1 万 5 千字的论文，其中中文摘要 300 字左右，并写出相应的英文摘要；3、完成与课题相关的外文资料翻译，字数不少于 5000 字。主要参考资料： 1、张志涌 杨祖樱等编著.MATLAB 教程M.北京:北京航空航天大学出版社.2006 年 8 月 2、谢中华 编著.MATLAB统计分析与应用：40个案例分析M.北京:北京航空航天大学出版社.2010年6月 3、韩中庚 著.数学建模方法及应用M.北京:高等教育出版社.2000年完成期限： 指导教师签名： 专业负责人签名： 年 月 日东北石油大学本科生毕业设计（论文） 摘 要 现实生活中, 由于客观事物内部规律的复杂性及人们认识程度的限制, 人们常搜集大量的数据, 基于数据的统计分析建立合乎机理规律的数学模型, 然后通过计算得到的模型结果来解释实际问题。而随着社会的不断发展，统计的知识越来越重要，运用抽样数据进行推断已经成为现代社会一种普遍适用并且强有力的思考方式。 本文在论述统计思想的广泛应用和简介 MATLAB 工具箱的基础上，通过介绍方差分析、判别分析、主成分分析三种常见统计分析方法的基础原理及其应用性实例，结合 MATLAB 编程，归纳总结了较为简便的计算机求解方法，并给出了具体答案。 最后，本文着重的探讨了典型数学模型应用统计方法的 MATLAB 求解，无论是 DNA 序列分类模型中运用的主成分分析法、FISHER 判别分析法，还是在解决葡萄酒评价问题中所用的单因素方差分析思想，都很好的体现了统计思想方法在解决典型数学模型问题时应用的广泛性和有效性。关键词：数学模型；方差分析；判别分析；主成分分析；MatlabAbstractIn real life, because of the complexity of the internal law of objective things and people know the degree of restriction, it is often to collect a large amount of data, based on statistical analysis of data create the mechanism mathematical model of the law, and then through the calculated model results to explain the practical problems.And with the continuous development of society, the knowledge of statistics is more and more important, using the sampling data to infer that has become a modern society widely used and powerful ways of thinking.This paper discusses the widely used of statistical thinking and on the basis of introduction to MATLAB toolbox , by introducing the variance analysis, discriminant analysis, principal component analysis based on analysis of three common statistical principle and its application examples, combining with the MATLAB programming, summarized the simple computer solving method, and presents a concrete answers.Finally, this paper discusses the typical statistical method of MATLAB to solve the mathematical model of the application, whether in the DNA sequence classification model using principal component analysis , Fisher discriminant analysis method, is used in solving for the evaluation of wine, the single factor analysis of variance of thoughts, all good embodies the statistical thinking method in solving mathematical model for a typical application of universality and effectiveness.Key words：Mathematical model; Analysis of variance; Discriminant analysis; Principal component analysis; MATLAB I目 录第第 1 章章 绪绪 论论.11.1 统计分析相关介绍统计分析相关介绍 .11.2 MATLAB 相关介绍 .11.3 本章小结 .3第第 2 章章 三类常见统计方法及其解决典型数学模型实例三类常见统计方法及其解决典型数学模型实例.42.1 方差分析 .42.2 判别分析 .122.3 主成分分析 .162.4 本章小结 .21第第 3 章章 典型数学模型应用统计方法求解典型数学模型应用统计方法求解.223.1 葡萄酒的评价 .223.2 DNA 序列分类.313.3 本章小结 .38结结 论论.39参考文献参考文献.40致致 谢谢.42附附 录录.431第1章 绪 论1.1 统计统计分析相关介绍分析相关介绍统计是随着概率论的发展而发展起来的。只有当人们认识到必须把数据视为来自具有一定概率分布的总体，所研究的对象是这个总体而不能局限于数据本身的时候，统计诞生了。统计分析作为一项长期的社会实践活动，在历史的发展中已有几千年。“统而计之”就是人们对数据统计分析的最朴素的认识。随着社会生产力的不断进步，社会经济的不断发展，当代的数据统计分析已不圄于“统而计之”的范畴。从高斯（Gauss）与勒让德（Legendre）关于观测数据的误差分析和最小二乘估计方法的研究成果，到克拉默的统计学的数学方法，统计才真正的发展为一门成熟的科学。回顾概率统计的发展历史，大致可分为如下三个阶段：古典时期(10 世纪末以前)：主要是描述性的统计学的形成和发展，也可以说是数理统计的萌芽时期。近代时期(19 世纪末至 1945 年)：小样本理论形成，数理统计的主要分支廷立，可以说是数理统计的形成时期。现代时期(1945 年以后)：这个时期，由于计算机的发明和使用，推动着数理统计在理论研究和应用方面都不断地向纵深发展，并产生一些新的分支和边缘性的新学科，如最优设计和非参数统计推断等。随着社会的不断发展，统计的知识越来越重要，运用抽样数据进行推断已经成为现代社会一种普遍适用并且强有力的思考方式。目前，统计的很多原理方法已被越来越多地应用到交通、经济、医学、气象等各种与人们生活息息相关的领域。本文将就统计的常用方法与思想，在经典数学模型的应用中展开一些讨论，推导出某些表面上并非直观的结论，从中可以看出统计思想方法在解决问题中的实用性、简捷性和高效性。1.2 MATLAB 相关介绍MATLAB （Matrix & Laboratory）是美国 MathWorks 公司自 20 世纪 80 年代中期推出的数学软件，优秀的数值计算能力和卓越的数据可视化能力使其很快在数学软件中脱颖而出.目前，MATLAB 已经发展成为多学科、多种工作平台的功能强大的大型软件，是涵括线性代数、自动控制理论、概率论及数理统计、数字信号处理、时间序列分析、动态系统仿真等高级课程的基本教学工具。2MATLAB 用法简易、灵活、程式结构强又兼具延展性.以下为其几个特色： 1.2.1 数值计算和符号计算功能（1）数值计算功能包括：矩阵运算、多项式和有理分式运算、数据统计分析、数值积分、优化处理等。（2）符号计算将得到问题的解析解。1.2.2 MATLAB 语言除了命令行的交互式操作以外，还可以程序方式工作。使用 MATLAB 很容易地实现 C 或 FORTRAN 语言的几乎全部功能。1.2.3 图形功能 两个层次的图形命令：一种是对图形句柄进行的低级图形命令，另一种是建立在低级图形命令之上的高级图形命令。利用 MATLAB 的高级图形命令可以轻而易举地绘制二维、三维乃至四维图形，并可进行图形和坐标的标识、视角和光照设计、色彩精细控制等等。如今，MATLAB 广泛应用于：数值计算、图形处理、符号运算、数学建模、系统辨识、小波分析、实时控制、动态仿真等领域.MATLAB 包括拥有数百个内部函数的主包和三十几种工具包.工具包又可以分为功能性工具包和学科工具包.功能型工具箱主要用来扩充 MATLAB 的符号计算功能、图形建模仿真功能、文字处理功能以及与硬件实时交互功能，能用于多种学科.而领域型工具箱是专业性很强的，如控制系统工具箱(Control System Toolbox)、财政金融工具箱(Financial Toolbox)等.在统计问题上，我们能用到的最基础的工具箱就是统计工具箱.它包括下面几个方面的内容：概率分布.给出了常见的 20 种概率分布类型的概率密度函数、累加分布函数、逆累加分布函数、参数估计函数、随机数生成函数和统计量计算函数.参数估计.提供了多种分布类型分布参数及其置信区间的估计方法. 样本描述.中心趋势和离中趋势的统计量函数，及统计量计算函数. 方差分析.包括单因子方差分析、双因子方差分析和多因子方差分析. 回归分析.包括线性回归、多项式拟合、非线性回归等.假设检验、分布的检验、非参数检验.判别分析、聚类分析、因子分析.统计过程控制.提供了常用的过程管理图和过程性能图.统计图.箱形图、误差条图、q-q 图、交互插值等值线图和威布尔图等.试验设计.包括完全析因设计、不完全析因设计和 D-优化设计.31.3 本章小结本章主要探讨了统计问题的背景及应用，并概述了统计思想在各个领域的广泛应用。此外，本章还简单介绍了 MATLAB 软件的特色和本文重点应用的MATLAB 统计工具箱。4第2 章 三类常见统计方法及其解决典型数学模型实例2.1 方差分析2.1.1 单因素一元方差分析2.1.1.1 Anova1 函数 MATLAB 统计工具箱中提供了 anova1 函数，用来作单因素一元方差分析，其调用格式为：（1）p = anova1（X）根据样本观测值矩阵 X 进行均衡试验的单因素一元方差分析，检验矩阵 X 的各列所对应的总体是否具有相同的均值，原假设是 X 的各列所对应的总体具有相同的均值。矩阵 X 的列数表示因素的水平数，X 的每一列对应因素的一个水平，矩阵 X 的行数表示因素的每个水平下重复试验的次数（即样本容量），所谓均衡试验是指因素的每个水平下重复试验次数相同的试验。anova1 函数的输出参数 p是检验的 p 值，对于给定的显著性水平 ，若 p，则拒绝原假设，认为 X 的各列所对应的总体具有不完全相同的均值，否则接受原假设，认为 X 的各列所对应的总体具有相同的均值。anova1 函数还生成两个图形。第 1 个图为方差分析表，它将 X 中数据的误差分成两部分：由于列值的差异导致的误差（组间差）；由于每一列数据与该列均值的差异导致的误差（组内差）。方差分析表中有 6 列：第 1 列显示误差的来源；第 2 列显示每一个误差来源的平方和（SS）；第 3 列显示与每一个误差来源相关的自由度（df）；第 4 列显示均值平方和（MS），它是误差来源平方和与自由度的比值，即SS/df；第 5 列显示 F 统计量，它是均值平方和的比值；第 6 列显示 p 值，p 值是 F 函数(fcdf)；当 F 增加时 p 值减小。第 2 个图显示 X 的每一列的箱形图。箱形图中心线上较大的差异对应于较大的 F 值和较小的 p 值。5（2）p = anova1（X，group） 当 X 是一个矩阵时，这种调用只适合于均衡试验，anova1 函数把 X 的每一列作为一个独立的组，检验各组多对应总体是否具有相同的均值。输入参数 group可以是字符数组或者字符串元胞数组，用来指定每组的组名，X 的每一列对应一个组名字符串，在箱线图中，组名字符串被作为箱线图的标签。如果不需要指定组名，可以输入空数组（）或者忽略 group 这个输入。当 X 是一个向量，这种调用不仅适用于均衡试验，还适用于非均衡试验。anova1 函数对 X 中的样本进行单因素方差分析，通过输入变量 group 来标识 X中的每个元素的水平，所以，group 与 X 的长度必须相等。group 中包含的标签同样用于箱形图的标注。anova1 函数的矢量输入形式不需要每个样本中的观测值个数相同，所以它适用于不平衡数据。（3）p = anova1（X，group，displayopt）当displayopt参数设置为on(默认设置)时，激活 ANOVA 表和箱形图的显示；displayopt参数设置为off时，不予显示。（4）p,table = anova1（.）返回单元数组表中的 ANOVA 表（包含列标签和行标签）。（使用“Edit”菜单中的“Copy Text”选项可以将 ANOVA 表以文本形式复制到记事本中。）（5）p,table,table = anova1（.）返回 stats 结构，用于进行多重比较检验。anova1 检验评价所有样本均值相等的零假设和均值不等的备择假设。有时进行检验，决定哪对均值差异显著，哪对均值差异不显著是很有效的。提供 stats 结构作为输入，使用 multcompare 函数可以进行此项检验。注意注意：方差分析要求样本数据满足下面的假设条件：所有样本数据满足正态分布条件；所有样本数据具有相等的方差；所有观测值相互独立。在基本满足前二个假设条件的情况下，一般认为 ANOVA 检验是稳健的。2.1.1.2 Multcompare 函数 MTALAB 统计工具箱中提供了 Multcompare 函数，用来作多重比较，其调用格式如下：（1）c = multcompare（stats） 根据结构体变量 stats 中的信息进行多重比较，返回两两比较的结果矩阵c。c 是一个 5 列的矩阵，它的每一行对应一次两两比较的检验，每一行上的元素包括作比较的两个组的组标号、两个组的均值差、均值差的置信区间。例如 c 的某行元素为：62.0000 5.0000 1.9442 8.2206 14.4971表示对第二组和第五组进行两两比较的检验，两组的均值（即第二组的均值减去第五组的均值）为 8.2206，均值差的 95%置信区间为1.9442,14.4971,这个区间不包含 0，说明在显著性 0.05 下，两组间均值的差异是显著的。（2）c = multcompare（stats，param1，val1，param2，val2，.）指定一个或多个成对出现的参数名与参数值来控制多重比较。（3）c,m = multcompare（.）还返回一个多行 2 列的矩阵 m，第一列为每一组组均值的估计值，第二列为相应的标准误差。（4）c,m,h = multcompare（.）还返回交互式多重比较的图形的句柄值 h，可通过 h 修改图形属性，如图形标题和 X 轴标等。（5）c,m,h,games = multcompare（.）还返回组名变量 gnames，它是一个元胞数组，每一行对应一个组名。2.1.1.3 应用实例例例 2-12-1：将四种工艺下生产的灯泡进行寿命测试，得到数据：表 2-1 灯泡寿命测试表 试验 工艺1A2A3A4A116201580146015002167016001540155031700164016201610417501720168051800试检验工艺对寿命有无显著影响。(=0.05)解答：程序：clear;X=1620,1670,1700,1750,1800,1580,1600,1640,1720,1460,1540,1620,1500,1550,1610,1680;group=1,1,1,1,1,2,2,2,2,3,3,3,4,4,4,4;p,tab,stats=anova1(X,group,on)运行结果：p = 0.03317tab = Source SS df MS F ProbF Groups 62820 3 20940 4.0608 0.0331 Error 61880 12 5.1567e+003 Total 124700 15 stats = gnames: 4x1 cell n: 5 4 3 4 source: anova1 means: 1708 1635 1540 1585 df: 12 s: 71.8099图 2-1 灯泡寿命测试运行结果图图 2-2 灯泡寿命测试水箱平衡图8例例 2-22-2：研究 6 种农药对杀虫效果的影响，试验所得数据如下表：表 2-2 农药杀虫量表农药编号杀虫量12345687.485.080.290.588.587.394.356.262.455.048.292.099.295.391.575.272.381.3设该试验的线性统计模型为：iijijiijnjiNdi iy, 1;6 , 1. ), 0(. .,2诸其中，。3, 4, 2, 2, 4, 3654321nnnnnn解答：程序：y1=87.4 85.0 80.2;y2=90.5 88.5 87.3 94.3;y3=56.2 62.4;y4=55.0 48.2;y5=92.0 99.2 95.3 91.5;y6=75.2 72.3 81.3;y=y1 y2 y3 y4 y5 y6;A1=ones(numel(y1),1);A2=2*ones(numel(y2),1);A3=3*ones(numel(y3),1);A4=4*ones(numel(y4),1);A5=5*ones(numel(y5),1);A6=6*ones(numel(y6),1);A=A1;A2;A3;A4;A5;A6;p,table,stats=anova1(y,A)B=multcompare(stats,0.05)运行结果：p = 1.0671e-0079table = Source SS df MS F ProbF Groups 3.8177e+003 5 763.5487 51.5921 1.0671e-007 Error 177.5967 12 14.7997 Total 3.9953e+003 17 stats = gnames: 6x1 cell n: 3 4 2 2 4 3 source: anova1 means: 84.2000 90.1500 59.3000 51.6000 94.5000 76.2667 df: 12 s: 3.8470B = 1.0000 2.0000 -15.8193 -5.9500 3.9193 1.0000 3.0000 13.1040 24.9000 36.6960 1.0000 4.0000 20.8040 32.6000 44.3960 1.0000 5.0000 -20.1693 -10.3000 -0.4307 1.0000 6.0000 -2.6174 7.9333 18.4840 2.0000 3.0000 19.6593 30.8500 42.0407 2.0000 4.0000 27.3593 38.5500 49.7407 2.0000 5.0000 -13.4872 -4.3500 4.7872 2.0000 6.0000 4.0141 13.8833 23.7526 3.0000 4.0000 -5.2219 7.7000 20.6219 3.0000 5.0000 -46.3907 -35.2000 -24.0093 3.0000 6.0000 -28.7627 -16.9667 -5.1706 4.0000 5.0000 -54.0907 -42.9000 -31.7093 4.0000 6.0000 -36.4627 -24.6667 -12.8706 5.0000 6.0000 8.3641 18.2333 28.1026：图 2-3 运行结果图10图 2-4 分布示意图2.1.2 双因素一元方差分析2.1.2.1 Anova2 函数MATLAB 统计工具箱中提供了 anova2 函数，用来作双因素一元方差分析，其调用格式如下：（1）p=anova2(X,reps) 进行平衡双因子方差分析，不同列中的数据代表一个因子 A 的变化。不同行中的数据代表另一因子 B 的变化。若在每一个行列匹配点上有一个以上的观测值，则变量 reps 指示每一个单元中观测值的个数。当 reps=1 (默认值可省略)时，anova2 函数返回两个 p 值到 p 矢量中：零假设 H0A 的 p 值。零假设为源于因子 A 的所有样本（如 X 中的所有列样本）取自相同的总体。零假设 H0B 的 p 值。零假设为源于因子 B 的所有样本（如 X 中的所有行样本）取自相同的总体。当 reps1 时，anova2 在矢量中返回第 3 个值：p零假设 H0AB 的 p 值。零假设为因子 A 和因子 B 之间没有交互效应。如果任意一个 p 值接近于 0，则认为相关的零假设不成立。对于零假设 H0A，一个足够小的 p 值表示至少有一个列样本均值明显地不同于其他列样本均值，即因子 A 存在主效应；对于零假设 H0B，一个足够小的 p 值表示至少有一个行样本均值明显地不同于其他行样本均值，即因子 B 存在主效应；对于零假设 H0AB，一个足够小的 p 值表示因子 A 与因子 B 之间存在交互效应。11为了决定结果是否是“统计上显著的”，需要确定 p 值。一般地，当 p 值小于 0.05 或 0.01 时，认为结果是显著的。anova2 函数还显示一个含方差分析表的图形。它将 X 中数据的误差根据reps 的值分为 3 部分或 4 部分：由于列均值差异引起的误差；由于行均值差异引起的误差；由于行列交互作用引起的误差（如果 reps 大于它的默认值 1）；剩下的误差为不能被任何系统因素解释的误差。该方差分析表中包含 6 列：第 1 列显示误差的来源；第 2 列显示来源于每一个误差来源的平方和（SS）；第 3 列为与每一个误差来源相关的自由度（df）；第 4 列为均值平方和（MS），它是误差平方和与自由度的比值，即 SS/df；第 5 列为 F 统计量，它是均值平方和的比值；第 6 列为 p 值，它是 F 函数(fcdf)；当 F 增加时 p 值减小。（2）p=anova2(X,group,displayopt) 当displayopt参数设置为on（默认设置）时，激活 ANOVA 表和箱形图的显示；displayopt参数设置为off时，不予显示。（3）p,table=anova2() 返回单元数组表中的 ANOVA 表（包含列标签和行标签）。（使用 Edit 菜单中的Copy Text 选项可以将 ANOVA 表以文本形式复制到剪贴板中。）（4）p,table,stats=anova2() 返回 stats 结构，用于进行列因子均值的多重比较检验。2.1.2.2 应用实例例例 2-32-3：为了考察某种电池的最大输出电压受板极材料与使用电池的环境温度的影响，材料类型(因子)取 3 个水平(即 3 种不同的材料)，温度(因子)也AB取 3 个水平，每个水平组合下重复 4 次试验，数据如下：表 2-3 电池最大输出电压受板极材料与环境温度的影响试验数据表温度B1525351130 155 174 18034 40 80 7520 70 82 582150 188 159 126136 122 106 11525 70 58 45A材料类型3138 110 168 160174 120 150 13996 104 82 60解答：程序：12%两因子等重复试验（方差分析）y1=130 155 174 18034 40 80 7520 70 82 58;y2=150 188 159 126136 122 106 11525 70 58 45;y3=138 110 168 160174 120 150 13996 104 82 60;y=y1 y2 y3;p=anova2(y,4)运行结果：p =0.0043 0.0000 0.0008图 2-5 运行结果图2.2 判别分析 判别分析是用于判别个体所属群体的一种统计方法，判别分析的特点是根据已掌握的、历史上每个类别的若干样本的数据信息，总结出客观事物分类的规律性，建立判别公式和判别准则。然后，当遇到新的样本点时，只要根据总结出来的判别公式和判别准则，就能判别该样本点所属的类别。判别分析是一种应用性很强的统计数据分析方法。判别分析的原理:为了能识别待判断的对象 x =是属于已知类中的哪Tm21) x, x,(x rAAA.,21一类？事先必须要有一个一般规则, 一旦知道了 x 的值, 便能根据这个规则立即作出判断, 称这样的一个规则为判别规则.判别规则往往通过的某个函数来表达, 把它称为判别函数, 记作 W(i; x).一旦知道了判别函数并确定了判别规则，最好将已知类别的对象代入检验，这一过程称为回代检验，以便检验你的判别函数和判别规则是否正确。132.2.1 相关函数2.2.1.1 classify 函数 用 classify 函数进行线性判别分析，其调用格式如下。（1）class=classify(sample,training,group)将 sample 数据中的每一行分到 training 参数中的某一个类中。sample 和training 必须为列数相同的矩阵。group 参数为 training 参数的分类变量。它的唯一值定义不同的组，每个元素定义 training 的对应行属于哪一类。group 可以是数值矢量、字符串数组成或字符串单元数组。group 和 training 必须具有相同的行数。classify 函数将 group 中的 NaN 或空字符串视为缺失值，并忽略training 的对应行。该函数返回 class，它是一个与 sample 具有相同行数的矢量。class 的每个元素指定 sample 中每行元素的类别。（2）class,err=classify(.)返回错判误差率的估计。该函数返回 trainin 中被错判的观测量的百分比。2.2.1.2 mahal 函数 用 mahal 函数计算马氏距离，其调用格式如下。（1）mahal(Y,X)计算 X 矩阵中的样本至 Y 矩阵中每个点（行）的马氏距离。Y 的列数必须等于 X 的列数，但它们的行数可以不同。X 的行数必须大于列数。马氏距离是一个将数据集与空间点分隔开来的多变量度量。在线性判别分析中，使它最小化。2.2.2 距离判别法 距离判别方法：判别函数 W(i; x) = d (x, ), 其中 d (x, )为待判iAiA别对象 x 到第 i 类的距离. 判别规则为若 W(k; x) = minW(i; x)| i =1, 2, iA , r , 则 x . kA2.2.3 贝叶斯判别法 判别函数 W(i; x ) = 其中为待判别对象 x的概率,如果没有）（xqpiiipiA任何这种附加的先验信息,通常取 pi = 1/r。为已知类别的分布密度判别)(xqiiA规则为若 W(k; x ) = maxW(i; x )| i =1, 2, , r ,则 x. kA142.2.4 Fisher 判别法 Fisher 判别法是一种线性判别的方法。它的工作思路是对原数据系统进行坐标变换，寻求能将总体尽可能分开的方向。具体的作法是先引入一个与样本同维的待定向量 u，再将 y 取为 x 坐标的线性组合 y=而 u 的选取，要使同一类别xuT产生的 y 尽量聚拢，不同类别产生的 y 尽量拉开。 这样，我们便可将样品 x 到某一类 G 的距离定义为 y=与 xuTcuycT 之间的欧氏距离：| )(|),(cxuycyGXLT 其中 c 为 G 的几何中心Fisher 分类的判据为： 1 若 ，则判定义为 A 类；),(),(BXLAXL 2 若 ，则判定又为 B 类；),(),(BXLAXL 3 若 ，则判定又为不可判类。),(),(BXLAXL2.2.5 应用实例例例 2-42-4：某地大气样品污染分类如下：表 2-4 大气样品污染分类表气体氯硫化氢二氧化硫碳 4环氧氯丙烷环己烷污染分类10.0560.0840.0310.0380.00810.022120.0400.0550.1000.1100.02200.0073130.0500.0740.0410.0480.00710.020140.0450.0590.1100.1000.02500.0063150.0380.1300.0790.1700.05800.043260.0300.1100.0700.1600.05000.046270.0340.0950.0580.1600.2000.029180.0300.0900.0680.1800.2200.039190.0840.0660.0290.3200.0120.0412100.0850.0760.0190.3000.0100.0402110.0640.0720.0200.2500.0280.0382120.0540.0650.0220.2800.0210.0402130.0480.0890.0620.2600.0380.036215140.0450.0920.0720.2000.0350.0322150.0690.0870.0270.0500.0890.0211在此地某大型化工厂的厂区及邻近地区挑选 4 个有代表性的大气样本取样，获得数据如下：表 2-5 大气样本取样污染分类表气体氯硫化氢二氧化硫碳 4环氧氯丙烷环己烷污染分类样品 10.0520.0840.0210.0370.00710.022样品 20.0410.0550.1100.1100.02100.0073样品 30.0300.1120.0720.1600.0560.021样品 40.0740.0830.1050.1900.0201.000求它们的污染分类。解答：程序：training=0.056 0.084 0.031 0.038 0.0081 0.022; 0.040 0.055 0.100 0.110 0.0220 0.0073; 0.050 0.074 0.041 0.048 0.0071 0.020; 0.045 0.050 0.110 0.100 0.0250 0.0063; 0.038 0.130 0.079 0.170 0.0580 0.043; 0.030 0.110 0.070 0.160 0.0500 0.046; 0.034 0.095 0.058 0.160 0.200 0.029; 0.030 0.090 0.068 0.180 0.220 0.039; 0.084 0.066 0.029 0.320 0.012 0.041; 0.085 0.076 0.019 0.300 0.010 0.040; 0.064 0.072 0.020 0.250 0.028 0.038; 0.054 0.065 0.022 0.280 0.021 0.040; 0.048 0.089 0.062 0.260 0.038 0.036; 0.045 0.092 0.072 0.200 0.035 0.032; 0.069 0.087 0.027 0.050 0.089 0.021; group=1;1 ;1 ;1 ;2 ;2 ;1; 1; 2 ;2 ;2 ;2 ;2 ;2 ;1; sample=0.052 0.084 0.021 0.037 0.0071 0.022; 0.041 0.055 0.110 0.110 0.0210 0.0073; 0.030 0.112 0.072 0.160 0.056 0.021; 0.074 0.083 0.105 0.190 0.020 1.000; 16class,err=classify(sample,training,group,mahalanobis)运行结果：class = 1 1 2 2err = 02.3 主成分分析2.3.1 主成分分析方法原理主成分分析是把原来多个变量化为少数几个综合指标的一种统计分析方法，从数学角度来看，这是一种降维处理技术。假定有 n 个地理样本，每个样本共有p 个变量描述，这样就构成了一个 np 阶的地理数据矩阵： （21）如何从这么多变量的数据中抓住地理事物的内在规律性呢？要解决这一问题，自然要在 p 维空间中加以考察，这是比较麻烦的。为了克服这一困难，就需要进行降维处理，即用较少的几个综合指标来代替原来较多的变量指标，而且使这些较少的综合指标既能尽量多地反映原来较多指标所反映的信息，同时它们之间又是彼此独立的。那么，这些综合指标（即新变量)应如何选取呢？显然，其最简单的形式就是取原来变量指标的线性组合，适当调整组合系数，使新的变量指标之间相互独立且代表性最好。17如果记原来的变量指标为，它们的综合指标新变量指标为pxxx21,，（mp)。则21,xxmz （22）111 11221221 122221 122ppppmmmmppxl xl xl xxl xl xlxzl xlxlx在（2)式中，系数 lij由下列原则来决定：（1)zi与 zj（ij；i，j=1，2，m)相互无关；（2)z1是 x1，x2，xp的一切线性组合中方差最大者；z2是与 z1不相关的x1，x2，xp的所有线性组合中方差最大者；zm是与 z1，z2，zm-1都不相关的 x1，x2，xp的所有线性组合中方差最大者。这样决定的新变量指标 z1，z2，分别称为原变量指标 x1，x2，xpmz的第一，第二，第 m 主成分。其中，z1在总方差中占的比例最大，z2，z3，zm的方差依次递减。在实际问题的分析中，常挑选前几个最大的主成分，这样既减少了变量的数目，又抓住了主要矛盾，简化了变量之间的关系。从以上分析可以看出，找主成分就是确定原来变量 xj（j=1，2，p)在诸主成分 zi（i=1，2，m)上的载荷 lij（i=1，2，m；j=1，2，p)，从数学上容易知道，它们分别是 x1，x2，xp的相关矩阵的 m 个较大的特征值所对应的特征向量。2.3.2 主成分分析的计算步骤（1)计算相关系数矩阵 （23）111212122212PPPPPPrrrrrrRMMMMrrr18 在公式（3)中，（i，j=1，2，p)为原来变量与的相关系数，其ijrixjx计算公式为： （24）12211()()()()nkiikjjkijnnkiikjjkkxxxxrxxxx因为 R 是实对称矩阵（即 rij=rji)，所以只需计算其上三角元素或下三角元素即可。（2)计算特征值与特征向量首先解特征方程I-R=0 求出特征值 i（i=1，2，p)，并使其按大小顺序排列，即 12，p0；然后分别求出对应于特征值 i的特征向量 ei（i=1，2，p)。（3)计算主成分贡献率及累计贡献率。pkkmkkpkkiipirz111), 2 , 1(/累计贡献率：贡献率：主成分一般取累计贡献率达 85-95的特征值，所对应的第一，第二，m,21，第 m（mp)个主成分。（4)计算主成分载荷 （2-5）), 2 , 1,(),(pkiexzpkikik由此可以进一步计算主成分得分：Z= （2-6）nmnnmmzzzMMMMzzzzzz212222111211192.3.3 主成分分析相关函数调用 MATLAB 工具箱 princomp 函数进行主成分分析：格式： PC = princomp（X） PC，SCORE，latent，tsquare = princomp（X）说明： PC,SCORE,latent,tsquare=princomp(X)对数据矩阵 X 进行主成分分析，给出各主成分 (PC)、所谓的 Z-得分(SCORE)、X 的方差矩阵的特征值(latent)和每个数据点的 Hotelling统计量(tsquare）。2T2.3.4 应用实例例例 2-52-5：城市环境生态化是城市发展的必然趁势，表现为社会、经济、环境与生态全方位的现代化水平，一个符合生态规律的生态城市应该是结构合理、功能高效和关系协调的城市生态系统。所谓结构合理是指适度的人口密度，合理的土地利用，良好的环境质量，充足的绿地系统，完善的基础设施，有效的自然保护；功能高效是指资源的优化配置、物力的经济投入、人力的充分发挥、物流的畅通有序、信息流的快捷；关系协调是指人和自然协调、社会关系协调、城乡协调、资源利用和更新协调。一个城市要实现生态城市的发展目标，关键是在市场经济的体制下逐步改善城市的生态环境质量，防止生态环境质量恶化，因此，对城市的生态环境水平调查评价很有必要。我们对江苏省十个城市的生态环境状况进行了调查，得到生态环境指标的指数值。表 2-6 生态环境水平分析和评价表20解答：利用 Matlab 中的 princomp 命令实现。具体程序如下：x = 0.7883 0 .7391 0.8111 0.6587 0.6543 0.8259 0.8486 0.6834 0.8495 0.7846 0.7633 0.7287 0.7629 0.8552 0.7564 0.7455 0.7800 0.9490 0.8918 0.8954 0.4745 0.5126 0.8810 0.8903 0.8288 0.7850 0.8032 0.8862 0.3987 0.3970 0.8246 0.7603 0.6888 0.8977 0.7926 0.7856 0.6509 0.8902 0.6799 0.9877 0.8791 0.8736 0.8183 0.9446 0.9202 0.9263 0.9185 0.9505 0.8620 0.8873 0.9538 0.9257 0.9285 0.9434 0.9154 0.8871 0.9357 0.8760 0.9579 0.9741 0.8785 0.8542 0.8537 0.9027 0.8729 0.8485 0.8473 0.9044 0.8866 0.9035 0.6305 0.6187 0.6313 0.7415 0.6398 0.6142 0.5734 0.8980 0.6186 0.7382 0.8928 0.7831 0.5608 0.8419 0.8464 0.7616 0.8234 0.6384 0.9604 0.8514x = ;xstdr = std(x); %求各变量标准差n,m = size(x);sddata = x./stdr(ones(n,1),:); %标准化变换p,princ,egenvalue = princomp(sddata) %调用主成分分析程序p3 = p(:,1:3) %输出前三个主成分系数sc = princ(:,1:3) %输出前三个主成分得分egenvalue %输出特征根per = 100*egenvalue/sum(egenvalue) %输出各个主成分贡献率执行后得到所要结果，这里是前三个主成分，主成分得分，特征根，即： 5933. 04589. 00643. 02804. 00267. 04815. 00958. 03116. 04236. 00270. 05273. 01318. 05664. 01601. 03355. 00582. 01812. 04048. 00498. 05299. 01364. 03535. 02313. 03702. 03282. 01442. 03677. 0p6469. 07034. 23736. 14498. 03442. 22134. 19765. 05024. 19332. 35432. 06907. 08984. 15837. 06636. 19266. 04326. 19002. 03634. 09043. 01623. 03754. 20552. 26267. 18408. 12729. 02159. 03364. 13946. 03897. 18301. 0princegenvalue = , per = 0597. 1 ,6407. 2 ,8811. 397.11,34.29,12.4321这样，前三个主成分为:98765432110643. 04815. 04236. 01318. 03355. 04048. 01364. 03702. 03677. 0 xxxxxxxxxz98765432124589. 00267. 03116. 05273. 01601. 01812. 05299. 02313. 01442. 0 xxxxxxxxxz98765432135933. 02804. 00958. 00270. 05664. 00582. 00498. 03535. 03282. 0 xxxxxxxxxz 第一主成分贡献率为 43.12，第二主成分贡献率为 29.34，第三主成分贡献率为 11.97，前三个主成分累计贡献率达 84.24。如果按 80 以上的信息量选取新因子，则可以选取前三个新因子。第一新因子 z1 包含的信息量最大为 43.12，它的主要代表变量为 X8(城市文明)、X7(生产效率)、X6 (城市绿化)，其权重系数分别为 0.4815、0.4236、0.4048，反映了这三个变量与生态环境水平密切相关，第二新因子 Z2包含的信息量次之为29.34，它的主要代表变量为 X3(地理结构)、X6(资源配置)、X9(可持续性)，其权重系数分别为 0.5299、0.5273、0.4589，第三新因子 Z3包含的信息量为11.97，代表总量为 X9(可持续性)、 X5(物质还原)，权重系数分别为 0.5933、0.5664。这些代表变量反映了各自对该新因子作用的大小，它们是生态环境系统中最重要的影响因素。根据前三个主成分得分，用其贡献率加权，即得十个城市各自的总得分:3077. 1 ,1107. 0 ,1383. 1 ,9562. 0,8178. 0,0640. 0 ,1801. 1 ,5170. 1,6069. 0,0970. 0)3(:,%97.11)2(:,%34.29) 1