数据处理专题 数据处理是指用简明而严格的方法把获得的实验数据所代表的事物内在的规律提炼出来，得出结果的加工过程，包括数据记录、描绘曲线，从带有误差的数据中提取参数，验证和寻找经验规律，外推实验数据等等。本章介绍一些最基本的数据处理方法。

数据处理的过程： 1、获得数据（标准化处理）。 2、将数据分类（聚类分析）。 3、提取主要影响因素（ 主成分分析）。 4、数据分析（相关性分析，回归分析）。

系统聚类法、有序样品聚类法、动态聚类法、模糊聚类法、图论聚类法、聚类预报法等。我们着重讲述系统聚类法。对样品分类成Q分类，对指标分类称R分类 聚类分析 聚类也就是分类，在社会经济领域中存在大量的分类问题，比如三十个省市自治区独立核算工业企业经济效益进行分析，一般不是逐个省市自治区去分析，而较好的做法是选取具有代表性的指标如，百元固定资产实现利税，资金利税率、产值利税率、百元销售收入实现利润、全员劳动生产率等等，根据这些指标对省市自治区分类，然后根据分类结果对企业经济效益进行综合评价。 聚类分析方法 系统聚类法、有序样品聚类法、动态聚类法、模糊聚类法、图论聚类法、聚类预报法等。我们着重讲述系统聚类法。对样品分类成Q分类，对指标分类称R分类

我们通过距离来分类。方法有:最短距离法、最长距离法、中间距离法、重心法等。我们用最短距离法来讲述，其它方法读者自己翻阅相关的多元统计教材。 聚类的三种尺度： 1、间隔尺度：变量是用连续量来表示，如长度、重量等 2、有序尺度：用一些等级来表示。如上中下三等。 3、名义尺度：既没有数量表示也没有次序表示。如红黄蓝三色等 我们通过距离来分类。方法有:最短距离法、最长距离法、中间距离法、重心法等。我们用最短距离法来讲述，其它方法读者自己翻阅相关的多元统计教材。

其中，Xi 为样品的p个指标组成的向量。 最短距离法步骤如下： 【1】定义样品之间的距离，计算样品两两距离，得一距离记为D(0) 开始每个样品自成一类，显然这时Dij =dij。其中D表示类G之间的距离，d表示样品之间的距离。 【2】找出D(0) 的非对角线最小元素，设为Dpq，则将Gp和Gq合并为一新类，记为Gr 。 【3】给出计算新类与其他的类的距离公式： 距离公式有：欧氏距离，马氏距离，兰氏距离等。我们一般用马氏距离，应为它即排除了各指标之间相关性的干扰，而且还不受各指标量纲的影响。 两个样本间的距离定义： 其中，Xi 为样品的p个指标组成的向量。 协方差阵的逆矩阵

协方差阵定义如下： 样品到总体的距离定义： 总体均值向量

Dkr=min{Dkp，Dkq}将D(0)中的第p、q行及p、q列用上面公式并成一个新行新列，新行新列对应Gr，所得到得矩阵记为D(1) 【4】对D(1)重复上述对D(0)的(2)(3)两步得D(2);如此下去，直到所有的元素并为一类。 注意：如果某一步中非对角线最小的元素不止一个，则对应这些最小元素的类可以同时合并。 为了大家便于掌握我们举例如下：

例：设抽取五个样品，每个样品只测一个指标，它们是1，2，3.5，7，9，试用最短距离法对这五个样品进行分类。 解：我们距离选用我们所熟悉的绝对值距离。 G1={X1} G2={X2} G3={X3} G4={X4} G5={X5} 1 2.5 1.5 6 5 3.5 8 7 5.5 2

G6={X1,x2} G3={X3} G4={X4} G5={X5} 1.5 5 3.5 7 5.5 2 G6={X1,x2,x3} G4={X4} G5={X5} 3.5 5.5 2 G6={X1,x2,x3} G7={x4,X5} 2 最终我们分为两类比较合适，{x1,x2,x3}与{x4,x5}

matlab做聚类分析 分步聚类：（1）找到数据集合中变量两两之间的相似性和非相似性，用pdist函数计算变量之间的距离；（2）用 linkage函数定义变量之间的连接；（3）用 cophenetic函数评价聚类信息；（4）用cluster函数创建聚类。 Step1 寻找变量之间的相似性 用pdist函数计算相似矩阵，有多种方法可以计算距离，进行计算之前最好先将数据用zscore函数进行标准化。 X=[1，2，3.5，7，9] X2=zscore(X); %标准化数据 Y2=pdist(X2); %计算距离 Step2 定义变量之间的连接 Z2=linkage(Y2); Step3 评价聚类信息 C2=cophenet(Z2,Y2); //0.94698 Step4 创建聚类，并作出谱系图 T=cluster(Z2,2); H=dendrogram(Z2);%画出聚类图

H=dendrogram(Z2);%画出聚类图 例 为了更深入了解我国人口的文化程度状况，1990年全国人口普查数据对全国30个省直辖市、自治区进行聚类分析。分析选用了三个指标：【1】大学以上文化程度的人口占全部人口的比例(DXBZ)；【2】初中以上文化程度的人口占全部人口的比例(CZBZ)；【3】文盲半文盲的人口占全部人口的比例(WMBZ);分别用来反映较高、中等、较低文化程度人口的状况，原始数据如附件： clear clc X=load('data1.txt') Y2=pdist(X);%计算距离 Z2=linkage(Y2); C2=cophenet(Z2,Y2); T=cluster(Z2,4); H=dendrogram(Z2);%画出聚类图

pdist函数 调用格式：Y=pdist(X,’metric’) 说明：用 ‘metric’指定的方法计算 X 数据矩阵中对象之间的距离。 X：一个m×n的矩阵，它是由m个对象组成的数据集，每个对象的大小为n。 metric’取值如下： ‘euclidean’：欧氏距离（默认）；‘seuclidean’：标准化欧氏距离； ‘mahalanobis’：马氏距离；‘cityblock’：布洛克距离； ‘minkowski’：明可夫斯基距离；‘cosine’： ‘chebychev’：Chebychev距离。 linkage函数 调用格式：Z=linkage(Y,’method’) 说 明：用‘method’参数指定的算法计算系统聚类树。 Y：pdist函数返回的距离向量； method：可取值如下： ‘single’：最短距离法（默认）； ‘complete’：最长距离法； ‘average’：未加权平均距离法； ‘weighted’： 加权平均法； ‘centroid’：质心距离法； ‘median’：加权质心距离法； ‘ward’：内平方距离法（最小方差算法）

练习题 根据信息基础设施的发展状况，对二十个国家的地区进行分类。

主成分分析 在实际问题中，研究多指标的问题是经常遇到的，然而在多数情况下，不同指标之间是有一定关系的。由于指标较多再加上指标之间有一定的相关性，势必增加了分析问题的复杂性。主成分分析就是设法将原来指标重新组合成一组新的互相无关的几个综合指标来代替原来指标，同时根据实际需要从中可取几个较少的综合指标尽可能多滴反映原来指标的信息。这种多个指标化为少数互不干扰的综合指标的统计方法叫做主成分分析法，如某人要做一件上衣要测量很多尺寸，如身长、袖长、胸围、腰围、肩宽、肩厚等十几项指标。但是某服装产生产一批新型服装绝不可能吧尺寸型号分的过多。而是从其中选取几个综合性的指标作为分类型号。1、反映胖瘦。2、反映特体。3反映长度。

设有n个样品，每个样品观测p个指标，将原始数据写成矩阵形式 计算步骤 设有n个样品，每个样品观测p个指标，将原始数据写成矩阵形式 特征值大的贡献大。 贡献率=特征值/所有特征值和 1、将原始数据标准化 2、建立变量的相关系数阵 3、求R的特征根及相应的单位特征向量a1,a2,.....ap 4、写出主成分 一般取累计贡献率达85—95%的特征值 所对应的第一、第二，…，第m（m≤p）个主成分。

例 中国大陆35个大城市某年的10项社会经济统计指标指标做主成分分析数据见下表。 pcacov 功能：运用协方差矩阵进行主成分分析 格式：PC=pcacov(X) [PC,latent,explained]=pcacov(X) 说明：[PC,latent,explained]=pcacov(X)通过协方差矩阵X进行主成分分析，返回主成分(PC)、协方差矩阵X的特征值(latent)和每个特征向量表征在观测量总方差中所占的百分数(explained)。 例 中国大陆35个大城市某年的10项社会经济统计指标指标做主成分分析数据见下表。 相关系数矩阵: std = 1.0000 -0.3444 0.8425 0.3603 0.7390 0.6215 0.4039 0.4967 0.6761 0.4689 -0.3444 1.0000 -0.4750 0.3096 -0.3539 0.1971 0.3571 0.2600 0.1570 0.3090 0.8425 -0.4750 1.0000 0.3358 0.5891 0.5056 0.3236 0.4456 0.5575 0.3742 0.3603 0.3096 0.3358 1.0000 0.1507 0.7664 0.9412 0.8480 0.7320 0.8614 0.7390 -0.3539 0.5891 0.1507 1.0000 0.4294 0.1971 0.3182 0.3893 0.2595 0.6215 0.1971 0.5056 0.7664 0.4294 1.0000 0.8316 0.8966 0.9302 0.9027 0.4039 0.3571 0.3236 0.9412 0.1971 0.8316 1.0000 0.9233 0.8376 0.9527 0.4967 0.2600 0.4456 0.8480 0.3182 0.8966 0.9233 1.0000 0.9201 0.9731 0.6761 0.1570 0.5575 0.7320 0.3893 0.9302 0.8376 0.9201 1.0000 0.9396 0.4689 0.3090 0.3742 0.8614 0.2595 0.9027 0.9527 0.9731 0.9396 1.0000

特征值(val) val = 0.0039 0 0 0 0 0 0 0 0 0 0 0.0240 0 0 0 0 0 0 0 0 0 0 0.0307 0 0 0 0 0 0 0 0 0 0 0.0991 0 0 0 0 0 0 0 0 0 0 0.1232 0 0 0 0 0 0 0 0 0 0 0.2566 0 0 0 0 0 0 0 0 0 0 0.3207 0 0 0 0 0 0 0 0 0 0 0.5300 0 0 0 0 0 0 0 0 0 0 2.3514 0 0 0 0 0 0 0 0 0 0 6.2602 特征根排序： 6.26022 2.35138 0.530047 0.320699 0.256639 0.123241 0.0990915 0.0307088 0.0240355 0.00393387

特征向量(vec)： -0.1367 0.2282 -0.2628 0.1939 0.6371 -0.2163 0.3176 -0.1312 -0.4191 0.2758 -0.0329 -0.0217 0.0009 0.0446 -0.1447 -0.4437 0.4058 -0.5562 0.5487 0.0593 -0.0522 -0.0280 0.2040 -0.0492 -0.5472 -0.4225 0.3440 0.3188 -0.4438 0.2401 0.0067 -0.4176 -0.2856 -0.2389 0.1926 -0.4915 -0.4189 0.2726 0.2065 0.3403 0.0404 0.1408 0.0896 0.0380 -0.1969 -0.0437 -0.4888 -0.6789 -0.4405 0.1861 -0.0343 0.2360 0.0640 -0.8294 0.0377 0.2662 0.1356 -0.1290 0.0278 0.3782 0.2981 0.4739 0.5685 0.2358 0.1465 -0.1502 -0.2631 0.1245 0.2152 0.3644 0.1567 0.3464 -0.6485 0.2489 -0.4043 0.2058 -0.0704 0.0462 0.1214 0.3812 0.4879 -0.5707 0.1217 0.1761 0.0987 0.3550 0.3280 -0.0139 0.0071 0.3832 -0.7894 -0.1628 0.1925 0.2510 -0.0422 0.2694 0.0396 0.0456 0.1668 0.3799 于是的三个指标为： Y1=-0.1312*x1-0.5562*x2+0.3188*x3+....+ 0.0456*x10 Y2=-0.4191*x1+0.5487*x2+...........+0.1668*x10 Y3=0.2758*x1+ 0.0593*x2+............+0.3799*x10

通过观察我们发现Y1当中x2,x5的系数比较大，即影响Y1比较明显因此我们可将Y1看做反映非农业人口比与客运总量的综合指标。 进一步还可做因子分析。 练习、我们给出了各地的企业的经济效益状况，通过相关的方法对各地的经济效益做分析。数据如下表：

相关性分析 在一元统计分析中，研究两入随机变量之间的线性相关关系、 可用相关系数(称为简单相关系数)；研究一个随机变量与多个随 机变量之间的线性相关关系，可用复相关系数(称为全相关系 数)将它推广到研究多个随机变量与多个随机变量之间的相关关系的讨论中， 提出了典型相关分析。 实际问题中，两组变量之间具有相关关系的问题很多，例如 几种主要产品如猪肉、牛肉、鸡蛋的价格(作为第一组变量)和 相应这些产品的销售量(作为第二组变量)有相关关系；投资性 变量(如劳动各人数、货物周转量、生产建设投资等)与国民收 入变量(如工农业国民收入、运输业国民收入、建筑业国民收入 等)只有相关关系；患某种疾病的病人的各种症状程度(第一组 变量)和用物理化学方法检验的结果(第二组变量)具有相关关 系；运动员的体力测试指标(如反复横向跳、纵跳、背力、握力 等)与运动能力测试指标(如耐力跑、跳远、投球等)之间具有 相关关系等等。

典型相关分析就是研究两组变量之间相关关系的一种多元统计方法，设两组变量用x1,x2,…xn和y1,y2…yn表示，要研究两组变量的相关关系，一种方法是分别研究X和Y之间的相关关系，然后列出相关系数表进行分析，当两组变量较多时，这样做法不仅烦琐也不易抓住问题的实际；另一种方法采用 类似主成分分析的做法在每一组变量中都选择若干个有代表性的综合指标(变量的线性组合)，通道研究两组的综合指标之间的关系来反映两组变量之间关系比如猪肉价格和牛肉价格用x1，X2表示，它们的销售售量用X，xl表示，研究它们之间的相又关系，从经济学观点就是希望构造一个X1、x2的线性函数入y1＝a11X1十a12x2称为价格指数及x3、x4的线性函数y2＝a21x3十a22X4称为销售指数，要求它们之间具有最大相关性，这就是一个典型相关分析问题。

1.插值拟合 2.线性回归 4.灰色分析 5.神经网络

在解决实际问题的生产（或工程）实践和科学实验过程中，通常需要通过研究某些变量之间的函数关系来帮助我们认识事物的内在规律和本质属性，而这些变量之间的未知函数关系又常常隐含在从试验、观测得到的一组数据之中。因此，能否根据一组试验观测数据找到变量之间相对准确的函数关系就成为解决实际问题的关键。 例如在工程实践和科学实验中，常常需要从一组试验观测数据(xi ,yi ) ，i = 0,1,....,n之中找到自变量x与因变量y 之间的函数关系，一般可用一个近似函数y = f (x)来表示。函数y = f (x)的产生办法因观测数据和要求不同而异，通常可采用数据拟合与函数插值两种办法来实现。

数据拟合主要是考虑到观测数据受随机观测误差的影响，进而寻求整体误差最小、能较好反映观测数据的近似函数y = f (x)，此时并不要求所得到的近似函数y = f (x)满足yi= f (xi) ， i = 0,1,…,n。 函数插值则要求近似函数y = f (x)在每一个观测点i x 处一定要满足y i= f (xi) ， i = 0,1,…,n ，在这种情况下，通常要求观测数据相对比较准确，即不考虑观测误差的影响。

在实际问题中，通过观测数据能否正确揭示某些变量之间的关系，进而正确认识事物的内在规律与本质属性，往往取决于两方面因素。其一是观测数据的准确性或准确程度，这是因为在获取观测数据的过程中一般存在随机测量误差，导致所讨论的变量成为随机变量。其二是对观测数据处理方法的选择，即到底是采用插值方法还是用拟合方法，插值方法之中、拟合方法之中又选用哪一种插值或拟合技巧来处理观测数据。插值问题忽略了观测误差的影响，而拟合问题则考虑了观测误差的影响。但由于观测数据客观上总是存在观测误差，而拟合函数大多数情况下是通过经验公式获得的，因此要正确揭示事物的内在规律，往往需要对大量的观测数据进行分析，尤为重要的是进行统计分析。统计分析的方法有许多，如方差分析、回归分析等。

数据拟合虽然较有效地克服了随机观测误差的影响，但从数理统计的角度看，根据一个样本计算出来的拟合函数（系数），只是拟合问题的一个点估计，还不能完全说明其整体性质。因此，还应该对拟合函数作区间估计或假设检验，如果置信区间太大或包含零点，则由计算得到的拟合函数系数的估计值就毫无意义。这里所采用的统计分析方法就是所谓的回归分析。另外还可用方差分析的方法对模型的误差作定量分析。 对于插值方法，本章简单介绍最常用的插值法的基本结论及其Matlab实现问题。由于数据拟合问题必须作区间估计或假设检验，所以除了在本章介绍最基本的数据拟合方法——最小二乘法的基本结论及其Matlab实现问题外，我们在专门介绍了对数值拟合问题进行区间估计或假设检验的统计方法，

即介绍回归分析方法及其Matlab实现。

插值方法 1、拉格朗日插值法 2、分段线性插值法 分段线性插值的Matlab实现 用Matlab实现分段线性插值不需要编制函数程序，Matlab中有现成的一维插值函数interp1。 y=interp1(x0,y0,x,'method') method指定插值的方法，默认为线性插值。其值可为： 'nearest' 最近项插值 'linear' 线性插值 'spline' 立方样条插值 'cubic' 立方插值。 2、分段线性插值法

对于三次样条插值，我们提倡使用函数csape，csape的返回值是pp形式，要求插值点的函数值，必须调用函数ppval。 Matlab中三次样条插值也有现成的函数： y=interp1(x0,y0,x,'spline')； y=spline(x0,y0,x)； pp=csape(x0,y0,conds), pp=csape(x0,y0,conds,valconds)，y=ppval(pp,x)。 其中x0,y0是已知数据点，x是插值点，y是插值点的函数值。 3、三次样条插值法 对于三次样条插值，我们提倡使用函数csape，csape的返回值是pp形式，要求插值点的函数值，必须调用函数ppval。

例1 机床加工 待加工零件的外形根据工艺要求由一组数据(x, y)给出（在平面情况下），用程控铣床加工时每一刀只能沿x方向和y 方向走非常小的一步，这就需要从已知数据得到加工所要求的步长很小的(x, y)坐标。 表中给出的x, y数据位于机翼断面的下轮廓线上，假设需要得到x坐标每改变0.1时的y坐标。试完成加工所需数据，画出曲线，并求出x = 0处的曲线斜率和13 ≤ x ≤ 15范围内y 的最小值。 x 0 3 5 7 9 11 12 13 14 15 y 0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6 要求用分段线性和三次样条two种插值方法计算。

x0=[0 3 5 7 9 11 12 13 14 15]; y0=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6]; x=0:0.1:15; y2=interp1(x0,y0,x,’linear’); y3=interp1(x0,y0,x,'spline'); pp1=csape(x0,y0); y4=ppval(pp1,x); pp2=csape(x0,y0,'second'); y5=ppval(pp2,x); subplot(2,2,2) plot(x0,y0,'+',x,y2) title('Piecewise linear') subplot(2,2,3) plot(x0,y0,'+',x,y3) title('Spline1') subplot(2,2,4) plot(x0,y0,'+',x,y4) title('Spline2') dx=diff(x); dy=diff(y3); dy_dx=dy./dx; dy_dx0=dy_dx(1) ytemp=y3(131:151); ymin=min(ytemp); index=find(y3==ymin); xmin=x(index); [xmin,ymin] 计算结果略。 可以看出，分段线性插值的光滑性较差（特别是在x = 14附近弯曲处），建议选用三次样条插值的结果。

五 一维插值总结 插值函数一般是已知函数的线性组合或者称为加权平均。在已知数据点较少时，插值技术在工程实践和科学实验中有着广泛而又十分重要的应用。例如在信息技术中的图像重建、图像放大过程中为避免图像失真、扭曲而增加的插值补点，建筑工程的外观设计，化学工程试验数据与模型分析，天文观测数据、地理信息数据的处理，社会经济现象的统计分析等方面，插值技术的应用是不可或缺的。 插值技术（或方法）远不止这里所介绍的这些，但在解决实际问题时，对于一维插值问题而言，前面介绍的插值方法已经足够了。剩下的问题关键在于什么情况下使用、怎样使用和使用何种插值方法的选择上。 拉格朗日插值函数在整个插值区间上有统一的解析表达式，其形式关于节点对称，光滑性好。但缺点同样明显，这主要体现在高次插值收敛性差（龙格现象）；增加节点时前期计算作废，导致计算量大；一个节点函数值的微小变化（观测误差存在）将导致整个区间上插值函数都发生改变，因而稳定性差等几个方面。因此拉格朗日插值法多用于理论分析，在采用拉格朗日插值方法进行插值计算时通常选取n < 7。分段线性插值函数（仅连续）与三次样条插值函数（二阶导数连续）虽然光滑性差，但他们都克服了拉格朗日插值函数的缺点，不仅收敛性、稳定性强，而且方法简单实用，计算量小。因而应用十分广泛。

数据拟合 在科学计算中经常要建立实验数据的数学模型。给定函数的实验数据，需要用比较简单和合适的函数来逼近（或拟合）实验数据。这种逼近的特点是： (a) 适度的精度是需要的； (b) 实验数据有小的误差； (c) 对于某些问题，可能有某些特殊的信息能够用来选择实验数据的数学模型。 逼近离散数据的基本方法就是曲线拟合，常采用最小二乘拟合 曲线拟合问题的数学描述是，已知一组（二维）数据(xi,yi ) ，i = 1,2,。。。,n（即平面上的n个点(xi, yi ) ，i = 1,2,。。,n）， x i 互不相同。寻求一个函数（曲线） y = f (x)，使f (x)在某种准则下与所有数据点最为接近，即曲线拟合得最好。 最小二乘拟合分为线性最小二乘拟合和非线性最小二乘拟合。

一 、线性最小二乘拟合 线性最小二乘法是解决曲线拟合问题最常用的方法，基本思路是，令 f(x)=a1*r1(x)+a2*r2(x)+…….+am*rm(x) 其中rk(x)是一组事先选定的线性无关的函数，(ak)是一组待定系数。寻求系数(ak)使得yi与f(xi)的距离 d i （i = 1,2,,n）的平方和最小。这种准则称为最小二乘准则，其求系数的方法称为线性最小二乘拟合方法。 1、 系数(ak)的求法 若记 则J 为a1‥am的二次函数。由高等数学的极值理论，J 达到最小的充分必要条件是a1‥am 满足dJ/dak =0(k = 1,。。。,m)。于是得到求使J 达到最小的a1‥am的方法是求解线性方程组（称为法方程组）

即求解线性方程组 若记 则以上方程组可表示为 RT RA = RTY 。 由于当{ r1(x ), ‥ , rm( x)} 线性无关时，R列满秩，RT R可逆，所以方程组（10.7）有唯一解 A = (RT R)-1RTY 。

用以上方法作线性最小二乘拟合的误差通常考虑以下两种形式： 最小平方误差： 最大偏差： 2、 函数组{ r1(x ), ‥ , rm( x)}的选取 面对一组数据 (x i , y i ), i = 1,2,。。。,n ，用线性最小二乘法作曲线拟合时，首要的、也是关键的一步是恰当地选取{ r1(x ), ‥ , rm( x)} 。

如果通过机理分析，能够知道y 与x之间应该有什么样的函数关系，则{ r1(x ), ‥ , rm( x)}容易确定。若无法知道y 与x之间的关系，通常可以将数据(x i , y i ), i = 1,2,。。,n,作图，直观地判断应该用什么样的曲线去作拟合。人们常用的曲线有： （i）直线 y = a1 x + a2 ； （ii）多项式 y = a1 x m + ‥ +am x+am+1, （一般m = 2,3，不宜太高） （iii）双曲线（一支）y=a1/x+a2拟合前作变量替换t = 1/x求解 a 1 ,a2 较简单 （iv）指数曲线 y =a1ea2x 拟合前作变量代换z = ln y，t = 1/x，则指数曲线y =a1ea2x 转化为关于lna1 ,a2 的线性函数z = ln a1+ a2t，这样做求 解 a1 ,a2 较简单。

f (x) = aebx ，f (x) = aeb/x ， f (x) = axb ， f (x) = a + b / x 在实际计算过程中，面对一组已知数据，到底用什么样的曲线拟合最好，可以在直观判断的基础上，选几种曲线分别拟合，然后比较，看哪条曲线的最小二乘指标J 最小。 二 、非线性最小二乘拟合 非线性最小二乘法是假设f (x)是待定系数{ak}的任意非线性函数，在最小二乘准则下求其系数{ak} 。 例如上述人们常用的双曲线和指数曲线就是非线性最小二乘拟合中最常用的非线性函数，只不过在上面使用中将它们转变成线性最小二乘拟合方法。 对于给定的实验数据，通常应根据实验数据的走向、趋势选择合适的数学模型，即拟合函数。例如当实验数据具有单调性和凸性时，可选择下述适当的数学模型y = f (x)来拟合实验数据。 f (x) = aebx ，f (x) = aeb/x ， f (x) = axb ， f (x) = a + b / x

在有可能的情况下，一般将非线性拟合函数转化为线性拟合函数求解，这一方面是如此求解简单，另一方面也是因为一般情况下求解法方程组dJ/dak=0得到的(a1, ‥, am) 通常仅是J 的驻点，不一定是极值点。也可以直接解J 极小化问题。 三、 最小二乘拟合法的Matlab实现 命令为 A = R \Y 例3 用最小二乘法求一个形如y = a1 + bx2的经验公式，使它与下表所示的数据拟合. x 19 25 31 38 44 y 19.0 32.3 49.0 73.3 97.8 解 编写程序如下 x=[19 25 31 38 44]'; y=[19.0 32.3 49.0 73.3 97.8]'; r=[ones(5,1),x.^2]; ab=inv(r’*r)*r’*y x0=19:0.1:44; y0=ab(1)+ab(2)*x0.^2; plot(x,y,'o',x0,y0,'r')

y=amxm+…+a1x+a0系数a=[ am, …, a1, a0]。 2 线性最小二乘拟合（多项式拟合）方法 在线性最小二乘拟合中，用的较多的是多项式拟合。如果取 { r1( x), ‥, rm+1( x)} ={1, ‥ ,xm } ，即用m 次多项式拟合给定数据，则Matlab中有现成的函数 a=polyfit(x0,y0,m)， 其中输入参数x0,y0为要拟合的数据，m为拟合多项式的次数，输出参数a为拟合多项式 y=amxm+…+a1x+a0系数a=[ am, …, a1, a0]。 多项式在x处的值y可用下面的函数计算 y=polyval(a,x)。 例4 某乡镇企业1990-1996年的生产利润如下表： 年份 1990 1991 1992 1993 1994 1995 1996 利润（万元） 70 122 144 152 174 196 202 试预测1997年和1998年的利润。

解 作已知数据的的散点图， x0=[1990 1991 1992 1993 1994 1995 1996]; y0=[70 122 144 152 174 196 202]; plot(x0,y0,'*') 发现该乡镇企业的年生产利润几乎直线上升。因此，我们可以用y = a1 x + a0 作为拟合函 数来预测该乡镇企业未来的年利润。编写程序如下： a=polyfit(x0,y0,1) y97=polyval(a,1997) y98=polyval(a,1998) 求得a1 = 20 ，a0 = -4.0705×104 ，1997年的生产利润y97=233.4286，1998年的生产利润y98=253.9286。

3 非线性最小二乘拟合 Matlab的优化工具箱中提供了两个求非线性最小二乘拟合的函数：curvefit和leastsq。使用这两个命令时，都要先建立M文件fun.m，但它们定义f (x)的方式是不同的。 1 curvefit 设已知xdata=(xdata1,xdata2,…,xdatan ),ydata=(ydata1,ydata2,…,ydatan ), curvefit用以求含参量x（向量）的向量值函数F(x,xdata)=(f(x,data1), …,f(x,xdata n )) T中的参变量x（向量），使得 Sum(F(x,xdatai)-ydatai)2最小 输入格式为： （1）x=curvefit('fun',x0,xdata,ydata)； （2）x=curvefit('fun',x0,xdata,ydata,options)； （3）x=curvefit('fun',x0,xdata,ydata,options, 'grad')； （4）[x，options]=curvefit('fun',x0,xdata,ydata,…)； （5）[x，options,funval]=curvefit('fun',x0,xdata,ydata,…)； （6）[x，options,funval,Jacob]=curvefit('fun',x0,xdata,ydata,…). 输出目标函数值格式：f=fun(x,xdata). 其中x0为迭代初值，options为控制参数。

c(t) = a + be-0.02kt 中的参数a,b, k 。 2 leastsq 设已xdata=(xdata1,xdata2,…,xdatan ),ydata=(ydata1,ydata2,…,ydatan ), leastsq 用以求含参量x（向量）的向量值函数 输入格式为： （1）x= leastsq ('fun',x0)； （2）x= leastsq ('fun',x0,options)； （3）x= leastsq ('fun',x0,options, 'grad')； （4）[x，options]= leastsq ('fun',x0,…)； （5）[x，options,funval]= leastsq ('fun',x0,…)； 例5 用下面一组数据拟合函数 c(t) = a + be-0.02kt 中的参数a,b, k 。 t 100 200 300 400 500 600 700 800 900 1000 cj×103 4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59

1 用命令curvefit。此时 F(x，tdata)=(a+b e-0.02kt1，…，a+be-0.02kt10)T，x=(a，b，k) （1） 编写M文件curvefun1.m function f=curvefun1(x,tdata) f=x(1)+x(2)*exp(-0.02*x(3)*tdata) %其中x(1)=a；x(2)=b；x(3)=k； （2） 输入命令 tdata=100:100:1000 cdata=1e03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59]; x0=[0.2,0.05,0.005]; x=curvefit(‘curvefun1’,x0,tdata,cdata) f=curvefun1(x,tdata) （3）运算结果为： x=0.0070 -0.0030 0.1012 f= Columns 1 through 7 0.0045 0.0050 0.0054 0.0057 0.0059 0.0061 0.0063 Columns 8 through 10 0.0064 0.0065 0.0066 即拟合得a=0.0070，b=-0.0030，k=0.0066

2 用命令leastsq。此时 f(x)=F(x,tdata,cdata)=(a+be-0.02kt1-c1，…，a+be-0.02kt10-c10)T，x=(a，b，k) （1） 编写M文件curvefun2.m function f=curvefun2(x) tdata=100:100:1000; cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59]; f=cdata-x(1)-x(2)*exp(-0.02*x(3)*tdata) %其中x(1)=a；x(2)=b；x(3)=k； （2） 输入命令 x0=[0.2,0.05,0.005]; x=leastsq(‘curvefun2’,x0) f=curvefun2(x) （3）运算结果为： x=0.0070 -0.0030 0.1012 f=1.0e-005* Columns 1 through 7 0.0221 0.2081 -0.3933 -0.2872 0.2973 0.3561 0.0693 Columns 8 through 10 -0.2327 -0.0970 0.0296 可以看出，两个命令的计算结果是相同的

回归分析 回归分析是处理很难用一种精确方法表示出来的变量之间关系的一种数学方法，它是最常用的数理统计方法，能解决预测、控制、生产工艺优化等问题。它在工农业生产和科学研究各个领域中均有广泛的应用。回归分析一般分为线性回归分析和非线性回归分析。本节着重介绍线性回归分析的基本结论及其在Matlab中的相应命令。线性回归分析是两类回归分析中较简单的一类，也是应用较多的一类。 一 、一元线性回归分析 针对一组（二维）数据( xi,yi )，i = 1,2,。。,n（其中xi 互不相同），其最简单的数据拟合形式为寻求直线y = b1 + b2*x ，使b1 + b2*x在最小二乘准则下与所有数据点最为接近。但由于随机观测误差的存在，满足上述数据点的直线应该是 y = b1 + b2* x +e， （1.1） 其中x, y是准确的，b1 ,b2 是两个未知参数，e 是均值为零的随机观测误差，具有不可观测性，可以合理地假设这种观测误差服从正态分布。于是我们得到一元线性回归模型为：

其中s 未知，固定的未知参数 b1 、 b2 称为回归系数，自变量x称为回归变量。 y = b1 + b2* x +e 正态分布 E(e) =0,D(e)=s2 其中s 未知，固定的未知参数 b1 、 b2 称为回归系数，自变量x称为回归变量。 式两边同时取期望得： Y = b1 + b2*x ，称为y 对x的回归直线方程。 在该模型下，第i个观测值可以看作样本Yi = b1 + b2*xi +ei （这些样本相互独立但不同分布， i = 1,2,,n）的实际抽样值，即样本值。 一元线性回归分析的主要任务是：用实验值（样本值）对b1 、b2和s 作点估计；对回归系数b1 、b2 作假设检验；在x = x0 处对y 作预测，并对y 作区间估计。

2 模型的假设、预测、控制 1回归方程的显著性检验 在实际问题中，因变量y 与自变量x之间是否有线性关系(1.1)只是一种假设，在求出回归方程之后，还必须对这种回归方程同实际观测数据拟合的效果进行检验。由(10.10)可知，| b2 |越大， y 随x变化的趋势就越明显；反之，| b2 |越小， y 随x变化的趋势就越不明显。特别当b2 =0时，则认为y 与x之间不存在线性关系，当 b2 ≠0时，则认为y与x之间有线性关(1.1)。因此，问题归结为对假设 H0:b2=0;H1:b2 ≠0 进行检验。假设: H0:b2=0被拒绝，则回归显著，认为y 与x之间存在线性关系，所求的线性回归方程有意义；否则回归不显著， y 与x的关系不能用一元线性回归模型来描述，所得的回归方程也无意义。此时，可能有如下几种情况： （i）x对y 没有显著影响，此时应丢掉变量x； （ii）x对y 有显著影响，但这种影响不能用线性关系来表示，应该用非线性 回归； （iii）除x之外，还有其他不可忽略的变量对y 有显著影响，从而削弱了x对 y 的影响。 此时应用多元线性回归模型。因此，在接受H 0 的同时，需要进一步查明原因以便分别处理。

比较F与Fa大小，来判断x，y是否存在线性关系。 如果F>Fa,则两者有显著的线性关系。反之没有。 下面介绍两种检验方法，分别是 Yi为根据回归公式计算的到值。 (a)F检验法 对样本方差 进行分解，有 上式中的 是由实际观测值没有落在回归直线上引起的（否则为零），U 是由回归直线引起的。因此，U 越大， 就越小，表示y 与x的线性关系就越显著；否则，U 越小， 就越大，表示y 与x的线性关系就越不显著。这样我们就找到了一种判别回归直线拟合程度好坏的方法： 如果 较大时，则对拟合效果感到满意。 Matlab计算公式 x=finv(1-α,n1,n2) 比较F与Fa大小，来判断x，y是否存在线性关系。 如果F>Fa,则两者有显著的线性关系。反之没有。

xi，yi 为实际数据 ，Yi为根据回归公式计算的到值。

（b）t检验法 判别指标 x=-tinv(α/2,n1 ) 当|t|>ta (n-2)时，x，y存在明显的线性关系，当|t|<ta (n-2)时, x，y不存在明显的线性关系。

|R|越大x，y线性关系越强，反之线性关系越弱。 检验指标 0≤|R| ≤1 |R|>Ra |R|越大x，y线性关系越强，反之线性关系越弱。

1．考察温度x对产量y的影响，测得下列10组数据： 求y关于x的线性回归方程，检验回归效果是否显著，并预测x=42℃时产量的估值及预测区间（置信度95%）. 解： clear clc x=[ 20 25 30 35 40 45 50 55 60 65 ] y=[ 13.2 15.1 16.4 17.1 17.9 18.7 19.6 21.2 22.5 24.3] plot(x,y,'r*') %y=a*x+b y1=[ones(10,1),x'] A=inv(y1'*y1)*y1'*y'%求的系数a，b y0=A(2).*x+A(1); %假设检验

假设H0:a=0，H1:a≠0 我们分别采用t检验和F检验来考察x，y的关系是否正确。

clear clc x=[ 20 25 30 35 40 45 50 55 60 65 ]; y=[ 13.2 15.1 16.4 17.1 17.9 18.7 19.6 21.2 22.5 24.3]; plot(x,y,'r*') y1=[ones(10,1),x']; A=inv(y1'*y1)*y1'*y'%求的系数a，b y0=A(2).*x+A(1); Lyy=sum((y-mean(y)).^2); Lxx=sum((x-mean(x)).^2); Lxy=sum((x-mean(x)).*(y-mean(y))); U=sum((y0-mean(y)).^2); Q=Lyy-U; %F检验 F=U*(10-2)/Q Fa=finv(1-0.05,1,8) %t检验 t=sqrt((10-2)*Lxx)*A(2)/sqrt(Q) ta=-tinv(0.025,8)

回归分析的Matlab实现 Matlab统计工具箱中提供了一些回归分析的命令，现介绍如下。 1 多元线性回归 多元线性回归的命令是regress，此命令也可用于一元线性回归。其格式为： （1）确定回归系数的点估计，用命令：b=regress(Y，X)。 （2）求回归系数的点估计和区间估计，并检验回归模型，用命令： [b，bint，r，rint，stats]=regress(Y，X，alpha)。 （3）画出残差及其置信区间，用命令： rcoplot（r，rint）。

在上述命令中，各符号的含义为： （i） b为回归方程的系数，Y，X的定义同本部分前面所述。对一元线性 回归，Y，X中取k=1即可； （ii）alpha为显著性水平（缺省时为0.05）； （iii）bint为回归系数的区间估计； （iv）r与rint分别为残差及其置信区间； （v）stats是用于检验回归模型的统计量，有三个数值，第一个是 R2，第二个是F值，第三个是与F对应的概率P。其中R2与F定义同前，值越大，说明回归方程越显著，P <a 时拒绝H0 ，回归模型成立。

clc,clear x1=[ 20 25 30 35 40 45 50 55 60 65 ]'; y=[ 13.2 15.1 16.4 17.1 17.9 18.7 19.6 21.2 22.5 24.3]'; x=[ones(10,1),x1]; [b,bint,r,rint,stats]=regress(y,x); b,bint,stats, rcoplot(r,rint)

例 合金的强度y 与其中的碳含量x有比较密切的关系，今从生产中收集了一批数据如下表。试先拟合一个函数y(x)，再用回归分析对它进行检验。 解 先画出散点图： x=0.10:0.01:0.18; y=[42.0,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0]; plot(x,y, ' +' ) 可知y与x大致为线性关系。 设回归模型为y= b1 + b2*x ，用regress和rcoplot编程如下： clc,clear x1=[0.10:0.01:0.18] ' ; y=[42.0,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0] ' ; x=[ones(9,1),x1]; [b,bint,r,rint,stats]=regress(y,x); b,bint,stats,rcoplot(r,rint)

得到 b=27.4722 137.5000 bint=18.6851 36.2594 75.7755 199.2245 stats=0.7985 27.7469 0.0012 即b1=27.4722， b2=137.5000， b1的置信区间是[18.6851，36.2594]， b2的置信区间是[75.7755，199.2245]； R2 =0.7985，F =27.7469， p = 0.0012。

可知所设回归模型成立。 观察命令rcoplot(r,rint)所画的残差分布，除第8个数据外其余残差的置信区间均包含零点， 第8个点应视为异常点，将其剔除后重新计算，可得 b=30.7280 109.3985 bint=26.2805 35.2834 76.9014 141.8955 stats=0.9188 67.8534 0.0002 应该用修改后的这个结果。

练习1. 某地区某商品的价格x(单位：元/公斤)与到货量y (单位：吨)的样本数据如下 xi 7 12 6 9 10 8 12 6 11 9 12 10 yi 57 72 51 57 60 55 70 55 70 53 76 56 求关于的线性回归方程.

3） 预测与控制 当检验结果拒绝了H0:b2=0，接下来的问题是如何利用回归方程Y = b1 + b2*x 进行预测和控制。预测就是对固定的x值预测相应的y 值，控制就是通过控制x的值，以便把y 的值控制在制定的范围内。 (a)预测 设y 与x满足模型(10.11)，并且通过了假设检验，即检验结果拒绝了: H0:b2=0 。令x0表示x的某个固定值，且y0 = b1 + b2 x0 +e ， e 服从正态分布。假设y0 , y1 ,…. , yn 相互独立， 则y0 的预测值和预测区间如下： y0=b1+b2x0+e

注意：上述区间适用于样本容量比较大的情况下，另外，在进行预测时，x0的取值最好能落在所给数据x的范围之内，而对超过范围的x0，预测范围可能会变得很大，置信度也会大大降低，预测也会失去意义。在工作中经常取a=0.05 (b)控制 控制问题是预测问题的反问题它是要以1-a置信度来确定，将x控制在什么范围时，才能使y落入预先给定的区间(y1,y2),也就是说，对于给定的置信度1-a，求相应的新x1，x2，使得当x属于(x1,x2)时，可使x所对应的y属于(y1,y2).

x的控制范围为： 注意：在实现控制时，必须限制：

二 可线性化的一元非线性回归（曲线回归） 在工程技术中，自变量x与因变量y 之间有时呈现出非线（或曲线）关系，这是通常出现两种情况：一种是呈现多项式的关系，这种情况通过变量替换可化为多元线性回归问题给予解决；另一种是呈现出其它非线性关系，通过变量替换可化为一元线性回归问题给予解决。 若匹配曲线（经验公式）为含参量a,b的非线性曲线，采用的办法是通过变量替换把非线性回归化为线性回归。 通常匹配的含参量a,b的非线性曲线有以下六类，具体的替换方法如下：

1 双曲线 1/y=a+b*1/x 对于这类曲线，只要令y’=1/y,x’=1/x,则有 y’=a+bx’ 2 幂函数曲线 y = axb ，其中x > 0,a > 0。 两边求对数则有 lny=lna+blnx，那么我们令 y’=lny,x’=lnx,则会有 y’=a+bx’

3 指数函数曲线 y = aebx，其中a > 0。 两边求对数则有 lny=lna+bx，那么我们令 y’=lny,a’=lna,则会有 y’=a‘+bx 4 倒指数函数曲线y = aeb/x，其中a > 0。 两边求对数则有 lny=lna+b/x，那么我们令 y’=lny,a’=lna,x’=1/x则会有 y’=a‘+bx’

5 对数函数曲线y = a + bln x ，其中x > 0。 6 S型函数曲线 y=1/(a+be-x) 那么我们令y’=1/y,x’=e-x则会有 y’=a‘+bx’ 注：对于非线性回归问题的Matlab实现问题，一种方法是化为相应的线性模型实现，另一种方法是直接应用Matlab中相应的命令，其结果是一致的。详见本节第五部分。

y= b0 + b1*x1 + ..+ bk *xk+e ，将观测数据带入的 y= b0 + b1*x1j + ..+ bk *xkj+ej 三 多元线性回归分析 一般地，在实际问题中影响应变量y 的自变量往往不止一个，不妨设有k 个为 x1 ， x2 ，…，xk 。通过观测得到一组（k +1维）相互独立的试验观测数据(x1j ,x2j , …,xkj ,yj ) y，i = 1,2,,n，其中n > k +1。假设变量y 与变量x1 ， x2 ，…， xk 之间有线性关系： y= b0 + b1*x1 + ..+ bk *xk+e ，将观测数据带入的 y= b0 + b1*x1j + ..+ bk *xkj+ej 注：b代表 e代表 Y=Xb+e

对线性模型y= b0 + b1*x1 + ..+ bk *xk+e 所要考虑的主要问题是： （i）用实验观测数据对未知参数b0 ，…， bk 和 做点估计和假设检验，从而建立因变量y 和自变量 x1 ，…， xk 之间的线性关系； （ii）在x1 = x10 ，…， xk = x1k 处对y 的值作预测和控制，并对y 作区间估计。本部分总是假设n > k +1。 1 未知参数bi估计 XTXb = XTY b=(XTX)-1 XTY

2 多元线性回归中的假设检验 在实际问题中，往往事先不知道或不能确定随机变量y 与自变量x1 ，…， xk 之间确有线性关系。因而(10.13)往往是一种假设，因此在求出线性回归方程之后，还必须对求出的线性回归方程同实际观测数据拟合效果进行检验。类似于一元线性回归，可提出以下原假设 H0:b1=b2=‥=bk =0。 当拒绝H0时表示线性关系成立，否则不成立。

检验指标 知F检验法的检验规则为： 如果F>Fa ，则拒绝H0 ，认为因变量y 与自变量1 x ，…， k x 之间的线 性关系显著；否则，认为y 与x1 ，… xk 之间的线性关系不显著。 需要注意的是， y 与x1 ，…， xk 之间的线性关系不显著，可能出现几种情况：如y 于其中某些自变量无关系，可以去掉这些自变量； y 与1 x ，…， k x 之间的存在非线性关系；还有其它变量与y 有关系等。当然还有其它检验方法。

例8 某厂生产的一种电器的销售量y 与竞争对手的价格x1和本厂的价格x2有关。下表是该厂商品在10个城市的销售记录。试根据这些数据建立y与x1和x2的关系式，对得到的模型和系数进行检验。若某市本厂产品售价160（元），竞争对手售价170（元），预测商品在该市的销售量。 x1元 120 140 190 130 155 175 125 145 180 150 x2元 100 110 90 150 210 150 250 270 300 250 Y个 102 100 120 77 46 93 26 69 65 85 解 分别画出y关于x1和y关于x2的散点图，可以看出y与x2有较明显的线性关系，而y与x1之间的关系则难以确定，我们将作几种尝试，用统计分析决定优劣。

设回归模型为 y = b0 + b1 x1 + b2 x2 。 编写如下程序： x1=[120 140 190 130 155 175 125 145 180 150] ' ; x2=[100 110 90 150 210 150 250 270 300 250] ' ; y=[102 100 120 77 46 93 26 69 65 85] ' ; x=[ones(10,1),x1,x2]; [b,bint,r,rint,stats]=regress(y,x); b,bint,stats 得到 b=66.5176 0.4139 -0.2698 bint=-32.5060 165.5411 -0.2018 1.0296 -0.4611 -0.0785 stats=0.6527 6.5786 0.0247

可以看出结果不是太好，p=0.0247，取a =0.05时所设回归模型可用，但取a =0.01时所设回归模型不能用； R2 =0.6527较小； b0,b1的置信区间包含了零点。后面将试图用x1，x2的二次函数改进它。

说明所有自变量与Y间的线性相关程度, 即观察值Y与估计值 之间的相关程度。 软件有关结果 反映了回归方程的精度，其值越小说明回归效果越好 Root MSE （残差标准差 反映了回归方程的精度，其值越小说明回归效果越好 R-Square （决定系数） 说明所有自变量能解释Y变化的百分比。取值（0，1），越接近1模型拟合越好 Adj R-Sq （校正决定系数） 说明所有自变量与Y间的线性相关程度, 即观察值Y与估计值 之间的相关程度。

根据对未来GNP及PI的估计，预测未来投资额 投资额与国民生产总值和物价指数 建立投资额模型，研究某地区实际投资额与国民生产总值 ( GNP ) 及物价指数 ( PI ) 的关系 问题 根据对未来GNP及PI的估计，预测未来投资额 该地区连续20年的统计数据 年份序号 投资额 国民生产总值 物价 指数 年份 序号 投资额 国民生产总值 物价 指数 1 90.9 596.7 0.7167 11 229.8 1326.4 1.0575 2 97.4 637.7 0.7277 12 228.7 1434.2 1.1508 3 113.5 691.1 0.7436 13 206.1 1549.2 1.2579 4 125.7 756.0 0.7676 14 257.9 1718.0 1.3234 5 122.8 799.0 0.7906 15 324.1 1918.3 1.4005 6 133.3 873.4 0.8254 16 386.6 2163.9 1.5042 7 149.3 944.0 0.8679 17 423.0 2417.8 1.6342 8 144.2 992.7 0.9145 18 401.9 2631.7 1.7842 9 166.4 1077.6 0.9601 19 474.9 2954.7 1.9514 10 195.0 1185.9 1.0000 20 424.5 3073.0 2.0688

四 逐步线性回归分析 从多元线性回归分析中我们知道，采用的自变量越多，则回归平方和越大，残差平方和越小。然而，采用较多的变量来拟合回归方程，得到的方程稳定性差，每个自变量的区间误差的积累将影响总体误差，用这样建立起来的回归方程作预测的可靠性差、精度低。另一方面，如果采用了对因变量影响小的自变量而遗漏了重要变量，可导致估计量产生偏倚和不一致性。因而希望得到最优的回归方程。逐步线性回归分析方法就是一种自动从大量可供选择的变量中选择那些对建立回归方程比较重要的变量的方法，它是在多元线性回归基础上派生的一种算法技巧，详可参阅相应的文献。其基本思路为：从一个自变量开始，视自变量对y 作用的显著程度，从大到小依次逐个引入

回归方程。当引入的自变量由于后面自变量的引入而变得不显著时，要将其剔除掉。引入一个自变量或从回归方程中剔除一个自变量，为逐步回归的一步。对于每一步，都要进行y 值检验，以确保每次引入新的显著性变量前回归方程中只包含对y 作用显著的变量。这个过程反复进行，直至即无不显著的变量从回归方程中剔除，又无显著变量可引入回归方程止。 使用逐步线性回归时要注意：要适当选择引入变量的显著性水平和剔除变量的显著性水平； 应尽量选择那些相互独立性强的变量。

4 逐步回归 逐步回归的命令是stepwise，它提供了一个交互式画面，通过此工具可以自由地选择变量，进行统计分析。通常用法是： stepwise(x，y，inmodel，alpha)， 其中x是自变量数据， y 是因变量数据，分别为n ´m和n ´1矩阵，inmodel是矩阵的列数指标，给出初始模型中包括的子集（缺省时设定为全部自变量），alpha为显著水平（缺省时为0.05）。 运行stepwise命令时产生三个图形窗口：Stepwise Plot，Stepwise Table Stepwise History。 所有这些图形界面都由热区，即当鼠标移到图形的某个区域时，鼠标的指针会变成一个小圆，点击后会产生交互作用。

在Stepwise Plot窗口，显示出各项的回归系数及其置信区间。其中：点表示回归系数的值，点两边的水平（实或虚）直线段表示其置信区间（虚线表示该变量的拟合与0无显著差异，实线表示有显著差异）；绿色的线表示当前在模型中的项，红色的线表示当前不在模型中的项。 点击一条线会改变其状态，即在模型中的项（绿线）会被移去（变为红线），不在模型中的项（红线）会被加入（变为绿线）。次窗口中的Export下拉式菜单可以向Matlab工作区传送各种数据。次窗口中的Scale Inputs可对输入数据的每列进行正态化处理，使其标准差为1。 在Stepwise Table窗口中列出了一个统计表，包括回归系数及其置信区间，以及模型的统计量剩余标准差(RMSE)、相关系数(R-square)、F值、与F值对应的概率P。

例11 水泥凝固时放出的热量y与水泥中4种化学成分x1,x2,x3,x4有关，今测得一组数据如下， 试用逐步回归来确定一个线性模型。 序号 x1 x2 x3 x4 y 1 7 26 6 60 78.5 2 1 29 15 52 74.3 3 11 56 8 20 104.3 4 11 31 8 47 87.6 5 7 52 6 33 95.9 6 11 55 9 22 109.2 7 3 71 17 6 102.7 8 1 31 22 44 72.5 9 2 54 18 22 93.1 10 21 47 4 26 115.9 11 1 40 23 34 83.8 12 11 66 9 12 113.3 13 10 68 8 12 109.4

解 编程如下： clc,clear x0=[1 7 26 6 60 78.5 2 1 29 15 52 74.3 3 11 56 8 20 104.3 4 11 31 8 47 87.6 5 7 52 6 33 95.9 6 11 55 9 22 109.2 7 3 71 17 6 102.7 8 1 31 22 44 72.5 9 2 54 18 22 93.1 10 21 47 4 26 115.9 11 1 40 23 34 83.8 12 11 66 9 12 113.3 13 10 68 8 12 109.4]; x0=(:,2:5); y=x0(:,6); stepwise(x,y)

得到图Stepwise Plot和Stepwise Table表。图Stepwise Plot中四条直线都是虚线，说明模型的显著性不好，从Stepwise Table表中可以看出变量x3和x4显著最差。移去这两个变量，图Stepwise Plot中点击直线3和直线4，这两条直线变为红色，同时直线1和直线2变为实线，说明移去变量x3和x4后模型具有显著性。从新的统计结果可以看出，虽然剩余标准差（RMSE）没有太大的变化，但是统计量F的值明显增大，因此新的回归模型更好一些。 对变量y和x1、x2作线性回归： stepwise(x,y,[1,2]) 得到结果： y =52.5773+1.4683 x1 +0.6623 x2 。

五 回归分析的Matlab实现 Matlab统计工具箱中提供了一些回归分析的命令，现介绍如下。 1 多元线性回归 多元线性回归的命令是regress，此命令也可用于一元线性回归。其格式为： （1）确定回归系数的点估计，用命令：b=regress(Y，X)。 （2）求回归系数的点估计和区间估计，并检验回归模型，用命令： [b，bint，r，rint，stats]=regress(Y，X，alpha)。 （3）画出残差及其置信区间，用命令： rcoplot（r，rint）。

在上述命令中，各符号的含义为： （i） b为回归方程的系数，Y，X的定义同本部分前面所述。对一元线性 回归，Y，X中取k=1即可； （ii）alpha为显著性水平（缺省时为0.05）； （iii）bint为回归系数的区间估计； （iv）r与rint分别为残差及其置信区间； （v）stats是用于检验回归模型的统计量，有三个数值，第一个是 R2，第二个是F值，第三个是与F对应的概率P。其中R2与F定义同前，值越大，说明回归方程越显著，P <a 时拒绝H0 ，回归模型成立。

2 多项式回归 （1）一元多项式的回归和预测 一元多项式的回归和预测可用命令polyfit或polytool和polyval或polyconf来实现。其命令格式如下： 令 x =( x1, x2....xn) y =( y1, y2... yn) ， P=(a1,a2.... am+1) 是多项式 y=a1*xm +a2*xm-1 - +amx+am+1 .的系数，S 是一个矩阵，用来估计预测误差。 回归可以用命令[P,S]=polyfit(x,y,m)或polytool(x,y,m)实现，其[P,S]=polyfit(x,y,m)是确定多项式系数的命令；命令polytool(x,y,m)命令产生一个交互式的画面，在画面中绿色曲线为拟合曲线，它两侧的红线是y 的置信区间。可以用鼠标移动图中的十字线来改变图下方的x值，

也可以在窗口内输入，左边就给出y的预测值与置信区间。通过左下方的Export下拉式菜单，可以输出回归系数等。 预测和预测误差估计的命令为polyval或者polyconf，其中Y=polyval(p,x)求polyfit所得的回归多项式在x处的预测值Y；[Y,DELTA]=polyconf(p,x,S,alpha)求polyfit所得的回归多项式在x处的预测值Y及预测值的显著性为1-alpha的置信区间Y± DELTA；alpha缺省时为0.05。 一元多项式回归也可化为多元线性回归来解。 一元多项式回归也可化为多元线性回归来解。

例9 观测物体降落的距离y 与时间x的关系，得到数据如下表，求y 关于x的回归方程 y = a + bx + cx2。 x 1/30 2/30 3/30 4/30 5/30 6/30 7/30 y 11.86 15.67 20.60 26.67 33.71 41.93 51.13 x 8/30 9/30 10/30 11/30 12/30 13/30 14/30 y 61.49 72.90 85.44 99.08 113.77 129.54 146.48 解 方法（一） 用一元多项式回归，编写程序如下： x=1/30:1/30:14/30; y=[11.86 15.67 20.60 26.67 33.71 41.93 51.1361.49 72.90 85.44 99.08 113.77 129.54 146.48]; [p,S]=polyfit(x,y,2); 得到p=489.2946 65.8896 9.1329 即a=9.1329 b=65.8896 c=489.2846

方法（二） 化为多元线性回归，其程序为： x=1/30:1/30:14/30; y=[11.86 15.67 20.60 26.67 33.71 41.93 51.1361.49 72.90 85.44 99.08 113.77 129.54 146.48]; T=[ones(14,1),t,t.^2] [b,bint,r,rint,stats]=regress(y,T); b,stats 得到结果：b=9.1329 65.8896 489.2946 stats=1.0e+007* 0.0000 1.0378 0 可以看出，两种方法的出的结果是一样的

（2）多元二项式回归 多元二项式回归可用命令：rstool(x,y,model,alpha)。其中，输入数据x、y分别为n×m矩阵和n维列向量；alpha为显著性水平（缺省时为0.05）；model由下列4个模型中选择1个（用字符串输入，缺省时为线性模型）： linear（线性）： y = b0 + b1x1 +....+ bmxm ； purequadratic（纯二次）： y = b0 + b1x1 +....+ bmxm + interaction（交叉）： y = b0 + b1x1 +....+ bmxm +

命令rstool产生一个交互式画面，画面中有m个图形，这m个图形分别给出了一个独立变量xi（另m-1个变量取固定值）与y的拟合曲线，以及y的置信区间。可以通过键入不同的xi的值来获得相应的y值。 图的左下方有两个下拉式菜单，一个菜单Export用以向Matlab工作区传送数据，包括beta（回归系数）、rmse（剩余标准差）、residuals（残差）。另一个菜单model用以在上述4个模型中选择。可以分别选4个模型，并比较它们的剩余标准差，其中最接近于0的模型是最好的。 我们再作一遍例8商品销售量与价格问题，选择纯二次模型，即 y = b0 + b1* x1 + b2*x2 + b3*x12 + b4* x22 。

编程如下： x1=[120 140 190 130 155 175 125 145 180 150] ' ; x2=[100 110 90 150 210 150 250 270 300 250] ' ; y=[102 100 120 77 46 93 26 69 65 85] ' ; x=[x1 x2]; rstool(x,y, ' purequadratic ' ) 得到一个交互式画面，给出两幅图形。左边图形是x1固定时的曲线y(x1)及其置信区间，右边图形是x2固定时的曲线y(x2)及其置信区间。用鼠标移动图中的十字线，或在图下方窗口内输入，可改变x1、x2。画面左边给出y的预测值即其置信区间，用这种画面可以回答例8提出的“若谋市本厂产品售价160（元），竞争对手售价170（元），预测商品在该市的销售量”问题。

在画面左下方的下拉式菜单Export中选择“all”，则beta、rmse和residuals都传送到Matlab工作区中。在Matlab工作区中输入命令： -0.0228 0.0037 rmse=16.6436 如果在另一菜单model选择其它多元二项式模型，比较它们的剩余标准差就会发现，本例的所选模型的 rmse=16.6436最小。 注：本例中的模型亦可化为多元线性回归来做。请读者自己编程并比较结果。

3 非线性回归 非线性回归可用命令nlinfit,nlintool,nlparci,nlpredci来实现。命令格式如下： 回归：回归可用命令[beta,r,J]=nlinfit(x,y,model,beta0)或者nlintool(x,y,model,beta0,alpha)来实现。其中命令[beta,r,J]=nlinfit(x,y,model,beta0)的作用为确定回归系数；而命令nlintool(x,y,model,beta0,alpha)产生一个交互式的画面，画面中有拟合曲线和y的置信区间。通过 左下方的Export下拉式菜单，可以输出回归系数等。

某些非线性回归也可化为多元线性回归来解。 例10 在研究化学动力学反应过程中，建立了一个反应速度和反应物含量的数学模型，形式为 其中 b1 ，…，b5 式未知系数， x1 , x2 , x3 是三种反应物（氢，n戊烷，异构戊烷）的含量， y 是反应速度。今测的一组数据如下表，试由此确定参数 b1 ，…， b5 ，并给出置信区间。 b1 ，…， b5 的参考值为（0.1,0.05,0.02,1,2）。

序号 反应速度y 氢x1 n戊烷x2 异构戊烷x3 1 8.55 470 300 10 2 3.79 285 80 10 3 4.82 470 300 120 4 0.02 470 80 120 5 2.75 470 80 10 6 14.39 100 190 10 7 2.54 100 80 65 8 4.35 470 190 65 9 13.00 100 300 54 10 8.50 100 300 120 11 0.05 100 80 120 12 11.32 285 300 10 13 3.12 285 190 120

解 首先，以回归系数和自变量为输入变量，并将要拟合的模型写成函数文件 huaxue.m: function yhat=huaxue(beta,x); yhat=(beta(4)*x(2)-x(3)/beta(5))./(1+beta(1)*x(1)+beta(2)*x(2)+beta(3)*x(3)); 然后，用nlinfit计算回归系数，用nlparci计算回归系数的置信区间，用nlpredci计算预测值 即其置信区间，编程如下： clc,clear x0=[1 8.55 470 300 10 2 3.79 285 80 10 3 4.82 470 300 120 4 0.02 470 80 120

5 2.75 470 80 10 6 14.39 100 190 10 7 2.54 100 80 65 8 4.35 470 190 65 9 13.00 100 300 54 10 8.50 100 300 120 11 0.05 100 80 120 12 11.32 285 300 10 13 3.12 285 190 120]; x=x0(:,3:5); y=x0(:,2); beta=[0.1,0.05,0.02,1,2]; %回归系数的初值 [betahat,f,j]=nlinfit(x,y, ' huaxue ' ,beta); % f,j是下面命令用的信息 betaci=nlparci(betahat,f,j); betaa=[betaa,betaci] %回归系数及其置信区间 [yhat,delta]=nlpredci( ' huaxue ' ,x,betahat,f,j) %y的预测值及其置信区间半径，置信区间为yhat ±delta。 用命令nlinfit(x,y, ' huaxue ' ,beta)可看到画面，并传出剩余标准差rmse=0.1933

建模案例：水塔水流量的估计 某居民区供水机构有一供居民用水的圆柱形水塔，由于该机构没有测量流入或流出该水塔的水量的装置，水塔的管理者只能通过隔一段时间测量一次水塔的水位来估计水的流量。但面临的困难是，当水塔水位下降到设定的最低水位时，水泵自动启动向水塔供水，到设定的最高水位时停止供水，这段时间无法测量水塔的水位和水泵的供水量。通常每天水泵供水一到二次，每次约2个小时。水塔是一个高为12.2米，直径为17.4米的正圆柱。按照设计，水塔水位降至约8.2米时，水泵自动启动，水位升到约10.8米是水泵停止工作。表1是某一天管理人员测量该水塔的水位测量记录（符号“//”表示水泵启动）。请您帮助该管理人员估计一天中任何时刻（包括水泵正在供水时刻）从水塔中流出的水流量及一天的居民用水总量。

表1 水位测量记录 时刻(h) 0 0.921 1.843 2.949 3.871 4.978 5.900 7.006 7.928 8.967 水位(m) 9.677 9.479 9.308 9.125 8.982 8.814 8.686 8.525 8.388 8.220 时刻(h) 9.9811 10.925 10.954 12.032 12.954 13.875 14.982 15.903 16.826 17.931 水位(m) // // 10.820 10.500 10.210 9.936 9.653 9.409 9.180 8.921 时刻(h) 19.037 19.959 20.839 22.015 22.958 23.880 24.986 25.908 水位(m) 8.662 8.433 8.220 // 10.820 10.591 10.354 10.

一、问题提出 某居民区供水机构有一供居民用水的圆柱形水塔，由于该机构没有测量流入或流出该水塔的水量的装置，水塔的管理者只能通过隔一段时间测量一次水塔的水位来估计水的流量。但面临的困难是，当水塔水位下降到设定的最低水位时，水泵自动启动向水塔供水，到设定的最高水位时停止供水，这段时间无法测量水塔的水位和水泵的供水量。通常每天水泵供水一到二次，每次约2个小时。水塔是一个高为12.2米，直径为17.4米的正圆柱。按照设计，水塔水位降至约8.2米时，水泵自动启动，水位升到约10.8米是水泵停止工作。表1是某一天管理人员测量该水塔的水位测量记录（符号“//”表示水泵启动）。请您帮助该管理人员估计一天中任何时刻（包括水泵正在供水时刻）从水塔中流出的水流量及一天的居民用水总量。

二 问题的分析 流量是单位时间流出的水的体积。由于水塔是正圆柱形，截面面积是常数，在水泵停止供水的时段，流量很容易从水位对时间的变换率算出，问题是如何估计水泵供水时段的水流量。水泵供水时段的流量与供水时段前后的流量有关。如果能够通过测量数据，产生若干个时刻的水流量，则计算水流量问题就转化为插值问题。这样水泵供水时段的流量就可以用供水时段前后的流量插值得到。作为用于函数插值的原始数据我们希望水泵不工作时段的流量越准确越好。 有了任何时刻的流量，就不难计算一天的总用水量。水泵不工作时段的用水量可以由测量记录直接得到，如由表1可知从t=0到t=8.967(h)水位下降了9.677-8.220=1.457(m)，乘以水塔的 截面积就是这一时段的用水量。这个数值可以用来检验插值的结果。 注：此处我们只考虑水流量的插值问题，当然亦可以考虑水位的插值问题；或者亦可以将此问题按拟合问题进行处理。

三 模型假设 1 流量只取决于水位差，与水位本身无关。按照Torricelli定律从小孔流出的流体的流速（流量）正比于水面高度的平方根，问题给出水塔的最低和最高水位分别约是8.22米和10.82米（设出水口的水位为零），因为10.82 / 8.22 » 1，所以可以忽略水位对流速（流量）的影响。 2 水塔中水流量是时间的连续光滑函数，与水泵是否工作无关。 3 水泵是否工作完全取决于水塔内的水位的高度，且每次加水的工作时间约为2小时，从表1中可知两次供水时间段分别为[8.967，10.954]和[20.839，22.958]。 4 水泵工作时单位时间的供水量大致是常数，此常数大于单位时间的平均流量。 5 一天之中可以从任何时刻开始，即从t=0(h)到t=24(h)结束与从t=a (h)到t=24+a (h)结束差别不大。每天同一时刻水流量保持相同，并且一天中开始时刻的不同对一天中的总用水量影响很小。 6 水塔内部底面直径为17.4米。 四 数据处理

表2 水塔中水的体积（单位：时刻（小时），体积（立方米）） 1 水体积计算 水塔是一个正圆柱，塔内水的体积为V =pi/4*D2*h，其中D为水塔直径，h为水面高度。近似地取pi =3.141592654，得到不同时刻水塔中水的体积如表2。 表2 水塔中水的体积（单位：时刻（小时），体积（立方米）） 时刻 0 0.921 1.843 2.949 3.871 4.978 5.900 7.006 7.928 8.967 水位 2294 2247 2206 2163 2129 2089 2059 2020 1988 1948 时刻9.9811 10.925 10.954 12.032 12.954 13.875 14.982 15.903 16.826 17.931 水位 // // 2564 2489 2420 2355 2288 2230 2176 2114 时刻 19.037 19.959 20.839 22.015 22.958 23.880 24.986 25.908 水位 2053 1999 1948 // 2564 2510 2454 2413

2 水流速(流量)的近似 水流速度（流量）是水塔中水的体积对时间的导数。由于没有水的体积关于时间的函数表达式，而只是一个离散的函数值表2，因此考虑常用的差商代替导数处理方法。为提高精度，采用二阶差商公式。具体地，表2的数据点被水泵两次工作分割成三组数据，对每组数据中间数据采用二阶中心差商，两边数据采用二阶向前或向后差商得到数据表3。其公式分别为： 中心差商公式 向前差商公式 向后差商公式

表3 水塔中水的流速（单位：时刻（小时），流速（立方米/小时）） 0 0.921 1.843 2.949 3.871 4.978 5.900 7.006 7.928 8.967 54.516 42.320 38.085 41.679 33.297 37.817 30.748 38.455 32.122 41.718 9.9811 10.925 10.954 12.032 12.954 13.875 14.982 15.903 16.826 17.931 // // 73.686 76.434 71.686 60.190 68.333 59.217 52.011 56.626 19.037 19.959 20.839 22.015 22.958 23.880 24.986 25.908 63.023 54.859 55.439 // 57.602 57.776 51.891 36.464

五 模型建立与求解 1 做出水流速的散点图 Matlab程序如下： r=[54.516 42.320 38.085 41.679 33.297 37.817 30.748 38.455 32.122 41.718 73.686 76.434 71.68660.190 68.333 59.217 52.011 56.626 63.023 54.859 55.439 57.602 57.776 51.891 36.464]; plot(t,r,' +' ); title('流速散点图'); xlabel('时间（小时）'); ylabel('流速（立方米/小时）') 做出的散点图略。 2 模型及计算结果 现在问题已转变为根据流速f (t)的一个函数值表3，产生函数f (t)在整个区间（二十四小时）上的函数或函数值，插值和拟合是最常用的两种方法。如果建立拟合模型，需要根据散点图的趋势，选择适当的拟合函数形式。若采用插值模型，可以考虑分段线性插值（光滑性较差）、分段三次多项式插值、三次样条插值等插值方法。我们采用插值方法。通过对不同插值方法的比较，结合假设2，即流速是时间的连续光滑函数，下面采用三次样条插值模型。

1）水流量及一天用水总量 首先由三次样条插值计算得到水流量函数f (t)。其Matlab程序为： r=[54.516 42.320 38.085 41.679 33.297 37.817 30.748 38.455 32.122 41.718 73.686 76.434 71.68660.190 68.333 59.217 52.011 56.626 63.023 54.859 55.439 57.602 57.776 51.891 36.464]; [1,n]=size(t); d1=t(n)-t(1); x=t(0):1/3600:t(n); %被插值点 ys=interp1(t,r,x,'spline'); %样条插值输出 plot(x,ys); title('样条插值下的水流速图'); xlabel('时间（小时）'); ylabel('流速（立方米/小时）') 做出的水流速曲线图略。

2）一天（24小时）的用水总量 用三次样条插值模型得到的函数f (t)在区间[0,24]上积分得到的结果为1257.3立方米。 3）误差分析 直接由表3的水流速数据，用梯形公式进行数值积分得到的结果为1250.3立方米，与上面得到的结果比较，则有：三次样条插值的绝对误差为7个立方米，相对误差为0.56%（不到百分之一）。另外，测量数据的误差可控制在0.5%；数值微分、数值积分的误差亦可以控制在一定的范围内。 六 模型的稳定性分析与检验 1 稳定性分析 用不同时刻作为起始点，使用三次样条插值函数得到的水流量函数f (t)，在长度为24小时的时间区间上积分，所得结果相差无几，详见表4。这说明我们所建立的数值微分模型和三次样条插值模型是稳定的。

表4 不同起点计算出的24小时用水量（单位：立方米） 起始点 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 用水量 1257.3 1258.5 1260.3 1262.6 1265.1 1267.6 1269.6 1270.7 2 检验 1） 分三段（水泵未工作）的实际用水量与由模型推算的用水量的差异很小，见表5，最大的不超过4%，三段总计不超过0.5%。其中由模型推算的用水量为水流量函数f (t)的积分；实际用水量为水塔内水的体积之差。 表5 分三段的实际用水量与模型计算出的24小时用水量比较（单位：立方米） 实际用水量 模型用水量 绝对误差 相对误差 第一段[0，8.967] 345.3682 338.0511 7.3171 2.1% 第二段[10.954，20.839] 616.3161 620.8341 4.5180 0.7% 第三段[22.958，25.908] 151.7308 156.6000 4.8692 3.2% 三段总计 1113.4151 1115.4852 2.0701 0.2%

2）两次充水期间，水泵注入量的差异大约为3个立方米，不到0.5%。 水泵充水量=充水后的水量+充水期间的流出量—充水前的水量。 第一次水泵充水量=2564+117.4463-1948=733.4463立方米； 第二次水泵充水量=2564+120.7052-1948=736.7052立方米。 由于水泵功率=充水量/充水时间，由上述计算知，水泵的功率大约为367立方米/小时，而且两次冲水期间计算出的功率大致相等，这与实际问题是一致的。 3）用水高峰的比较 实际与模型之间相差无几。实际用水高峰可近似地用差商最大值点表示为t=11，即上午11点钟左右；模型得到的用水高峰可由模型得到的水流量函数依次在单位区间上积分或者找出水流量函数的最大值点得出，它们都在11点左右。 七 模型的优缺点及推广 1 优点 （1）模型灵活性好、稳定性强，可用于那些拥有地方性的竖直圆柱型水塔的小城镇和乡镇。 模型中的输入数据可以是任何近似均匀的时间间隔时的水位，时间间隔大约2小时。 （2）模型中的数学概念简单，并且容易理解。只用到数值计算知识。 （3）模型容易实现，且给出了一天里水流速度和用水量的精确估计。

2 缺点 （1）本模型受水塔的几何形状限制。 （2）光滑曲线的逼近方法不能模拟真实流动中流速的微小变化，实际流动中流速可能会有一种高程度的“噪音”，即激波出现。 3 改进与推广 （1）我们可以在模型中用一个参数来限定不同几何形状的水塔。 （2）可以通过对流速数据进行回归分析等一系列处理，以便得到一些随机变化的特征。