1、普通高等院校计算机课程规划教材,MATLAB数据分析方法,李柏年 吴礼斌 主编 张孔生 丁 华 参编,回归分析是最常用的数据分析方法之一。它是根据已得的试验结果以及以往的经验来建立统计模型，并研究变量间的相关关系，建立起变量之间关系的近似表达式即经验公式，并由此对相应的变量进行预测和控制等.,3.1一元回归模型,3.1.1一元线性回归模型,1.一元线性回归的基本概念,通常，我们对总体(x,Y)进行n次的独立观测，获得n组数据（称为样本观测值）,(x1,y1),(x2,y2),(xn,yn),利用最小二乘法可以得到回归模型参数0，1的最小二乘估计,设Y是一个可观测的随机变量，它受到一个非随机变量
2、因素x和随机误差的影响。若Y与x有如下线性关系：,(3.1.1),且E=0,D=2,则称(3.1.1)为一元线性回归模型. 其中0，1为回归系数，x为自变量，Y为因变量.,(3.1.2),其中,于是建立经验公式模型：,(3.1.3),一元线性回归分析的主要任务：一是利用样本观测值对回归系数0，1和作点估计；二是对方程的线性关系即1作显著性检验；三是在x=x0处对Y作预测等. 以下举例说明建立经验公式（3.1.3）的方法。,例3.1.1 近10年来，某市社会商品零售总额与职工工资总额（单位：亿元）数据如下表3.1。 表3.1 商品零售总额与职工工资表 （单位：亿元）,建立社会商品零售总额与职工工
3、资总额数据的回归模型,解：% 首先输入数据x=23.80,27.60,31.60,32.40,33.70,34.90,43.20,52.80,63.80,73.40; y=41.4,51.8,61.70,67.90,68.70,77.50,95.90,137.40,155.0,175.0;,% 然后作散点图 plot(x,y,*) %作散点图 xlabel(x(职工工资总额) %横坐标名 ylabel(y(商品零售总额) %纵坐标名,图3.1商品零售总额与职工工资总额数据散点图,% 计算最佳参数 Lxx=sum(x-mean(x).2); Lxy=sum(x-mean(x).* (y-mean
4、(y); b1=Lxy/Lxx; b0=mean(y)-b1*mean(x);,运行后得到： b1 = 2.7991，b0 = -23.5493 所以，回归模型为,问题1：当x=0，得到y=-23.5493亿元如何理解？,问题2：如何检验E=0? D=2?,2. 一元多项式回归模型,在一元回归模型中，如果变量y与x的关系是n次多项式，即,其中是随机误差，服从正态分布N(0,2),a0,a1,an为回归系数，则称(3.1.4)为多项式回归模型.,(3.1.4),(1)多项式曲线拟合,在MATLAB7的统计工具箱中，有多项式曲线拟合的命令polyfit，其调用格式有以下三种：,p=polyfit(
5、x,y,n) p,S=polyfit(x,y,n) p,S,mu=polyfit(x,y,n),其中，输入x,y分别为自变量与因变量的样本观测数据向量；n是多项式的阶数，对于一元线性回归则取n=1;输出p是按照降幂排列的多项式的系数向量，S是一个矩阵，用于估计预测误差或供MATLAB的其它函数的调用 。,例3.1.2 某种合金中的主要成分为A,B两种金属，经过试验发现：这两种金属成分之和x与合金的膨胀系数y有如下关系，建立描述这种关系的数学表达式.,表3.2 合金的膨胀系数表,解：%首先输入数据 x=37:0.5:43; y=3.4,3,3,2.27,2.1,1.83,1.53,1.7,1.8
6、,1.9,2.35,2.54,2.9;%其次做散点图 plot(x,y,*) xlabel(x(两种合金之和) %横坐标名 ylabel(y(合金膨胀系数) %纵坐标名%然后根据散点图猜测曲线类别,（2.1.7）,由于散点图呈抛物线，故选择二次函数曲线进行拟合.,p = polyfit(x,y,2) %注意取n=2 运行得到回归系数： p=0.1660 -13.3866 271.6231 即二次回归模型为：,多项式曲线拟合预测的命令polyval，其调用格式有以下两种：,Y=polyval(p,x0) Y,Delta=polyconf(p,x0,S,alpha) 其中，输入p,S是由多项式拟合
7、命p,S=polyfit(x,y,n) 的输出，x0是要预测的自变量的值.输出Y是polyfit所得的回归多项式在x处的预测值。,(2) 多项式回归的预测与置信区间,如果输入数据的误差相互独立，且方差为常数，则YDelta至少包含95%的预测值；alpha缺省时为0.05。 (Y-Delta, Y+Delta)即95%的置信区间,(3) 多项式回归的GUI界面命令,多项式回归的GUI界面命令polytool，其典型调用格式polytool(x,y,n,alpha) 其中，输入x,y分别为自变量与因变量的样本观测数据向量；n是多项式的阶数；置信度为(1-alpha)%，alpha缺省时为0.05
8、。,该命令可以绘出总体拟合图形以及(1-alpha)上、下置信区间的直线（屏幕上显示为红色）.此外，用鼠标拖动图中纵向虚线，就可以显示出对于不同的自变量数值所对应的预测状况，与此同时图形左端数值框中会随着自变量的变化而得到的预报数值以及(1-alpha) 置信区间长度一半的数值。,例3.1.3为了分析X射线的杀菌作用，用200千伏的X射线来照射细菌，每次照射6分钟用平板计数法估计尚存活的细菌数，照射次数记为t，照射后的细菌数y如表3.3所示。,表3.3 X射线照射次数与残留细菌数,试求： 给出y与t的二次函数回归模型； 在同一坐标系内做出原始数据与拟合结果的散点图 预测t=16时残留的细菌数；
9、 根据问题实际意义选择多项式函数是否合适？,数据来源：http/www.ilr.cornell.edu/hadi/RABE,解：% 输入原始数据 t=1:15; y=352,211,197,160,142,106,104,60,56,38,36,32,21,19,15; p=polyfit(t,y,2); % 作二次多项式回归 y1= polyval(p,t); % 模型估计与作图 plot(t,y,-*,t,y1,-o); legend(原始数据,二次函数) xlabel(t(照射次数) ylabel(y(残留细菌数) t0=16; yc1= polyconf(p,t0) % 预测t0=16
10、时残留的细菌数,运行结果为 p =1.9897 -51.1394 347.8967，yc1 =39.0396,即二次回归模型为,yc1 =39.0396,表明照射16次后，用二次函数计算出细菌残留数为39.0396，显然与实际不相符合。 调用多项式回归的GUI界面命令polytool,如图3.4,原始数据与拟合结果的散点图如图3.3所示，从图形可知拟合效果较好.,图 3.3 原始数据与拟合结果的散点图,根据实际问题的意义可知：尽管二次多项式拟合效果较好，但是用于预测并不理想。因此如何根据原始数据散点图的规律，选择适当的回归曲线是非常重要的，因此有必要研究非线性回归分析.,图 3.4 二次函数预
11、测交互图,3.1.2一元非线性回归模型,为了便于正确地选择合适的函数进行回归分析建模，我们给出通常选择的六类曲线如下所示：,1. 非线性曲线选择,（1）双曲线1/y=a+b/x(见图3.5)。,图3.5双曲线,图3.5双曲线,(2) 幂函数曲线y=axb, 其中x0,a0(图3.6)。,图3.6 幂函数曲线,（3）指数曲线y=aebx，其中参数a0(见图3.7)。,图3.7 指数曲线,（4）倒指数曲线 ,其中a0(图3.8)。,图3.8 倒指数曲线,（5）y=a+blnx (见图3.9)。,图3.9 对数曲线,（6）S型曲线 (见图3.10)。,图3.10 S型曲线,对于非线性回归建模通常有两
12、种方法：一是通过适当的变换转化为线性回归模型，例如双曲线模型(图3.5)。如果无法实现线性化，可以利用最小二乘法直接建立非线性回归模型，求解最佳参数。,2.非线性回归的MATLAB命令,MATLAB统计工具箱中实现非线性回归的命令有nlinfit、nlparci、lpredci和nlintool。下面逐一介绍调用格式。,非线性拟合命令nlinfit，调用格式： beta,r,J = nlinfit(x,y,model,beta0) 其中，输人数据x，y分别为nm矩阵和n维列向量，对一元非线性回归，x为n维列向量，model是事先用M文件定义的非线性函数，beta0是回归系数的初值(需要通过解方
13、程组得到)，beta是估计出的最佳回归系数，r是残差，J是Jacobian矩阵，它们是估计预测误差需要的数据。,非线性回归预测命令nlpredci，调用格式：ypred = nlpredci(FUN,inputs,beta,r,J) 其中，输入参数beta,r,J是非线性回归命令nlinfit的输出结果, FUN 是拟合函数，inputs是需要预测的自变量；输出量ypred是inputs的预测值。,非线性回归置信区间命令nlparci，调用格式：ci = nlparci(beta,r,J,alpha) 输入参数beta,r,J就是非线性回归命令nlinfit输出的结果，输出ci是一个矩阵，每一
14、行分别为每个参数的(1-alpha)% 的置信区间，alpha缺省时默认为0.05.,非线性回归的GUI界面命令nlintool，典型调用格式nlintool(x,y,fun,beta0) 其中参数x,y,fun,beta0与命令nlinfit中的参数含义相同.,例3.1.4. 在M文件中建立函数y=a(1-be-cx)，其中a,b,c为待定的参数。,解：fun=inline(b(1)*(1-b(2)*exp(-b(3)*x),b,x); 此处，将b看成参变量，b(1),b(2),b(3)为其分量.,例3.1.5 炼钢厂出钢时所用盛钢水的钢包，由于钢水对耐火材料的侵蚀，容积不断增大，我们希望找
15、出使用次数与增大容积之间的函数关系.实验数据如表3.4。,表3.4 钢包使用次数与增大容积,（1）建立非线性回归模型1/y=a+b/x； （2）预测钢包使用x0=17次后增大的容积y0； （3）计算回归模型参数的95%的置信区间。,MATLAB脚本程序如下: x=2:16; y=6.42,8.2,9.58,9.5,9.7,10,9.93,9.99,10.49,10.59,10.6,10.8,10.6,10.9,10.76;,%建立非线性双曲线回归模型 b0=0.084,0.1436; % 初始参数值 fun=inline(x./(b(1)*x+b(2),b,x); beta,r,J=nlinf
16、it(x,y,fun,b0); beta % 输出最佳参数 y1=x./(0.0845*x+0.1152); % 拟合曲线 plot(x,y,*,x,y1,-or) legend(原始数据,拟合曲线),注意：初始值要先计算后，才能得到上面程序中的b0，由于确定两个参数值，因此我们选择已知数据中的两点（2,6.42）和（16,10.76）代入设定方程，得到方程组,上述方程组有两种解法：手工方法与Matlab方法。下面用Matlab方法解方程组： a,b=solve(6.42*(2*a+b)=2,10.76*(16*a+b)=16) 输出为 a = .83961597702347450462657
17、355615004e-1 b = .14360328434608391527406223581049,图3.11钢包使用次数与增大容积的非线性拟合图,在例3.1.5中，预测钢包使用17次后增大的容积，可在执行上面的程序中，继续输入命令 ypred=nlpredci(fun,17,beta,r,J) 得到： ypred =10.9599 即钢包使用17次后增大的容积10.9599。,求回归模型参数的95%的置信区间，只要继续添加程序 ci = nlparci(beta,r,J) 运行后得到 ci =0.0814 0.08760.0934 0.1370,即回归模型 中参数a,b的95%的置信区间分
18、别为 (0.0814 ,0.0876) 与（0.0934，0.1370）. 我们求出的最佳参数分别为 a=0.0845,b=0.1152均属于上述置信区间。,图3.12给出钢包使用次数与增大容积的非线性拟合的交互图形，图中的的圆圈是实验的原始数据点，两条虚线为95%上、下置信区间的曲线（屏幕上显示为红色），中间的实线（屏幕上显示为绿色）是回归模型曲线，纵向的蓝色虚线显示了在自变量为8.9502，横向的虚线给出了对应的预测值为10.2734.,图3.12 钢包使用次数与增大容积的非线性拟合交互图,例3.1.6 对例题3.1.3进行非线性回归，并预测照射16次后细菌残留数目，给出模型参数的95%的
19、置信区间，绘出模型交互图形.,解：我们选取函数y=aebt进行非线性回归，该方程的两个参数具有简单的物理解释，a表示实验开始时的细菌数目，b表示细菌死亡（或衰变）的速率。 MATLAB脚本程序如下：,t=1:15; y=352 211 197 160 142 106 104 60 56 38 36 32 21 19 15; fun=inline(b(1)*exp(b(2)*t),b,t) % 非线性函数 beta0=148,-0.2; % 参数初始值 beta,r,J=nlinfit(t,y,fun,beta0); % 非线性拟合 beta % 输出最佳参数 y1=nlpredci(fun,t
20、,beta,r,J); % 模型数值计算,plot(t,y,*,t,y1,-or), legend(原始数据,非线性回归) xlabel(t(照射次数) ylabel(y(残留细菌数) ypred = nlpredci(fun,16,beta,r,J) % 预测残留细菌数 ci = nlparci(beta,r,J) % 参数95%区间估计 nlintool(t,y,fun,beta0) % 作出交互图形,运行后结果如下： beta =400.0904 -0.2240 即，最佳参数为：a=400.0904，b=-0.2240 故非线性回归模型为,预测为：ypred =11.1014 即，照射1
21、6次后细菌残留数目为11.1014，该预测符合实际，显然比例3.1.3中多项式回归的结果合理。 ci =355.2481 444.9326-0.2561 -0.1919 即参数a置信度为95%的置信区间 (ci的第一行)为：355.2481 , 444.9326 参数b的置信度为95%的置信区间 (ci的第二行)为-0.2561 -0.1919 显然，最佳参数a=400.0904，b=-0.2240，均属于各自置信度为95%的置信区间。,图3.13原始数据与非线性回归图形,图3.14 原始数据与非线性回归GUI图形,从交互图形3.14可以看出：圆圈为原始数据，两条虚线（屏幕上显示红色）是置信区
22、间曲线；两条虚线内的实线（屏幕上显示绿色）是回归模型曲线；纵向虚线指示照射8次，此时对应的水平虚线表示模型得到的残留细菌数为：66.6451。,图3.14 原始数据与非线性回归GUI图形,3.1.3一元回归建模实例,例3.1.7在四川白鹅的生产性能研究中，得到如下一组关于雏鹅重（g）与70日龄重(g)的数据，试建立70日龄重(y)与雏鹅重(x)的直线回归方程，计算模型误差平方和以及可决系数，当雏鹅重分别为：85,95 ,115时预测其70日龄重，以及置信区间。,表3.5 四川白鹅重与70日龄重测定结果 （单位：g）,解：（1）作散点图。以雏鹅重（x）为横坐标，70日龄重（y）为纵坐标作散点图，
23、如图2-14。 在MATLAB命令窗口中输入: x=80 86 98 90 120 102 95 83 113 105 110 100; % 雏鹅重 y=2350 2400 2720 2500 3150 2680 2630 2400 3080 2920 2960 2860; %70日龄重 plot(x,y,*) %作散点图 xlabel(x(雏鹅重) %横坐标名 ylabel(y(70日龄重) %纵坐标名,图3.15 四川白鹅的雏鹅重与70日龄重散点图和回归直线图,由图形3.15可见白鹅的70日龄重与雏鹅重间存在直线关系，且70日龄重随雏鹅重的增大而增大。因此，可认为y与x符合一元线性回归模型
24、。,（2）建立直线回归方程。在MATLAB中调用命令 polyfit，从而求出参数0，1的最小二乘估计. 在MATLAB命令窗口中继续输入:,n= size(x,1) % 计算样本容量 p,s=polyfit(x,y,1); % 调用命令polyfit计算回归参数 y1=polyval(p,x); % 计算回归模型的函数值 hold on plot(x,y1) % 作回归方程的图形，结果如图3.15 p % 显示参数的最小二乘估计结果 p=582.185021.7122,即参数 的最小二乘估计为,所以70日龄重(y)与雏鹅重(x)的直线回归经验方程为,（3）误差估计与决定系数。在MATLAB命
25、令窗口中继续输入: TSS=sum(y-mean(y).2) %计算总离差平方和 RSS=sum(y1-mean(y).2) %计算回归平方和 ESS=sum(y-y1).2) %计算残差平方和 R2=RSS/TSS; %计算样本决定系数R2.,输出: TSS =8.314917e+005 RSS =7.943396e+005 ESS =3.715217e+004 R2= 0.9553 TSS =8.314917e+005 RSS =7.943396e+005 ESS =3.715217e+004 R2=0.9553,由于样本决定系数R2=0.9553接近于1，因此模型的拟合的效果较好。,（4
26、）回归方程关系显著性的F检验。在MATLAB命令窗口中继续输入: F=(n-2)*RSS/ESS %计算的F统计量 F1=finv(0.95,1,n-2) %查F统计量0.05的分位数 F2=finv(0.99,1,n-2) %查F统计量0.01的分位数,输出结果: F=2.138e+002，F1 =4.9646，F2 =10.0442 为了方便，将以上的计算结果列成表3.6。,表3.6 四川白鹅70日龄重与雏鹅重回归关系方差分析表,因为 表明四川白鹅70日龄重与雏鹅重间存在显著的线性关系。,（5）回归关系显著性的t检验。在MATLAB命令窗口中继续输入: T=p(2)/sqrt(ESS/(n
27、-2)*sqrt(sum(x-mean(x).2)%计算T统计量 T1=tinv(0.975,n-2) %t统计量0.05的分位数 T2=tinv(0.995,n-2) %t统计量0.01的分位数,输出: T =14.622, T1 =2.228, T2 =3.169,因为T=14.62t0.01(10)，否定H0，接受H1即四川白鹅70日龄重(y)与雏鹅重(x)的线性回归系数是显著的,可用所建立的回归方程进行预测和控制。,(6)预测 x1=85,95,115; % 输入自变量 yc=polyval(p,x1) % 计算预测值 Y,Delta=polyconf(p,x1,s); I1=Y-De
28、lta,Y+Delta % 置信区间,输出: yc = 2427.722644.843079.08 I1 =2279.47 2575.962503.01 2786.672927.55 3230.62,所以当雏鹅重分别为85,95,115时，白鹅70日龄重分别为 2427.72, 2644.84, 3079.08; 且95%的置信区间分别为： 2279.47 ,2575.96， 2503.01,2786.67, 2927.55,3230.62.,在程序中加入： polytool(x,y) % 交互功能 bar(x,y-y1), % 残差图 legend(残差) h=lillietest(y-y1
29、) % 残差正态性检验 输出h = 0 得到交互图形如图3.16所示，可以看出当雏鹅重为100时，模型给出70日龄鹅重为2753.4016.,图3.16 四川白鹅70日龄重与雏鹅重线性模型交互图,3.2多元线性回归模型,3.2.1多元线性回归模型及其表示,对于总体,的n组观测值,它应满足式（3.2.1），即,其中i (i=1,2,n)相互独立，且设 记,则模型(3.2.2)可用矩阵形式表示为Y=X+ (3.2.3) 其中Y称为观测向量，X称为设计矩阵，称为待估计向量，是不可观测的n维随机向量，它的分量相互独立，假定 .,2. 多元线性回归建模的基本步骤 (1)对问题进行直观分析，选择因变量与解
30、释变量，作出与因变量与各解释变量的散点图，初步设定多元线性回归模型的参数个数； (2)输入因变量与自变量的观测数据(y,X)调用命令b,bint,r,rint,s=regress(y,X,alpha)， 计算参数的估计。 (3)调用命令rcoplot(r,rint)，分析数据的异常点情况。,(4)作显著性检验，若检验通过，则用模型作预测。 (5)对模型进一步研究：如残差的正态性检验，残差的异方差检验，残差进行自相关性的检验等。,3.2.2 MATLAB的回归分析命令 在MATLAB7.0的统计工具箱中，与多元回归模型有关的命令有多个，下面逐一介绍。 1.多元回归建模命令regeress，其调用
31、格式有以下三种： （1）b = regress(y,X) （2）b,bint,r,rint,stats = regress(Y,X) （3）b,bint,r,rint,stats = regress(Y,X,alpha),三种方式的主要区别在输出项参数多少上，第3种方式可称为全参数方式。以第3种为例来说明regeress命令的输入与输出参数的含义。,输入参数：输入量Y表示模型（3.1.1）中因变量的观测向量；X是一个的矩阵，其中第一列元全部是数“1”，第j列是自变量Xj的观测向量，即,对一元线性回归，取p=1即可；alpha为显著性水平 输出参数：输出向量b为回归系数估计值，bint为回归系数
32、的(1-alpha)置信区间；输出向量r表示残差列向量,输出rint为模型的残差的 (1- )的置信区间；输出stats是用于检验回归模型的统计量，有4个分量值：第一个是R2，其中R是相关系数，第二个是F统计量值，第三个是与统计量F对应的概率P，当P时拒绝H0，即认为线性回归模型有意义，第四个是方差2的无偏估计.,例3.2.1某销售公司将库存占用资金情况、广告投入的费用、员工薪酬以及销售额等方面的数据作了汇总,该公司试图根据这些数据找到销售额与其他变量之间的关系，以便进行销售额预测并为工作决策提供参考依据。(1)建立销售额的回归模型；(2)如果未来某月库存资金额为150万元，广告投入预算为45
33、万元，员工薪酬总额为27万元，试根据建立的回归模型预测该月的销售额。,表3.7 占用资金、广告投入、员工薪酬、销售额（单位：万元）,解：为了确定销售额与库存占用资金、广告投入、员工薪酬之间的关系，分别作出y与x1,x2,x3的散点图，若散点图显示它们之间近似线性关系，则可设定y与x1,x2,x3的关系为三元线性回归模型,%输入数据并作散点图（图3.18） A=75.2 30.6 21.1 1090.4;77.6 31.3 21.4 1133 80.7 33.9 22.9 1242.1;76 29.6 21.4 1003.2 79.5 32.5 21.5 1283.2;81.8 27.9 21.
34、7 1012.2 98.3 24.8 21.5 1098.8;67.7 23.6 21 826.3 74 33.9 22.4 1003.3;151 27.7 24.7 1554.6 90.8 45.5 23.2 1199;102.3 42.6 24.3 1483.1 115.6 40 23.1 1407.1;125 45.8 29.1 1551.3 137.8 51.7 24.6 1601.2;175.6 67.2 27.5 2311.7 155.2 65 26.5 2126.7;174.3 65.4 26.8 2256.5;,m,n=size(A); subplot(3,1,1),plot(
35、A(:,1),A(:,4),+), xlabel(x1(库存资金额) ylabel(y(销售额) subplot(3,1,2),plot(A(:,2),A(:,4),*), xlabel(x2(广告投入) ylabel(y(销售额) subplot(3,1,3),plot(A(:,3),A(:,4),x), xlabel(x3(员工薪酬) ylabel(y(销售额),所得图形如图3.18所示，可见销售额y与库存资金、广告投入、员工薪酬具有线性关系，因此可以建立三元线性回归模型.,图3.18销售额与库存、广告、薪酬散点图,% 调用命令regress建立三元线性回归模型 x=ones(m,1),
36、A(:,1), A(:,2), A(:,3); y=A(:,4) b,bint,r,rint,stats=regress(y,x); b,bint,stats, % 输出结果,程序运行结果 b =162.06327.273913.9575-4.3996,bint =-580.3603 904.48674.3734 10.17437.1649 20.7501-46.7796 37.9805,stats =0.9574804050 105.0866520891 0.0000000008 10077.9867891125,输出结果说明：,b就是模型中的参数0 ,1 ,2 ，因此回归模型为,b就是模型
37、中的参数0 ,1 ,2 ，因此回归模型为,bint的各行分别为参数0 ,1 ,2的95%的置信区间。,stats的第一列为模型可决系数，第二列为F统计量的观测值，第三列得到概率p，最后一列为模型的残差平方和,2.多元回归辅助图形命令 （1）残差图命令rcoplot，其调用格式 rcoplot(r,rint) 其中，输入参数r,rint 是多元回归建模命令regress 输出的结果，运行该命令后展示了残差与置信区间的图形。该命令有助于对建立的模型进行分析，如果图形中出现红色的点，则可以认作异常点，此时可删除异常点，重新建模，最终得到改进的回归模型。,在上面的程序中加入rcoplot(r,rint
38、) 得到如下图形,图3.19 残差与置信区间图,从图形中可以看到第五个点为异常点，实际上从表3.7可以发现第5个月库存占用资金、广告投入、员工薪酬均比3月份少，为何销售额反而增加？这就可以促使该公司的经理找出原因，寻找对策。下面的例题介绍如何删除异常点，对模型进行改进的方法。,例3.2.2 葛洲坝机组发电耗水率的主要影响因素为库水位,出库流量。数据如表3.8所示，利用多元线性回归分析方法建立耗水率与出库流量、库水位的模型。,表3.8 某天耗水率与出库流量、库水位的数据,解：% 输入原始数据 A=65.08 15607 60.46 65.10 15565 60.28 65.12 15540 60
39、.10 65.17 15507 59.78 65.21 15432 59.44 65.37 15619 59.25 65.38 15536 58.91 65.39 15514 58.76 65.40 15519 58.73 65.43 15510 58.63 65.47 15489 58.48 65.53 15437 58.31 65.62 16355 57.96 65.58 14708 57.06 65.70 14393 56.43 65.84 14296 55.83;,% 做散点图 subplot(1,2,1),plot(A(:,1),A(:,3),+) xlabel(x1(库水位) yl
40、abel(y(耗水率) subplot(1,2,2),plot(A(:,2),A(:,3),o) xlabel(x2(出库流量) ylabel(y(耗水率),运行后得到的图形如图3.20所示，从图中可以看到无论是库水位还是出库流量都与机组发电耗水率具有线性关系，因此，可以建立机组发电耗水率与库水位和出库流量的二元线性回归模型。,图3.20 库水位、出库流量与耗水率的散点图,% 建立模型 m,n=size(A); y=A(:,3); x=A(:,1:2); b,bint,r,rint,stats=regress(y,ones(m,1),x); b,bint,stats 输出回归模型的系数、系数置
41、信区间与统计量如表3.9所示,表3.9回归模型的系数、系数置信区间与统计量,由此可得模型为：,% 模型改进rcoplot(r,rint); 得到图形如图3.21所示，发现有一个异常点，下面给出删除异常点后，重新建模的程序。,由此可得模型为：,图3.21 残差示意图,% 删除异常点程序并建模 b1,bint1,r1,rint1,stats1=regress(y(1:12);y(14:m),ones(m-1,1),x(1:12,:);x(14:m,:) rcoplot(r1,rint1);,删除异常点后，残差示意图如图2-21所示，此时没有异常点，改进回归模型的系数、系数置信区间与统计量参见表3.
42、10,表3.10改进回归模型的系数、系数置信区间与统计量,我们将表3.9与表3.10加以比较，可以发现：可决系数从0.9863提高到0.9931，F统计量从468.4118提高到858.5846，由此可知改进后的模型显著性提高。,图 3-22删除异常点后残差示意图,图3.21 残差示意图,3.2.3 多元线性回归实例,例3.2.3现代服务业是社会分工不断深化的产物，随着经济的发展，科学技术的进步，现代服务业的发展受到多种因素和条件的影响。不仅受到经济总体发展水平的影响，还受到第二产业、就业、投入等因素的影响，从这几个主要方面出发，利用江苏省统计年鉴的有关数据，通过建立多元线性回归模型对1990
43、-2008年各种因素对现代服务业的影响进行回归分析。假如构建如下江苏省服务业增长模型：,Y代表江苏省服务业的增加值(单位：亿元)，反映了江苏省服务业发展的总体水平。x1x4表示影响江苏省服务业发展的四种主要因素和影响，其中 x1代表江苏省人均GDP(单位：元)，说明江苏省总体经济发展水平对服务业的影响； x2代表江苏省第二产业的增加值(单位：亿元)，说明了工业发展对服务业的影响，体现了生产性服务业的需求规模； x3表示江苏省服务业的就业人数(单位：万人)； x4表示江苏省服务业资本形成总额(单位：亿元)，主要体现服务业投资的经济效应。,表3.11 江苏省关于服务业发展及各影响因素相关数据,解：
44、%输入各影响因素的数据 x0=2038 70.24 589.74 252.01 2109 35.53 623.19 275.82 33928 2055.56 1713.33 5287.91; y=37.76,28.13,93.58,160.62,286.58,277.12,387.11,367.16,291.77,280.01,227.61,329.16,385.44,437.02,601.39,704.72,1291.11,1360.09,1769.28; %Y服务业增加值列向量 n,p=size(x0); % 矩阵的行数即样本容量 x=ones(n,1),x0; % 构造设计矩阵 db,d
45、bint,dr,drint,dstats=regress(y,x); % 调用多元回归分析命令,（1）回归参数的估计 输出：db =345.24930.16720.1962-0.7012-0.6537,所以，服务业增加值Y对4个自变量的线性回归方程为,所以，服务业增加值Y对4个自变量的线性回归方程为,输出： dstats =1.0e+003 * 0.00010 0.1727 0.0000 5.7926 其中dstats的第4项是残差的方差估计值。所以，残差方差2的无偏估计值为,下面对例3.2.3的回归模型进行显著性检验。,（1）F检验 接上面的程序，在MATLAB命令窗口中继续输入：,TSS=
46、y*(eye(n)-1/n*ones(n,n)*y; % 计算TSS H=x*inv(x*x)*x; % 计算对称幂等矩阵 ESS=y*(eye(n)-H)*y; % 计算ESS RSS=y*(H-1/n*ones(n,n)*y; % 计算RSS MSR=RSS/p; % 计算MSR MSE=ESS/(n-p-1); %计算MSE,%F检验 F0=(RSS/p)/(ESS/(n-p-1); % 计算F0 Fa=finv(p,n-p-1,0.95); % F分布时的临界值,%t检验 S=MSE*inv(x*x); % 计算回归参数的协方差矩阵 T0=db./sqrt(diag(S); % 每个回归参数的T统计量 Ta=tinv(n-p-1,0.975); % t分布的分位数 pp=tpdf(T0,n-p-1); % 每个回归参数的T统计量对应的概率,
