数学建模方法大全

数学建模方法大全

二○一二年九月九日星期日

9时59分32秒

目录

一、主成分分析法 ...................................................................................................................... 3 二、因子分析法 .......................................................................................................................... 6 三、聚类分析 ............................................................................................................................ 10 四、最小二乘法与多项式拟合 ............................................................................................. 17 五、回归分析（略） ............................................................................................................... 23 六、概率分布方法（略） ...................................................................................................... 23 七、插值与拟合（略） ........................................................................................................... 23 八、方差分析法 ....................................................................................................................... 24 九、逼近理想点排序法 ........................................................................................................... 29 十、动态加权法 ........................................................................................................................ 30 十一、灰色关联分析法 ........................................................................................................... 32 十二、灰色预测法 .................................................................................................................... 34 十三、模糊综合评价 ............................................................................................................... 36 十四、隶属函数的刻画（略） ............................................................................................. 38 十五、时间序列分析法 ........................................................................................................... 39 十六、蒙特卡罗(MC)仿真模型 ............................................................................................. 43 十七、BP神经网络方法 .......................................................................................................... 45 十八、数据包络分析法（DEA） ........................................................................................... 52 十九、多因素方差分析法（）基于SPSS） ..................................................................... 55 二十、拉格朗日插值 ........................................................................................................... 71 二十一、回归分析（略） ...................................................................................................... 76 二十二、概率分布方法（略） ............................................................................................. 76 二十三、插值与拟合（略） .................................................................................................. 76 二十四、隶属函数的刻画（参考《数学建模及其方法应用》） ............................... 76 二十五、0-1整数规划模型（参看书籍） ....................................................................... 76 二十六、Board评价法（略） .............................................................................................. 76 二十七、纳什均衡（参看书籍） ........................................................................................ 76 二十八、微分方程方法与差分方程方法（参看书籍） ............................................... 76 二十九、莱斯利离散人口模型（参看数据） .................................................................. 76 三十、一次指数平滑预测法（主要是软件的使用） .................................................... 76 三十一、二次曲线回归方程（主要是软件的使用） .................................................... 76 三十二、成本-效用分析（略） ........................................................................................... 76 三十三、逐步回归法（主要是软件的使用） .................................................................. 76 三十四、双因子方差分析（略） ........................................................................................ 76

一、主成分分析法 一）、主成分分析法介绍：

主成分分析（principal components analysis，PCA）又称：主分量分析，主成分回归分析法。旨在利用降维的思想，把多指标转化为少数几个综合指标。它是一个线性变换。这个变换把数据变换到一个新的坐标系统中，使得任何数据投影的第一大方差在第一个坐标(称为第一主成分)上，第二大方差在第二个坐标(第二主成分)上，依次类推。主成分分析经常用减少数据集的维数，同时保持数据集的对方差贡献最大的特征。这是通过保留低阶主成分，忽略高阶主成分做到的。这样低阶成分往往能够保留住数据的最重要方面。但是，这也不是一定的，要视具体应用而定。

二）、主成分分析法的基本思想：

在实证问题研究中，为了全面、系统地分析问题，我们必须考虑众多影响因素。这些涉及的因素一般称为指标，在多元统计分析中也称为变量。因为每个变量都在不同程度上反映了所研究问题的某些信息，并且指标之间彼此有一定的相关性，因而所得的统计数据反映的信息在一定程度上有重叠。在用统计方法研究多变量问题时，变量太 多会增加计算量和增加分析问题的复杂性，人们希望在进行定量分析的过程中，涉及的变量较少，得到的信息量较多。主成分分析正是适应这一要求产生的，是解决这类题的理想工具。

同样，在科普效果评估的过程中也存在着这样的问题。科普效果是很难具体量化的。在实际评估工作中，我们常常会选用几个有代表性的综合指标，采用打分的方法来进行评估，故综合指标的选取是个重点和难点。如上所述，主成分分析法正是解决这一问题的理想工具。因为评估所涉及的众多变量之间既然有一定的相关性，就必然存在着起支配作用的因素。根据这一点，通过对原始变量相关矩阵内部结构 的关系研究，找出影响科普效果某一要素的几个综合指标，使综合指标为原来变量的线 性拟合。这样，综合指标不仅保留了原始变量的主要信息，且彼此间不相关，又比原始 变量具有某些更优越的性质，就使我们在研究复杂的科普效果评估问题时，容易抓住主 要矛盾。 上述想法可进一步概述为：设某科普效果评估要素涉及个指标，这指标构 成的维随机向量为。对作正交变换，令，其中为正交阵，的各分量是不相关的，使得的各分量在某个评估要素中的作用容易解释，这就使得我们有可能从主分量中选择主要成分，削除对这一要素影响微弱的部分，通过 对主分量的重点分析，达到对原始变量进行分析的目的。的各分量是原始变量线性组合，不同的分量表示原始变量之间不同的影响关系。由于这些基本关系很可能与特定的作用过程相联系，主成分分析使我们能从错综复杂的科普评估要素的众多指标中，找出一些主要成分，以便有效地利用大量统计数据，进行科普效果评估分析，使我们在研究科普效果评估问题中，可能得到深层次的一些启发，把科普效果评估研究引向深入。

例如，在对科普产品开发和利用这一要素的评估中，涉及科普创作人数百万人、科 普作品发行量百万人、科普产业化（科普示范基地数百万人）等多项指标。经过主成分分析计算，最后确定个或个主成分作为综合评价科普产品利用和开发的综合指标，变量数减少，并达到一定的可信度，就容易进行科普效果的评估。

三）、主成分分析法的数学模型：

其中：

为第 j个指标对应于第 个主成分的初始因子载荷， 为第 l个主成分对应的特征值 根据主成分表达式得出综合得分模型：

四）、主成分分析法的基本原理：

主成分分析法是一种降维的统计方法，它借助于一个正交变换，将其分量相关的原随机向量转化成其分量不相关的新随机向量，这在代数上表现为将原随机向量的协方差阵变换成对角形阵，在几何上表现为将原坐标系变换成新的正交坐标系，使之指向样本点散布最开的p 个正交方向，然后对多维变量系统进行降维处理，使之能以一个较高的精度转换成低维变量系统，再通过构造适当的价值函数，进一步把低维系统转化成一维系统。 五）、主成分分析法的作用：

概括起来说，主成分分析主要由以下几个方面的作用。

1．主成分分析能降低所研究的数据空间的维数。即用研究m维的Y空间代替p维的X空间(m＜p)，而低维的Y空间代替 高维的x空间所损失的信息很少。即：使只有一个主成分Yl(即 m＝1)时，这个Yl仍是使用全部X变量(p个)得到的。例如要计算Yl的均值也得使用全部x的均值。在所选的前m个主成分中，如果某个Xi的系数全部近似于零的话，就可以把这个Xi删除，这也是一种删除多余变量的方法。

2．有时可通过因子负荷aij的结论，弄清X变量间的某些关系。

3．多维数据的一种图形表示方法。我们知道当维数大于3时便不能画出几何图形，多元统计研究的问题大都多于3个变量。要把研究的问题用图形表示出来是不可能的。然而，经过主成分分析后，我们可以选取前两个主成分或其中某两个主成分，根据主成分的得分，画出n个样品在二维平面上的分布况，由图形可直观地看出各样品在主分量中的地位，进而还可以对样本进行分类处理，可以由图形发现远离大多数样本点的离群点。

4．由主成分分析法构造回归模型。即把各主成分作为新自变量代替原来自变量x做回归分析。

5．用主成分分析筛选回归变量。回归变量的选择有着重的实际意义，为了使模型本身易于做结构分析、控制和预报，好从原始变量所构成的子集合中选择最佳变量，构成最佳变量集合。用主成分分析筛选变量，可以用较少的计算量来选择量，获得选择最佳变量子集合的效果。 六）、主成分分析法的计算步骤：

1、原始指标数据的标准化采集p 维随机向量x = (x1,X2,...,Up)T)n 个样品xi = (xi1,xi2,...,dip)T ，I=1,2,…,n，

n＞p，构造样本阵，对样本阵元进行如下标准化变换：

其中

2、对标准化阵Z 求相关系数矩阵

其中,

3、解样本相关矩阵R 的特征方程按

。

得p 个特征根,确定主成分 ，得标准化阵Z。

确定m 值，使信息的利用率达85%以上，对每个job,

j=1,2,...,m, 解方程组Rib = job得单位特征向量 。 4、将标准化后的指标变量转换为主成分

U1称为第一主成分,U2 称为第二主成分,…,Up 称为第p 主成分。 5 、对m 个主成分进行综合评价

对m 个主成分进行加权求和，即得最终评价值，权数为每个主成分的方差贡献率。

PS另一种易于理解的步骤： 1、数据标准化;

2、求相关系数矩阵;

3、一系列正交变换，使非对角线上的数置0，加到主对角上; 得特征根xi（即相应那个主成分引起变异的方差),并按照从大到小的顺序把特征根排列;

4、求各个特征根对应的特征向量; 用下式计算每个特征根的贡献率Vi; VI=xi/(x1+x2+........)

5、根据特征根及其特征向量解释主成分物理意义 七）、主成分分析法的案例： 参见：基于主成分分析的力量结构指标的权重的计算、基于主成分析的江苏省地方高校创新力研究

二、因子分析法

一）因子分析法介绍：

主成分分析通过线性组合将原变量综合成几个主成分，用较少的综合指标来代替原来较多的指标(变量)。在多变量分析中，某些变量间往往存在相关性。是什么原因使变量间有关联呢？是否存在不能直接观测到的、但影响可观测变量变化的公共因子？因子分析法(Factor Analysis)就是寻找这些公共因子的模型分析方法，它是在主成分的基础上构筑若干意义较为明确的公因子，以它们为框架分解原变量，以此考察原变量间的联系与区别。

例：随着年龄的增长，儿童的身高、体重会随着变化，具有一定的相关性，身高和体重之间为何会有相关性呢？因为存在着一个同时支配或影响着身高与体重的生长因子。那么，我们能否通过对多个变量的相关系数矩阵的研究，找出同时影响或支配所有变量的共性因子呢？因子分析就是从大量的数据中“由表及里”、“去粗取精”，寻找影响或支配变量的多变量统计方法。 因此，可以说因子分析是主成分分析的推广，也是一种把多个变量化为少数几个综合变量的多变量分析方法，其目的是用有限个不可观测的隐变量来解释原始变量之间的相关关系。

因子分析主要用于：1、减少分析变量个数；2、通过对变量间相关关系探测，将

原始变量进行分类。即将相关性高的变量分为一组，用共性因子代替该组变量。 二）、因子分析法的基本模型：

因子分析法是从研究变量内部相关的依赖关系出发，把一些具有错综复杂关系的变量归结为少数几个综合因子的一种多变量统计分析方法。它的基本思想是将观测变量进行分类，将相关性较高，即联系比较紧密的分在同一类中，而不同类变量之间的相关性则较低，那么每一类变量实际上就代表了一个基本结构，即公共因子。对于所研究的问题就是试图用最少个数的不可测的所谓公共因子的线性函数与特殊因子之和来描述原来观测的每一分量。 因子分析模型描述如下：

1、X=(x1，x2，…，xp)是可观测随机向量，均值向量E(X)=0，协方差阵Cov(X)=∑，且协方差阵∑与相关矩阵R相等（只要将变量标准化即可实现）。 2、F=(F1，F2，…，Fm)（m3、e=(e1，e2，…，ep)与F相互独立，且E(e)=0，e的协方差阵∑是对角阵，即各分量e之间是相互独立的，则模型：

x1=a11F1+a12F2+…+a1mFm+e1 x2=a21F1+a22F2+…+a2mFm+e2 xp=ap1F1+ap2F2+…+apmFm+ep

称为因子分析模型，由于该模型是针对变量进行的，各因子又是正交的，所以也称为R型正交因子模型。其矩阵形式为：

x=AF+e

其中：

x=，A=，F=，e=

这里

（1）m￡p；

（2）Cov(F，e)=0，即F和e是不相关的；

（3）D(F)=Im，即F1，F2，…，Fm不相关且方差均为1； （4）D(e)=，即e1，e2，…，ep不相关，且方差不同。

我们把F称为X的公共因子或潜因子，矩阵A称为因子载荷矩阵，e称为X的特殊因子。

A=(aij)，aij为因子载荷。数学上可以证明，因子载荷aij就是第i变量与第j因子的相关系数，反映了第i变量在第j因子上的重要性。 三）、模型的统计意义：

模型中F1，F2，…，Fm叫做主因子或公共因子，它们是在各个原观测变量的表达式中都共同出现的因子，是相互独立的不可观测的理论变量。公共因子的含义，必须结合具体问题的实际意义而定。e1，e2，…，ep叫做特殊因子，是向量x的分量xi(i=1，2，…，p)所特有的因子，各特殊因子之间以及特殊因子与所有公共因子之间都是相互独立的。模型中载荷矩阵A中的元素(aij)是为因子载荷。因子载荷aij是xi与Fj的协方差，也是xi与Fj的相关系数，它表示xi依赖Fj的程度。可将aij看作第i个变量在第j公共因子上的权，aij的绝对值越大(|aij|￡1)，表明xi与Fj的相依程度越大，或称公共因子Fj对于xi的载荷量越大。为了得到因子分析结果的经济解释，因子载荷矩阵A中有两个统计量十分重要，即变量共同度和公共因子的方差贡献。

因子载荷矩阵A中第i行元素之平方和记为hi2，称为变量xi的共同度。它是全部公共因子对xi的方差所做出的贡献，反映了全部公共因子对变量xi的影响。hi2大表明x的第i个分量xi对于F的每一分量F1，F2，…，Fm的共同依赖程度大。

将因子载荷矩阵A的第j列( j =1，2，…，m)的各元素的平方和记为gj2，称为公共因子Fj对x的方差贡献。gj2就表示第j个公共因子Fj对于x的每一分量xi(i= 1，2，…，p)所提供方差的总和，它是衡量公共因子相对重要性的指标。gj2越大，表明公共因子Fj对x的贡献越大，或者说对x的影响和作用就越大。如果将因子载荷矩阵A的所有gj2(j=1，2，…，m)都计算出来，使其按照大小排序，就可以依此提炼出最有影响力的公共因子。 四）、因子旋转：

建立因子分析模型的目的不仅是找出主因子，更重要的是知道每个主因子的意义，以便对实际问题进行分析。如果求出主因子解后，各个主因子的典型代表变量不很突出，还需要进行因子旋转，通过适当的旋转得到比较满意的主因子。

旋转的方法有很多，正交旋转(orthogonal rotation)和斜交旋转(oblique rotation)是因子旋转的两类方法。最常用的方法是最大方差正交旋转法(Varimax)。进行因子旋转，就是要使因子载荷矩阵中因子载荷的平方值向0和1两个方向分化，使大的载荷更大，小的载荷更小。因子旋转过程中，如果因子对应轴相互正交，则称为正交旋转；如果因子对应轴相互间不是正交的，则称为斜交旋转。常用的斜交旋转方法有Promax法等。 五）、因子得分：

因子分析模型建立后，还有一个重要的作用是应用因子分析模型去评价每个样品在整个模型中的地位，即进行综合评价。例如地区经济发展的因子分析

模型建立后，我们希望知道每个地区经济发展的情况，把区域经济划分归类，哪些地区发展较快，哪些中等发达，哪些较慢等。这时需要将公共因子用变量的线性组合来表示，也即由地区经济的各项指标值来估计它的因子得分。 设公共因子F由变量x表示的线性组合为：

Fj=uj1xj1+uj2xj2+…+ujpxjpj=1，2，…，m

该式称为因子得分函数，由它来计算每个样品的公共因子得分。若取m=2，则将每个样品的p个变量代入上式即可算出每个样品的因子得分F1和F2，并将其在平面上做因子得分散点图，进而对样品进行分类或对原始数据进行更深入的研究。

但因子得分函数中方程的个数m小于变量的个数p，所以并不能精确计算出因子得分，只能对因子得分进行估计。估计因子得分的方法较多，常用的有回归估计法，Bartlett估计法，Thomson估计法。 具体方法为： （1）回归估计法

F=Xb=X(X￠X)-1A￠=XR-1A￠(这里R为相关阵，且R=X￠X)。 （2）Bartlett估计法

Bartlett估计因子得分可由最小二乘法或极大似然法导出。 F=(W-1/2A)￠W-1/2A]-1(W-1/2A)￠W-1/2X=(A￠W-1A)-1A￠W-1X （3）Thomson估计法

在回归估计法中，实际上是忽略特殊因子的作用，取R = X ￠X，若考虑特殊因子的作用，此时R = X ￠X＋W，于是有： F=XR-1A￠=X(X￠X＋W)-1A￠

这就是Thomson估计的因子得分，使用矩阵求逆算法(参考线性代数文献)可以将其转换为：

F=XR-1A￠=X(I+A￠W-1A)-1W-1A￠ 六）、因子分析的步骤：

因子分析的核心问题有两个：一是如何构造因子变量；二是如何对因子变量进行命名解释。因此，因子分析的基本步骤和解决思路就是围绕这两个核心问题展开的。

因子分析常常有以下四个基本步骤：

1、确认待分析的原变量是否适合作因子分析。 2、构造因子变量。

3、利用旋转方法使因子变量更具有可解释性。 4、计算因子变量得分。 因子分析的计算过程：

1、将原始数据标准化，以消除变量间在数量级和量纲上的不同。 2、求标准化数据的相关矩阵；

3、求相关矩阵的特征值和特征向量； 4、计算方差贡献率与累积方差贡献率；

5、确定因子：设F1，F2，…，Fp为p个因子，其中前m个因子包含的数据信息总量（即其累积贡献率）不低于80%时，可取前m个因子来反映原评价指标；

6、因子旋转：

若所得的m个因子无法确定或其实际意义不是很明显，这时需将因子进行旋转以获得较为明显的实际含义。

7、用原指标的线性组合来求各因子得分：采用回归估计法，Bartlett估计法或Thomson估计法计算因子得分。

8、综合得分：以各因子的方差贡献率为权，由各因子的线性组合得到综合评价指标函数。

F=(w1F1+w2F2+…+wmFm)／(w1+w2+…+wm)

此处wi为旋转前或旋转后因子的方差贡献率。 9、得分排序：利用综合得分可以得到得分名次。 七）、主成分分析法的使用范围：

1、简化系统结构，探讨系统内核。可采用主成分分析、因子分析、对应分析等方法，在众多因素中找出各个变量最佳的子集合，从子集合所包含的信息描述多变量的系统结果及各个因子对系统的影响。“从树木看森林”，抓住主要矛盾，把握主要矛盾的主要方面，舍弃次要因素，以简化系统的结构，认识系统的内核。

2、构造预测模型，进行预报控制。在自然和社会科学领域的科研与生产中，探索多变量系统运动的客观规律及其与外部环境的关系，进行预测预报，以实现对系统的最优控制，是应用多元统计分析技术的主要目的。在多元分析中，用于预报控制的模型有两大类。一类是预测预报模型，通常采用多元线性回归或逐步回归分析、判别分析、双重筛选逐步回归分析等建模技术。另一类是描述性模型，通常采用聚类分析的建模技术。

3、进行数值分类，构造分类模式。在多变量系统的分析中，往往需要将系统性质相似的事物或现象归为一类。以便找出它们之间的联系和内在规律性。过去许多研究多是按单因素进行定性处理，以致处理结果反映不出系统的总的特征。进行数值分类，构造分类模式一般采用聚类分析和判别分析技术。

如何选择适当的方法来解决实际问题，需要对问题进行综合考虑。对一个问题可以综合运用多种统计方法进行分析。例如一个预报模型的建立，可先根据有关生物学、生态学原理，确定理论模型和试验设计；根据试验结果，收集试验资料；对资料进行初步提炼；然后应用统计分析方法(如相关分析、逐步回归分析、主成分分析等)研究各个变量之间的相关性，选择最佳的变量子集合；在此基础上构造预报模型，最后对模型进行诊断和优化处理，并应用于生产实际。

三、聚类分析

一）聚类分析的概念：

聚类分析指将物理或抽象对象的集合分组成为由类似的对象组成的多个类的分析过程。它是一种重要的人类行为。

聚类与分类的不同在于，聚类所要求划分的类是未知的。

聚类是将数据分类到不同的类或者簇这样的一个过程，所以同一个簇中的对象有很大的相似性，而不同簇间的对象有很大的相异性。

从统计学的观点看，聚类分析是通过数据建模简化数据的一种方法。传统的统计聚类分析方法包括系统聚类法、分解法、加入法、动态聚类法、有序样品聚类、有重叠聚类和模糊聚类等。采用k-均值、k-中心点等算法的聚类分析工具已被加入到许多著名的统计分析软件包中，如SPSS、SAS等。 二）、聚类分析的主要应用： 在商业上

聚类分析被用来发现不同的客户群，并且通过购买模式刻画不同的客户群的特征；

在生物上

聚类分析被用来动植物分类和对基因进行分类，获取对种群固有结构的认识 在地理上

聚类能够帮助在地球中被观察的数据库商趋于的相似性 在保险行业上

聚类分析通过一个高的平均消费来鉴定汽车保险单持有者的分组，同时根据住宅类型，价值，地理位置来鉴定一个城市的房产分组 在因特网应用上

聚类分析被用来在网上进行文档归类来修复信息 在电子商务上

聚类分析在电子商务中网站建设数据挖掘中也是很重要的一个方面，通过分组聚类出具有相似浏览行为的客户，并分析客户的共同特征，可以更好的帮助电子商务的用户了解自己的客户，向客户提供更合适的服务。 三）聚类分析的主要步骤： 1、数据预处理，

2、为衡量数据点间的相似度定义一个距离函数， 3、聚类或分组， 4、评估输出。

数据预处理包括选择数量，类型和特征的标度，它依靠特征选择和特征抽取，特征选择选择重要的特征，特征抽取把输入的特征转化为一个新的显著特征，它们经常被用来获取一个合适的特征集来为避免“维数灾”进行聚类，数据预处理还包括将孤立点移出数据，孤立点是不依附于一般数据行为或模型的数据，因此孤立点经常会导致有偏差的聚类结果，因此为了得到正确的聚类，我们必须将它们剔除。

既然相类似性是定义一个类的基础，那么不同数据之间在同一个特征空间相似度的衡量对于聚类步骤是很重要的，由于特征类型和特征标度的多样性，距离度量必须谨慎，它经常依赖于应用，例如，通常通过定义在特征空间的距离度量来评估不同对象的相异性，很多距离度都应用在一些不同的领域，

一个简单的距离度量，如Euclidean距离，经常被用作反映不同数据间的相异性，一些有关相似性的度量，例如PMC和SMC，能够被用来特征化不同数据的概念相似性，在图像聚类上，子图图像的误差更正能够被用来衡量两个图形的相似性。

将数据对象分到不同的类中是一个很重要的步骤，数据基于不同的方法被分到不同的类中，划分方法和层次方法是聚类分析的两个主要方法，划分方法一般从初始划分和最优化一个聚类标准开始。Crisp Clustering，它的每一个数据都属于单独的类；Fuzzy Clustering，它的每个数据可能在任何一个类中，Crisp Clustering和Fuzzy Clusterin是划分方法的两个主要技术，划分方法聚类是基于某个标准产生一个嵌套的划分系列，它可以度量不同类之间的相似性或一个类的可分离性用来合并和分裂类，其他的聚类方法还包括基于密度的聚类，基于模型的聚类，基于网格的聚类。

评估聚类结果的质量是另一个重要的阶段，聚类是一个无管理的程序，也没有客观的标准来评价聚类结果，它是通过一个类有效索引来评价，一般来说，几何性质，包括类间的分离和类内部的耦合，一般都用来评价聚类结果的质量，类有效索引在决定类的数目时经常扮演了一个重要角色，类有效索引的最佳值被期望从真实的类数目中获取，一个通常的决定类数目的方法是选择一个特定的类有效索引的最佳值，这个索引能否真实的得出类的数目是判断该索引是否有效的标准，很多已经存在的标准对于相互分离的类数据集合都能得出很好的结果，但是对于复杂的数据集，却通常行不通，例如，对于交叠类的集合。

四）聚类分析的计算方法： 1、划分法(partitioning methods)：给定一个有N个元组或者纪录的数据集，分裂法将构造K个分组，每一个分组就代表一个聚类，K录的个数无关的，它只与把数据空间分为多少个单元有关。代表算法有：STING算法、CLIQUE算法、WAVE-CLUSTER算法； 5、基于模型的方法(model-based methods):基于模型的方法给每一个聚类假定一个模型，然后去寻找能个很好的满足这个模型的数据集。这样一个模型可能是数据点在空间中的密度分布函数或者其它。它的一个潜在的假定就是：目标数据集是由一系列的概率分布所决定的。通常有两种尝试方向：统计的方案和神经网络的方案。 具体的有：

1、K-MEANS算法

k-means 算法接受输入量 k ；然后将n个数据对象划分为 k个聚类以便使得所获得的聚类满足：同一聚类中的对象相似度较高；而不同聚类中的对象相似度较小。聚类相似度是利用各聚类中对象的均值所获得一个“中心对象”（引力中心）来进行计算的。

k-means 算法的工作过程说明如下：首先从n个数据对象任意选择 k 个对象作为初始聚类中心；而对于所剩下其它对象，则根据它们与这些聚类中心的相似度（距离），分别将它们分配给与其最相似的（聚类中心所代表的）聚类；然后再计算每个所获新聚类的聚类中心（该聚类中所有对象的均值）；不断重复这一过程直到标准测度函数开始收敛为止。一般都采用均方差作为标准测度函数. k个聚类具有以下特点：各聚类本身尽可能的紧凑，而各聚类之间尽可能的分开。 2、K-MEDOIDS算法

K-MEANS有其缺点：产生类的大小相差不会很大，对于脏数据很敏感。 改进的算法：k—medoids 方法。这儿选取一个对象叫做mediod来代替上面的中心的作用，这样的一个medoid就标识了这个类。步骤：

（1）、任意选取K个对象作为medoids（O1,O2,…Oi…Ok）。 以下是循环的：

（2）、将余下的对象分到各个类中去（根据与medoid最相近的原则）； （3）、对于每个类（Oi）中，顺序选取一个Or，计算用Or代替Oi后的消耗—E（Or）。选择E最小的那个Or来代替Oi。这样K个medoids就改变了，下面就再转到2。

（4）、这样循环直到K个medoids固定下来。

这种算法对于脏数据和异常数据不敏感，但计算量显然要比K均值要大，一般只适合小数据量。 3、Clara算法

上面提到K-medoids算法不适合于大数据量的计算。现在介绍Clara算法，这是一种基于采用的方法，它能够处理大量的数据。

Clara算法的思想就是用实际数据的抽样来代替整个数据，然后再在这些抽样的数据上利用K-medoids算法得到最佳的medoids。Clara算法从实际数据中抽取多个采样，在每个采样上都用K-medoids算法得到相应的（O1,O2…Oi…Ok），然后在这当中选取E最小的一个作为最终的结果。 4、Clarans算法

Clara算法的效率取决于采样的大小，一般不太可能得到最佳的结果。

在Clara算法的基础上，又提出了Clarans的算法，与Clara算法不同的是：在Clara算法寻找最佳的medoids的过程中，采样都是不变的。而Clarans

算法在每一次循环的过程中所采用的采样都是不一样的。与上次课所讲的寻找最佳medoids的过程不同的是，必须人为地来限定循环的次数。

模糊聚类分析方法

聚类分析方法形成思路 变量的数据预处理

分类前，对原始数据进行预处理，使其所有变量尺度均匀化。方法有以下几种： 变量的标准化

设有n个样品，m个特征变量，设第i个样品，第j个变量的观测值为xij(i1,2,,n;j1,2,,m)，由此可构成一个nm阶矩阵为

x1mx11x12xxx2m （1） X(xij)nm2122xxxnmn1n2将式（1）中每个变量xij根据以下公式变换，称为标准化。 对每个变量的标准化计算公式为

xxj(i1,2,ijxij(j1,2,Sj,n) （2） ,m)1n1nSj[(xijxj)2]1/2 式中，xjxijni1ni1标准化后变量的平均值为0，标准离差为1。

变量的正规化

对每个变量施行以下变换，称为正规化。

xijxj(min)(i1,2,xijxj(max)xj(min)(j1,2,,n) （3） ,m)1。式中，xj(max)和xj(min)分别为第j个变量的最大值和最小值。显然，0xij

变量的规格化

对每个变量施行以下变换，称为规格化。

xij(i1,2,,n) （4） xijxj(max)(j1,2,,m)1。 式中，，xj(max)为第j个变量的最大值。显然，0xij注：数据的预处理以不丢失原有信息为前提。三种预处理方法的选择应根据现有

数据的特点来考虑。

分类统计量的确定及其聚类方法的选择

分类统计量的确定

一般是把相似程度大的并成一类，把相似程度小的分为不同的类，因此要定量地表示样品间的相似程度。设论域U{x1,x2,,xn}，xi{xi1,xi2,,xim}(i1,2,,n)，即数据矩阵为

Axijnm，如果xi与xj的相似程度为rijR(xi,xj)(i,j1,2,,n)，则称之为相似系

~数，确定相似系数rij有多种不同的方法。常用的方法如下：

(1) 数量积法

对

于

xi{xi1,xi2,,xim}U，令

mMmaxxikxjkijk1，则取

ij1,rij1mr0r，显然rij[0,1]。若出现有某些ij，可令ij，rij1xx,ij2Mikjkk1则有rij[0,1]。也可以用平移-极差变换将其压缩到[0,1]上，即可以得到模糊相似矩阵Rrijnm。

m(2) 夹角余弦法（相似系数统计量）: 令

xrijk1mk1ikxjkm(i,j1,2,,n)

22xxikjkk1则Rrijnn。

(3) 相关系数法（相关系数统计量）: 令

xrijk1mikxixjkxjxk1mikxix2k1mjkxjnn(i,j1,2,,n)

21m1m其中xixik，xjxjk，则Rrijmk1mk1。

注意：xi{xi1,xi2,,xim}中的样本xik属于同一个样本空间Xi(k1,2,,m)。 (4) 指数相似系数法: 令

21m3(xikxjk)rijexp 2mk14sk21n1n其中skxikxk,xkxik(k1,2,,m)。则Rrijnn。

ni1ni1注意：xi{xi1,xi2,,xim}中的样本xik属于不同的样本空间Xk，即

xikXk(k1,2,,m)。

(5) 最大最小值法: 令

rijxxk1k1mmikxjkxjkik(xij0;i,j1,2,,n)

则Rrijnn。

(6) 算术平均值法: 令

rij则Rrijxk1mmikxjk1xikxjk2k1(xij0;i,j1,2,,n)

nn。

(7) 几何平均值法：令

rij则Rrijxk1mmikxjk(xij0;i,j1,2,,n)

k1xikxjknn。

(8) 绝对值倒数法：令

ij1,1rijm Mxx,ijikjkk1其中M为使得所有rij[0,1](i,j1,2,,n)的确定常数，则Rrijnn。

(9) 绝对值指数法：令

则Rrijmrijexpxikxjk(i,j1,2,,n)

k1。

nn (10) 海明距离法（距离系数统计量。如果变量的量纲不同，原始数据变异

范围相差悬殊时，建议首先进行数据的标准化处理，然后再计算距离）：令

rij1Hd(xi,xj)m(i,j1,2,,n) d(xi,xj)xikxjkk1其中H为使得所有rij[0,1](i,j1,2,,n)的确定常数。则Rrijnn。

(11) 欧氏距离法（最常用）：令

rij1Ed(xi,xj)m(i,j1,2,,n) 2d(xi,xj)xikxjkk1其中E为使得所有rij[0,1](i,j1,2,,n)的确定常数。则Rrijnn。 (12) 契比雪夫距离法：令

rij1Qd(xi,xj)md(xi,xj)xikxjkk1(i,j1,2,,n)

其中Q为使得所有rij[0,1](i,j1,2,,n)的确定常数。则Rrijnn。

(13) 主观评分法：设有N个专家组成专家组{p1,p2,,pN}，让每一位专家对所研究的对象xi与xj相似程度给出评价，并对自己的自信度作出评估。如果第k位专家pk关于对象xi与xj的相似度评价为rij(k)，对自己的自信度评估为

aij(k)(i,j1,2,,n)，则相关系数定义为

rij则Rrijnn。

ak1Nij(k)rij(k)(i,j1,2,,n)

ijak1N(k)综上所述，以上给出了实际中能够使用的一些方法，具体地选择要根据具体问题的性质和使用的方便来确定。

在实际工作中，当需要研究样品与样品之间关系时，一般用距离系数统计量或者相似系数统计量作为分类计算依据，这种方法又称为Q型聚类法；当需要研究变量与变量之间的关系时，常用相关系数统计量作为分类计算依据，这种方法又称R型聚类法。

选择适当的聚类方法

聚合法

开始把每个样品看成自成一类，计算各类之间的相似程度的统计量，把最相似的两类合并为一类，再计算各类相似程度统计量，把最相似的两类合并，照此继续下去，一直到所有样品都聚合成一类为止，最后人为确定合适的分类数，得到分类结果。 分解法

它的聚类过程恰好和聚合法相反，开始把全体样品看成一类，然后分成二类，……，一直到每个样品为一类或分到不能再分时为止，通常要设计一个分类函数（目标函数）来控制整个分类过程。 调优法

开始人为将样品作初始分类，在一定准则下判断这个分类是否最优，如果不是最优，则对分类进行修改，再判断修改后的分类是否最优，若仍不是最优，再作修改，不断重复上述步骤，一直到分类方案最优为止。 *动态聚类法

步骤：

1、按照一定的原则选择一批凝聚点（聚核），

2、让样品向最近的凝聚点凝聚，这样就由点凝聚成类，得到初始分类。 3、初始分类不一定合理，可按最近距离原则进行修改，直到分类合理得到最终的分类为止。

四、最小二乘法与多项式拟合 一）、最小二乘法的基本原理

从整体上考虑近似函数p(x)同所给数据点(xi,yi)(i=0,1,…,m)误差

rip(xi)yi(i=0,1,…,m)

的大小，常用的方法有以下三种：一是误差

ririp(xi)yi(i=0,1,…,m)绝对值的最大值max0im，即误差 向量

r(r0,r1,rm)的∞—范数；二是误差绝对值的和i0Tmri，即误差向量r的1—

范数；三是误差平方和的算术平方根，即误差向量r的2—范数；前两种方法简单、自然，但不便于微分运算 ，后一种方法相当于考虑 2—范数的平方，

ii0rm2因此在曲线拟合中常采用误差平方和体大小。

ri0m2i来 度量误差ri(i=0，1，…，m)的整

数据拟合的具体作法是：对给定数据 (xi,yi) (i=0,1,…，m)，在取定的函数类中,求p(x),使误差rip(xi)yi(i=0,1,…,m)的平方和最小，即

i0rim2=i02p(x)yiiminm

从几何意义上讲，就是寻求与给定点(xi,yi)(i=0,1,…,m)的距离平方和为最

小的曲线yp(x)（图6-1）。函数p(x)称为拟合 函数或最小二乘解，求拟合函数p(x)的方法称为曲线拟合的最小二乘法。

在曲线拟合中，函数类可有不同的选取方法.

6—1

二）、多项式拟合

假设给定数据点(xi,yi)(i=0,1,…,m)，为所有次数不超过n(nm)的多项式构成的函数类，现求一

mpn(x)akxkk0n,使得

2Ipn(xi)yii02nakxikyimini0k0 (1)

m当拟合函数为多项式时，称为多项式拟合，满足式（1）的pn(x)称为最小二乘

拟合多项式。特别地，当n=1时，称为线性拟合或直线拟合。 显然

为a0,a1,an的多元函数，因此上述问题即为求II(a0,a1,an)的极值 问题。由多元函数求极值的必要条件，得

mnI2(akxikyi)xij0,j0,1,,naji0k0 (2)

i0k0I(akxikyi)2mn即

(xk0i0nmjki)akxijyi,i0mj0,1,,n (3)

（3）是关于a0,a1,an的线性方程组，用矩阵表示为

mmmnym1xixiiai0i0i00mmmm2n1xxixia1xiyiii0i0i0i0mammmnnnn12nxixixiyixii0 (4) i0i0i0式（3）或式（4）称为正规方程组或法方程组。

可以证明，方程组（4）的系数矩阵是一个对称正定矩阵，故存在唯一解。从式（4）中解出ak(k=0,1,…，n)，从而可得多项式

(5)

可以证明，式（5）中的pn(x)满足式（1），即pn(x)为所求的拟合多项式。我

k0pn(x)akxkn们把i0pmn(xi)yi2称为最小二乘拟合多项式pn(x)的平方误差，记作

r22

pn(xi)yii0nmm2

由式(2)可得

r22yak(xikyi)2ii0k0i0m (6)

多项式拟合的一般方法可归纳为以下几步：

(1) 由已知数据画出函数粗略的图形——散点图，确定拟合多项式的次数n；

(2) 列表计算i0和i0(3) 写出正规方程组，求出a0,a1,an；

pn(x)akxknxmji(j0,1,,2n)xmjiyi(j0,1,,2n)；

k0(4) 写出拟合多项式。

在实际应用中，nm或nm；当nm时所得的拟合多项式就是拉格朗日或牛

顿插值多项式。

例1 测得铜导线在温度Ti(℃)时的电阻Ri()如表6-1，求电阻R与温度 T的近似函数关系。 i 0 1 2 3 4 5 6 Ti(℃) 19.1 25.0 30.1 36.0 40.0 45.1 50.0 Ri() 76.30 77.80 79.25 80.80 82.35 83.90 85.10 解 画出散点图（图6-2），可见测得的数据接近一条直线，故取n=1，拟合函数为

Ra0a1T

列表如下

i 0 1 2 3 4 5 6 Ti Ri Ti2 TiRi 19.1 25.0 30.1 36.0 40.0 45.1 50.0 76.30 77.80 79.25 80.80 82.35 83.90 85.10 364.81 625.00 906.01 1296.00 1600.00 2034.01 2500.00 1457.330 1945.000 2385.425 2908.800 3294.000 3783.890 4255.000 正规方程组为 245.3 565.5 9325.83 20029.445 245.3a0565.57245.39325.83a20029.4451

解方程组得

a070.572,a10.921

故得R与T的拟合直线为

R70.5720.921T

利用上述关系式，可以预测不同温度时铜导线的电阻值。例如，由R=0得T=-242.5，即预测温度T=-242.5℃时，铜导线无电阻。

例2 已知实验数据如下表

6-2

i 0 1 2 3 4 5 6 7 8 xi 1 3 4 5 6 7 8 9 10 10 5 4 2 1 1 2 3 4 试用最小二乘法求它的二次拟合多项式。 解 设拟合曲线方程为

ya0a1xa2x2

列表如下

I 0 1 2 3 4 5 6 7 8 xi yi xi2 xi3 xi4 yixiyi xi2yi 1 3 4 5 6 7 8 9 10 10 5 4 2 1 1 2 3 4 1 9 16 25 36 49 64 81 100 1 27 64 125 216 343 512 729 1000 1 81 256 625 1296 2401 4096 6561 10000 10 15 16 10 6 7 16 27 40 10 45 64 50 36 49 128 243 400 得正规方程组 53 32 381 3017 25317 147 1025 52381a0329523813017a1471381301725317a21025

解得

a013.4597,a13.6053a20.2676

故拟合多项式为

y13.45973.60530.2676x2

*三 最小二乘拟合多项式的存在唯一性

定理1 设节点x0,x1,,xn互异，则法方程组（4）的解存在唯一。

证 由克莱姆法则，只需证明方程组（4）的系数矩阵非奇异即可。用反证法，设方程组（4）的系数矩阵奇异，则其所对应的齐次方程组

m1mxii0mxini0nmxxi0i0mmi2in1ixi0jkimmxyiai0i00mmn1axi1xiyii0i0ammn2nnxiyixii0 （7）i0nim

有非零解。式(7)可写为

j0,1,,n(xk0i0)ak0,aj （8）

n

将式（8）中第j个方程乘以(j=0,1,…，n)，然后将新得到的n+1个方程左

nmjkaj(xi)ak00右两端分别 相加，得j0k0i0

因为

mnnmnmjkmnn2jkjkaj(xi)akakajxi(ajxi)(akxi)pn(xi)j0i0j0k0i0k0i0i0j0k0 其中 npn(x)akxkk0n

所以

pn(xi)0 (i=0,1,…,m)

pn(x)是次数不超过n的多项式，它有m+1＞n个相异零点，由代数基本定理，必

须有a0a1an0，与齐次方程组有非零解的假设矛盾。因此正规方程组（4）

aa,,an必有唯一解 。定理2 设0,1是正规方程组（4）的解，则是满足式（1）的最小二乘拟合多项式。

pn(x)akxkk0n证 只需证明，对任意一组数

b0,b1,,bnm组成的多项式

2Qn(x)bkxkk0n，恒有

Qi0mn(xi)yipn(xi)yi2i0

即可。

2Q(x)yp(x)yniinii2i0i0mmQn(xi)pn(xi)2Qn(xi)pn(xi)pn(xi)yi2i0i0mm02i0j0mnnmnnk(bjaj)xiakxiyi2bjajakxikyixijj0i0k0k0

j因为ak(k=0,1,…，n)是正规方程组（4）的解，所以满足式（2），因此有

Qi0mn(xi)yipn(xi)yi022i0m

故pn(x)为最小二乘拟合多项式。

*四 多项式拟合中克服正规方程组的病态

在多项式拟合中，当拟合多项式的次数较高时，其正规方程组往往是病态的。而且

①正规方程组系数矩阵的阶数越高，病态越严重；

②拟合节点分布的区间x0,xm偏离原点越远，病态越严重； ③xi(i=0,1,…，m)的数量级相差越大，病态越严重。 为了克服以上缺点，一般采用以下措施：

①尽量少作高次拟合多项式，而作不同的分段低次拟合；

②不使用原始节点作拟合，将节点分布区间作平移，使新的节点xi关于原 点对称，可大大降低正规方程组的条件数，从而减低病态程度。 平移公式为：

xxmxixi0,i0,1,,m2 (9) ③对平移后的节点xi(i=0,1,…，m),再作压缩或扩张处理：

xipxi,p2r(m1)其中

i0,1,,m （10） ,（r是拟合次数） （11）

(x)ii0m2rxi经过这样调整可以使的数量级不太大也不太小，特别对于等距节点

xix0ih(i0,1,,m)，作式（10）和式（11）两项变换后，其正规方程

组的系数矩阵设 为A，则对1～4次多项式拟合，条件数都不太大，都可以得到满意的结果。

变换后的条件数上限表如下： 拟合次数 1 2 3 4 =1 <9.9 <50.3 <435 cond2(A) ④在实际应用中还可以利用正交多项式求拟合多项式。一种方法是构造离散正交多项式；另一种方法是利用切比雪夫节点求出函数值后再使用正交多项式。这两种方法都使正规方程 组的系数矩阵为对角矩阵，从而避免了正规方程组的病态。我们只介绍第一种，见第三节。 例如 m=19,x0=328,h=1, x1=x0+ih，i=0,1,…，19，即节点 分布在［328,347］，作二次多项式拟合时

① 直接用xi构造正规方程组系数矩阵A0，计算可得

cond2(A0)2.251016

严重病态，拟合结果完全不能用。 ② 作平移变换

328347xixi,i0,1,,192

用xi构造正规方程组系数矩阵A1，计算可得

cond2(A1)4.4838681016

比cond2(A0)降低了13个数量级，病态显著改善，拟合效果较好。③ 取压缩因子

p420(x)ii0190.14984

i0,1,,19 作压缩变换

xi用构造正规方程组系数矩阵A2，计算可得

xipxi,cond2(A2)6.839

又比cond2(A1)降低了3个数量级，是良态的方程组，拟合效果十分理想。

p(x)中使用原来节点所对应的变量x，可写为n如有必要，在得到的拟合多项式

x0xm))2

仍为一个关于x的n次多项式，正是我们要求的拟合多项式。

五、回归分析（略） Qn(x)pn(p(x六、概率分布方法（略） 七、插值与拟合（略）

八、方差分析法

一）、方差分析的意义

前述的t检验和u检验适用于两个样本均数的比较，对于k个样本均数的比较，如果仍用

t检验或u检验，需比较

次，如四个样本均数需比较次。假设每次比较所确定的检验水准=0.05，则每次检验拒绝H0不犯第一类错误的概率为1-0.05=0.95；那么6次检验都不犯第一类错误的概率为(1-0.05)6=0.7351，而犯第一类错误的概率为0.2649，因而t检验和u检验不适用于多个样本均数的比较。用方差分析比较多个样本均数,可有效地控制第一类错误。方差分析(analysis of variance,ANOVA)由英国统计学家R.A.Fisher首先提出，以F命名其统计量，故方差分析又称F检验。

二）、方差分析的基本思想

下面通过表5.1资料介绍方差分析的基本思想。

例如，有4组进食高脂饮食的家兔，接受不同处理后，测定其血清肾素血管紧张素转化酶（ACE）浓度（表5.1），试比较四组家兔的血清ACE浓度。

表5.1对照组及各实验组家兔血清ACE浓度（u/ml） 实验组 对照组 A降脂药 B降脂药 C降脂药 61.24 82.35 26.23 25.46 58.65 56.47 46.87 38.79 46.79 61.57 24.36 13.55 37.43 48.79 38.54 19.45 66.54 62.54 42.16 34.56 59.27 60.87 30.33 10.96 20.68 48.23 ) 329.92 372.59 229.17 191.00 1122.68 (6 54.99 6 7 7 26 （N ） 62.10 32.74 27.29 43.18 （） 18720.97 23758.12 8088.59 6355.43 56923.11 ( ) 由表5.1可见，26只家兔的血清ACE浓度各不相同，称为总变异；四组家兔的血清ACE浓度均数也各不相同，称为组间变异；即使同一组内部的家兔血清ACE浓度相互间也不相同，称为组内变异。该例的总变异包括组间变异和组内变异两部分，或者说可把总变异分解为组间变异和组内变异。组内变异是由于家兔间的个体差异所致。组间变异可能由两种原因所致，一是抽样误差；二是由于各组家兔所接受的处理不同。正如第四章所述，在抽样研究中抽样误差是不可避免的，故导致组间变异的第一种原因肯定存在；第二种原因是否存在，需通过假设检验

作出推断。假设检验的方法很多，由于该例为多个样本均数的比较，应选用方差分析。

方差分析的检验假设H0为各样本来自均数相等的总体，H1为各总体均数不等或不全相等。若不拒绝H0时，可认为各样本均数间的差异是由于抽样误差所致，而不是由于处理因素的作用所致。理论上，此时的组间变异与组内变异应相等，两者的比值即统计量F为1；由于存在抽样误差，两者往往不恰好相等，但相差不会太大，统计量F应接近于1。若拒绝H0，接受H1时，可认为各样本均数间的差异，不仅是由抽样误差所致，还有处理因素的作用。此时的组间变异远大于组内变异，两者的比值即统计量F明显大于1。在实际应用中，当统计量F值远大于1且大于某界值时，拒绝H0，接受H1，即意味着各样本均数间的差异，不仅是由抽样误差所致，还有处理因素的作用。

（5.1）

方差分析的基本思想是根据研究目的和设计类型，将总变异中的离均差平方和SS及其自由度分别分解成相应的若干部分，然后求各相应部分的变异；再用各部分的变异与组内（或误差）变异进行比较，得出统计量F值；最后根据F值的大小确定P值，作出统计推断。 例如，完全随机设计的方差分析，是将总变异中的离均差平方和SS及其自由度分别分解成组间和组内两部分，SS组间/组间和SS组内/组内分别为组间变异（MS组间）和组内变异（MS组内），两者之比即为统计量F（MS组间/MS组内）。 又如，随机区组设计的方差分析，是将总变异中的离均差平方和SS及其自由度分别分解成处理间、区组间和误差3部分，然后分别求得以上各部分的变异（MS处理、MS 区组和MS误差），进而得出统计量F值（MS处理/MS误差、MS区组/MS误差）。

3、方差分析的计算方法

下面以完全随机设计资料为例，说明各部分变异的计算方法。将N个受试对象随机分为k组，分别接受不同的处理。归纳整理数据的格式、符号见下表：

处理组（i）

1

…

2 …

3 … k … …

… … … …

合计

… …

1）总离均差平方和（sum of squares,SS）及自由度（freedom,ν）

总变异的离均差平方和为各变量值与总均数（）差值的平方和，离均差平方和和自由度分别为：

（5.2）

=N-1（5.3）

2）组间离均差平方和、自由度和均方 组间离均差平方和为各组样本均数(

)与总均数()差值的平方和

（5.4）

（5.5）

（5.6）

3）组内离均差平方和、自由度和均方

组内离均差平方和为各处理组内部观察值与其均数（）差值的平方和之和，

数理统计证明，总离均差平方和等于各部分离均差平方和之和，因此，

（5.7）

（5.8）

（5.9）

4）三种变异的关系：

= N-1= (k-1)+(N-k) =

可见，完全随机设计的单因素方差分析时，总的离均差平方和（SS总）可分解为组间离均差平方和（SS组间）与组内离均差平方和（SS组内）两部分；相应的总自由度（分解为组间自由度（

）和组内自由度（

）两部分。

）也

5）方差分析的统计量：

（5.10）

4、方差分析的应用条件与用途

方差分析的应用条件为①各样本须是相互独立的随机样本；②各样本来自正态分布总体；③各总体方差相等，即方差齐。

方差分析的用途①两个或多个样本均数间的比较；②分析两个或多个因素间的交互作用；③回归方程的线性假设检验；④多元线性回归分析中偏回归系数的假设检验；⑤两样本的方差齐性检验等。

九、逼近理想点排序法

原理：通过测度各个被测评对象的指标评价值向量与评价的理想解和负理想解的相对距离进行测评排序，同时计算各评价对象的综合评价指数。 确定规范化决策矩阵

无量纲化处理

bijaijamin(i)amax(i)amin(i) -----------------→ 规范化决策矩阵 B(bij)ij

（第j个被测评对象的第i个指标的无量纲化处理公式） 确定指标的权重系数（以变异系数法为例）  先求不同指标下指标评价的均值ai和标准差Si  再计算各指标的变异系数，取其绝对值为Vi  对作归一化处理，得各指标的权重WiVi/Vi

 再由规范化决策矩阵B和权重构成加权规范阵RWB(rij)

确定理想解x和负理想解x

xri/i1,xri/i1,,mmaxrij/j1,jjij,n ,mminr/j1,,n

计算各被测评对象到理想解距离d与负理想解的距离d

dj(ri1nijri) d2j(ri1nijri)2 (j=1,…,n)

计算被测评对象与理想解的相对接近度，作为其综合评价指数

cjdjddjj100(j1,,n)

cj值越大，则顾客满意程度越高

十、动态加权法 动态加权：

关于不同的指标可以取相同的权函数，也可以取不同的权函数。

举例：长江水质…… 数据：

求解：

十一、灰色关联分析法

灰色关联度是两个系统或两个因素间关联性大小的量度，它描述系统发展过程中因素间相对变化的情况，也就是变化大小、方向与速度等的相对性。

如果两因素在发展过程中相对变化态势一致性高，则两者的灰色关联度大；反之，灰色关联度就小。

所谓灰色关联分析，就是系统的因素分析，是对一个系统发展变化态势的定量比较和反映。灰色关联分析是通过灰色关联度来分析和确定系统因素间的影响程度或因素对系统主行为的贡献测度的一种方法。

灰色关联分析的基本思想是根据序列曲线几何形状的相似程度来判断其联系是否紧密。曲线越接近,相应序列之间的关联度就越大,反之就越小。灰色关联分析方法弥补了用数理统计作系统分析所导致的缺憾。它对样本量的多少和样本有无规律都同样适用,而且计算量小,十分方便,更不会出现量化结果与定性分析结果不符的情况。

具体步骤：灰色系统关联分析的具体计算步骤如下 ：

（1）确定反映系统行为特征的参考数列和影响系统行为的比较数列 反映系统行为特征的数据序列，称为参考数列。影响系统行为的因素组成的数据序列，称比较数列。

（2）对参考数列和比较数列进行无量纲化处理

由于系统中各因素的物理意义不同，导致数据的量纲也不一定相同，不便于比较，或在比较时难以得到正确的结论。因此在进行灰色关联度分析时，一般都要进行无量纲化的数据处理。

（3）求参考数列与比较数列的灰色关联系数ξ（Xi）

所谓关联程度，实质上是曲线间几何形状的差别程度。因此曲线间差值大小，可作为关联程度的衡量尺度。对于一个参考数列X0有若干个比较数列X1, X2,…, Xn，各比较数列与参考数列在各个时刻（即曲线中的各点）的关联系数ξ（Xi）可由下列公式算出：

(k)ˆ0kX0kmaxmaxXˆ0kX0kminminXˆ0kX0kmaxmaxXˆ0kX0kX

称为关联系数，其中称为分辨系数，（0，1），常取0.5.实数 第二级最小差，记为Δmin。 两级最大差，记为Δmax。 为各比较数列Xi曲线上的每一个点与参考数列X0曲线上的每一个点的绝对差值。记为Δoi(k)。所以关联系数ξ（Xi）也可简化如下列公式： r(xo(k),xi(k))(minminoi(k)maxmaxoi(k))(oi(k)maxmaxoi(k)) ikikik（4）求关联度ri 因为关联系数是比较数列与参考数列在各个时刻（即曲线中的各点）的关联程度值，所以它的数不止一个，而信息过于分散不便于进行整体性比较。因此有必要将各个时刻（即曲线中的各点）的关联系数集中为一个值，即求其平均值，作为比较数列与参考数列间关联程度的数量表示，关联度ri公式如下：

1nˆ0k的关联度 rk称为X0k与Xnk1（5）排关联序

因素间的关联程度，主要是用关联度的大小次序描述，而不仅是关联度的大小。将m个子序列对同一母序列的关联度按大小顺序排列起来，便组成了关联序，记为{x}，它反映了对于母序列来说各子序列的“优劣”关系。若r0i>r0j，则称{xi}对于同一母序列{x0}优于{xj}，记为{xi}>{xj} ；若r0i表1 代表旗县参考数列、比较数列特征值。

十二、灰色预测法 灰色预测

注：参考人口预测论文<纪江版>（灰色预测+时间序列的一次平滑指数预测法） 1、灰色预测一般有四种类型：

（1）、数列预测。对某现象随时间的顺延而发生的变化所做的预测定义为数列预测。例如对消费物价指数的预测，需要确定两个变量，一个是消费物价指数的水平。另一个是这一水平所发生的时间。

（2）、灾变预测。对发生灾害或异常突变时间可能发生的时间预测称为灾变预测。例如对地震时间的预测。

（3）、系统预测。对系统中众多变量间相互协调关系的发展变化所进行的预测称为系统预测。例如市场中替代商品、相互关联商品销售量互相制约的预测。 （4）、拓扑预测。将原始数据作曲线，在曲线上按定值寻找该定值发生的所有时点，并以该定值为框架构成时点数列，然后建立模型预测未来该定值所发生的时点。

2、使用方法前一定要在段前作一个引子，连接问题分析和数据特点，以下便是： 通过对已知数据的分析，随着时间的变化，排污量一直呈增长趋势，并且增长的很快。在这里利用灰色预测模型对（ ）进行预测。 通过对数据的分析，传统的数理统计预测方法往往需要足够多的数据，而本问题的数据给出的数据偏小，如果采用传统的方法误差太大。根据上述的特点可采用灰色预测模型。

3、灰色预测具体步骤：

（1）、首先是数据的检验处理，要求级比

x(0)(i1)(i)(0)(en1,en1（)i2,3,,n)

x(i)2n12n122A、如果不全属于(e,e)，则要做必要的变换处理（如取适当的常数C，作平移变换），使其落入区域中。

B、若A不成立，则建立GM（1，1）模型 （2）、建立GM（1，1）模型

步骤一：一次累加生成数列AGO，（目的是弱化原始时间序列的随机性，增加其稳定程度）

x(1)(x(1)(1),x(1)(2),x(1)(n) 步骤二：求均值数列

z(1)(k)z(1)(k)(1)z(1)(k1)(k2,3,,n) z(1)(z(1)(2),z(1)(3),,z(1)(n)

步骤三：建立GM（1，1）模型相应的白化微分方程 dX1aX1 dt其中：α称为发展灰数；μ称为内生控制灰数。

步骤四：求的参数估计a、b（最小二乘法）

ˆBTBBTYn 1步骤五： 给出累加时间数列预测模型

ˆ1k1X01eak，k0,1,2...,n Xaa步骤六：做差得到原始预测值

bˆ(0)(k1)xˆ(1)(k1)xˆ(1)(k)(x(0)(1))(eatea(k1)) xa4、检验预测值 A、残差检验

x(0)(k)x(0)(k)(k)（i1,2,,n)（若(k)<0.2,则达到一般要求;若(k)<0.1,(0)x(k)则效果好

B\\级比偏差值检验

步骤一;首先有参考数据

x(0)(x(0)(1),x(0)(2),x(0)(3),,x(0)(n))计算出级比0(k)，再由发展系数a，求出相应级比偏差

10.5ak10(k)

10.5a若(k)<0.2,则达到一般要求;若(k)<0.1,则效果好 程序实现：

采用EXCEl的方法实现灰色预测。

十三、模糊综合评价 1. 模糊综合评判的一般提法

设U{u1,u2,,un}为研究对象的n种因素（或指标），称之为因素集（或指标集）．V{v1,v2,,vm}为诸因素（或指标）的m种评判所构成的评判集（或称语集、评价集、决策集等），它们的元素个数和名称均可根据实际问题的需要和决策人主观确定．实际中，很多问题的因素评判集都是模糊的，因此，综合评判应该是

V上的一个模糊子集

B(b1,b2,,bm)F(V)

其中bk为评判vk对模糊子集B的隶属度：B(vk)bk(k1,2,,,m)，即反映了第k种评判vk在综合评价中所起的作用．综合评判B依赖于各因素的权重，即它应该是U上的模糊子集A(a1,a2,,an)F(U)，且ai1，其中ai表示第i种因素

i1n的权重．于是，当权重A给定以后，则相应地就可以给定一个综合评判B． 2. 模糊综合评判的一般步骤

(1) 确定因素集U{u1,u2,,un}； (2) 确定评判集V{v1,v2,,vm}；

(3) 确定模糊评判矩阵R(rij)nm：

首先，对每一个因素ui做一个评判f(ui)(i1,2,,n)，则可以得U到V的一个模糊映射f，即

f:UF(U)

uif(ui)(ri1,ri2,,rim)F(V)然后，由模糊映射f可以诱导出模糊关系RfF(UV)，即

Rf(ui,vj)f(ui)(vj)rij(i1,2,,n;j1,2,,m)

因此，可以确定出模糊评判矩阵R(rij)nm．而且称(U,V,R)为模糊综合评判模型，U,V,R称为该模型的三要素．

(4) 综合评判：对于权重A(a1,a2,,an)F(U)，用模型M(,)取最大－最小合成运算，可以得到综合评判

BAR(bj(airij),j1,2,,m)

i1n注意到：关于评判集V的权重A(a1,a2,,an)的确定在综合评判中起重要的作用，通常情况下可以由决策人凭经验给出，但往往带有一定的主观性．要从实际出发，或更客观地反映实际情况可采用专家评估法、加权统计法和频数统计法，或更一般的模糊协调决策法、模糊关系方法等来确定． 综合评判模型的构成

如果模糊综合评判模型为(U,V,R),对于权重A(a1,a2,,an)F(U)，模糊评判矩阵为R(rij)nm，则用模型M(,)运算得综合评判为

BAR(b1,b2,,bm)F(V)，其中bj(airij)(j1,2,,m)．

n事实上，由于ai1，对于某些情况可能会出现airij，即airijai．这

i1ni1样可能导致模糊评判矩阵R中的许多信息的丢失，即人们对某些因素ui所作的评判信息在决策中未得到充分的利用．从而导致综合评判结果失真．为此，实际中可以对模型M(,)进行改进．

(1) 模型M(•,)法：对于A(a1,a2,,an)F(U)和R(rij)nm，则用模型

M(•,)运算得BAR，即bj(ai•rij)(j1,2,,m)．

n(2) 模型M(,)法：对于A(a1,a2,,an)F(U)和R(rij)nm，则用模型

M(,)运算得BAR，即bj(airij)(j1,2,,m)．

ni1(3) 模型M(•,)法：对于A(a1,a2,,an)F(U)和R(rij)nm，则用模型

M(•,)运算得BAR，即bj(ai•rij)(j1,2,,m)．

i1在实际应用时，主因素（即权重最大的因素）在综合中起主导作用时，则可首选“主因素决定型”模型M(,)；当模型M(,)失效时，再来选用“主因素突出型”模型M(•,)和M(,)；当需要对所有因素的权重均衡时，可选用加权平均模型M(•,)．在模型的选择时，还要特别注意实际问题的需求．

ni1多层次模糊综合评判

对于实际中的许多问题往往都是涉及因素多，各因素的权重分配较为均衡的情况，此时，可采用将诸因素分为若干个层次进行研究．即首先分别对单层次的各因素进行评判，然后再对所有的各层次因素作综合评判．这里仅就两个层次的情况进行说明，具体方法如下：

将因素集U{u1,u2,,un}分成若干个组U1,U2,,Uk(1kn)使得UUi，

i1k且UiUj(ij)，称U{U1,U2,,Uk}为一级因素集。

不妨设Ui{u,u,,u}(i1,2,,k;(i)1(i)2(i)nini1kin)，称之为二级因素集．

(i)(i)设评判集V{v1,v2,,vm}，对二级因素集Ui{u1(i),u2,,un}的ni个因素进行i单因素评判，即建立模糊映射

fi:UiF(V)i)(i)u(ji)fi(u(ji))(rj(1i),rj(2,,rjm)(j1,2,,ni)

于是得到评判矩阵为

(i)r11(i)r21Ri(i)rni1(i)r12(i)r22)rn(ii2i)r1(mi)r2(m

)rn(iim(i)(i)(i)(i)(i)不妨设Ui{u1(i),u2的权重为则可以求得综合评,,un}A(a,a,,a}，i12nii判为

(i)(i)BiAiRi(b1(i),b2,,bm)(i1,2,,k)

其中b(ji)由模型M(,)，或M(•,)、M(,)、M(•,)确定．

对于一级因素集U{U1,U2,,Uk}作综合评判，不妨设其权重A(a1,a2,,ak)，总评判矩阵为R[B1,B2,,Bk]T．按模型M(,)，或M(•,)、M(,)，、M(•,)运算得到综合评判BAR(b1,b2,,bm)F(V)．

十四、隶属函数的刻画（略）

十五、时间序列分析法

ARIMA(autoregressive integrated moving average models)时间序列模型 一般概念;

系统中某一变量的观测值按时间序列（时间间隔相同）排列成一个数值序列，展示研究对象在一定时期内的变动过程，从中寻找和分析事物的变化特征、发展趋势和规律。他是系统中某一变量受其他各种因素影响的总结果。

变动特点： 趋势性：某个变量随时间进展或自变量变化，呈现一种比较缓慢而长期的持续上升、下降、停留的同性质变动趋势，但变动幅度可能不等。

周期性：某因素由于外部影响随着自然季节的交替出现高峰与低谷的规律。 随机性：个别为随机变动，整体呈统计规律 综合性：实际变化情况一般是几种变动的叠加或组合。预测时一般设法过滤去不规则变动，突出反映趋势性和周期性变动。

特征识别：认识时间序列所具有的变动特征，以便在系统预测时选择采用不同的方法

随机性：均匀分布、无规则分布，可能符合某统计分布（用因变量的散点图和直方图及其包含的正态分布检验随机性，大多服从正态分布） 平稳性：样本序列的自相关函数在某一固定水平线附近摆动，即方差和数学期望稳定为常数

特征识别利用自相关函数ACF:kk/0，其中k是yt的k阶自协方差，且

01，-1<k<1

平稳过程的自相关系数和偏自相关系数都会以某种方式衰减趋于0，前者测度当前序列与先前序列之间简单和常规的相关程度，后者是在控制其它先前序列的影响后，测度当前序列与某一先前序列之间的相关程度。

实际上，预测模型大都难以满足这些条件，现实的经济、金融、商业等序列都是非稳定的，但通过数据处理可以变换为平稳的。

基本步骤：

分析数据序列的变化特征 选择模型形式和参数检验 利用模型进行趋势预测 评估预测结果并修正模型

自回归AR(p)模型（自己影响自己，但可能存在误差，误差即没有考虑到的因素）

模型意义

仅通过时间序列变量的自身历史观测值来反映有关因素对预测目标的影响和作用，不受模型变量互相独立的假设条件约束，所构成的模型可以消除普通回归预测方法中由于自变量选择、多重共线性的比你更造成的困难

用PACF函数判别（从p阶开始的所有偏自相关系数均为0）

移动平均MA(q)模型

模型含义

用过去各个时期的随机干扰或预测误差的线性组合来表达当前预测值。AR(q)的假设条件不满足时可以考虑用此形式。

用ACF函数判别（从q阶开始的所有自相关系数均为0）

自回归移动平均ARMA(p,q)模型

识别条件

平稳时间序列的偏相关系数k和自相关系数rk均不截尾，但较快收敛到0，

则该时间序列可能是ARMA(p,q)模型。实际问题中，多数要用此模型。因此建模解模的主要工作时求解p,q和、的值，检验t和yt的值。 模型阶数

实际应用中p,q一般不超过2.

自回归综合移动平均ARIMA(p,d,q)模型 模型含义

模型形式类似ARMA(p,q)模型，但数据必须经过特殊处理。特别当线性时间序列非平稳时，不能直接利用ARMA(p,q)模型，但可以利用有限阶差分使非平稳时间序列平稳化，实际应用中d（差分次数）一般不超过2. 模型识别

平稳时间序列的偏相关系数k和自相关系数rk均不截尾，且缓慢衰减收敛，则该时间序列可能是ARIMA(p,d,q)模型。

若时间序列存在周期性波动，则可按时间周期进行差分，目的是将随机误差有长久影响的时间序列变成仅有暂时影响的时间序列。即差分处理后新序列符合ARMA(p,q)模型，元序列符合ARIMA(p,d,q)模型。

一个平稳的随机过程有以下要求：均数不随时间变化，方差不随时间变化，自相关系数只与时间间隔有关，而与所处的时间无关。

偏自相关函数（PACF）解决如下问题： 高阶的自相关是否真的非常重要？

是他的确有意义，还是因为低阶自相关系数较大才引起高阶自相关系数也大？ 如果建立一个以前值预测现在值的回归模型，需要包括多少个以前值？

指数平滑法用序列过去值的加权均数来预测将来的值，并且给序列中近期的数据以较大的权重，远期的数据给以较小的权重。理由是随着时间流逝，过去值的影响逐渐减小。

指数平滑法应用时存在以下问题：

指数平滑法只适合于影响时间的消逝呈指数下降的数据、 指数平滑法的每次预测都是根据上一个数来的，一般来说，用序列的第一个数作为初始值。如果数据点较多，那么经过指数衰减后，初始值的影响就不明显了。但是如果数据点少，则初始值的影响会很大，甚至大于近期的数据点，这就违背指数平滑影响呈指数衰减的假设了。所以，如果数据点少时应该考虑初始值的问题，一般来说，数据点大于40初始值的影响就不太明显。

需要指出的是，时间序列模型的预测一般不能太超前，对过于遥远的时间预测结果大多是不准确的。

十六、蒙特卡罗(MC)仿真模型 模型介绍：

蒙特卡罗（Monte Carlo）方法，又称随机抽样或统计试验方法，属于计算数学的一个分支，它是在本世纪四十年代中期为了适应当时原子能事业的发展而发展起来的。传统的经验方法由于不能逼近真实的物理过程，很难得到满意的结果，而蒙特卡罗方法由于能够真实地模拟实际物理过程，故解决问题实际非常符合，可以得到很圆满的结果。这也是我们采用该方的原因。

蒙特卡罗的基本原理及思想：

当所要求解的问题是某种事件出现的概率，或者是某个随机变量的期望值时，它们可以通过某种“试验”的方法，得到这种事件出现的频率，或者这个随机变数的平均值，并用它们作为问题的解。这就是蒙特卡罗方法的基本思想。蒙特卡罗方法通过抓住事物运动的几何数量和几何特征，利用数学方法来加以模拟，即进行一种数字模拟实验。它是以一个概率模型为基础，按照这个模型所描绘的过程，通过模拟实验的结果，作为问题的近似解。可以把蒙特卡罗解题归结为三个主要步骤：构造或描述概率过程；实现从已知概率分布抽样；建立各种估计量。

蒙特卡罗解题的三个主要步骤：

（1）、构造或描述概率过程：

对于本身就具有随机性质的问题，如粒子输运问题，主要是正确描述和模拟这个概率过程，对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。即要将不具有随机性质的问题转化为随机性质的问题。 （2）、实现从已知概率分布抽样： 构造了概率模型以后，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量（或随机向量），就成为实现蒙特卡罗方法模拟实验的基本手段，这也是蒙特卡罗方法被称为随机抽样的原因。最简单、最基本、最重要的一个概率分布是(0,1)上的均匀分布（或称矩形分布）。随机数就是具有这种均匀分布的随机变量。随机数序列就是具有这种分布的总体的一个简单子样，也就是一个具有这种分布的相互独立的随机变数序列。产生随机数的问题，就是从这个分布的抽样问题。在计算机上，可以用物理方法产生随机数，但价格昂贵，不能重复，使用不便。另一种方法是用数学递推公式产生。这样产生的序列，与真正的随机数序列不同，所以称为伪随机数，或伪随机数序列。不过，经过多种统计检验表明，它与真正的随机数，或随机数序列具有相近的性质，因此可把它作为真正的随机数来使用。由已知分布随机抽样有各种方法，与从(0,1)上均匀分布抽样不同，这些方法都是借助于随机序列来实现的，也就是说，都是以产生随机数为前提的。由此可见，随机数是我们实现蒙特卡罗模拟的基本工具。 （3）、建立各种估计量：

一般说来，构造了概率模型并能从中抽样后，即实现模拟实验后，我们就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计。建立各种估计量，相当于对模拟实验的结果进行考察和登记，从中得到问题的解。 蒙特卡罗的特点及优缺点：

蒙特卡罗方法与一般计算方法有很大区别，一般计算方法对于解决多维或因素复杂的问题非常困难，而蒙特卡罗方法对于解决这方面的问题却比较简单。其特点如下：

· 直接追踪粒子，物理思路清晰，易于理解。

· 采用随机抽样的方法，较真切的模拟粒子输运的过程，反映了统计涨落的规律。

· 不受系统多维、多因素等复杂性的限制，是解决复杂系统粒子输运问题的好方法。

· MC程序结构清晰简单。

· 研究人员采用MC方法编写程序来解决粒子输运问题，比较容易得到自己想得到的任意中间结果，应用灵活性强。

十七、BP神经网络方法

摘 要 人工神经网络是一种新的数学建模方式，它具有通过学习逼近任意非线性映射的能力。本文提出了一种基于动态BP神经网络的预测方法，阐述了其基本原理，并以典型实例验证。 关键字 神经网络，BP模型，预测 1 引言

在系统建模、辨识和预测中，对于线性系统，在频域，传递函数矩阵可以很好地表达系统的黑箱式输入输出模型；在时域，Box-Jenkins方法、回归分析方法、ARMA模型等，通过各种参数估计方法也可以给出描述。对于非线性时间序列预测系统，双线性模型、门限自回归模型、ARCH模型都需要在对数据的内在规律知道不多的情况下对序列间关系进行假定。可以说传统的非线性系统预测，在理论研究和实际应用方面，都存在极大的困难。相比之下，神经网络可以在不了解输入或输出变量间关系的前提下完成非线性建模[4，6]。神经元、神经网络都有非线性、非局域性、非定常性、非凸性和混沌等特性，与各种预测方法有机结合具有很好的发展前景，也给预测系统带来了新的方向与突破。建模算法和预测系统的稳定性、动态性等研究成为当今热点问题。目前在系统建模与预测中，应用最多的是静态的多层前向神经网络，这主要是因为这种网络具有通过学习逼近任意非线性映射的能力。利用静态的多层前向神经网络建立系统的输入/输出模型，本质上就是基于网络逼近能力，通过学习获知系统差分方程中的非线性函数。但在实际应用中，需要建模和预测的多为非线性动态系统，利用静态的多层前向神经网络必须事先给定模型的阶次，即预先确定系统的模型，这一点非常难做到。近来，有关基于动态网络的建模和预测的研究，代表了神经网络建模和预测新的发展方向。 2 BP神经网络模型

BP网络是采用Widrow-Hoff学习算法和非线性可微转移函数的多层网络。典型的BP算法采用梯度下降法，也就是Widrow-Hoff算法。现在有许多基本的优化算法，例如变尺度算法和牛顿算法。如图1所示，BP神经网络包括以下单元：①处理单元(神经元)（图中用圆圈表示），即神经网络的基本组成部分。输入层的处理单元只是将输入值转入相邻的联接权重，隐层和输出层的处理单元将它们的输入值求和并根据转移函数计算输出值。②联接权重(图中如V,W)。它将神经网络中的处理单元联系起来，其值随各处理单元的联接程度而变化。③层。神经网络一般具有输入层x、隐层y和输出层o。④阈值。其值可为恒值或可变值，它可使网络能更自由地获取所要描述的函数关系。⑤转移函数F。它是将输入的数据转化为输出的处理单元，通常为非线性函数。

图1 BP神经网络结构

2.1 基本算法

BP算法主要包含4步，分为向前传播和向后传播两个阶段： 1）向前传播阶段

（1）从样本集中取一个样本（Xp，Yp），将Xp输入网络； （2）计算相应的实际输出Op

在此阶段，信息从输入层经过逐级的变换，传送到输出层。这个过程也是网络在完成训练后正常运行时的执行过程。 2）向后传播阶段

（1）计算实际输出Op与相应的理想输出Yp的差； （2）按极小化误差的方式调整权矩阵。

这两个阶段的工作受到精度要求的控制，在这里取

作为网络关于第p个样本的误差测度，而将网络关于整个

样本集的误差测度定义为

。图2是基本BP算法的流程图。

图2 BP基本算法流程

2.2 动态BP神经网络预测算法

在经典的BP算法以及其他的训练算法中都有很多变量，这些训练算法可以确定一个ANN结构，它们只训练固定结构的ANN权值（包括联接权值和结点转换函数）。在自动设计ANN结构方面，也已有较多的尝试，比如构造性算法和剪枝算法。前一种是先随机化网络，然后在训练过程中有必要地增加新的层和结点；而剪枝法则正好相反。文献[2]中提出了演化神经网络的理念，并把EP算法与BP进行了组合演化；也有很多学者把遗传算法和BP进行结合，但这些算法都以时间复杂度以及空间复杂度的增加为代价。根据Kolmogorov定理,对于任意给定的L2型连续函数f: [ 0, 1 ]n →Rm , f可以精确地用一个三层前向神经网络来实现，因而可以只考虑演化网络的权值和结点数而不影响演化结果。基于此，在BP原有算法的基础上，增加结点数演化因子，然后记录每层因子各异时演化出的结构，最后选取最优的因子及其网络结构，这样就可以避免由于增加或剪枝得到的局部最优。根据实验得知，不同的预测精度也影响网络层神经元的结点数，所以可根据要求动态地建立预测系统。具体步骤如下： （1）将输入向量和目标向量进行归一化处理。

（2）读取输入向量、目标向量，记录输入维数m、输出层结点数n。 （3）当训练集确定之后，输入层结点数和输出层结点数随之而确定，首先遇到的一个十分重要而又困难的问题是如何优化隐层结点数和隐层数。实验表明，如果隐层结点数过少，网络不能具有必要的学习能力和信息处理能力。反之，若过多，不仅会大大增加网络结构的复杂性（这一点对硬件实现的网络尤其重要），网络在学习过程中更易陷入局部极小点，而且会使网络的学习速度变得很慢。隐层结点数的选择问题一直受到神经网络研究工作者的高度重视。Gorman指出隐层结点数s与模式数N的关系是：s＝log2N；Kolmogorov定理表明，隐层结点数s＝2n＋1（n为输入层结点数）；而根据文献[7]：s＝sqrt（0.43mn＋0.12nn＋2.54m＋0.77n＋0.35）＋0.51[7]。

（4）设置结点数演化因子a。为了快速建立网络，可以对其向量初始化， 并从小到大排序[4，7]。

（5）建立BP神经网络。隐含层传递函数用tansig，输出层用logsig，训练函数采用动态自适应BP算法，并制订停止准则：目标误差精度以及训练代数。

（6）初始化网络。

（7）训练网络直到满足停止判断准则。

（8）用测试向量对网络进行预测，并记录误差和逼近曲线，评估其网络的适应性。其适应度函数采取规则化均方误差函数。

（9）转到（5），选取下一个演化因子，动态增加隐含层结点数，直到最后得到最佳预测网络。 3 基于神经网络的预测原理 3.1 正向建模

正向建模是指训练一个神经网络表达系统正向动态的过程，这一过程建立的神经网络模型称为正向模型，其结构如图3所示。其中，神经网络与待辨识的系统并联，两者的输出误差用做网络的训练信号。显然，这是一个典型的有导师学习问题，实际系统作为教师，向神经网络提供算法所需要的期望输出。当系统是被控对象或传统控制器时，神经网络多采用多层前向网络的形式，可直接选用BP网络或它的各种变形。而当系统为性能评价器时，则可选择再励学习算法，这时网络既可以采用具有全局逼近能力的网络（如多层感知器），也可选用具有局部逼近能力的网络（如小脑模型控制器等）。

[4]

图3 正向建模结构

3.2 逆向建模

建立动态系统的逆模型，在神经网络中起着关键作用，并且得到了广泛的应用。其中，比较简单的是直接逆建模法，也称为广义逆学习。其结构如图4所示，拟预报的系统输出作为网络的输入，网络输出与系统输入比较，相应的输入误差用于训练，因而网络将通过学习建立系统的逆模型。但是，如果所辨识的非线性系统是不可逆的，利用上述方法将得到一个不正确的逆模型。因此，在建立系统时，可逆性应该先有所保证。

图4 直接逆建模结构

4 应用实例分析

以我国西南某地震常发地区的地震资料作为样本来源，实现基于动态神经网络的地震预报。根据资料，提取出7个预报因子和实际发生的震级M作为输入和目标向量。预报因子为半年内M>=3的地震累计频度、半年内能量释放积累值、b值、异常地震群个数、地震条带个数、是否处于活动期内以及相关地震区地震级。在训练前，对数据进行归一化处理。由于输入样本为7维的输入向量，一般情况下输入层设7个神经元。根据实际情况，输出层神经元个数为1。隐含层神经元的传递函数为S型正切函数，输出层也可以动态选择传递函数。实例数据来自文献[4]，将数据集分为训练集、测试集和确定集。表1中的7×7数组表示归一化后的训练向量，第一个7表示预报因子数，第二个7表示样本数。

表1 归一化后的训练向量

在不同神经元数情况下，对网络进行训练和仿真，得到如图5所示的一组预测误差曲线。其中，曲线A表示隐层结点数为6时的预测误差曲线，曲线B表示隐含层结点数为3时的预测误差曲线，曲线C表示隐含层结点数为5时的预测误差曲线，曲线D表示隐含层结点数为4时的预测误差曲线。将五种情况下的误差进行对比，曲线C表示的网络预测性能最好，其隐含层神经元数为5，图中曲线E表示的是隐含层结点数为15时的预测误差曲线（文献[4]中的最好结果）。同时也证明，在设计BP网络时，不能无限制地增加层神经元的个数。若过多，不仅会大大增加网络结构的复杂性，网络在学习过程中更易陷入局部极小点，而且会使网络的学习速度、预测速度变得很慢。

图5 不同神经元数预测误差对比曲线

5 结论

本文针对基本的BP神经网络，提出了可动态改变神经元数（与精度相关）的BP神经网络预测方法，可以根据实际情况建立预测系统。用此种方法可以建立最好的神经网络，不会有多余的神经元，也不会让网络在学习过程中过早陷于局部极小点。

神经网络的优缺点：

多层前向BP网络是目前应用最多的一种神经网络形式, 但它也不是非常完美的, 为了更好的理解应用神经网络进行问题求解, 这里对它的优缺点展开讨论:

多层前向BP网络的优点：

①网络实质上实现了一个从输入到输出的映射功能，而数学理论已证明它具有实

现任何复杂非线性映射的功能。这使得它特别适合于求解内部机制复杂的问题； ②网络能通过学习带正确答案的实例集自动提取“合理的”求解规则，即具有自学习能力；

③网络具有一定的推广、概括能力。

多层前向BP网络的问题：

①BP算法的学习速度很慢，其原因主要有：

a 由于BP算法本质上为梯度下降法，而它所要优化的目标函数又非常复杂，因此，必然会出现“锯齿形现象”，这使得BP算法低效；

b 存在麻痹现象，由于优化的目标函数很复杂，它必然会在神经元输出接近0或1的情况下，出现一些平坦区，在这些区域内，权值误差改变很小，使训练过程几乎停顿；

c 为了使网络执行BP算法，不能用传统的一维搜索法求每次迭代的步长，而必须把步长的更新规则预先赋予网络，这种方法将引起算法低效。 ②网络训练失败的可能性较大，其原因有：

a 从数学角度看，BP算法为一种局部搜索的优化方法，但它要解决的问题为求解复杂非线性函数的全局极值，因此，算法很有可能陷入局部极值，使训练失败； b 网络的逼近、推广能力同学习样本的典型性密切相关，而从问题中选取典型样本实例组成训练集是一个很困难的问题。

③难以解决应用问题的实例规模和网络规模间的矛盾。这涉及到网络容量的可能性与可行性的关系问题，即学习复杂性问题；

④网络结构的选择尚无一种统一而完整的理论指导，一般只能由经验选定。为此，有人称神经网络的结构选择为一种艺术。而网络的结构直接影响网络的逼近能力及推广性质。因此，应用中如何选择合适的网络结构是一个重要的问题； ⑤新加入的样本要影响已学习成功的网络，而且刻画每个输入样本的特征的数目也必须相同；

⑥网络的预测能力（也称泛化能力、推广能力）与训练能力（也称逼近能力、学习能力）的矛盾。一般情况下，训练能力差时，预测能力也差，并且一定程度上，随训练能力地提高，预测能力也提高。但这种趋势有一个极限，当达到此极限时，随训练能力的提高，预测能力反而下降，即出现所谓“过拟合”现象。此时，网络学习了过多的样本细节，而不能反映样本内含的规律。

十八、数据包络分析法（DEA） 数据包络分析法介绍：

通过明确地考虑多种投入（即资源）的运用和多种产出（即服务）的产生，它能够用来比较提供相似服务的多个服务单位之间的效率，这项技术被称为数据包络线分析（DEA）。它避开了计算每项服务的标准成本，因为它可以把多种投入和多种产出转化为效率比率的分子和分母，而不需要转换成相同的货币单位。因此，用DEA衡量效率可以清晰地说明投入和产出的组合，从而，它比一套经营比率或利润指标更具有综合性并且更值得信赖。

DEA是一个线形规划模型，表示为产出对投入的比率。通过对一个特定单位的效率和一组提供相同服务的类似单位的绩效的比较，它试图使服务单位的效率最大化。在这个过程中，获得100%效率的一些单位被称为相对有效率单位，而另外的效率评分低于100%的单位本称为无效率单位。 这样，企业管理者就能运用DEA来比较一组服务单位，识别相对无效率单位，衡量无效率的严重性，并通过对无效率和有效率单位的比较，发现降低无效率的方法。

DEA线形规划模型的建立： 1) 定义变量

设Ek（k=1，2，……， K）为第k个单位的效率比率，这里K代表评估单位的总数。

设uj（j=1，2，……， M）为第j种产出的系数，这里M代表所考虑的产出种类的总数。变量uj用来衡量产出价值降低一个单位所带来的相对的效率下降。 设vI（I=1，2，……，N）为第I种投入的系数，这里N代表所考虑的投入种类的综合素。变量vI用来衡量投入价值降低一个单位带来的相对的效率下降。 设Ojk为一定时期内由第k个服务单位所创造的第j种产出的观察到的单位的数量。

设Iik为一定时期内由第k个服务单位所使用的第i种投入的实际的单位的数量。

2) 目标函数

目标是找出一组伴随每种产出的系数u和一组伴随每种投入的系数ν，从而给被评估的服务单位最高的可能效率。

（*）

式中，e是被评估单位的代码。这个函数满足这样一个约束条件，当同一组投入和产出的系数（uj和vi）用于所有其他对比服务单位时，没有一个服务单位将超过100%的效率或超过1.0的比率。 3) 约束条件

（**） k=1,2,……,K

式中所有系数值都是正的且非零。

为了用标准线性规划软件求解这个有分数的线性规划，需要进行变形。要注意，目标函数和所有约束条件都是比率而不是线性函数。通过把所评估单位的投入人为地调整为总和1.0，这样等式（*）的目标函数可以重新表述为：

满足以下约束条件：

对于个服务单位，等式（**）的约束条件可类似转化为：

k=1,2,…,K

式中 uj≥0 j=1，2，…，M vi≥0 i=1，2，…，N

关于服务单位的样本数量问题是由在分析种比较所挑选的投入和产出变量的数量所决定的。下列关系式把分析中所使用的服务单位数量K和所考虑的投入种类数N与产出种类数M联系出来，它是基于实证发现和DEA实践的经验：

3、C2R模型 设

hjuTYjvXjT，j1，2，，n

为第j个决策单元的评价指数，总可以选择适当的权系数u、v，使得

hj1，j1，2，，n

第j个决策单元的评价指数hj的意义是：在全系数u、v下，投入为vTXj，产出为uTYj的投入产出比。

如果我们需要考虑某个决策单元j0的效率评价指数hj0为目标，在约束

hj1，j1，2，，n的最大值，即分式线性规划

MAX 4、C2R模型的等价模型

MAX VPuTYjvXjTuTYj0vTXj01

u0，v0VC2RTYj0TXjTYj1Xj10，0Tj1,2,,n

（注：详细参见《优化建模与LINDO、LINGO软件》）

5、DEA使用条件: 1)、所选择的投入指标体系与产出指标体系具有明确的因果关系，以及指标体系满足某些一般公认的要求，如系统性、层次性、全面性等要求。 2）、DEA模型是针对同质的生产部门，对投入规模相当的评价单元的生产有效性进行数量分析，它不直接具备考察不同的投入规模的评价单元之间的比较分析功能

十九、多因素方差分析法（）基于SPSS）

多因素方差分析是对一个独立变量是否受一个或多个因素或变量影响而进行的方差分析。SPSS调用“Univariate”过程，检验不同水平组合之间因变量均数，由于受不同因素影响是否有差异的问题。在这个过程中可以分析每一个因素的作用，也可以分析因素之间的交互作用，以及分析协方差，以及各因素变量与协变量之间的交互作用。该过程要求因变量是从多元正态总体随机采样得来，且总体中各单元的方差相同。但也可以通过方差齐次性检验选择均值比较结果。因变量和协变量必须是数值型变量，协变量与因变量不彼此独立。因素变量是分类变量，可以是数值型也可以是长度不超过8的字符型变量。固定因素变量(Fixed Factor)是反应处理的因素;随机因素是随机地从总体中抽取的因素。 例子:

研究不同温度与不同湿度对粘虫发育历期的影响，得试验数据如表5-7。分析不同温度和湿度对粘虫发育历期的影响是否存在着显著性差异。

表5-7 不同温度与不同湿度粘虫发育历期表 重 复 温度℃ 相对湿度（%） 1 2 3 4 25 91.2 95.0 93.8 93.0 27 87.6 84.7 81.2 82.4 100 29 79.2 67.0 75.7 70.6 31 65.2 63.3 63.6 63.3 25 93.2 89.3 95.1 95.5 27 85.8 81.6 81.0 84.4 80 29 79.0 70.8 67.7 78.8 31 70.7 86.5 66.9 64.9 25 100.2 103.3 98.3 103.8 27 90.6 91.7 94.5 92.2 40 29 77.2 85.8 81.7 79.7 31 73.6 73.2 76.4 72.5 准备分析数据:

在数据编辑窗口中输入数据。建立因变量历期“历期”变量，因素变量温度“A”，湿度为“B”变量，重复变量“重复”。然后输入对应的数值，如图5-6所示。

图5-6 数据输入格式

启动分析过程

点击主菜单“Analyze”项，在下拉菜单中点击“General Linear Model”项，在右拉式菜单中点击“Univariate”项，系统打开单因变量多因素方差分析设置窗口如图5-7。

图5-7 多因素方差分析窗口

设置分析变量:

设置因变量: 在左边变量列表中选“历期”，用“Dependent Variable：”框中。

设置因素变量: 在左边变量列表中选“a”和“b”变量，用

向右拉按钮移向右拉按钮选入到

到“Fixed Factor(s):”框中。可以选择多个因素变量。由于内存容量的限制，选择的因素水平组合数(单元数)应该尽量少。

设置随机因素变量: 在左边变量列表中选“重复”变量，用向右拉按钮移到“到Random Factor(s)”框中。可以选择多个随机变量。

设置协变量：如果需要去除某个变量对因素变量的影响，可将这个变量移到“Covariate(s)”框中。

设置权重变量：如果需要分析权重变量的影响，将权重变量移到“WLS Weight”框中。 选择分析模型

在主对话框中单击“Model”按钮，打开“Univariate Model”对话框。见图5-8。

图5-8 “Univariate Model” 定义分析模型对话框

在Specify Model栏中，指定分析模型类型。 ① Full Factorial选项

此项为系统默认的模型类型。该项选择建立全模型。全模型包括所有因素变量的主效应和所有的交互效应。例如有三个因素变量，全模型包括三个因素变量的主效应、两两的交互效应和三个因素的交互效应。选择该项后无需进行进一步的操作，即可单击“Continue”按钮返回主对话框。此项是系统缺省项。 ② Custom选项

建立自定义的分析模型。选择了“Custom”后，原被屏蔽的“Factors & Covariates”、“Model”和“Build Term(s)”栏被激活。在“Factors & Covariates”框中自动列出可以作为因素变量的变量名，其变量名后面的括号中标有字母“F”；和可以作为协变量的变量名，其变量名后面的括号中标有字母

“C”。这些变量都是由用户在主对话框中定义过的。根据表中列出的变量名建立模型，其方法如下：

在“Build Term(s)”栏右面的有一向下箭头按钮（下拉按钮），单击该按钮可以展开一小菜单，在下拉菜单中用鼠标单击某一项，下拉菜单收回，选中的交互类型占据矩形框。有如下几项选择：

• • • • • •

Interaction 选中此项可以指定任意的交互效应； Main effects 选中此项可以指定主效应； All 2-way 指定所有2维交互效应； All 3-way 指定所有3维交互效应； All 4-way 指定所有4维交互效应 All 5-way 指定所有5维交互效应。

③ 建立分析模型中的主效应：

在“Build Term(s)”栏用下拉按钮选中主效应“Main effects”。 在变量列表栏用鼠标键单击某一个单个的因素变量名，该变量名背景将改变颜色(一般变为蓝色)，单击“Build Term(s)”栏中的右拉箭头按钮，该变量出现在“Model”框中。一个变量名占一行称为主效应项。欲在模型中包括几个主效应项，就进行几次如上的操作。也可以在标有“F”变量名中标记多个变量同时送到“Model”框中。

本例将“a”和“b”变量作为主效应，按上面的方法选送到“Model”框中。 ④ 建立模型中的交互项

要求在分析模型中包括哪些变量的交互效应，可以通过如下的操作建立交互项。

例如，因素变量有“a(F)”和“b(F)”，建立它们之间的相互效应。

•

连续在“Factors &”框的变量表中单击“a(F)”和“b(F)”变量使其选中。

•

单击“Build Term(s)”栏内下拉按钮，选中交互效应“Interaction”项。

•

单击“Build Term(s)”栏内的右拉按钮，“a*b”交互效应就出现在“Model”框中，模型增加了一个交互效应项：a*b

⑤ Sum of squares 栏分解平方和的选择项

•

Type I项，分层处理平方和。仅对模型主效应之前的每项进行调整。一般适用于：平衡的AN0VA模型，在这个模型中一阶交互

效应前指定主效应，二阶交互效应前指定一阶交互效应，依次类推；多项式回归模型。嵌套模型是指第一效应嵌套在第二 效应里，第二效应嵌套在第三效应里，嵌套的形式可使用语句指定。

•

Type II项，对其他所有效应进行调整。一般适用于：平衡的AN0VA模型、主因子效应模型、回归模型、嵌套设计。

•

Type III项，是系统默认的处理方法。对其他任何效应均进行调整。它的优势是把所估计剩余常量也考虑到单元频数中。对没

有缺失单元格的不平衡模型也适用，一般适用于：Type I、Type II所列的模型：没有空单元格的平衡和不平衡模型。

•

Type IV顶，没有缺失单元的设计使用此方法对任何效应F计算平方和。如果F不包含在其他效应里，Type IV = Type IIIl =

TypeII。如果F包含在其他效应里，Type IV只对F的较高水平效应参数作对比。一般适用于：Type I、Type lI所列模型； 没有空单元的平衡和不平衡模型。

⑥ Include intercept in model栏选项

系统默认选项。通常截距包括在模型中。如果能假设数据通过原点，可以不包括截距，即不选择此项。

选择比较方法

在主对话框中单击“Contrasts”按钮，打开“Contrasts”比较设置对话框，如图5-9所示。

如图5-9 Contrasts对比设置框

在“Factors”框中显示出所有在主对话框中选中的因素变量。因素变量名后的括号中是当前的比较方法。 ① 选择因子

在“Factors”框中选择想要改变比较方法的因子，即鼠标单击选中的因子。这一操作使“Change Contrast”栏中的各项被激活。 ② 选择比较方法

单击“Contrast”参数框中的向下箭头，展开比较方法表。用鼠标单击选中的对照方法。可供选择的对照方法有：

• •

None，不进行均数比较。

Deviation，除被忽略的水平外，比较预测变量或因素变量的每个水平的效应。可以选择“Last”(最后一个水平)或

“First”(第一个水平)作为忽略的水平。

•

Simple，除了作为参考的水平外，对预测变量或因素变量的每一水平都与参考水平进行比较。选择“Last”或“First”作为 参考水平。

•

Difference，对预测变量或因素每一水平的效应，除第一水平以外，都与其前面各水平的平均效应进行比较。与Helmert对照 方法相反。

•

Helmert，对预测变量或因素的效应，除最后一个以外，都与后续的各水平的平均效应相比较。

•

Repeated，对相邻的水平进行比较。对预测变量或因素的效应，除第一水平以外，对每一水平都与它前面的水平进行比较。

•

Polynomial，多项式比较。第一级自由度包括线性效应与预测变量或因素水平的交叉。第二级包括二次效应等。各水平彼此 的间隔被假设是均匀的。

③ 修改比较方法

先按步骤①选中因子变量，再选比较方法，然后单击“Change”按钮，选中的(或改变的)比较方法显示在步骤①选中的因子变量后面的括号中。 ④设置比较的参考类

在“Reference Category”栏比较的参考类有两个，只有选择了

“Deviation”或“Simple”方法时才需要选择参考水平。共有两种可能的选择，最后一个水平“Last”选项和第一水平“First”项。系统默认的参考水平是“Last”。 选择均值图

在主对话框中单击“Plot”按钮，打开“Profile Plots”对话框，如图5-10所示。在该对话框中设置均值轮廓图。

如图5-10 “Profile Plots”对话框

均值轮廓图(Profile Plots)用于比较边际均值。轮廓图是线图，图中每个点表明因变量在因素变量每个水平上的边际均值的估计值。如果指定了协变量，该均值则是经过协变量调整的均值。因变量做轮廓图的纵轴；一个因素变量做横轴。 做单因素方差分析时，轮廓图表明该因素各水平的因变量均值。

双因素方差分析时，指定一个因素做横轴变量，另一个因素变量的每个水平产生不同的线。如果是三因素方差分析，可以指定第三个因素变量，该因素每个水平产生一个轮廓图。双因素或多因素轮廓图中的相互平行的线表明在因素间无交互效应；不平行的线表明有交互效应。

• •

Factors 框中为因素变量列表。

Horlzontal Axis 横坐标框，选择选择“Factors”框中一个因素变量做横坐标变量。被选的变量名反向显示，单击向右拉箭

头按钮，将变量名送入相应的横坐标轴框中。 如果只想看该因素变量各水平的，因变量均值分布，单击“Add”按钮，将所选因素变量移入下面的“Plots”框中。否 则，不点击“Add”按钮，接着做下步。

•

Separate Lines 分线框。如果想看两个因素变量组合的各单元格中因变量均值分布，或想看两个因变量间是否存在交互效应，

选择“Factors”框中另一个因素变量，单击右拉按钮将变量名送入“Separate Lines”框中。单击“Add”按钮，将自动生成

的图形表达式送入到“Plots”栏中。分线框中的变量的每个水平将在图中是一条线。图形表达式是用“*”连接的两个因素变 量名。

•

Separate Plots 分图框。如果在“Factors”栏中还有因素变量，可以按上述方法，将其送入“Separate Plot”框中，单击

“Add”按钮，将自动生成的图形表达式送入到“Plots”栏中。图形表达式是用“*’连接的三个因素变量名。分图变量的每个 水平生成一张线图。

•

将图形表达式送到“Plots”框后发现有错误，单击选错的变量，单击“Remove”按钮，将其取消，再重新输入正确内容。

在检查无误后，按“Continue”按钮确认，返回到主对话框。如果取消做的设置单击“Cancel”按钮 选择多重比较

在主对话框中单击“Post Hoc”选项，打开“Post Hoc Multiple Comparisons for Observed Means”对话框，从“Factor(s)”框选择变量，单击向右拉按钮，使被选变量进入“Post Hoc test for”框。本例子选择了“a”和“b”。 然后选择多重比较方法。在对话框中选择多重比较方法。本例子选择了“Duncan”和“Tamhane's T2”。 选择保存运算值

图5-11 Save对话框

在主对话框中，单击“Save”按钮，打开“Save”设置对话框，如图5-11所示。通过在对话框中的选择，可以将所计算的预测值、残差和检测值作为新的变量保存在编辑数据文件中。以便于在其他统计分析中使用这些值。 ① Predicted Values 预测值 Unstsndardized，非标准化预测值。

Weighted，如果在主对话框中选择了WLS变量，选中该复选项，将保存加权非标准化预测值。

Standard error，预测值标准误。 ② Diagnostics 诊断值

Cook’s distance，Cook 距离。

Leverage values，非中心化 Leverage 值。 ③ Residuals 残差

Unstsndardized，非标准化残差值，观测值与预测值之差。

Weighted，如果在主对话框中选择了WLS变量，选中该复选项，将保存加权非标准化残差。

Standardized，标准化残差，又称Pearson残差。 Studentized，学生化残差。

Deleted，剔除残差，自变量值与校正预测值之差。 ④ Save to New File 保存协方差矩阵

选中”Coefficient statistics”项，将参数协方差矩阵保存到一个新文件中。单击“File”按钮，打开相应的对话框将文件保存。 选择输出项

在主对话框中单击“Options”按钮，打开“Options”输出设置对话框，见图5-12。

图5-12 “Options”输出设置对话框

① Estimated Marginal Means 估测边际均值设置

•

在“Factor(s) and Factor Interactions”框中列出“Model”对话框中指定的效应项，在该框中选定因素变量的各种效应项，

单击右拉按钮就将其复制到“Display Means for”框中。选择主

效应，则产生估计的边际均值表；选择二维交互效应产生的估计 边际均值表实际上是典型的单元格均值表。选择三维交互效应也是单元格均值表。

•

在“Display Means for”框中有主效应时激活此框下面的“Compare main effects”复选项，对主效应的边际均值进行组间的配 对比较。

•

Confidence interval adjustment参数框，进行多重组间比较。打开下拉菜单，共有三个选项：

LSD(none)、Bonferroni、Sidak.。

② 在“Display”栏中指定要求输出的统计量

Descriptive statistics项，输出描述统计量：观测量的均值、标准差和每个单元格中的观测量数。

Estimates of effect size项，效应量估计。选择此项，给出η2(eta-Square)值。它反应了每个效应与每个参数估计值可以归于 因素的总变异的大小。

Observed power复选项，选中此项给出在假设是基于观测值时各种检验假设的功效。计算功效的显著性水平，系统默认的临界值 是0.05。

Parameter estimates项。选择此项给出了各因素变量的模型参数估计、标准误、t检验的t值、显著性概率和95％的置信区间。

Contrast coefficient matrix项，显示协方差矩阵。 Homogeneity test项，方差齐次性检验。本例子选中该项。

Spread vs.level plot项，绘制观测量均值对标准差和观测量均值对方差的图形。

Residual plot项，绘制残差图。给出观测值、预测值散点图和观测量数目，观测量数目对标准化残差的散点图，加上正态和标准化 残差的正态概率图。

Lack of fit项，检查独立变量和非独立变量间的关系是否被充分描述。 General estimable function项，可以根据一般估计函数自定义假设检验。对比系数矩阵的行与一般估计函数是线性组合的。 ③ Significance level 框设置

改变“Confidence intervals”框内多重比较的显著性水平 10）提交执行

设置完成后，在多因素方差分析窗口框中点击“OK”按钮，SPSS就会根据设置进行运算，并将结算结果输出到SPSS结果输出窗口中。

11）结果与分析 主要输出结果：

结果分析:

方差不齐次性检验显著

表5-8 方差齐次性检验表明：方差不齐次性显著，p<0.05。 方差分析:

表5-9 主效应方差分析表：在表的左上方标明研究的对象是粘虫历期。

偏差来源和偏差平方和:

•

Source 列是偏差的来源。其次列是“Type III Sum of Squares”偏差平方和。

•

Corrected Model 校正模型，其偏差平方和等于两个主效应a、b平方和加上交互a*b的平方和之和。

• •

Intercept 截距。

a 温度主效应，其偏差平方和反应的是不同温度造成对粘虫历期的差异。与b偏差平方相同均属于组间偏差平方和。

•

b 湿度主效应，其偏差平方和反应的是不同湿度计量造成的粘虫历期之差异。

•

a*b 温度和湿度交互效应，其偏差平方和反应的是不同温度和湿度共同造成的粘虫历期的差异。

• •

Error 误差。其偏差平方和反应的是组内差异。也称组内偏差平方和。 Total 是偏差平方和在数值上等于截距、主效应、次效应和误差偏差平方和之总和。

•

Corrected Total 校正总和。其偏差平方和等于校正模型与误差之偏差平方和之总和。

• • • •

df 自由度

Mean Square 均方，数值上等于偏差平方和除以相应的自由度。 F 值，是各效应项与误差项的均方之比值

Sig 进行F检验的p值。p≤0.05，由此得出“温度”和“湿度”对因变量“粘虫历期”在0.05水平上是有显著性差异的。

根据方差分析表明：

•

不同温度（a）对粘虫历期的偏差均方是1575.434，F值为90.882，显著性水平是0.000，即p<0.05存在显著性差异；

•

不同湿度（b）对粘虫历期的偏差均方是322.000，F值为18.575，显著性水平是0.000，即p<0.05存在显著性差异；

•

不同温度和不同湿度（a*b）共同对粘虫历期的偏差均方是19.809，F值为1.143，显著性水平是0.358，即p>0.05存在不显著性 差异。

多重比较

由于方差不齐次性，应选择方差不具有齐次性时的“Tamhane's T2”t检验进行配对比较。表5-10 多重比较表就是“温度”各水平“Tamhane's T2”方法比较的结果。表中的各项说明参见表5-6（5.2.2节）。

•

温度25℃与27℃、29℃和31℃之间都有显著性差异；

• • •

温度27℃与25℃、29℃和31℃之间都有显著性差异；

温度29℃与26℃和27℃之间都有显著性差异；与31℃无显著性差异； 温度31℃与25℃和27℃之间都有显著性差异；与29℃无显著性差异。

不同湿度水平之间无显著性差异存在，这里没有列出多重比较表。

二十、拉格朗日插值

注:参考<数学建模与数学实验>

5.2.1 基函数

由上一节的证明可以看到，要求插值多项式pn(x)，可以通过求方程组（5.1.22）的解a0,a1,...an得到，但这样不但计算复杂，且难于得到pn(x)的简单表达式。

考虑简单的插值问题：设函数在区间［a，b］上n＋1个互异节点x0,x1,...xn上的函数值为

1,jiyjij0,ji, j=0,1,…,n

求插值多项式li(x)，满足条件

j=0,1,…n, i=0,1,…,n

由上式知，x0,x1,...,xi1,xi1,...,xn是li(x)的根，且li(x)∈Hn，可令

li(x)ij

li(x)Ai(xx0)(xx1)...(xxi1)(xxi1)...(xxn) 再由li(xi)1得

1Ai(xix0)(xix1)...(xixi1)(xixi1)...(xixn)

于是

(xx0)(xx1)...(xxi1)(xxi1)...(xxn)li(x)(xix0)(xix1)...(xixi1)(xixi1)...(xixn)

n+1个n次多项式l0(x),l1(x),...,ln(x)称为以为x0,x1,...xn节点的n次插值基函数。

n=1时的一次基函数为（图5-2）：

xx0xx1l0(x),l1(x)x0x1 x1x0 .

n=2时的二次基函数为（图5-3）：

(xx1)(xx2)l0(x)(x0x1)(x0x2)(xx0)(xx2)l1(x)(x1x0)(x1x2)(xx0)(xx1)l2(x)(x2x0)(x2x1)

5-2

5-3

5.2.2 拉格朗日插值多项式

现在考虑一般的插值问题：设函数在区间［a,b］上n+1个互异节点x0,x1,...xn上的函数值分别为,y0,y1,...yn，求n次插值多项式pn(x)，满足条件 令

pn(xj)yj, j=0,1,…n

n（5.2.3）

其中l0(x),l1(x),...,ln(x)为以x0,x1,...xn为节点的n次插值基函数，则Ln(x)是一次

i0Ln(x)y0l0(x)y1l1(x)...ynln(x)yili(x)数不超过n的多项式，且满足

Ln(xj)yj， j=0，1，…，n

再由插值多项式的唯一性，得

pn(x)Ln(x)

式（5.2.3）表示的插值多项式称为拉格朗日（Lagrange）插值多项式。特别地，n=1 时称为线性插值（图5-4（a）），n=2时称为抛物插值或二次插值（图5-4（b））。

值得注意的是，插值基函数l0(x),l1(x),...,ln(x)仅由插值节点x0,x1,...xn确定，与被插函数f(x)无关。因此，若以 x0,x1,...xn为插值节点对函数f(x)≡1作 插值多项式，则由式(5.2.3)立即得到基函数的一个性质

≡1

还应注意，对于插值节点x0,x1,...xn，只要求它们互异，与大小次序无关。

i0l(x)in

5-4

例1 已知y=x，x0=4，x1=9，用线性插值求7的近似值。解 y0=2，y1=3，基函数分别为

l0(x)

x41x41(x9),l1(x)(x4)495945

插值多项式为

11L1(x)y0l0(x)y1l1(x)2(x9)3(x4)551(x6) 5 所以

137L1(7)2.65

例2 求过点(-1，-2)，(1，0)，(3，-6)，(4，3)的三次插值多项式。 解 以x0＝-1，x1＝1，x2＝3，x3＝4为节点的基函数分别为

(x1)()x3)(x4)1(x1)(x3)(x4)(11)(13)(14)40(x1)(x3)(x4)1l1(x)(x1)(x3)(x4)(11)(13)(14)12(x1)(x1)(x4)1l2(x)(x1)(x1)(x4)(31)(31)(34)8(x1)(x1)(x3)1l3(x)(x1)(x1)(x3)(41)(41)(43)15

插值多项式为

l0(x)L3(x)yili(x)i1311(x1)(x3)(x4)0(x1)(x3)(x4)401211(6)(x1)(x1)(x4)3(x1)(x1)(x3)81532 x4x3 5.2.3 插值余项

(2)插值多项式的余项Rn(x)=f(x)-Ln(x)，也就是插值的截断误差或方法误差。关于余项有 如下的余项定理：

(n1)(x)在开定理 设被插函数f(x)在闭区间［a，b］上n阶导数连续，f区间(a,b)

内存在，x0,x1,...xn是［a，b］上n+1个互异节点，记

ni0

n1(x)(xxi)(xx0)(xx1)...(xxn)

则插值多项式Ln(x)的余项为

f(n1)()Rn(x)f(x)Ln(x)n1(x),x[a,b](n1)!其中(x)(a,b). （5.2.4）

证明 由插值条件和n1(x)的定义，当x=xk时式（5.2.4）显然成立， 并且有

Rn(xk)0 k=0,1,…，n （5.2.5） 这表明x0,x1,...xn都是函数Rn(x)的零点，从而Rn(x)可表示为 Rn(x)=f(x)- Ln(x)=K(x) n1(x) (5.2.6) 其中K(x)是待定函数。

对于任意固定的x［a,b］，x≠ xk (k=0,1,…，n)，构造自变量t的辅助函数

(t)f(t)Ln(t)K(x)n1(t) (5.2.7) 由式（5.2.5）和式（5.2.6）可知x0,x1,...xn和x是(t)在区间［a,b］上的n+2个互异 零点，因此，根据罗尔（Rolle）定理，至少存在一点(x)∈(a,b)，使得

(n1)()0于是，由式（5.2.7）得到

f(n1)()K(x)(n1)!

代入式（5.2.6）即得式（5.2.4）。

由于(x)一般无法确定，因此式(5.2.4)只能用作余项估计。如果

f(n1)(x)在区 间(a,b)上有界，即存在常数Mn1＞0，使得

(n1)f(x)|≤ Mn1 x∈(a,b) |

|Rn(x)|Mn1|n1(x)|(n1)!

则有余项估计

(5.2.8)

max|f(n1)(x)|x[a,b]当f(n1)(x)在闭区间［a,b］上连续时，可取

Mn1。

推论 设节点x0＜x1, f″(x)在闭区间［x0,x1］上连续，记

M2maxx(a,b)｜f″(x)｜，则过点(x0,f(x0))，(x1,f(x1))的线性插 值余项为f''(x)R1(x)(xx0)(xx1),(x)(x0,x1)2 (5.2.9)

(x1x0)(x1x0)242由于在［x0, x1］上，｜(x-x0)(x-x1)｜在x=达到 最大值，

可得余项的一个上界估计：x∈［x0,x1］有

|R1(x)|M2(x1x0)283121在本节例1中f''(x)x,max|f\"(x)|,由上式432x[4,9]11|R1(7)||7L1(7)|..(94)20.097656258321例三设f(x),节点x02,x12.5,x34,求f(x)的抛物插值多项式L2(x)，x

且计算f(3)的近似值并估计误差。

解 由于y0f(2)0.5,y1f(2.5)0.4,y2f(4)0.25，插值多项式为

(x2.5)(x4)(x2)(x4)(x2)(x2.5)L2(x)0.50.40.25(22.5)(24)(2.52)(2.54)(42.5)(42)0.05x20.425x1.15

于是

f(3)≈L2(x)=0

因为

f\"(x)325

代入式(5.2.8)得

|R2(3)||f(3)L2(3)|63,|f'''(x)||f'''(2)|8 x4maxx[2,4]例4

字。

(1)用线性插值求sin0.33的近似值；

(2)证明在区间［0.32，0.34］上用线性插值计算sinx时至少有4位有效数字。

解 (1)用线性插值

0.330.340.330.32sin0.33L1(0.33)0.3145670.3334870.320.340.340.321(0.3145670.333487)0.3240272 (2)由式(5.2.10)在区间［0.32，0.34］用线性插值计算sinx时的余项

13.|(32)(32.5)(34)|0.0312568

已知sin0.32=0.314567，sin0.34=0.333487有六位有效数

满足

|R1(x)|0.333487(0.340.32)20.5104,x[0.32,0.34]8 因此结果至少有4位有效数字。

二十一、回归分析（略）

二十二、概率分布方法（略） 二十三、插值与拟合（略）

二十四、隶属函数的刻画（参考《数学建模及其方法应用》） 二十五、0-1整数规划模型（参看书籍） 二十六、Board评价法（略） 二十七、纳什均衡（参看书籍）

二十八、微分方程方法与差分方程方法（参看数据） 二十九、莱斯利离散人口模型（参看数据）

三十、一次指数平滑预测法（主要是软件的使用） 三十一、二次曲线回归方程（主要是软件的使用） 三十二、成本-效用分析（略）

三十三、逐步回归法（主要是软件的使用） 三十四、双因子方差分析（略）