实验11_统计回归模型(4学时)

1、数学建模实验王平实验11 统计回归模型（4学时）（第10章 统计回归模型）1. 牙膏的销售量p325332下面给出一组数据，其中：第1列 销售周期；第2列 某公司牙膏销售价格（元）x4；第3列 其它厂家平均价格（元）x3；第4列 广告费用（百万元）x2；第5列 价格差（元）x1(x3-x4)；第6列 销售量（百万支）y。存放在一个名为p325.txt的文件中。1 3.85 3.80 5.50 -0.05 7.382 3.75 4.00 6.75 0.25 8.513 3.70 4.30 7.25 0.60 9.524 3.70 3.70 5.50 0 7.505 3.60 3.85 7.00

2、0.25 9.336 3.60 3.80 6.50 0.20 8.287 3.60 3.75 6.75 0.15 8.758 3.80 3.85 5.25 0.05 7.879 3.80 3.65 5.25 -0.15 7.1010 3.85 4.00 6.00 0.15 8.0011 3.90 4.10 6.50 0.20 7.8912 3.90 4.00 6.25 0.10 8.1513 3.70 4.10 7.00 0.40 9.1014 3.75 4.20 6.90 0.45 8.8615 3.75 4.10 6.80 0.35 8.9016 3.80 4.10 6.80 0.30 8

3、.8717 3.70 4.20 7.10 0.50 9.2618 3.80 4.30 7.00 0.50 9.0019 3.70 4.10 6.80 0.40 8.7520 3.80 3.75 6.50 -0.05 7.9521 3.80 3.75 6.25 -0.05 7.6522 3.75 3.65 6.00 -0.10 7.2723 3.70 3.90 6.50 0.20 8.0024 3.55 3.65 7.00 0.10 8.5025 3.60 4.10 6.80 0.50 8.7526 3.65 4.25 6.80 0.60 9.2127 3.70 3.65 6.50 -0.05

4、8.2728 3.75 3.75 5.75 0 7.6729 3.80 3.85 5.80 0.05 7.9330 3.70 4.25 6.80 0.55 9.261.1（验证）基本模型p325329先保存上面的p325.txt文件。(1) 绘制y对x1的散点图程序如下：M=dlmread('p325.txt');%读取ASCII码文件x1=M(:,5); y=M(:,6);plot(x1,y, 'bo');提示：dlmread将以ASCII码分隔的数值数据文件读入到矩阵中dlmread：读取ASCII码文件的MATLAB函数M=dlmread('fun

5、.txt');fun.m是一个数据文件，存放一个数据矩阵，将文件内容写入M。(1) 运行程序并给出结果（比较327图1）：(2) 确定y对x1的拟合，绘制散点图与拟合曲线组合图形从y对x1的散点图可以发现，可用线性模型（直线）来拟合（其中是随机误差）。程序如下：clc; format short g;M=dlmread('p325.txt');%读取ASCII码文件x1=M(:,5); y=M(:,6);plot(x1,y, 'bo');b=regress(y,ones(size(x1),x1); % b=0 1 '，列向量x1=sort(x1)

6、; %按升序排序，用于画图y=ones(size(x1),x1*b;%使用矩阵乘法hold on;plot(x1,y, '-r');hold off;提示：regress多元线性回归函数调用格式b,bint,r,rint,stats=regress(y,x,alpha)例，多元回归模型为：输入：y为n（=30）维列向量数据。x为对应于回归系数 ( 0, 1, 2, 3 )' 的数据矩阵 1 x1 x2 x22（30×4矩阵，第1列全1）。alpha为置信水平（缺省时为0.05）。输出：b为=( 0, 1, 2, 3 )'估计值，4维列向量。bint为b

7、的置信区间，4×2矩阵。r为残差n(=30)维列向量y-x。rint为r的置信区间，30×2矩阵。stats为回归模型的检验统计量，含4个值：R2回归方程的决定系数（R是相关系数）F统计值P与F统计量对应的概率值s2剩余方差(2) 运行程序并给出结果（比较327图1）：(3) 绘制y对x2的散点图程序如下：clc; format short g;M=dlmread('p325.txt');%读取ASCII码文件x2=M(:,4); y=M(:,6);plot(x2,y,'bo');(3) 运行程序并给出结果（比较327图2）：(4) 确定y对

8、x2的的拟合，绘制散点图与拟合曲线组合图形从y对x2的散点图可以发现，可用二次函数模型来拟合。程序如下：clc;format short g;M=dlmread('p325.txt');%读取ASCII码文件x2=M(:,4); y=M(:,6);plot(x2,y,'bo');b=regress(y,ones(size(x2),x2,x2.2); % b=0 1 2'，列向量x2=sort(x2);y=ones(size(x2),x2,x2.2*b; %使用矩阵乘法hold on;plot(x2,y,'-r');hold off;(4)

9、 运行程序并给出结果（比较327图2）：(5) y对x1, x2的回归模型及其求解，销售量预测综上得回归模型变量x1, x2为回归变量，参数b0, b1, b2, b3为回归系数。程序如下：clc; format compact; format short g;M=dlmread('p325.txt');%读取ASCII码文件x1=M(:,5);x2=M(:,4); y=M(:,6);b,bint,r,rint,stats=regress(y,ones(size(x1),x1,x2,x2.2,0.05);fprintf('%2s%5s%11sn','参数

10、','估计值','置信区间');%1个汉字算1个字符for i=1:length(b) fprintf ('%1d%9.4f %7.4f, %7.4fn',i-1,b(i,:),bint(i,:);end % %d将i当整数输出，%7.4f按实数格式输出，区域宽7个字符，4位小数fprintf('nR2=%.4f F=%.4f p<%.4e s2=%.4fn',stats); x1=0.2; x2=6.5; y=1 x1,x2,x22*b; %使用矩阵乘法fprintf('n销售量预测：x1=%.1f, x2

11、=%.1f, y=%.4fn',x1,x2,y);提示：fprintf输出到命令窗口或写数据到文本文件见参考资料：MATLAB函数和命令的用法。(5) 运行程序并给出结果（比较328表2,329的预测结果）：1.2（验证，编程）模型改进p329332仍使用题1的数据。(1)（编程）y对x1, x2的回归模型的改进和求解，销售量预测改进的模型参考题1(5)的程序，编写一个类似的程序，运行结果与教材p329330的表3及相关结果相比较。(1) 给出程序和运行结果（比较329表3）：clc; format compact; format short g;M=dlmread('p325

12、.txt');%读取ASCII码文件x1=M(:,5);x2=M(:,4); y=M(:,6);b,bint,r,rint,stats=regress(y,ones(size(x1),x1,x2,x2.2,x1.*x2,0.05);fprintf('%2s%5s%11sn','参数','估计值','置信区间');%1个汉字算1个字符for i=1:length(b) fprintf ('%1d%9.4f %7.4f, %7.4fn',i-1,b(i,:),bint(i,:);end % %d将i当整数输出，

13、%7.4f按实数格式输出，区域宽7个字符，4位小数fprintf('nR2=%.4f F=%.4f p<%.4e s2=%.4fn',stats); x1=0.2; x2=6.5; y=1 x1,x2,x22,x1.*x2*b; %使用矩阵乘法fprintf('n销售量预测：x1=%.1f, x2=%.1f, y=%.4fn',x1,x2,y);(2)（验证）完全二次多项式模型运行以下程序（参考教材p331332）：clear; clc; format compact; format short g;M=dlmread('p325.txt'

14、); %读取ASCII码文件x1=M(:,5); x2=M(:,4); y=M(:,6);rstool(x1,x2,y,'quadratic')得以下的交互画面。画面中的两个座标系给出y的估计值和预测区间。用鼠标移动交互式画面中的十字线，或在图下方的窗口内输入，可改变x1和x2的数值。改变x1=0.2，x2=6.5，观察窗口左边的y估计值和预测区间。点击所得交互画面左下方的输出按钮“Export”，所得画面（导出到工作空间）第1个复选框是“将拟合参数存到一个名为beta的MATLAB变量中”，点击OK。在命令窗口提示符键入变量名beta将得到参数( 0, 1, 2, 3, 4,

15、 5 )' 的值。(2) 运行程序并给出结果：参数(0,1,2,3,4,5)' 的值（比较331）：2. 软件开发人员的薪金p332338在下面给出的数据中：（存入文件p333.txt）第1列 编号第2列 薪金y第3列 资历x1（从事专业工作的年数）第4列 管理x2（1表示管理人员，0表示非管理人员）第5列 教育x3，x4（1表示中学程度x3x4=10，2为大学x3x4=01，3为更高程度x3x4=00）01 13876 1 1 102 11608 1 0 303 18701 1 1 304 11283 1 0 205 11767 1 0 306 20872 2 1 207 1

16、1772 2 0 208 10535 2 0 109 12195 2 0 310 12313 3 0 211 14975 3 1 112 21371 3 1 213 19800 3 1 314 11417 4 0 115 20263 4 1 316 13231 4 0 317 12884 4 0 218 13245 5 0 219 13677 5 0 320 15965 5 1 121 12366 6 0 122 21352 6 1 323 13839 6 0 224 22884 6 1 225 16978 7 1 126 14803 8 0 227 17404 8 1 128 22184 8

17、 1 329 13548 8 0 130 14467 10 0 131 15942 10 0 232 23174 10 1 333 23780 10 1 234 25410 11 1 235 14861 11 0 136 16882 12 0 237 24170 12 1 338 15990 13 0 139 26330 13 1 240 17949 14 0 241 25685 15 1 342 27837 16 1 243 18838 16 0 244 17483 16 0 145 19207 17 0 246 19346 20 0 1假设满足多元线性回归模型y = 0 + 1 x1 +

18、2 x2 + 3 x3 + 4 x4 + 2.1（验证）基本模型p332335求回归系数及其置信区间（置信水平a = 0.05）、检验统计量R2、F、p、s2，有关散点图的MATLAB程序如下：%10.2 软件开发人员的薪金基本模型%模型：y=0+1*x1+2*x2+3*x3+4*x4+clear;clc;format compact;format short g;M=dlmread('p333.txt'); %读取ASCII码文件y=M(:,2); x1=M(:,3); x2=M(:,4);x3=M(:,5)=1; x4=M(:,5)=2; %教育程度a,aint,r,rin

19、t,stats=regress(y,ones(size(M,1),1) x1 x2 x3 x4 );fprintf('%2s%4s%9sn','参数','估计值','置信区间');%1个汉字算1个字符 for i=1:length(a) fprintf('%1d%7.0f %5.0f, %5.0fn',i-1,a(i,:),aint(i,:);endfprintf('nR2=%.3f F=%.0f p<%.4e s2=%.3en',stats); subplot(121);plot(x1,r,

20、'+'); %模型的残差与资历x1的关系 subplot(122);x234=(1+M(:,4).*(M(:,5)=1)+(3+M(:,4).*(M(:,5)=2).+ (5+M(:,4).*(M(:,5)=3);%见p336表3plot(x234,r,'+'); %模型的残差与管理-教育x2-x3,x4组合x234的关系 给出程序的运行结果（数据和图形）（比较334，335）：数据结果（比较334表2）：图形结果（比较335图1、图2）：2.2（编程）更好的模型p335336在题2.1的模型中增加x2与x3，x4的交互项后，多元线性回归模型为y = 0 + 1

21、 x1 + 2 x2 + 3 x3 + 4 x4 + 5 x2 x3 + 6 x2 x4 + 要求：同题2.1，区别在于模型不同，所以要根据新模型修改题2.1的程序，仍后运行。 给出程序和运行结果（程序、数据和图形）（比较336表4、图3、图4）：程序：%10.2 软件开发人员的薪金基本模型%模型：y=0+1*x1+2*x2+3*x3+4*x4+clear;clc;format compact;format short g;M=dlmread('p333.txt'); %读取ASCII码文件y=M(:,2); x1=M(:,3); x2=M(:,4);x3=M(:,5)=1;

22、x4=M(:,5)=2; %教育程度a,aint,r,rint,stats=regress(y,ones(size(M,1),1) x1 x2 x3 x4 x2.*x3 x2.*x4);fprintf('%2s%4s%9sn','参数','估计值','置信区间');%1个汉字算1个字符 for i=1:length(a) fprintf('%1d%7.0f %5.0f, %5.0fn',i-1,a(i,:),aint(i,:);endfprintf('nR2=%.3f F=%.0f p<%.4e s2

23、=%.3en',stats); subplot(121);plot(x1,r,'+'); %模型的残差与资历x1的关系 subplot(122);x234=(1+M(:,4).*(M(:,5)=1)+(3+M(:,4).*(M(:,5)=2).+ (5+M(:,4).*(M(:,5)=3);%见p336表3plot(x234,r,'+'); %模型的残差与管理-教育x2-x3,x4组合x234的关系数据结果（比较p336表4）：图形结果（比较p336图3、图4）：2.3（编程，验证）模型应用p337338在题2.1的模型中增加x2与x3，x4的交互项后，

24、多元线性回归模型为y =0 + 1 x1 + 2 x2 + 3 x3 + 4 x4 + 5 x2 x3 + 6 x2 x4 + (1)（编程）同题2.2，但要去除题2.1中给出数据中编号为33的数据（异常点）(1) 给出程序和运行结果（数据和图形）（比较337表5、图5、图6）：程序：%10.2 软件开发人员的薪金基本模型%模型：y=0+1*x1+2*x2+3*x3+4*x4+clear;clc;format compact;format short g;M=dlmread('p333.txt'); %读取ASCII码文件M(33,:)=;y=M(:,2); x1=M(:,3)

25、; x2=M(:,4);x3=M(:,5)=1; x4=M(:,5)=2; %教育程度a,aint,r,rint,stats=regress(y,ones(size(M,1),1) x1 x2 x3 x4 x2.*x3 x2.*x4);fprintf('%2s%4s%9sn','参数','估计值','置信区间');%1个汉字算1个字符 for i=1:length(a) fprintf('%1d%7.0f %5.0f, %5.0fn',i-1,a(i,:),aint(i,:);endfprintf('nR2

26、=%.3f F=%.0f p<%.4e s2=%.3en',stats); subplot(121);plot(x1,r,'+'); %模型的残差与资历x1的关系 subplot(122);x234=(1+M(:,4).*(M(:,5)=1)+(3+M(:,4).*(M(:,5)=2).+ (5+M(:,4).*(M(:,5)=3);%见p336表3plot(x234,r,'+'); %模型的残差与管理-教育x2-x3,x4组合x234的关系数据结果（比较p337表5）：图形结果（比较p337图5、图6）：(2)（验证） 6种管理教育组合人员的“基

27、础”薪金（资历为0）。组合管 理 教 育 “基础”薪金资历10 1 ？021 1 ？030 2 ？041 2 ？050 3 ？061 3 ？0利用(1)所求得的模型求解，程序如下（与教材p337表6比较）： clear; clc; format compact; format short g;M=dlmread('p333.txt'); %读取ASCII码文件M(33,:)=; %修改y=M(:,2); x1=M(:,3); x2=M(:,4);x3=M(:,5)=1; x4=M(:,5)=2; %教育程度a,aint,r,rint,stats=regress(y,ones(s

28、ize(M,1),1) x1 x2 x3 x4 x2.*x3 x2.*x4);x1=zeros(6,1); %资历为0x2=0 1 0 1 0 1' %管理x3=1 1 0 0 0 0' %教育x4=0 0 1 1 0 0'y=ones(6,1),x1,x2,x3,x4,x2.*x3,x2.*x4*a;fprintf('%2s%4s%4s%8sn', '组合', '管理', '教育', '“基础”薪金');%1个汉字算1个字符for k=1:6fprintf('%3d%6d%6d%1

29、1.0fn',k,x2(k),fix(k+1)/2),y(k);end(2) 运行程序并给出结果（比较337表6）：3. 酶促反应p338345嘌呤霉素实验中的反应速度与底物浓度数据底物浓度(ppm)0.020.060.110.220.561.10反应速度处理764797107123139159152191201207200未处理6751848698115131124144158160/3.1（编程）线性化模型p338340(1) 根据给出的数据，试编写一个程序绘制出如下图所示的散点图（比较：p339图1、图2）(1) 给出程序和运行结果（比较339图1、图2）：程序：clear; c

30、lc; format compact; format short g;M=dlmread('p338.txt'); %读取ASCII码文件y=M(2,1:2:11);y2=M(2,2:2:12);x1=M(1,1:6);x2=M(1,1:5);y4=M(3,1:2:11) ;y3=M(3,2:2:10);subplot(1,2,1)plot(x1,y, 'bo');hold on;plot(x1,y2, 'bo');title('y对x经处理的散点图');subplot(1,2,2)plot(x1,y4, '+')

31、;hold on;plot(x2,y3, '+');title('y对x未经处理的散点图');图形（比较p339图1、图2）：(2) 线性化模型Michaelis-Menten模型通过代换或可化为线性模型试编写一个程序，对经过处理的实验数据，求出：(a) 参数的估计；(b) 1/y与1/x的散点图和回归直线；(c) 用线性化得到的原始数据拟合图。程序运行的结果如下：p340表2p339图3、p340图4(2) 给出程序和运行结果（比较339，340）：程序：数据结果（比较p340表2）：图形(比较p339图3、p340图4)：3.2（验证）非线性模型及求解p34

32、0341Michaelis-Menten模型提示：求解非线性模型用到的操作(a) 非线性最小二乘拟合命令beta, R, J=nlinfit(x, y, 'model', beta0)输入：x为自变量数据矩阵，每列一个变量；y为因变量数据列向量；model为模型的m函数文件，形式为function y=f(beta,x)，beta为待估计参数；beta0为给定的参数初值；输出：beta为参数的估计值；R为残差；J为用于估计预测误差的Jacobi矩阵。(b) 参数的置信区间命令betaci=nlparci(beta,R,J)(c) 非线性模型的GUI工具预测命令nlintool(

33、x,y, 'model',beta)功能：给出一个交互式画面，拖动画面中的十字线可改变自变量x的取值，直接得到因变量y的预测值和预测区间，同时在左下方的Export中，可向MATLAB工作区传送统计数据。(1) 运行以下程序( p340 )：clc; clear; format compact; %命令m文件x=0.02 0.02 0.06 0.06 0.11 0.11 0.22 0.220.56 0.56 1.10 1.10'y=76 47 97 107 123 139 159 152 191 201 207 200'beta0=195.802 0.0484;

34、 %用线性化模型得到的作为非线性模型参数估计的初始迭代值。beta,R,J=nlinfit(x,y,'fun',beta0); %beta是行向量betaci=nlparci(beta,R,J); %betaci是2行2列矩阵fprintf('%2s%6s%14sn','参数','估计值','置信区间');%1个汉字算1个字符for i=1:length(beta) fprintf('%1d%11.5f %9.5f, %9.5fn',i,beta(i),betaci(i,:);endyy=beta(

35、1)*x./(beta(2)+x);plot(x,y,'o',x,yy,'+');nlintool(x,y,'fun',beta)function yhat=fun(beta,x) %必须作为单独一个函数式m文件yhat=beta(1)*x./(beta(2)+x);(1) 给出运行结果（数据和图形）（比较341）：数据结果（比较p341表3）：（Figure 1，比较p341图5）：拟合图（比较p341图6）：(2) 对上面程序运行中nlintool给出的交互式画面进行操作。单击画面左下方的导出按钮”Export”，出现下面画面，其中只选中Pa

36、rameters,Parameter CI和RMSE复选框，单击OK，则变量beta1,betaci1和rmse被送到工作空间中，在命令窗口键入这三个变量名显示其值。注：beta1为参数估计值，betaci1为参数的置信区间，rmse为剩余标准差。(2) 给出结果（beta1,betaci1和rmse的值）：3.3（验证）混合反应模型p342345模型其中：x1为底物浓度；x2为一示性变量（0-1变量），用来表示是否嘌呤霉素处理，令x2=1表示经过处理，x2=0表示未经处理；1是未经处理的最终反应速度；1是经处理后最终反应速度的增长值；2是未经处理的反应的半速度；2是经处理后反应的半速度点的增

37、长值；(1) 运行以下程序：clc; clear;x=0.02 0.02 0.06 0.06 0.11 0.11 0.22 0.220.56 0.56 1.10 1.10'y1=76 47 97 107 123 139 159 152 191 201 207 200'%已处理y2= 67 51 84 86 98 115 131 124 144 158 160' %未处理xx=x;x(1:11),ones(12,1);zeros(11,1);%第2列为示性变量值（取0, 1）yy=y1;y2;beta0=170 0.05 60 0.01; %参数初值beta,R,J=nl

38、infit(xx,yy,'fun',beta0);betaci=nlparci(beta,R,J);fprintf('%2s%5s%14sn','参数','估计值','置信区间');%1个汉字算1个字符for i=1:length(beta)if i<=2fprintf('%1d%10.4f %8.4f, %8.4fn',i,beta(i),betaci(i,:);elsefprintf('%1d%10.4f %8.4f, %8.4fn',i-2,beta(i),betaci(

39、i,:);endendsubplot(1,2,1);y=fun(beta,xx);plot(xx(:,1),yy,'o',xx(:,1),y,'+');xlabel('图7 预测图 o表示原始数据，+表示拟合结果');subplot(1,2,2);plot(xx(1:12,1),R(1:12),'+',xx(13:23,1),R(13:23),'o',0,1.4,0,0,'r-');xlabel('图8 残差图 +表示经过处理，o表示未经处理');nlintool(xx,yy,

40、9;fun',beta)function y=fun(beta,x)y=(beta(1)+beta(3)*x(:,2).*x(:,1)./(beta(2)+beta(4)*x(:,2)+x(:,1);(1) 给出运行结果（数据和图形）（比较342，343）：数据结果（比较p342表4）：图形（比较p343图7、图8)：拟合图（见p343图9）：(2) 对上面程序运行中nlintool给出的交互式画面进行操作。单击画面左下方的导出按钮”Export”，出现下面画面，其中只选中Parameters,Parameter CI和RMSE复选框，单击OK，则变量beta1,betaci1和rms

41、e被送到工作空间中，在命令窗口键入这三个变量名显示其值。注：beta1为参数估计值，betaci1为参数的置信区间，rmse为剩余标准差。(2) 给出结果（beta1,betaci1和rmse的值）：(3) 将上面的模型简化为按(1)和(2)完成本题。(3) 给出结果（程序、数据和图形）（比较343，345）：程序：数据结果（比较p343表5）：图形（比较p345图10、图11)：拟合图（见p345图12）：数据结果：4. 投资额与生产总值和物价指数p346352下在给出一组数据，其中：第1列 年份序号t第2列 投资额yt第3列 国民生产总值x1t第4列 物价指数x2t存放在一个名为p346.

42、txt的文件中。1 90.9 596.7 0.71672 97.4 637.7 0.72773 113.5 691.1 0.74364 125.7 756.0 0.76765 122.8 799.0 0.79066 133.3 873.4 0.82547 149.3 944.0 0.86798 144.2 992.7 0.91459 166.4 1077.6 0.960110 195.0 1185.9 1.000011 229.8 1326.4 1.057512 228.7 1434.2 1.150813 206.1 1549.2 1.257914 257.9 1718.0 1.323415

43、324.1 1918.3 1.400516 386.6 2163.9 1.504217 423.0 2417.8 1.634218 401.9 2631.7 1.784219 474.9 2954.7 1.951420 424.5 3073.0 2.06884.1（编程）基本的回归模型p346348(1) 编写一个程序，绘制yt对x1t的散点图和yt对x2t的散点图。(1) 给出程序和运行结果（比较347图1、2）：程序：图形（比较p347图1、图2)：(2) 假设满足多元线性回归模型yt = 0 + 1 x1t + 2 x2t + t编写一个程序求回归系数及其置信区间（置信水平a = 0.0

44、5）、检验统计量R2、F、p、s2。(2) 给出要求的程序和运行结果（比较347表2）：程序：运行结果（比较p347表2）：(3) 用下面程序求参数的估计值和剩余标准差。clc; clear; M=dlmread('p346.txt');yt=M(:,2); x1t=M(:,3); x2t=M(:,4);rstool(x1t,x2t,yt,'linear')点击“Expor”，选择第1和2复选框，OK，变量beta和rmse导出到工作空间，在命令窗口分别键入变量名。(3) 给出结果（beta和rmse值）：4.2（验证）加入自相关后的模型p348352(1) 给

45、出模型yt=0+1x1t+2x2t+t的残差et（见p348表3）和残差etet-1的散点图（见p349图4）。程序如下：clc; clear;M=dlmread('p346.txt');yt=M(:,2); x1t=M(:,3); x2t=M(:,4);b,bint,r,rint,stats=regress(yt,ones(size(M,1),1) x1t x2t);disp('表3 模型yt=0+1 x1t+2 x2t+t的残差et');for k=1:5:16 fprintf('%10d',k:k+4); fprintf('n

46、9;); fprintf('%10.4f',r(k:k+4); fprintf('nn');endplot(r(1:19),r(2:20),'r+',-25,20,0,0,'b-',0,0,-25,20,'b-');xlabel('et-1'); ylabel('et');title('图4 残差etet-1的散点图');(1) 给出程序运行结果（比较348349）：数据结果（见348表3）：图形结果（见349图4）：(2) 加入自相关后的模型求解。模型为yt = 0

47、 + 1 x1t + 2 x2t + t，t = t-1 + ut其中是自相关系数，|1，相互独立且服从均值为零的正态分布，t=1,2,n。作变换yt ( 0 + 1 x1t + 2 x2t )= yt-1 -( 0 + 1 x1t-1 + 2 x2t-1 )+ utyt yt-1 =0(1 ) + 1(x1t x1t-1) + 2(x1t x1t-1)+ utyt*= yt yt-1，xt*= xt xt-1，0*=0(1 )模型化为yt* = 0* + 1 x1t* + 2 x2t* + ut自相关系数的估计值为其中et为模型yt = 0 + 1 x1t + 2 x2t + t的残差。程序

48、如下：clc; clear; format compact;M=dlmread('p346.txt');yt=M(:,2); x1t=M(:,3); x2t=M(:,4);b,bint,r,rint,stats=regress(yt,ones(size(M,1),1) x1t x2t);format short g;DW=(r(2:20)-r(1:19)'*(r(2:20)-r(1:19)/(r(2:20)'*r(2:20)lou=1-DW/2yyt=yt(2:20)-lou*yt(1:19);xx1t=x1t(2:20)-lou*x1t(1:19);xx2t=x2t(2:20)-lou*x2t(1:19);b,bint,r,rint,stats=regress(yyt,