17实验十七 回归分析

实验十七回归分析一实验目的学习用MATLAB求解一元线性回归问题.学会正确使用命令Regress,并从输出表中读懂线性回归模型中各参数的估计，回归方程，线性假设的显著性检验结果,因变量Y在观察点％的预测区间等.二实验准备多元线性回归的数学模型为Y=XB+8，其中X=8~N(0qY=XB+8，其中X=8~N(0q2I)「1XX・・・X「y]「8一1121p1111XXX.・.，Y=y81222p22，8=•2・・.・・・・::1XX・・・XLy」8L8」1—1n2npnnn，记X，残差的每一分量定义pipa*=P0P=(P0'吧,…'P)T为参数为r=▲-*.，i=1'2''n.残差平方和Q=£r=£(y-V)2iiii=1i=1P=(P0,P；..,P「T使残差平方和最小.S=S=£(i=1111最小二乘估计就是选择参数七-y)2称为总的偏差平方和；U=£(g-y)2称为回归平方和；i=1可以证明平方和的分解公式：S=Q+U.决定系数R2=US是反映模型是否有效的指标之一.0vR2v1，R2越接近于1，模型的有效性越好.yyF=Q洋p-1)~F(P'"-PT)，是检验Ho：P0=Bi=-=Pp=0的检验统计量，在F的值偏大时拒绝Hq.寸二r2QiS2==~T是回归模型中。2的无偏估计.n-p-1n-p-1一元线性回归是多元线性回归的特例，只要令p=1.一元线性回归模型中个值y0的预测区间当H0被拒绝，即判断模型有效后，就可以从自变量x的一个给定值君预测因变量的理论值y0，预测值实际上是一个点估计，记作y0.显然可取y°=B°+6]x0(i)这里的y有如下的统计意义：y是无偏估计，即Ey=y.在给定显著性水平a下，y预00000

测区间(a=0.05时Vo一"2侦-2)邛,y有95%的可能落在此区间)为01+—+—0(n-2)sxx其中W(y-yIii=-4=1n—2s=Wxx(x0测区间(a=0.05时Vo一"2侦-2)邛,y有95%的可能落在此区间)为01+—+—0(n-2)sxx其中W(y-yIii=-4=1n—2s=Wxx(x0—x).虽然这个结果比较复杂，但是当n很_i=1大且x0接近x时，可以忽视上式根号内的后两项，且t(n—2)接近于N(0,1)的1—a2分1—a2如果因为F<F1a(1,n—2)或p1的置信区间包含零点而接受H0，那也只是说明y与x之间没有合适的线性模型，但是两者之间可能存在其它关系.°多元线性回归模型中个值y0的预测区间当模型通过有效性检验后，可由自变量的任一给定值x0=1。「…,x。,)预测因变量的理论值y，记作y,显然°y=p+&x+...+Bx00101p0p在给定显著性水平a下y的预测区间为[y—6Cx)y+8Cx)]，§(x)=t.(n—p—1)s^1+xrTXTXT1寸乙r2…Qi其中s2=n_,—1=n-1p—1•当n很大且x0接近x时,上述预测区间也可简化为y—z..s,y+zq1-01—a201—a2」三学习MATLAB命令1.MATLAB用于一元或多元线性回归的命令是regress，其格式为b=regress(y,X)%简单形式[b,bint,r,rint,stats]=regress(y,X,alpha)%完整形式输入参数及意义：因变量y(列向量)，矩阵X(n行p+1列，左边第一列为1，这时回归模型中有常数项)，alpha是显著性水平a(默认时设定为0.05).输出参数及意义3=(P0,料,…,P/为p+1个模型参数的最小二乘估计，bint是参数(Po’P],…,P)t的置信区间，r是残差(列向量)，rint是残差的置信区间，stats包含4个统计量：第1个是决定系数R2；第2个是作回归方程的显著性检验的F(p,n—p—1)统计量的值;第3个是F(p,n—p—1)分布大于F值的概率p(检验的P-值)，p<a时拒绝H°，回归模型有效；第4个是s2(剩余方差，MATLAB7.0以前版本没有第4个统计量，s2可由sum(r."2)/(n-pT)计算.2.多项式回归命令polyfit，其格式是a=polyfit(x,y,m)输入参数及意义：x,y是要拟合的数据组（长度相同的向量），m为拟合多项式的次数.输出a为拟合多项式y=axm+axm-ihfax+a的系数a=a1,a2,…,a,a+1]（降幕排列），注意：顺序与regress的输出相反.下面的命令常与polyfit连用，计算上述多项式在x处的值y：y=polyval（a,x）3.非线性回归命令nlinfit，其格式是[b,R,J]=nlinfit（x,y,’model’,b0）输入参数及意义：乂是自变量数据矩阵，每列一个变量；y是因变量数据向量；model是模型的函数名（M文件或inline形式的函数），形式为y=f（b,x）,b为参数，b0为参数迭代初始值.输出参数及意义：b是参数的估计值；R是残差，J返回用于估计预测误差的Jacobi矩阵.四实验内容一元线性回归例1表17.1列出了18个5-8岁儿童的重量（这是容易测得的）和体积（这是难以测量的）.表17.1重量x（千克）17.110.513.815.711.910.415.016.017.8体积y（立方分米）16.710.413.515.711.610.214.515.817.6重量x（千克）15.815.112.118.417.116.716.515.115.1体积y（立方分米）15.214.811.918.316.716.615.915.114.5（1）画出散点图.（2）求y关于x的线性回归方程y=a+bx,并作回归分析.（3）求x=14.0时y的置信水平为0.95的预测区间.解（1）输入数据，并输入作散点图命令x=[17.110.513.815.711.910.415.016.017.815.815.112.118.417.116.716.515.1];%输入自变量x的观察值y=[16.710.413.515.711.610.214.515.817.615.214.811.918.316.716.615.914.5];%输入因变量y的观察值plot（x,y,'*'）%作散点图输出结果是图17.1191「17'TOC\o"1-5"\h\z16'丁'15■++14-+1312■+10111213141516171819图17.1(2)作一元回归分析，输入n=length(y);X=[ones(n,1),x'];%1与自变量X组成的输入矩阵[b,bint,r,rint,s]=regress(y',X);b,bint,srcoplot(r,rint)%残差及其置信区间作图执行后得到输出结果：b=-0.10400.9881bint=-0.76550.55740.94451.0316s=1.0e+003*0.00102.31190这个结果可整理如表17.2表17.2例1(重量与体积)的计算结果回归系数回归系数估计值回归系数置信区间P0-0.1040[-0.7655,0.5574]&10.9881[0.9445,1.0316]R2r1F=2311.9p<0.001元回归方程为y=—0.1040+0.9881尤;从几个方面都可以检验模型是有效的:F检验的P-值接近于零；P1的置信区间不含零点;p<a；用MATLAB命令finv(0.95,1,16)计算得到F(1,16)=4.4940<F.残差及置信区间如图17.2^9524681012141618CaseNumber图17.2(3)将上面的回归系数估计值P=—0.1040,P=0.9881代入回归方程，对重量为0114.0kg的儿童的体积(个值)进行预测，得%=13.7287.在a=0.05下，体积预测区间简化为「y—uy+u..s］，输入计算指令1-01—a201—a2」t1=13.7287-norminv(0.975,0,1)*sqrt(sum(r.”2)/16)

t2=13.7287+norminv（0.975,0,1）*sqrt（sum（r.”2）/16）输出t1=13.3333t2=14.1241即14.0kg的儿童的体积（个值）的置信度为0.95的预测区间为[13.3340,14.1248].也可以用命令polytool（x,y,1,0.05）作出散点图及拟合直线，并对x=14.0时的y进行预报（输出结果略，详见例4）.iIMATLAB将它们作图，如图17.3.从图形直观地看，y与x大致呈线性关系，即iIMATLAB将它们作图，如图17.3.从图形直观地看，y与x大致呈线性关系，即y二。。+。产根据试验数据从统计推断的角度讨论P0,61的置信区间和假设检验，进而对任意的年龄x给出y的预测区间，属于一元线性回归分析.序号血压年龄序号血压年龄序号血压年龄1144391116264211363622154712150562214250313845131405923120394145471411034241202151626515128422516044614246161304826158537170671713545271446381244218114182813029915867191162029125251015456201241930,、17569表17.330个人的血压与年龄分析记血压（因变量）），'年龄（自变量）工表17.3的数据为C.,y*=1,2,MATLAB编程如下：y=[144215138145162142170124158154162150140110128130135114116124136142120120160158144130125175];x=[394745476546674267566456593442484518201936503921445363292569];n=30;%1与自变量组成的输入矩阵X=[ones(n,1),x'];%1与自变量组成的输入矩阵[b,bint,r,rint,s]=regress(y',X);%血压与年龄的散点图%残差及其置信区间作图%血压与年龄的散点图%残差及其置信区间作图b,bint,srcoplot(r,rint)b=98.40840.9732bint=78.7484118.06830.56011.3864s=0.454023.28340.0000这个结果可以整理如表17.4表17.4年龄与血压的计算结果回归系数回归系数估计值回归系数置信区间P098.4084[78.7484118.0683]Pi0.9732[0.56011.3864]R2=0.4540F=23.2834p<0.0001从几个方面都可以检验模型是有效的：F检验的P-值接近于零；P1的置信区间不含零点；p<a；用MATLAB命令finv（0.95,1,28）计算得到F095G,28）=4.1960<F.但是P1的置信区间较长，R2较小，说明模型精度不高..残差及置信区间如图17.4,图中第二个点仁2,七）残差的置信区间不包括零点，可以认为这个数据是异常的（残差应服从均值为0的正态分布），它偏离数据整体的变化趋势，给模型的有效性和精度带来不利的影响，称异常点或离群点，应予以剔除将原始数据中的第二个剔除后重新计算得到表17.5表17.5剔除第二个数据后的计算结果回归系数回归系数估计值回归系数置信区间P096.8665[85.4771108.2559]P10.9533[0.71401.1925]R2=0.7123F=66.8358p<0.0001可以看出6,6变化不大，但置信区间变短，R2和F变大，说明模型精度提高.残差01及其置信区间见图17.5.从图17.5又发现两个新的异常点（为什么会出现？），剔除它们后模型会有什么改进？

图17.4例2的残差及其置信区间图17.5例2剔除^2,>2）的残差及其置信区间将上面得到的回归系数的估计值P=96.8665,p=0.9533代入回归方程，对50岁

01（x0=50）人的血压进行预测，得y0=144.5298.在侦=0.05下按照（2）式计算的预测区间为[124.5406,164.5190],按照简化的（3）式计算的预测区间为[125.7887,163.2708].多元线性回归例3世界卫生颁布的“体重指数”的定义是体重除以身高的平方，显然它比体重本身更能反映人的胖瘦.对表17.3给出的三个人又测量了他（她）们的体重指数如表17.6所示.试建立血压与年龄（岁）和体重指数之间的模型，作回归分析.如果还记录了他（她）们的吸烟习惯（表17.6中0表示不吸烟，1表示吸烟），怎样在模型中考虑这个因素，吸烟会使血压升高吗？对50岁、体重指数25的吸烟者的血压作预测.表17.630个人的血压与年龄、体重指数、吸烟习惯模型分析记血压J，年龄尤模型分析记血压J，年龄尤，体重指数x，吸烟习惯x，用MATLAB将J与尤的1232数据作散点图（略），看出大致也呈线性关系，建立模型y=。+Px+px+px,由数0112233据估计系数p,p,p,P也可以看作曲线（曲面）拟合0123仍然用命令regress（y,X）,只需注意输入矩阵X的构成，编程如下：y=[144116124136x1=[3944536329215138145142120120474547652569];162142170160158144130125175];46674267124158154162150140110128130135114566456593442484518201936503921序号血压年龄体重指数吸烟习惯序号血压年龄体重指数吸烟习惯序号血压年龄体重指数吸烟习惯11443924.20111626428.01211363625.0022154731.11121505625.80221425026.2131384522.60131405927.30231203923.5041454724.01141103420.10241202120.3051626525.91151284221.70251604427.1161424625.10161304822.21261585328.6171706729.51171354527.40271446328.3081244219.70181141818.80281302922.0191586727.21191162022.60291252525.30101545619.30201241921.50301756927.41解x2=[24.231.122.624.025.925.129.519.727.219.328.025.827.320.121.727.418.822.621.525.026.223.520.327.128.628.322.025.327.4];x3=[010110101010000100000100110101];n=length(y);%输入矩阵X=[ones(n,1),x1',x2',x3'];%输入矩阵[bbintrrints]=regress(y',X);b,bint,s,rcoplot(r,rint)输出b=45.36360.36043.090611.8246bint=3.553787.1736-0.07580.79651.05305.1281-0.148223.7973s=0.685518.89060.0000计算得到的数据如表17.7表17.7例3的计算结果回归系数回归系数估计值回归系数置信区间P045.3636[3.553787.1736]P10.3604[-0.07580.7965]P23.0906[1.05305.1281]P311.8246[-0.148223.7973]R2=0.6855F=18.8906p<0.0001从残差及其置信区间的图(略)发现第2和第10个点为异常点，剔除它们重新计算得到表17.8.表17.8例3剔除第2和第10个数据后的计算结果回归系数回归系数估计值回归系数置信区间P058.5101[29.906487.1138]P10.4303[0.12730.7332]P22.3449[0.85093.8389]P310.3065[3.387817.2253]R2=0.8462F=44.0087p<0.0001由剔除异常点后数据得到的预测模型为宁=58.5101+0.4303x+2.3499x+10.3065x.123由这个结果可知，年龄和体重指数相同的人，吸烟者比不吸烟者的血压(平均)高10.3.对50岁、体重指数25的吸烟者的血压作预测：将气=50,x2=25,x3=1代入上面的预测模型得V=148.9525，在a=0.05下按照(5)式计算得预测区间为［1；3.1716,164.7334］，按照简化得(6)式计算得预测区间为［134.5951,163.3099］.多项式回归例4一种合金在某种添加剂的不同浓度下,各做三次试验，得到数据如下：浓度x10.015.020.025.030.025.229.831.231.729.4抗压强度y27.331.132.630.130.828.727.829.732.332.8(1)作散点图.(2)以模型Y=«+«x+b2x2+s,s~N(0q2)拟合数据，其中b,b,b舟2与x无关.求回归方程Y=b+bx+bx2,并作回归分析.012012解(1)输入数据与作散点图命令x=[10.010.010.015.015.015.020.020.020.025.025.025.030.030.030.0];y=[25.227.328.729.831.127.831.232.629.731.730.132.329.430.832.8];plot(x,y,'*')输出为图17.6

图17.6从图17.6可以看出y与x没有明显的线性关系.(2)依题意，建立回归模型Y=«+七x+b2x2+s(7)如记气=x,x2=x2，上回归模型仍属于多元线性模型，可以用MATLAB命令regress求解，输入n=length(x);X=[ones(n,1),x',(x."2)'];[b,bint,r,rint,s]=regress(y',X);b,bint,s输出为19.03330086-0.020489220.2320-0.039626.17451.7852-0.00120.61409.5449bints=bbints=0.0033将回归系数的估计值直接代入(7)式，得到J的预测方程为y=19.0333+1.0086x-0.0204x2回归系数置信区间和统计量等输出从略(检验结果模型是有效的).(7)式又称为二项式回归，一元多项式回归模型的一般形式为y=g+gx++gxm+8可以直接使用MATLAB多项式回归命令polyfit计逾.输入a=polyfit(x,y,2)%2是多项式次数输出a=-0.02041.008619.0333可见多项式回归和多元线性回归得到同样的预测方程.用MATLAB求解一元多项式回归，除了polyfit(x,y,m)外，还有更方便的命令:polytool(x,y,2,0.05)

输入参数x,y,m同polyfit；alpha是显著性水平侦（默认时设定为0.05），输出一个如图17.7的交互式画面，实线为多项式的回归曲线y，它两侧的虚线表出〉的置信区间.x的数值可以用鼠标移动来改变，也可在图下的窗口内输入，图左边相应地给出y的预测值及其置信区间.通过图左下方的Export下拉菜单，还可以输出回归系数估计值及其置信区间、残差等.3426EnportDegree2362410121416183426EnportDegree236241012141618202224262830I~~XValues22图17.7添加剂浓度和抗压强度回归模型输出的交互式画面4.非线性回归例5酶是一种高效生物催化剂，催化条件温和，经过酶催化的化学反应称为酶促反应酶促反应中的反应速度主要取决于反应物（称为底物）的浓度，浓度较低时反应物速度大致与底物浓度成正比（称为一级反应），浓度较高、渐进饱和时反应速度趋向于常数（称为零级反应），二者之间有一过渡，根据酶促反应的这种性质、描述反应速度与底物浓度关系的一类模型是Michaelis-Menten模型：_。,x。+x其中y是反应速度，x是底物浓度，p2,p为待定参数.容易知道P是饱和浓度下的速度，121但是可以通过下面的变量代换称为最终反应速度，而p2是达到最终反应速度一半时的底物浓度，称为半速度点试估计模型中的参数p,p.酶促反应试验中反应速度与底物浓度数据如下：但是可以通过下面的变量代换底物浓度1-~~2—10.020.060.110.220.561.10反应速度6751848698115131124144158160解⑴Michaelis-Menten模型对参数P『P2是非线性的，化为线性模型1=-1+土1=。+。1y吧％xi2x画出反应速度的倒数Ly与底0.561.10];此模型中1y对新的参数@,气是线性的.根据试验数据，画出反应速度的倒数Ly与底0.561.10];clear;x=[0.020.020.060.060.110.110.220.220.56y=[6751848698115131124144158160];x1=1./x;y1=1./y;

n=length(x);X=[ones(n,1),x1'];[b,bint,r,rint,s]=regress(y1',X);formatlongb,bint,s输出b=0.006972133705670.00021500118582bint=0.005687613673840.00015891279429s=0.89310360382203%1与自变量组成的输入矩阵%长格式（16位）0.00825665373749n=length(x);X=[ones(n,1),x1'];[b,bint,r,rint,s]=regress(y1',X);formatlongb,bint,s输出b=0.006972133705670.00021500118582bint=0.005687613673840.00015891279429s=0.89310360382203%1与自变量组成的输入矩阵%长格式（16位）0.008256653737490.0002710895773575.193670898092550.00001156213074yy=b(1)+b(2)*xx;plot(x1,y1,'+',xx,yy)输出散点与拟合的直线图(图17.8)，可以发现1x较小时拟合较好，1x较大时出现较c1C大的分散.c1C大的分散.由0]=厂,P21的拟合图，输入xx=0:0.01:1.2;yy2=1/b(1)*xx./(b(2)/b(1)+xx);将&1，&2代入Michaelis-Menten模型，得到与原始数据比较%旗=1他,&,=0,他，&,，&,的估计值分别为143.43和乙乙乙0.0308plot（x,y,'+',xx,yy2）输出拟合图（图17.9）,可以发现在x较大时y的计算值比实际值要小，这是因为在对线性化模型作参数估计时，底物浓度x较小的数据在很大程度上控制了参数的确定，从而对底物浓度较大数据的拟合，出现较大的偏差.000.20.40.60.81.21.4000.20.40.60.81.21.4为解决线性化模型中拟合欠佳的问题，需要直接考虑拟合非线性模型(Michaelis-Menten模型).(2)用非线性拟合求解，编程如下：x=[0.020.020.060.060.110.110.220.220.560.561.10];y=[6751848698115131124144158160];b0=[1430.03];fun=inline(，b(1)*x./(b(2)+x)','b',，x，);[b,R,J]=nlinfit(x,y,fun,b0);bi=nlparci(b,R,J);xx=0:0.01:1.2;yy=b(1)*xx./(b(2)+xx);plot(x,y,，+，,xx,yy),pause%数据散点图和预测曲线nlintool(x,y,fun,b)%交互式画面b,bi得到b=160.27810.0477bi=145.6191174.93720.03010.0653数据散点图和预测(拟合)曲线见图17.10，可与图17.9比较.交互式画面见图17.11.利用交互式画面也可以进行预报.000.20.40.6000.20.40.60.811.21.4图17.10非线性模型的预测（拟合曲线）顷-50图17.11非线性模型的交互式画面顷-50图17.11非线性模型的交互式画面五实验作业1.某乡镇企业办的小加工厂其产品年销售额x与所获纯利润y从1984年到1994年的数据如下（单位：百万元）：年度8485868788899091929394销售额x6.17.59.410.714.617.421.124.429.832.934.3纯利润y4.56.48.38