清华数学实验12-回归分析

实验12回归分析化学工程系化71李骥聪2007011861【实验目的】1．了解回归分析的基本原理，掌握MATLAB的实现方法；2．练习用回归分析方法解决实际问题。【实验内容】1.题目1.用切削机床加工时，为实时地调整机床需测定刀具的磨损程度，每隔一小时测量刀具的厚度得到以下的数据（见表1），建立刀具厚度对于切削时间的回归模型，对模型和回归系数进行检验，并预测7.5h和15h后刀具的厚度，用（30）和（31）式两种办法计算预测区间，解释计算结果。表1.时间/h012345678910刀具厚度/cm30.629.128.428.128.027.727.527.227.026.826.5模型及其求解利用MATLAB将刀具厚度和时间的关系绘图。输入程序如下：x1=0:10;y=[30.629.128.428.128.027.727.527.227.026.826.5]';plot(x1,y,'+')绘图得：从图中可以看出，刀具厚度y与时间x呈线性关系。建立模型为01yx。利用MATLAB计算模型和模型参数，并对回归系数进行检验。输入程序如下：alpha=0.05;n=11;X=[ones(n,1),x1'];[b,bint,r,rint,s]=regress(y,X,alpha);state=sum(r.^2)/(n-2);%剩余方差b,bint,s,state,rcoplot(r,rint)输出结果为：b=29.5455-0.3291bint=28.976930.1140-0.4252-0.2330s=0.869660.00180.0000state=0.1985整理得：表1.2.原始数据回归分析结果回归系数回归系数估计值回归系数置信区间029.5455[28.976930.1140]1-0.3291[-0.4252-0.2330]2R=0.8696，F=60.0018，p0.0001，2s=0.1985从以下几点可检验模型是有效的：1的置信区间不含零点；p；用MATLAB命令finv(0.95,1,11-2)计算得到(1,2),1nF=5.1174F。故拒绝0H，回归模型成立。残差及其置信区间如下图：图中第1个点残差的置信区间不包含零点，可以认为这个数据是异常的（残差应服从均值为0的正态分布），将第一个数据剔除后，重新计算，输入程序如下：x1=1:10;y=[29.128.428.128.027.727.527.227.026.826.5]';alpha=0.05;n=10X=[ones(n,1),x1'];[b,bint,r,rint,s]=regress(y,X,alpha);state=sum(r.^2)/(n-2);%剩余方差b,bint,s,state,rcoplot(r,rint)得到如下结果：表1.3.剔除第1个数据后的计算结果回归系数回归系数估计值回归系数置信区间029.0533[28.833429.2732]1-0.2588[-0.2942-0.2233]2R=0.9726，F=283.5599，p0.0001，2s=0.01952残差及其置信区间如下图：图中第2个数据（1，29.1）残差的置信区间也不包括零点，仍该视为异常点。剔除再次计算，得到如下结果：表1.4.剔除第2个数据后的计算结果回归系数回归系数估计值回归系数置信区间028.8667[28.779628.9537]1-0.2333[-0.2467-0.2200]2R=0.99593，F=1715.0，p0.0001，2s=0.0019残差及其置信区间如下图：这次的数据残差的置信区间全部包括零点，无异常点。对比分析：对比剔除前后的变化，发现置信区间明显缩小，决定系数和F都明显变大。表明异常点的剔除有利于更好的建立模型。回归模型为：/cm28.86670.2333/hyx。输出最终计算结果和图形（其中b为最后输出是的b）：x1=2:10;y=[28.428.128.027.727.527.227.026.826.5]';alpha=0.05;n=9X=[ones(n,1),x1'];[b,bint,r,rint,s]=regress(y,X,alpha);x2=0:10;y2=[30.629.128.428.128.027.727.527.227.026.826.5]';x0=0:0.1:12;y0=b(1)+b(2)*x0;plot(x2,y2,'+',x0,y0)y3=b(1)+b(2)*7.5y4=b(1)+b(2)*15解得：y(7.5)=27.1167y(15)=25.3667原始数据和回归分析结果（剔除第1、2个点后）比较显著性水平下的预测区间为：22000(2),1/20(2),1/2()()11ˆˆ[1,1]nnxxxxxxxxytsytssnsn……课本(30)式01/201/2ˆˆ[,]yusyus……课本(31)式利用上述两式计算预测区间，由于在回归分析中剔除了第1、2点，所以在预测时也不考虑这两点。输入程序如下：x1=2:10;y=[28.428.128.027.727.527.227.026.826.5]';alpha=0.05;n=9;X=[ones(n,1),x1'];[b,bint,r,rint,s]=regress(y,X,alpha);state=sum(r.^2)/(n-2);%剩余方差t=tinv(1-alpha/2,n-2);u=norminv(1-alpha/2,0,1);S=state^0.5;%剩余标准差xbar=mean(x1);Sxx=var(x1);x0=[7.5,15]y0=b(1)+b(2).*x0Y3=y0-t*S*((x0-xbar).^2/Sxx+1/n+1).^0.2;Y4=y0+t*S*((x0-xbar).^2/Sxx+1/n+1).^0.2;Y5=y0-u*S;Y6=y0+u*S;yt1=[Y3(1),Y4(1)]yt2=[Y3(2),Y4(2)]yu1=[Y5(1),Y6(1)]yu2=[Y5(2),Y6(2)]输出结果为：yt1=[27.006127.2272]yt2=[25.197325.5361]yu1=[27.031127.2022]yu2=[25.281125.4522]整理得到：预测方式点估计ˆy20(2),1/2()1()1nxxxxxtssn1/2()xusx=7.527.1167[27.006127.2272][27.031127.2022]x=1525.3667[25.197325.5361][25.281125.4522]可以看到用（31）式求得预测区间比用（30）式求得预测区间小。这是因为用（31）式预测时，1/2()xus的要求是n，显然这里n有限，用（30）式预测不能达到（31）式理想条件下预测的精确，所以用1/2()xus会得到较小的预测区间。结果分析：（1）模型建立前应对原始数据有一定认识：在建立模型之前应作图，直观的了解数据之间的关系。通过对数据关系的分析，可以建立线性模型。客观规律来讲：在刀具没有大幅度磨损之前，磨损速度应该是一定的，可以建立起线性模型。（2）模型检验与诊断异常点的存在不同：当1的置信区间不含零点；p；(1,2),1nF=5.1174F，几项指标都达到要求时，拒绝0H，回归模型成立，模型检验结束，不管原始数据点有没有异常，线性模型都成立。也就是说：异常点的存在并不能从根本上动摇和改变合理的模型，而不合理的模型也不会通过几个异常点的剔除就变得合理的。从显著性水平也可以证明这一点，即残差置信区间的显著性水平和模型检验的显著性水平并没有关系，是相互独立的。（3）从题中可以看到，两个异常点不是一次剔除的。在剔除第1个点之前，第2个点还是合理的，但剔除后就不合理了。从原始数据会出的图上可以看出，前两个点都有明显偏离，但由于它们的存在，影响回归模型发生的变化，使得第二个点恰好在回归模型上；而第1个点被剔除后，第2个点的异常性就显露出来。这也就是为什么在原始数据的残差区间图中第2点并不是异常点的原因。（4）考虑到实际情况，在刚开始工作时，刀具还没有进入稳定工作状态，致使开始时磨损过快，进入稳定工作状态后，磨损快慢应该达到设计值，而这个设计值应该与我们建立的线性模型相近。2.题目2.电影剧院调电视广告费用和报纸广告费用对每周收入的影响，得到下面的数据（见表2），建立回归分析模型并进行检验，诊断异常点的存在并进行处理。表2.每周收入9690959295959494电视广告费用1.52.01.52.53.32.34.22.5报纸广告费用5.02.04.02.53.03.52.53.0模型及其求解利用MATLAB将每周收入和电视广告费用、报纸广告费用的关系绘图。输入程序如下：y=[9690959295959494]';x1=[1.52.01.52.53.32.34.22.5]';x2=[5.02.04.02.53.03.52.53.0]';plot3(x1,x2,y,'+')绘图得：图1.绘图直接得到图2.转换一点角度后得到从图2可以看出每周收入y和电视广告费用1x、报纸广告费用2x呈一定线性关系。假定其满足线性关系，建立二元线性回归模型：01122yxx。利用MATLAB计算模型和模型参数，并对回归系数进行检验。输入程序如下：y=[9690959295959494]';x1=[1.52.01.52.53.32.34.22.5]';x2=[5.02.04.02.53.03.52.53.0]';n=8;m=3X=[ones(n,1)x1x2];[b,bint,r,rint,s]=regress(y,X);state=r'*r/(n-m);%剩余方差b,bint,s,state,rcoplot(r,rint)输出结果为：b=83.21161.29852.3372bint=78.805887.61740.40072.19621.48603.1883s=0.908924.94080.0025state=0.4897整理得：回归系数回归系数估计值回归系数置信区间083.2116[78.805887.6174]11.2985[0.40072.1962]22.3372[1.48603.1883]2R=0.9089，F=24.9408，p=0.0025，2s=0.4897从以下几点可检验模型是有效的：1、2的置信区间不含零点；p；用MATLAB命令finv(0.95,3-1,8-3)计算得到(2,3),1nF=5.7861F。故拒绝0H，回归模型成立。残差及其置信区间如下图：图中第1个点残差的置信区间不包含零点，可以认为这个数据是异常的（残差应服从均值为0的正态分布），将第一个数据剔除后，重新计算得到：b=81.48811.28772.9766bint=78.787884.18830.79641.77902.32813.6250s=0.976884.38420.0005state=0.1257整理得：回归系数回归系数估计值回归系数置信区间081.4881[78.787884.1883]11.2877[0.79641.7790]22.9766[2.32813.6250]2R=0.9768，F=84.3842，p=0.0005，2s=0.1257剔除第1个数据后的残差及其置信区间如下图：这次的数据残差的置信区间全部包括零点，无异常点。对比分析：对比剔除前后的变化，发现置信区间明显缩小，决定系数2R和F明显增大，p，2s明显减小，表明异常点的剔除有利于更好的建立模型。二元线性模型为：1281.48811.28772.9766yxx。结果分析：（1）将（1(1)x，2(1)x）代入计算可得到ˆ(1)y=98.3027，原始值为96，所以在原始残差图中，96为异常点且偏小。（