實驗11_統(tǒng)計回歸模型(4學時)

1、數(shù)學建模實驗王平實驗11 統(tǒng)計回歸模型（4學時）（第10章 統(tǒng)計回歸模型）1. 牙膏的銷售量p325332下面給出一組數(shù)據(jù)，其中：第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碼分隔的數(shù)值數(shù)據(jù)文件讀入到矩陣中dlmread：讀取ASCII碼文件的MATLAB函數(shù)M=dlmread('fun

5、.txt');fun.m是一個數(shù)據(jù)文件，存放一個數(shù)據(jù)矩陣，將文件內(nèi)容寫入M。(1) 運行程序并給出結(jié)果（比較327圖1）：(2) 確定y對x1的擬合，繪制散點圖與擬合曲線組合圖形從y對x1的散點圖可以發(fā)現(xiàn)，可用線性模型（直線）來擬合（其中是隨機誤差）。程序如下：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多元線性回歸函數(shù)調(diào)用格式b,bint,r,rint,stats=regress(y,x,alpha)例，多元回歸模型為：輸入：y為n（=30）維列向量數(shù)據(jù)。x為對應于回歸系數(shù) ( 0, 1, 2, 3 )' 的數(shù)據(jù)矩陣 1 x1 x2 x22（30×4矩陣，第1列全1）。alpha為置信水平（缺省時為0.05）。輸出：b為=( 0, 1, 2, 3 )'估計值，4維列向量。bint為b

7、的置信區(qū)間，4×2矩陣。r為殘差n(=30)維列向量y-x。rint為r的置信區(qū)間，30×2矩陣。stats為回歸模型的檢驗統(tǒng)計量，含4個值：R2回歸方程的決定系數(shù)（R是相關(guān)系數(shù)）F統(tǒng)計值P與F統(tǒng)計量對應的概率值s2剩余方差(2) 運行程序并給出結(jié)果（比較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) 運行程序并給出結(jié)果（比較327圖2）：(4) 確定y對

8、x2的的擬合，繪制散點圖與擬合曲線組合圖形從y對x2的散點圖可以發(fā)現(xiàn)，可用二次函數(shù)模型來擬合。程序如下：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、 運行程序并給出結(jié)果（比較327圖2）：(5) y對x1, x2的回歸模型及其求解，銷售量預測綜上得回歸模型變量x1, x2為回歸變量，參數(shù)b0, b1, b2, b3為回歸系數(shù)。程序如下：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','參數(shù)

10、','估計值','置信區(qū)間');%1個漢字算1個字符for i=1:length(b) fprintf ('%1d%9.4f %7.4f, %7.4fn',i-1,b(i,:),bint(i,:);end % %d將i當整數(shù)輸出，%7.4f按實數(shù)格式輸出，區(qū)域?qū)?個字符，4位小數(shù)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輸出到命令窗口或?qū)憯?shù)據(jù)到文本文件見參考資料：MATLAB函數(shù)和命令的用法。(5) 運行程序并給出結(jié)果（比較328表2,329的預測結(jié)果）：1.2（驗證，編程）模型改進p329332仍使用題1的數(shù)據(jù)。(1)（編程）y對x1, x2的回歸模型的改進和求解，銷售量預測改進的模型參考題1(5)的程序，編寫一個類似的程序，運行結(jié)果與教材p329330的表3及相關(guān)結(jié)果相比較。(1) 給出程序和運行結(jié)果（比較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','參數(shù)','估計值','置信區(qū)間');%1個漢字算1個字符for i=1:length(b) fprintf ('%1d%9.4f %7.4f, %7.4fn',i-1,b(i,:),bint(i,:);end % %d將i當整數(shù)輸出，

13、%7.4f按實數(shù)格式輸出，區(qū)域?qū)?個字符，4位小數(shù)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的估計值和預測區(qū)間。用鼠標移動交互式畫面中的十字線，或在圖下方的窗口內(nèi)輸入，可改變x1和x2的數(shù)值。改變x1=0.2，x2=6.5，觀察窗口左邊的y估計值和預測區(qū)間。點擊所得交互畫面左下方的輸出按鈕“Export”，所得畫面（導出到工作空間）第1個復選框是“將擬合參數(shù)存到一個名為beta的MATLAB變量中”，點擊OK。在命令窗口提示符鍵入變量名beta將得到參數(shù)( 0, 1, 2, 3, 4,

15、 5 )' 的值。(2) 運行程序并給出結(jié)果：參數(shù)(0,1,2,3,4,5)' 的值（比較331）：2. 軟件開發(fā)人員的薪金p332338在下面給出的數(shù)據(jù)中：（存入文件p333.txt）第1列 編號第2列 薪金y第3列 資歷x1（從事專業(yè)工作的年數(shù)）第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假設(shè)滿足多元線性回歸模型y = 0 + 1 x1 +

18、2 x2 + 3 x3 + 4 x4 + 2.1（驗證）基本模型p332335求回歸系數(shù)及其置信區(qū)間（置信水平a = 0.05）、檢驗統(tǒng)計量R2、F、p、s2，有關(guān)散點圖的MATLAB程序如下：%10.2 軟件開發(fā)人員的薪金基本模型%模型：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','參數(shù)','估計值','置信區(qū)間');%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的關(guān)系 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的關(guān)系 給出程序的運行結(jié)果（數(shù)據(jù)和圖形）（比較334，335）：數(shù)據(jù)結(jié)果（比較334表2）：圖形結(jié)果（比較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，區(qū)別在于模型不同，所以要根據(jù)新模型修改題2.1的程序，仍后運行。 給出程序和運行結(jié)果（程序、數(shù)據(jù)和圖形）（比較336表4、圖3、圖4）：程序：%10.2 軟件開發(fā)人員的薪金基本模型%模型：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','參數(shù)','估計值','置信區(qū)間');%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的關(guān)系 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的關(guān)系數(shù)據(jù)結(jié)果（比較p336表4）：圖形結(jié)果（比較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中給出數(shù)據(jù)中編號為33的數(shù)據(jù)（異常點）(1) 給出程序和運行結(jié)果（數(shù)據(jù)和圖形）（比較337表5、圖5、圖6）：程序：%10.2 軟件開發(fā)人員的薪金基本模型%模型：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','參數(shù)','估計值','置信區(qū)間');%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的關(guān)系 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的關(guān)系數(shù)據(jù)結(jié)果（比較p337表5）：圖形結(jié)果（比較p337圖5、圖6）：(2)（驗證） 6種管理教育組合人員的“基

27、礎(chǔ)”薪金（資歷為0）。組合管 理 教 育 “基礎(chǔ)”薪金資歷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', '組合', '管理', '教育', '“基礎(chǔ)”薪金');%1個漢字算1個字符for k=1:6fprintf('%3d%6d%6d%1

29、1.0fn',k,x2(k),fix(k+1)/2),y(k);end(2) 運行程序并給出結(jié)果（比較337表6）：3. 酶促反應p338345嘌呤霉素實驗中的反應速度與底物濃度數(shù)據(jù)底物濃度(ppm)0.020.060.110.220.561.10反應速度處理764797107123139159152191201207200未處理6751848698115131124144158160/3.1（編程）線性化模型p338340(1) 根據(jù)給出的數(shù)據(jù)，試編寫一個程序繪制出如下圖所示的散點圖（比較：p339圖1、圖2）(1) 給出程序和運行結(jié)果（比較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經(jīng)處理的散點圖');subplot(1,2,2)plot(x1,y4, '+')

31、;hold on;plot(x2,y3, '+');title('y對x未經(jīng)處理的散點圖');圖形（比較p339圖1、圖2）：(2) 線性化模型Michaelis-Menten模型通過代換或可化為線性模型試編寫一個程序，對經(jīng)過處理的實驗數(shù)據(jù)，求出：(a) 參數(shù)的估計；(b) 1/y與1/x的散點圖和回歸直線；(c) 用線性化得到的原始數(shù)據(jù)擬合圖。程序運行的結(jié)果如下：p340表2p339圖3、p340圖4(2) 給出程序和運行結(jié)果（比較339，340）：程序：數(shù)據(jù)結(jié)果（比較p340表2）：圖形(比較p339圖3、p340圖4)：3.2（驗證）非線性模型及求解p34

32、0341Michaelis-Menten模型提示：求解非線性模型用到的操作(a) 非線性最小二乘擬合命令beta, R, J=nlinfit(x, y, 'model', beta0)輸入：x為自變量數(shù)據(jù)矩陣，每列一個變量；y為因變量數(shù)據(jù)列向量；model為模型的m函數(shù)文件，形式為function y=f(beta,x)，beta為待估計參數(shù)；beta0為給定的參數(shù)初值；輸出：beta為參數(shù)的估計值；R為殘差；J為用于估計預測誤差的Jacobi矩陣。(b) 參數(shù)的置信區(qū)間命令betaci=nlparci(beta,R,J)(c) 非線性模型的GUI工具預測命令nlintool(

33、x,y, 'model',beta)功能：給出一個交互式畫面，拖動畫面中的十字線可改變自變量x的取值，直接得到因變量y的預測值和預測區(qū)間，同時在左下方的Export中，可向MATLAB工作區(qū)傳送統(tǒng)計數(shù)據(jù)。(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、 %用線性化模型得到的作為非線性模型參數(shù)估計的初始迭代值。beta,R,J=nlinfit(x,y,'fun',beta0); %beta是行向量betaci=nlparci(beta,R,J); %betaci是2行2列矩陣fprintf('%2s%6s%14sn','參數(shù)','估計值','置信區(qū)間');%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) %必須作為單獨一個函數(shù)式m文件yhat=beta(1)*x./(beta(2)+x);(1) 給出運行結(jié)果（數(shù)據(jù)和圖形）（比較341）：數(shù)據(jù)結(jié)果（比較p341表3）：（Figure 1，比較p341圖5）：擬合圖（比較p341圖6）：(2) 對上面程序運行中nlintool給出的交互式畫面進行操作。單擊畫面左下方的導出按鈕”Export”，出現(xiàn)下面畫面，其中只選中Pa

36、rameters,Parameter CI和RMSE復選框，單擊OK，則變量beta1,betaci1和rmse被送到工作空間中，在命令窗口鍵入這三個變量名顯示其值。注：beta1為參數(shù)估計值，betaci1為參數(shù)的置信區(qū)間，rmse為剩余標準差。(2) 給出結(jié)果（beta1,betaci1和rmse的值）：3.3（驗證）混合反應模型p342345模型其中：x1為底物濃度；x2為一示性變量（0-1變量），用來表示是否嘌呤霉素處理，令x2=1表示經(jīng)過處理，x2=0表示未經(jīng)處理；1是未經(jīng)處理的最終反應速度；1是經(jīng)處理后最終反應速度的增長值；2是未經(jīng)處理的反應的半速度；2是經(jīng)處理后反應的半速度點的增

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; %參數(shù)初值beta,R,J=nl

38、infit(xx,yy,'fun',beta0);betaci=nlparci(beta,R,J);fprintf('%2s%5s%14sn','參數(shù)','估計值','置信區(qū)間');%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表示原始數(shù)據(jù)，+表示擬合結(jié)果');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 殘差圖 +表示經(jīng)過處理，o表示未經(jīng)處理');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) 給出運行結(jié)果（數(shù)據(jù)和圖形）（比較342，343）：數(shù)據(jù)結(jié)果（比較p342表4）：圖形（比較p343圖7、圖8)：擬合圖（見p343圖9）：(2) 對上面程序運行中nlintool給出的交互式畫面進行操作。單擊畫面左下方的導出按鈕”Export”，出現(xiàn)下面畫面，其中只選中Parameters,Parameter CI和RMSE復選框，單擊OK，則變量beta1,betaci1和rms

41、e被送到工作空間中，在命令窗口鍵入這三個變量名顯示其值。注：beta1為參數(shù)估計值，betaci1為參數(shù)的置信區(qū)間，rmse為剩余標準差。(2) 給出結(jié)果（beta1,betaci1和rmse的值）：(3) 將上面的模型簡化為按(1)和(2)完成本題。(3) 給出結(jié)果（程序、數(shù)據(jù)和圖形）（比較343，345）：程序：數(shù)據(jù)結(jié)果（比較p343表5）：圖形（比較p345圖10、圖11)：擬合圖（見p345圖12）：數(shù)據(jù)結(jié)果：4. 投資額與生產(chǎn)總值和物價指數(shù)p346352下在給出一組數(shù)據(jù)，其中：第1列 年份序號t第2列 投資額yt第3列 國民生產(chǎn)總值x1t第4列 物價指數(shù)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) 給出程序和運行結(jié)果（比較347圖1、2）：程序：圖形（比較p347圖1、圖2)：(2) 假設(shè)滿足多元線性回歸模型yt = 0 + 1 x1t + 2 x2t + t編寫一個程序求回歸系數(shù)及其置信區(qū)間（置信水平a = 0.0

44、5）、檢驗統(tǒng)計量R2、F、p、s2。(2) 給出要求的程序和運行結(jié)果（比較347表2）：程序：運行結(jié)果（比較p347表2）：(3) 用下面程序求參數(shù)的估計值和剩余標準差。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) 給出結(jié)果（beta和rmse值）：4.2（驗證）加入自相關(guān)后的模型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) 給出程序運行結(jié)果（比較348349）：數(shù)據(jù)結(jié)果（見348表3）：圖形結(jié)果（見349圖4）：(2) 加入自相關(guān)后的模型求解。模型為yt = 0

47、 + 1 x1t + 2 x2t + t，t = t-1 + ut其中是自相關(guān)系數(shù)，|1，相互獨立且服從均值為零的正態(tài)分布，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自相關(guān)系數(shù)的估計值為其中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,