MATLAB数据分析方法-(2)

MATLAB数据分析方法-(2)回归分析是最常用的数据分析方法之一。它是根据已得的试验结果以及以往的经验来建立统计模型，并研究变量间的相关关系，建立起变量之间关系的近似达式即经验公式，并由此对相应的变量进行预测和控制等.3.1一元回归模型3.1.1一元线性回归模型1.一元线性回归的基本概念通常，我们对总体(x,Y)进行n次的独立观测，获得n组数据（称为样本观测值）(x1,y1),(x2,y2),…,(xn,yn)利用最小二乘法可以得到回归模型参数0，1的最小二乘估计设Y是一个可观测的随机变量，它受到一个非随机变量因素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商品零售总额与职工工资表（单位：亿元）建立社会商品零售总额与职工工资总额数据的回归模型工资总额23.827.631.632.433.734.943.252.863.873.4零售总额41.451.861.767.968.777.595.9137.4155.0175.0解：%首先输入数据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(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(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,1.9,2.35,2.54,2.9];%其次做散点图plot(x,y,‘*’)xlabel('x(两种合金之和)')%横坐标名ylabel(‘y(合金膨胀系数)’)%纵坐标名%然后根据散点图猜测曲线类别（2.1.7）x3737.53838.53939.54040.54141.54242.543y3.4332.272.11.831.531.71.81.92.352.542.9由于散点图呈抛物线，故选择二次函数曲线进行拟合.p=polyfit(x,y,2)%注意取n=2运行得到回归系数：p=0.1660-13.3866271.6231 即二次回归模型为：多项式曲线拟合预测的命令polyval，其调用格式有以下两种：Y=polyval(p,x0)[Y,Delta]=polyconf(p,x0,S,alpha)其中，输入p,S是由多项式拟合命[p,S]=polyfit(x,y,n)的输出，x0是要预测的自变量的值.输出Y是polyfit所得的回归多项式在x处的预测值。(2)多项式回归的预测与置信区间如果输入数据的误差相互独立，且方差为常数，则Y±Delta至少包含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。该命令可以绘出总体拟合图形以及(1-alpha)上、下置信区间的直线（屏幕上显示为红色）.此外，用鼠标拖动图中纵向虚线，就可以显示出对于不同的自变量数值所对应的预测状况，与此同时图形左端数值框中会随着自变量的变化而得到的预报数值以及(1-alpha)置信区间长度一半的数值。例3.1.3为了分析X射线的杀菌作用，用200千伏的X射线来照射细菌，每次照射6分钟用平板计数法估计尚存活的细菌数，照射次数记为t，照射后的细菌数y如表3.3所示。t123456789101112131415y3522111971601421061046056383632211915表3.3X射线照射次数与残留细菌数试求：①给出y与t的二次函数回归模型；②在同一坐标系内做出原始数据与拟合结果的散点图③预测t=16时残留的细菌数；④根据问题实际意义选择多项式函数是否合适？数据来源：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时残留的细菌数运行结果为p=1.9897-51.1394347.8967，yc1=39.0396即二次回归模型为yc1=39.0396,表明照射16次后，用二次函数计算出细菌残留数为39.0396，显然与实际不相符合。调用多项式回归的GUI界面命令polytool,如图3.4原始数据与拟合结果的散点图如图3.3所示，从图形可知拟合效果较好.图3.3原始数据与拟合结果的散点图根据实际问题的意义可知：尽管二次多项式拟合效果较好，但是用于预测并不理想。因此如何根据原始数据散点图的规律，选择适当的回归曲线是非常重要的，因此有必要研究非线性回归分析.图3.4二次函数预测交互图3.1.2一元非线性回归模型为了便于正确地选择合适的函数进行回归分析建模，我们给出通常选择的六类曲线如下所示：1.非线性曲线选择（1）双曲线1/y=a+b/x(见图3.5)。图3.5双曲线图3.5双曲线(2)幂函数曲线y=axb,其中x>0,a>0(图3.6)。图3.6幂函数曲线（3）指数曲线y=aebx，其中参数a>0(见图3.7)。图3.7指数曲线（4）倒指数曲线,其中a>0(图3.8)。图3.8倒指数曲线（5）y=a+blnx(见图3.9)。图3.9对数曲线（6）S型曲线(见图3.10)。图3.10S型曲线对于非线性回归建模通常有两种方法：一是通过适当的变换转化为线性回归模型，例如双曲线模型(图3.5)。如果无法实现线性化，可以利用最小二乘法直接建立非线性回归模型，求解最佳参数。2.非线性回归的MATLAB命令MATLAB统计工具箱中实现非线性回归的命令有nlinfit、nlparci、lpredci和nlintool。下面逐一介绍调用格式。非线性拟合命令nlinfit，调用格式：[beta,r,J]=nlinfit(x,y,'model',beta0)其中，输人数据x，y分别为n×m矩阵和n维列向量，对一元非线性回归，x为n维列向量，model是事先用M文件定义的非线性函数，beta0是回归系数的初值(需要通过解方程组得到)，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是一个矩阵，每一行分别为每个参数的(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炼钢厂出钢时所用盛钢水的钢包，由于钢水对耐火材料的侵蚀，容积不断增大，我们希望找出使用次数与增大容积之间的函数关系.实验数据如表3.4。使用次数(x)23456789增大容积(y)6.428.29.589.59.7109.939.99使用次数(x)10111213141516增大容积(y)10.4910.5910.610.810.610.910.76表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]=nlinfit(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=.83961597702347450462657355615004e-1b=.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.08140.08760.09340.1370即回归模型中参数a,b的95%的置信区间分别为(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%的置信区间，绘出模型交互图形.解：我们选取函数y=aebt进行非线性回归，该方程的两个参数具有简单的物理解释，a表示实验开始时的细菌数目，b表示细菌死亡（或衰变）的速率。MATLAB脚本程序如下：t=1:15;y=[3522111971601421061046056383632211915];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,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即，照射16次后细菌残留数目为11.1014，该预测符合实际，显然比例3.1.3中多项式回归的结果合理。ci=355.2481444.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可以看出：圆圈为原始数据，两条虚线（屏幕上显示红色）是置信区间曲线；两条虚线内的实线（屏幕上显示绿色）是回归模型曲线；纵向虚线指示照射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）编号123456789101112雏鹅重(x)80869890120102958311310511010070日龄重(Y)235024002720250031502680263024003080292029602860解：（1）作散点图。以雏鹅重（x）为横坐标，70日龄重（y）为纵坐标作散点图，如图2-14。在MATLAB命令窗口中输入:x=[808698901201029583113105110100]';%雏鹅重y=[235024002720250031502680263024003080292029602860]';%70日龄重plot(x,y,'*')%作散点图xlabel('x(雏鹅重)')%横坐标名ylabel('y(70日龄重)')%纵坐标名图3.15四川白鹅的雏鹅重与70日龄重散点图和回归直线图由图形3.15可见白鹅的70日龄重与雏鹅重间存在直线关系，且70日龄重随雏鹅重的增大而增大。因此，可认为y与x符合一元线性回归模型。（2）建立直线回归方程。在MATLAB中调用命令polyfit，从而求出参数0，1的最小二乘估计.在MATLAB命令窗口中继续输入:n=size(x,1)%计算样本容量[p,s]=polyfit(x,y,1);%调用命令polyfit计算回归参数y1=polyval(p,x);%计算回归模型的函数值holdonplot(x,y1)%作回归方程的图形，结果如图3.15p%显示参数的最小二乘估计结果p=582.185021.7122即参数的最小二乘估计为所以70日龄重(y)与雏鹅重(x)的直线回归经验方程为（3）误差估计与决定系数。在MATLAB命令窗口中继续输入:TSS=sum((y-mean(y)).^2)%计算总离差平方和RSS=sum((y1-mean(y)).^2)%计算回归平方和ESS=sum((y-y1).^2)%计算残差平方和R2=RSS/TSS;%计算样本决定系数R2.输出:TSS=8.314917e+005RSS=7.943396e+005ESS=3.715217e+004R2=0.9553TSS=8.314917e+005RSS=7.943396e+005ESS=3.715217e+004R2=0.9553由于样本决定系数R2=0.9553接近于1，因此模型的拟合的效果较好。（4）回归方程关系显著性的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日龄重与雏鹅重回归关系方差分析表自由度（df）平方和(SS)均方和（MS）F值F0.05F0.01回归1794339.60794339.60213.81**4.9610.04残差1037152.073715.21总离差11831491.67因为表明四川白鹅70日龄重与雏鹅重间存在显著的线性关系。（5）回归关系显著性的t检验。在MATLAB命令窗口中继续输入:T=p(2)/sqrt(ESS/(n-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.62>t0.01(10)，否定H0，接受H1即四川白鹅70日龄重(y)与雏鹅重(x)的线性回归系数是显著的,可用所建立的回归方程进行预测和控制。(6)预测x1=[85,95,115]';%输入自变量yc=polyval(p,x1)%计算预测值[Y,Delta]=polyconf(p,x1,s);I1=[Y-Delta,Y+Delta]%置信区间输出:yc=2427.722644.843079.08I1=2279.472575.962503.012786.672927.553230.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)%残差正态性检验输出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)对问题进行直观分析，选择因变量与解释变量，作出与因变量与各解释变量的散点图，初步设定多元线性回归模型的参数个数；(2)输入因变量与自变量的观测数据(y,X)调用命令[b,bint,r,rint,s]=regress(y,X,alpha)，计算参数的估计。(3)调用命令rcoplot(r,rint)，分析数据的异常点情况。(4)作显著性检验，若检验通过，则用模型作预测。(5)对模型进一步研究：如残差的正态性检验，残差的异方差检验，残差进行自相关性的检验等。3.2.2MATLAB的回归分析命令在MATLAB7.0的统计工具箱中，与多元回归模型有关的命令有多个，下面逐一介绍。1.多元回归建模命令regeress，其调用格式有以下三种：（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为回归系数的(1-alpha)置信区间；输出向量r表示残差列向量输出rint为模型的残差的(1-)的置信区间；输出stats是用于检验回归模型的统计量，有4个分量值：第一个是R2，其中R是相关系数，第二个是F统计量值，第三个是与统计量F对应的概率P，当P<时拒绝H0，即认为线性回归模型有意义，第四个是方差2的无偏估计.例3.2.1某销售公司将库存占用资金情况、广告投入的费用、员工薪酬以及销售额等方面的数据作了汇总,该公司试图根据这些数据找到销售额与其他变量之间的关系，以便进行销售额预测并为工作决策提供参考依据。(1)建立销售额的回归模型；(2)如果未来某月库存资金额为150万元，广告投入预算为45万元，员工薪酬总额为27万元，试根据建立的回归模型预测该月的销售额。表3.7占用资金、广告投入、员工薪酬、销售额（单位：万元）月份库存资金额(x1)广告投入(x2)员工薪酬总额(x3)销售额(y)175.230.621.11090.4277.631.321.41133380.733.922.91242.147629.621.41003.2579.532.521.51283.2681.827.921.71012.2798.324.821.51098.8867.723.621826.397433.922.41003.31015127.724.71554.61190.845.523.2119912102.342.624.31483.113115.64023.11407.11412545.829.11551.315137.851.724.61601.216175.667.227.52311.717155.26526.52126.718174.365.426.82256.5解：为了确定销售额与库存占用资金、广告投入、员工薪酬之间的关系，分别作出y与x1,x2,x3的散点图，若散点图显示它们之间近似线性关系，则可设定y与x1,x2,x3的关系为三元线性回归模型%输入数据并作散点图（图3.18）A=[75.230.621.11090.4;77.631.321.4113380.733.922.91242.1;7629.621.41003.279.532.521.51283.2;81.827.921.71012.298.324.821.51098.8;67.723.621826.37433.922.41003.3;15127.724.71554.690.845.523.21199;102.342.624.31483.1115.64023.11407.1;12545.829.11551.3137.851.724.61601.2;175.667.227.52311.7155.26526.52126.7;174.365.426.82256.5];[m,n]=size(A);subplot(3,1,1),plot(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),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.3996bint=-580.3603904.48674.373410.17437.164920.7501-46.779637.9805stats=0.9574804050105.08665208910.000000000810077.9867891125输出结果说明：b就是模型中的参数0,1,2，因此回归模型为b就是模型中的参数0,1,2，因此回归模型为bint的各行分别为参数0,1,2的95%的置信区间。stats的第一列为模型可决系数，第二列为F统计量的观测值，第三列得到概率p，最后一列为模型的残差平方和2.多元回归辅助图形命令（1）残差图命令rcoplot，其调用格式rcoplot(r,rint)其中，输入参数r,rint是多元回归建模命令regress输出的结果，运行该命令后展示了残差与置信区间的图形。该命令有助于对建立的模型进行分析，如果图形中出现红色的点，则可以认作异常点，此时可删除异常点，重新建模，最终得到改进的回归模型。在上面的程序中加入rcoplot(r,rint)得到如下图形图3.19残差与置信区间图从图形中可以看到第五个点为异常点，实际上从表3.7可以发现第5个月库存占用资金、广告投入、员工薪酬均比3月份少，为何销售额反而增加？这就可以促使该公司的经理找出原因，寻找对策。下面的例题介绍如何删除异常点，对模型进行改进的方法。例3.2.2葛洲坝机组发电耗水率的主要影响因素为库水位,出库流量。数据如表3.8所示，利用多元线性回归分析方法建立耗水率与出库流量、库水位的模型。表3.8某天耗水率与出库流量、库水位的数据时间年-月-天-时库水位（米)出库流量（立方米）机组发电耗水率(立方米/万千瓦)2005-10-15:0065.081560760.462005-10-15:0265.101556560.282005-10-15:0465.121554060.102005-10-15:0665.171550759.782005-10-15:0865.211543259.442005-10-15:1065.371561959.252005-10-15:1265.381553658.912005-10-15:1465.391551458.762005-10-15:1665.401551958.732005-10-15:1865.431551058.632005-10-15:2065.471548958.482005-10-15:2265.531543758.312005-10-16:0065.621635557.962005-10-16:0265.581470857.062005-10-16:0465.701439356.432005-10-16:0665.841429655.83解：%输入原始数据A=[65.081560760.4665.101556560.2865.121554060.1065.171550759.7865.211543259.4465.371561959.2565.381553658.9165.391551458.7665.401551958.7365.431551058.6365.471548958.4865.531543758.3165.621635557.9665.581470857.0665.701439356.4365.841429655.83];%做散点图subplot(1,2,1),plot(A(:,1),A(:,3),'+')xlabel('x1(库水位)')ylabel('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输出回归模型的系数、系数置信区间与统计量如表3.9所示表3.9回归模型的系数、系数置信区间与统计量回归系数回归系数估计值回归系数置信区间0373.8698[340.082,407.6577]1-4.9759[-5.4642,-4.4875]20.0007[0.0004,0.0009]R2=0.9863，F=468.4118，p<0.0001，s2=0.0278由此可得模型为：%模型改进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.10表3.10改进回归模型的系数、系数置信区间与统计量回归系数回归系数估计值回归系数置信区间0328.4616[290.6145,366.3087]1-4.3594[-4.8880,-3.8308]20.0010[0.00073,0.0012]R2=0.9931，F=858.5846，p<0.0001，s2=0.0150我们将表3.9与表3.10加以比较，可以发现：可决系数从0.9863提高到0.9931，F统计量从468.4118提高到858.5846，由此可知改进后的模型显著性提高。图3-22删除异常点后残差示意图图3.21残差示意图3.2.3多元线性回归实例例3.2.3现代服务业是社会分工不断深化的产物，随着经济的发展，科学技术的进步，现代服务业的发展受到多种因素和条件的影响。不仅受到经济总体发展水平的影响，还受到第二产业、就业、投入等因素的影响，从这几个主要方面出发，利用江苏省统计年鉴的有关数据，通过建立多元线性回归模型对1990-2008年各种因素对现代服务业的影响进行回归分析。假如构建如下江苏省服务业增长模型：Y代表江苏省服务业的增加值(单位：亿元)，反映了江苏省服务业发展的总体水平。x1~x4表示影响江苏省服务业发展的四种主要因素和影响，其中x1代表江苏省人均GDP(单位：元)，说明江苏省总体经济发展水平对服务业的影响；x2代表江苏省第二产业的增加值(单位：亿元)，说明了工业发展对服务业的影响，体现了生产性服务业的需求规模；x3表示江苏省服务业的就业人数(单位：万人)；x4表示江苏省服务业资本形成总额(单位：亿元)，主要体现服务业投资的经济效应。表3.11江苏省关于服务业发展及各影响因素相关数据年份Y服务业增加值省人均GDP第二产业增加值服务业就业人数服务业资本形成总额198937.76203870.24589.74252.01199028.13210935.53623.19275.82199193.582353101.33640.95330.711992160.623106325.34706.39439.321993286.584321478.79786.37620.971994277.125801588.72855.97858.911995387.117319528.49920.451102.711996367.168471358.86975.661293.431997291.779371337.741025.221370.211998280.0110049228.241102.311624.741999227.6110695280.051151.681773.372000329.1611765515.741192.021903.372001385.4412882471.571263.772131.872002437.0214396697.031341.862189.782003601.39168301182.621407.632686.572004704.72202231650.881443.373362.1920051291.11245601917.051542.463930.5620061360.09288141895.81625.064628.5920071769.28339282055.561713.335287.91解：%输入各影响因素的数据x0=[203870.24589.74252.01210935.53623.19275.82………………………………………………………………………339282055.561713.335287.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,dbint,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.000100.17270.00005.7926其中dstats的第4项是残差的方差估计值。所以，残差方差2的无偏估计值为下面对例3.2.3的回归模型进行显著性检验。（1）F检验接上面的程序，在MATLAB命令窗口中继续输入：TSS=y'*(eye(n)-1/n*ones(n,n))*y;%计算TSSH=x*inv((x'*x))*x';%计算对称幂等矩阵ESS=y'*(eye(n)-H)*y;%计算ESSRSS=y'*(H-1/n*ones(n,n))*y;%计算RSSMSR=RSS/p;%计算MSRMSE=ESS/(n-p-1);%计算MSE%F检验F0=(RSS/p)/(ESS/(n-p-1));%计算F0Fa=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统计量对应的概率%F检验F0=(RSS/p)/(ESS/(n-p-1));%计算F0Fa=finv(p,n-p-1,0.95);%F分布时的临界值%F检验F0=(RSS/p)/(ESS/(n-p-1));%计算F0Fa=finv(p,n-p-1,0.95);%F分布时的临界值%可决系数检验R2=RSS/TSS;%计算样本决定系数程序的输出结果列在表3.12,3.13中表3.12方差分析表方差来源平方和自由度均方和F值p值回归400051341000128.161172.6560误差81096.389145792.599　　总计408160918　　　表3.13回归系数　变量β值差t值p值常数项345.25150.3222.2970.038省人均GDP0.1670.0443.8120.002第二产业增加值0.1960.0822.390.031服务业就业人数-0.7010.216-3.2420.006服务业资本形成总额-0.6540.295-2.2150.044该方程的拟合优度判定系数调整后的拟合优度判定系数由此说明该多元线性回归方程的拟合程度比较理想。F检验：；从方方差分析表可知统计量：F0=172.656，给定=0.05，查分布表，得到一个临界值F=3.1122，因为F0>F,或者由F0的p值为0小于0.05，所以拒绝H0，接受备择假设H1，说明总体回归系数i不全为零，即表明模型的线性关系在95%的置信水平下显著成立。t检验：对于=0.05，从表3.13最后一列的概率可以看出均小于0.05，所以拒绝H0，接受备择假设，即回归系数i（i=0,1,2,3,4），显著不为零。方法2：从包含全部变量的回归方程中逐次剔除不显著因子。首先建立包含全部变量的回归方程，然后对每一个因子作显著性检验，剔除不显著因子中偏回归平方和最小的一个因子，重新建立包含全部变量（剔除的除外）的回归方程。然后重复上面的过程，对新建立回归方程的每一个因子作显著性检验，剔除不显著因子中偏回归平方和最小的因子，再重新建立回归方程。如此，当新建立回归方程中所有因子都显著时，回归方程就是“最优”的了。这种方法在因子，特别是不显著因子不多时，可以采用。但计算的工作量仍然可能较大。3.3逐步回归在建立经济预测问题的数学模型时，常常从可能影响预测量Y的许多因素中挑选一批因素作为自变量，应用回归分析的方法建立回归方程作预报或控制用。问题是如何在为数众多的因素中挑选变量，以建立我们称为是这批观测数据“最优”的回归方程。3.3.1最优回归方程的选择选择“最优”回归方程有以下几种方法：方法1：从所有可能的变量组合的回归方程中挑选最优者，即把所有包含1个，2个，…直至所有变量的线性回归方程全部计算，对每个方程及自变量作显著性检验，然后从中挑选一个方程，所有的变量全部显著，且剩余均方和MSE较小。这种方法只在变量较少时可用。方法3：从一个变量开始，把变量逐个引入回归方程。这一方法首先计算各因子与Y的相关系数，将绝对值最大的一个因子引入方程，并对回归平方和进行检验，若显著则引入。然后找出余下的因子中与Y的偏相关系数最大的那个因子，将其引入方程，检验显著性，等等，当引入的因子建立的方程检验不显著时，该因子就不再引入。这种方法尽管工作量较小，但并不保证最后所得到的方程是“最优”的，还得进一步作检验，剔除不显著因子。同时这种方法每一步要计算偏相关系数，也较麻烦。结合方法２与３，产生了一种建立“最优”回归方程的方法——逐步回归分析。逐步回归的基本思想是，将变量一个一个引入，引入变量的条件是偏回归平方和经检验是显著的，同时每引入一个新变量后，对已选入的变量要进行逐个检验，将不显著变量剔除，具体做法是将变量一个一个引入，当每引入一个自变量后，对已选的变量要进行逐个检验，当原引入的变量由于后面变量的引入而变得不再显著时，要将其剔除。引入一个变量或从回归方程中剔除一个变量，为逐步回归的一步，每一步都要进行F检验，以确保每次引入新的变量之前回归方程中只包含显著的变量。这个过程反复进行，直到既无显著的自变量选入回归方程，也无不显著自变量从回归方程中剔除为止。这样保证最后所得的变量子集中的所有变量都是显著的。这样经若干步以后便得“最优”变量子集。3.3.2逐步回归的MATLAB方法统计工具箱中用作逐步回归的是命令stepwise，它提供了一个交互式画面，通过这个工具你可以自由地选择变量，进行统计分析，其通常用法是stepwise(X,Y,in,penter,premove)其中X是自变量数据，Y是因变量数据，分别为矩阵，in是矩阵X的列数的指标，给出初始模型中包括的子集，缺省时设定为全部自变量不在模型中，penter为变量进入时显著性水平，缺省时=0.05，premove为变量剔除时显著性水平，缺省=0.10在应用stepwise命令进行运算时，程序不断提醒将某个变量加入（Movein）回归方程，或者提醒将某个变量从回归方程中剔除（Moveout）。注意：应用stepwise命令做逐步回归，数据矩阵X的第一列不需要人工加一个全1向量，程序会自动求出回归方程的常数项（intercept）。下面通过一个例子说明stepwise的用法。例3.3.1（Hald,1960）Hald数据是关于水泥生产的数据。某种水泥在凝固时放出的热量(单位：卡／克)与水泥中4种化学成分x1~x4所占的百分比有关：在生产中测得13组数据，见表3.14，试建立y关于这些因子的“最优”回归方程。表3.14水泥生产的Hald数据序号12345678910111213x17111117113122111110x226295631525571315447406668x3615886917221842398x46052204733226442226341212y78.574.3104.387.695.9109.2102.772.593.1115.983.8113.3109.4解：在Matlab软件包中写一个M文件“liti3_3_1.m”X=[7,26,6,60;1,29,15,52;11,56,8,20;11,31,8,47;7,52,6,33;11,55,9,22;3,71,17,6;1,31,22,44;2,54,18,22;21,47,4,26;1,40,23,34;11,66,9,12;10,68,8,12];%自变量数据Y=[78.5,74.3,104.3,87.6,95.9,109.2,102.7,72.5,93.1,115.9,83.8,113.3,109.4]';%因变量数据stepwise(X,Y,[1,2,3,4],0.05,0.10)%in=[1,2,3,4]表示x1、x2、x3、x4均保留在模型中程序执行后得到下列逐步回归的窗口(如图3.23所示)：图3.23逐步回归窗口在图3.23中，用蓝色行显示变量x1、x2、x3、x4均保留在模型中，窗口的右侧按钮上方提示：将变量x3剔除回归方程(Movex3out)，点击NextStep按钮，即，进行下一步运算，将第3列数据对应的变量x3剔除回归方程。点击NextStep按健后，剔除的变量x3所对应的行用红色表示，同时又得到提示：将变量x4剔除回归方程(Movex4out)，点击NextStep按钮，即，进行下一步运算，将第4列数据对应的变量x4剔除回归方程。点击NextStep按健后，x4所对应的行用红色表示，同时提示：MoveNoterms，即，没有需要加入(也没有需要剔除)的变量了(如图3.24所示)。在Matlab7.0软件包中，可以直接点击“AllSteps”按钮，直接求出结果(省略中间过程)。图3.24逐步回归结果由图3.24，最后得到回归方程（蓝色行是被保留的有效行，红色行表示被剔除的变量）：由图3.24，最后得到回归方程（蓝色行是被保留的有效行，红色行表示被剔除的变量）：图3.24中显示了模型参数R2=0.978698，修正的R2=0.972282，F_检验值＝229.504，与显著性概率相关的p值＝4.40658e-009<0.05，残差均方RMSE＝2.40634（这个值越小越好）。以上指标值都很好，说明回归效果比较理想。另外，截距intercept=52.5773，这就是回归方程的常数项。逐步回归窗口中对已建模型给出了在线与超链接的显示功能，当将光标指向“ModelHistory”框中的均方残差RMSE的第一个蓝色点时，光标在线显示“inmodel:x1、x2、x3、x4”，若双击光标，则超链接到图3.23所示的逐步回归窗口。从“ModelHistory”框中可以观察不同模型的均方残差变化。历史考题：某厂5年间工业增加值与劳动生产率的资料如下表所示(1)计算(2)求y对x的线性回归方程：(3)F0.05(1,3)=10.13对建立的方程进行显著性检验年份12345工业增加值y（万元）203570105120劳动生产率x（万元/人）1.52.03.55.06.0所以y对x的线性关系显著由于F=531.3055>10.13=F0.05(1,3)(3)故回归方程为：解(1)谢谢！THANKYOU！结束