MATLAB數(shù)據(jù)分析方法

回歸分析是最常用的數(shù)據(jù)分析方法之一。它是根據(jù)已得的試驗結(jié)果以及以往的經(jīng)驗來建立統(tǒng)計模型，并研究變量間的相關(guān)關(guān)系，建立起變量之間關(guān)系的近似表達式即經(jīng)驗公式，并由此對相應(yīng)的變量進行預(yù)測和控制等.3.1一元回歸模型

3.1.1一元線性回歸模型1.一元線性回歸的基本概念通常，我們對總體(x,Y)進行n次的獨立觀測，獲得n組數(shù)據(jù)（稱為樣本觀測值）(x1,y1),(x2,y2),…,(xn,yn)利用最小二乘法可以得到回歸模型參數(shù)

0，1的最小二乘估計

設(shè)Y是一個可觀測的隨機變量，它受到一個非隨機變量因素x和隨機誤差

的影響。若Y與x有如下線性關(guān)系：(3.1.1)且E=0,D=2,則稱(3.1.1)為一元線性回歸模型.其中0，1為回歸系數(shù)，x為自變量，Y為因變量.

(3.1.2)

其中于是建立經(jīng)驗公式模型：(3.1.3)一元線性回歸分析的主要任務(wù)：一是利用樣本觀測值對回歸系數(shù)

0，

1和

作點估計；二是對方程的線性關(guān)系即

1作顯著性檢驗；三是在x=x0處對Y作預(yù)測等.以下舉例說明建立經(jīng)驗公式（3.1.3）的方法。例3.1.1近10年來，某市社會商品零售總額與職工工資總額（單位：億元）數(shù)據(jù)如下表3.1。表3.1商品零售總額與職工工資表（單位：億元）建立社會商品零售總額與職工工資總額數(shù)據(jù)的回歸模型工資總額23.827.631.632.433.734.943.252.863.873.4零售總額41.451.861.767.968.777.595.9137.4155.0175.0解：%首先輸入數(shù)據(jù)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(職工工資總額)')%橫坐標(biāo)名ylabel('y(商品零售總額)')%縱坐標(biāo)名圖3.1商品零售總額與職工工資總額數(shù)據(jù)散點圖%計算最佳參數(shù)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：當(dāng)x=0，得到y(tǒng)=-23.5493億元如何理解？問題2：如何檢驗E=0?D=2?2.一元多項式回歸模型在一元回歸模型中，如果變量y與x的關(guān)系是n次多項式，即其中

是隨機誤差，服從正態(tài)分布N(0,

2)a0,a1,…,an為回歸系數(shù)，則稱(3.1.4)為多項式回歸模型.

(3.1.4)(1)多項式曲線擬合在MATLAB7的統(tǒng)計工具箱中，有多項式曲線擬合的命令polyfit，其調(diào)用格式有以下三種：p=polyfit(x,y,n)[p,S]=polyfit(x,y,n)[p,S,mu]=polyfit(x,y,n)其中，輸入x,y分別為自變量與因變量的樣本觀測數(shù)據(jù)向量；n是多項式的階數(shù)，對于一元線性回歸則取n=1;輸出p是按照降冪排列的多項式的系數(shù)向量，S是一個矩陣，用于估計預(yù)測誤差或供MATLAB的其它函數(shù)的調(diào)用。例3.1.2某種合金中的主要成分為A,B兩種金屬，經(jīng)過試驗發(fā)現(xiàn)：這兩種金屬成分之和x與合金的膨脹系數(shù)y有如下關(guān)系，建立描述這種關(guān)系的數(shù)學(xué)表達式.表3.2合金的膨脹系數(shù)表解：%首先輸入數(shù)據(jù)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(兩種合金之和)')%橫坐標(biāo)名ylabel(‘y(合金膨脹系數(shù))’)%縱坐標(biāo)名

%然后根據(jù)散點圖猜測曲線類別（2.1.7）

x3737.53838.53939.54040.54141.54242.543y3.4332.272.11.831.531.71.81.92.352.542.9由于散點圖呈拋物線，故選擇二次函數(shù)曲線進行擬合.p=polyfit(x,y,2)%注意取n=2運行得到回歸系數(shù)：p=0.1660-13.3866271.6231

即二次回歸模型為：多項式曲線擬合預(yù)測的命令polyval，其調(diào)用格式有以下兩種：Y=polyval(p,x0)[Y,Delta]=polyconf(p,x0,S,alpha)其中，輸入p,S是由多項式擬合命[p,S]=polyfit(x,y,n)的輸出，x0是要預(yù)測的自變量的值.輸出Y是polyfit所得的回歸多項式在x處的預(yù)測值。(2)多項式回歸的預(yù)測與置信區(qū)間如果輸入數(shù)據(jù)的誤差相互獨立，且方差為常數(shù)，則Y±Delta至少包含95%的預(yù)測值；alpha缺省時為0.05。(Y-Delta,Y+Delta)即95%的置信區(qū)間

(3)多項式回歸的GUI界面命令多項式回歸的GUI界面命令polytool，其典型調(diào)用格式

polytool(x,y,n,alpha)其中，輸入x,y分別為自變量與因變量的樣本觀測數(shù)據(jù)向量；n是多項式的階數(shù)；置信度為(1-alpha)%，alpha缺省時為0.05。該命令可以繪出總體擬合圖形以及(1-alpha)上、下置信區(qū)間的直線（屏幕上顯示為紅色）.此外，用鼠標(biāo)拖動圖中縱向虛線，就可以顯示出對于不同的自變量數(shù)值所對應(yīng)的預(yù)測狀況，與此同時圖形左端數(shù)值框中會隨著自變量的變化而得到的預(yù)報數(shù)值以及(1-alpha)置信區(qū)間長度一半的數(shù)值。例3.1.3為了分析X射線的殺菌作用，用200千伏的X射線來照射細菌，每次照射6分鐘用平板計數(shù)法估計尚存活的細菌數(shù)，照射次數(shù)記為t，照射后的細菌數(shù)y如表3.3所示。t123456789101112131415y3522111971601421061046056383632211915表3.3X射線照射次數(shù)與殘留細菌數(shù)試求：①給出y與t的二次函數(shù)回歸模型；②在同一坐標(biāo)系內(nèi)做出原始數(shù)據(jù)與擬合結(jié)果的散點圖③預(yù)測t=16時殘留的細菌數(shù)；④根據(jù)問題實際意義選擇多項式函數(shù)是否合適？數(shù)據(jù)來源：http///~hadi/RABE解：%輸入原始數(shù)據(jù)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('原始數(shù)據(jù)','二次函數(shù)')xlabel('t(照射次數(shù))')ylabel('y(殘留細菌數(shù))')t0=16;yc1=polyconf(p,t0)%預(yù)測t0=16時殘留的細菌數(shù)運行結(jié)果為p=1.9897-51.1394347.8967，yc1=39.0396即二次回歸模型為yc1=39.0396,表明照射16次后，用二次函數(shù)計算出細菌殘留數(shù)為39.0396，顯然與實際不相符合。調(diào)用多項式回歸的GUI界面命令polytool,如圖3.4原始數(shù)據(jù)與擬合結(jié)果的散點圖如圖3.3所示，從圖形可知擬合效果較好.圖3.3原始數(shù)據(jù)與擬合結(jié)果的散點圖根據(jù)實際問題的意義可知：盡管二次多項式擬合效果較好，但是用于預(yù)測并不理想。因此如何根據(jù)原始數(shù)據(jù)散點圖的規(guī)律，選擇適當(dāng)?shù)幕貧w曲線是非常重要的，因此有必要研究非線性回歸分析.

圖3.4二次函數(shù)預(yù)測交互圖3.1.2一元非線性回歸模型

為了便于正確地選擇合適的函數(shù)進行回歸分析建模，我們給出通常選擇的六類曲線如下所示：1.非線性曲線選擇（1）雙曲線1/y=a+b/x(見圖3.5)。圖3.5雙曲線圖3.5雙曲線(2)冪函數(shù)曲線y=axb,其中x>0,a>0(圖3.6)。圖3.6冪函數(shù)曲線（3）指數(shù)曲線y=aebx，其中參數(shù)a>0(見圖3.7)。圖3.7指數(shù)曲線（4）倒指數(shù)曲線

,其中a>0(圖3.8)。圖3.8倒指數(shù)曲線（5）y=a+blnx(見圖3.9)。圖3.9對數(shù)曲線（6）S型曲線(見圖3.10)。圖3.10S型曲線

對于非線性回歸建模通常有兩種方法：一是通過適當(dāng)?shù)淖儞Q轉(zhuǎn)化為線性回歸模型，例如雙曲線模型(圖3.5)。如果無法實現(xiàn)線性化，可以利用最小二乘法直接建立非線性回歸模型，求解最佳參數(shù)。2.非線性回歸的MATLAB命令MATLAB統(tǒng)計工具箱中實現(xiàn)非線性回歸的命令有nlinfit、nlparci、lpredci和nlintool。下面逐一介紹調(diào)用格式。非線性擬合命令nlinfit，調(diào)用格式：[beta,r,J]=nlinfit(x,y,'model',beta0)其中，輸人數(shù)據(jù)x，y分別為n×m矩陣和n維列向量，對一元非線性回歸，x為n維列向量，model是事先用M文件定義的非線性函數(shù)，beta0是回歸系數(shù)的初值(需要通過解方程組得到)，beta是估計出的最佳回歸系數(shù)，r是殘差，J是Jacobian矩陣，它們是估計預(yù)測誤差需要的數(shù)據(jù)。非線性回歸預(yù)測命令nlpredci，調(diào)用格式：

ypred=nlpredci(FUN,inputs,beta,r,J)其中，輸入?yún)?shù)beta,r,J是非線性回歸命令nlinfit的輸出結(jié)果,FUN是擬合函數(shù)，inputs是需要預(yù)測的自變量；輸出量ypred是inputs的預(yù)測值。非線性回歸置信區(qū)間命令nlparci，調(diào)用格式：ci=nlparci(beta,r,J,alpha)輸入?yún)?shù)beta,r,J就是非線性回歸命令nlinfit輸出的結(jié)果，輸出ci是一個矩陣，每一行分別為每個參數(shù)的(1-alpha)%的置信區(qū)間，alpha缺省時默認(rèn)為0.05.非線性回歸的GUI界面命令nlintool，典型調(diào)用格式nlintool(x,y,fun,beta0)其中參數(shù)x,y,fun,beta0與命令nlinfit中的參數(shù)含義相同.例3.1.4.在M文件中建立函數(shù)y=a(1-be-cx)，其中a,b,c為待定的參數(shù)。解：fun=inline('b(1)*(1-b(2)*exp(-b(3)*x))','b','x');此處，將b看成參變量，b(1),b(2),b(3)為其分量.例3.1.5煉鋼廠出鋼時所用盛鋼水的鋼包，由于鋼水對耐火材料的侵蝕，容積不斷增大，我們希望找出使用次數(shù)與增大容積之間的函數(shù)關(guān)系.實驗數(shù)據(jù)如表3.4。使用次數(shù)(x)23456789增大容積(y)6.428.29.589.59.7109.939.99使用次數(shù)(x)10111213141516增大容積(y)10.4910.5910.610.810.610.910.76表3.4鋼包使用次數(shù)與增大容積（1）建立非線性回歸模型1/y=a+b/x；（2）預(yù)測鋼包使用x0=17次后增大的容積y0；（3）計算回歸模型參數(shù)的95%的置信區(qū)間。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];%初始參數(shù)值fun=inline('x./(b(1)*x+b(2))','b','x');[beta,r,J]=nlinfit(x,y,fun,b0);beta%輸出最佳參數(shù)y1=x./(0.0845*x+0.1152);%擬合曲線plot(x,y,'*',x,y1,'-or')legend('原始數(shù)據(jù)','擬合曲線')注意：初始值要先計算后，才能得到上面程序中的b0，由于確定兩個參數(shù)值，因此我們選擇已知數(shù)據(jù)中的兩點（2,6.42）和（16,10.76）代入設(shè)定方程，得到方程組上述方程組有兩種解法：手工方法與Matlab方法。下面用Matlab方法解方程組：[a,b]=solve('6.42*(2*a+b)=2','10.76*(16*a+b)=16')輸出為a=.83961597702347450462657355615004e-1b=.14360328434608391527406223581049圖3.11鋼包使用次數(shù)與增大容積的非線性擬合圖在例3.1.5中，預(yù)測鋼包使用17次后增大的容積，可在執(zhí)行上面的程序中，繼續(xù)輸入命令ypred=nlpredci(fun,17,beta,r,J)得到：ypred=10.9599即鋼包使用17次后增大的容積10.9599。求回歸模型參數(shù)的95%的置信區(qū)間，只要繼續(xù)添加程序ci=nlparci(beta,r,J)運行后得到ci=0.08140.08760.09340.1370即回歸模型中參數(shù)a,b的95%的置信區(qū)間分別為(0.0814,0.0876)與（0.0934，0.1370）.我們求出的最佳參數(shù)分別為a=0.0845,b=0.1152均屬于上述置信區(qū)間。圖3.12給出鋼包使用次數(shù)與增大容積的非線性擬合的交互圖形，圖中的的圓圈是實驗的原始數(shù)據(jù)點，兩條虛線為95%上、下置信區(qū)間的曲線（屏幕上顯示為紅色），中間的實線（屏幕上顯示為綠色）是回歸模型曲線，縱向的藍色虛線顯示了在自變量為8.9502，橫向的虛線給出了對應(yīng)的預(yù)測值為10.2734.圖3.12鋼包使用次數(shù)與增大容積的非線性擬合交互圖例3.1.6對例題3.1.3進行非線性回歸，并預(yù)測照射16次后細菌殘留數(shù)目，給出模型參數(shù)的95%的置信區(qū)間，繪出模型交互圖形.解：我們選取函數(shù)y=aebt進行非線性回歸，該方程的兩個參數(shù)具有簡單的物理解釋，a表示實驗開始時的細菌數(shù)目，b表示細菌死亡（或衰變）的速率。MATLAB腳本程序如下：t=1:15;y=[3522111971601421061046056383632211915];fun=inline('b(1)*exp(b(2)*t)','b','t')%非線性函數(shù)beta0=[148,-0.2];%參數(shù)初始值[beta,r,J]=nlinfit(t,y,fun,beta0);%非線性擬合beta%輸出最佳參數(shù)y1=nlpredci(fun,t,beta,r,J);%模型數(shù)值計算plot(t,y,'*',t,y1,'-or'),legend('原始數(shù)據(jù)','非線性回歸')xlabel('t(照射次數(shù))')ylabel('y(殘留細菌數(shù))')ypred=nlpredci(fun,16,beta,r,J)%預(yù)測殘留細菌數(shù)ci=nlparci(beta,r,J)%參數(shù)95%區(qū)間估計nlintool(t,y,fun,beta0)%作出交互圖形運行后結(jié)果如下：beta=400.0904-0.2240即，最佳參數(shù)為：a=400.0904，b=-0.2240故非線性回歸模型為

預(yù)測為：ypred=11.1014即，照射16次后細菌殘留數(shù)目為11.1014，該預(yù)測符合實際，顯然比例3.1.3中多項式回歸的結(jié)果合理。ci=355.2481444.9326-0.2561-0.1919即參數(shù)a置信度為95%的置信區(qū)間(ci的第一行)為：[355.2481,444.9326]參數(shù)b的置信度為95%的置信區(qū)間(ci的第二行)為[-0.2561-0.1919]顯然，最佳參數(shù)a=400.0904，b=-0.2240，均屬于各自置信度為95%的置信區(qū)間。圖3.13原始數(shù)據(jù)與非線性回歸圖形圖3.14原始數(shù)據(jù)與非線性回歸GUI圖形從交互圖形3.14可以看出：圓圈為原始數(shù)據(jù)，兩條虛線（屏幕上顯示紅色）是置信區(qū)間曲線；兩條虛線內(nèi)的實線（屏幕上顯示綠色）是回歸模型曲線；縱向虛線指示照射8次，此時對應(yīng)的水平虛線表示模型得到的殘留細菌數(shù)為：66.6451。圖3.14原始數(shù)據(jù)與非線性回歸GUI圖形3.1.3一元回歸建模實例例3.1.7在四川白鵝的生產(chǎn)性能研究中，得到如下一組關(guān)于雛鵝重（g）與70日齡重(g)的數(shù)據(jù)，試建立70日齡重(y)與雛鵝重(x)的直線回歸方程，計算模型誤差平方和以及可決系數(shù)，當(dāng)雛鵝重分別為：85,95

,115時預(yù)測其70日齡重，以及置信區(qū)間。

表3.5四川白鵝重與70日齡重測定結(jié)果（單位：g）編號123456789101112雛鵝重(x)80869890120102958311310511010070日齡重(Y)235024002720250031502680263024003080292029602860解：（1）作散點圖。以雛鵝重（x）為橫坐標(biāo)，70日齡重（y）為縱坐標(biāo)作散點圖，如圖2-14。在MATLAB命令窗口中輸入:x=[808698901201029583113105110100]';%雛鵝重y=[235024002720250031502680263024003080292029602860]';

%70日齡重plot(x,y,'*')%作散點圖xlabel('x(雛鵝重)')%橫坐標(biāo)名ylabel('y(70日齡重)')%縱坐標(biāo)名圖3.15四川白鵝的雛鵝重與70日齡重散點圖和回歸直線圖由圖形3.15可見白鵝的70日齡重與雛鵝重間存在直線關(guān)系，且70日齡重隨雛鵝重的增大而增大。因此，可認(rèn)為y與x符合一元線性回歸模型。（2）建立直線回歸方程。在MATLAB中調(diào)用命令polyfit，從而求出參數(shù)

0，1的最小二乘估計.在MATLAB命令窗口中繼續(xù)輸入:n=size(x,1)%計算樣本容量[p,s]=polyfit(x,y,1);%調(diào)用命令polyfit計算回歸參數(shù)y1=polyval(p,x);%計算回歸模型的函數(shù)值holdonplot(x,y1)%作回歸方程的圖形，結(jié)果如圖3.15p%顯示參數(shù)的最小二乘估計結(jié)果p=582.185021.7122即參數(shù)的最小二乘估計為所以70日齡重(y)與雛鵝重(x)的直線回歸經(jīng)驗方程為（3）誤差估計與決定系數(shù)。在MATLAB命令窗口中繼續(xù)輸入:TSS=sum((y-mean(y)).^2)%計算總離差平方和RSS=sum((y1-mean(y)).^2)%計算回歸平方和ESS=sum((y-y1).^2)%計算殘差平方和R2=RSS/TSS;%計算樣本決定系數(shù)R2.輸出:TSS=8.314917e+005RSS=7.943396e+005ESS=3.715217e+004R2=0.9553TSS=8.314917e+005RSS=7.943396e+005ESS=3.715217e+004R2=0.9553由于樣本決定系數(shù)R2=0.9553接近于1，因此模型的擬合的效果較好。（4）回歸方程關(guān)系顯著性的F檢驗。在MATLAB命令窗口中繼續(xù)輸入:F=(n-2)*RSS/ESS%計算的F統(tǒng)計量F1=finv(0.95,1,n-2)%查F統(tǒng)計量0.05的分位數(shù)F2=finv(0.99,1,n-2)%查F統(tǒng)計量0.01的分位數(shù)輸出結(jié)果:F=2.138e+002，F(xiàn)1=4.9646，F(xiàn)2=10.0442為了方便，將以上的計算結(jié)果列成表3.6。表3.6

四川白鵝70日齡重與雛鵝重回歸關(guān)系方差分析表自由度（df）平方和(SS)均方和（MS）F值F0.05F0.01回歸1794339.60794339.60213.81**4.9610.04殘差1037152.073715.21總離差11831491.67因為

表明四川白鵝70日齡重與雛鵝重間存在顯著的線性關(guān)系。（5）回歸關(guān)系顯著性的t檢驗。在MATLAB命令窗口中繼續(xù)輸入:T=p(2)/sqrt(ESS/(n-2))*sqrt(sum((x-mean(x)).^2))%計算T統(tǒng)計量T1=tinv(0.975,n-2)%t統(tǒng)計量0.05的分位數(shù)T2=tinv(0.995,n-2)%t統(tǒng)計量0.01的分位數(shù)輸出:T=14.622,T1=2.228,T2=3.169因為T=14.62>t0.01(10)，否定H0，接受H1即四川白鵝70日齡重(y)與雛鵝重(x)的線性回歸系數(shù)是顯著的,可用所建立的回歸方程進行預(yù)測和控制。(6)預(yù)測x1=[85,95,115]';%輸入自變量yc=polyval(p,x1)%計算預(yù)測值[Y,Delta]=polyconf(p,x1,s);I1=[Y-Delta,Y+Delta]%置信區(qū)間輸出:yc=2427.722644.843079.08I1=2279.472575.962503.012786.672927.553230.62所以當(dāng)雛鵝重分別為85,95,115時，白鵝70日齡重分別為2427.72,2644.84,3079.08;且95%的置信區(qū)間分別為：[2279.47,2575.96]，[2503.01,2786.67],[2927.55,3230.62].在程序中加入：polytool(x,y)%交互功能bar(x,y-y1),%殘差圖legend('殘差')h=lillietest(y-y1)%殘差正態(tài)性檢驗輸出h=0得到交互圖形如圖3.16所示，可以看出當(dāng)雛鵝重為100時，模型給出70日齡鵝重為2753.4016.圖3.16四川白鵝70日齡重與雛鵝重線性模型交互圖3.2多元線性回歸模型3.2.1多元線性回歸模型及其表示對于總體的n組觀測值它應(yīng)滿足式（3.2.1），即其中

i(i=1,2,…,n)相互獨立，且設(shè)記,,,則模型(3.2.2)可用矩陣形式表示為

Y=X+(3.2.3)其中Y稱為觀測向量，X稱為設(shè)計矩陣，

稱為待估計向量，

是不可觀測的n維隨機向量，它的分量相互獨立，假定.2.多元線性回歸建模的基本步驟(1)對問題進行直觀分析，選擇因變量與解釋變量，作出與因變量與各解釋變量的散點圖，初步設(shè)定多元線性回歸模型的參數(shù)個數(shù)；(2)輸入因變量與自變量的觀測數(shù)據(jù)(y,X)調(diào)用命令

[b,bint,r,rint,s]=regress(y,X,alpha)，計算參數(shù)的估計。(3)調(diào)用命令rcoplot(r,rint)，分析數(shù)據(jù)的異常點情況。(4)作顯著性檢驗，若檢驗通過，則用模型作預(yù)測。(5)對模型進一步研究：如殘差的正態(tài)性檢驗，殘差的異方差檢驗，殘差進行自相關(guān)性的檢驗等。3.2.2MATLAB的回歸分析命令在MATLAB7.0的統(tǒng)計工具箱中，與多元回歸模型有關(guān)的命令有多個，下面逐一介紹。1.多元回歸建模命令regeress，其調(diào)用格式有以下三種：（1）b=regress(y,X)（2）[b,bint,r,rint,stats]=regress(Y,X)（3）[b,bint,r,rint,stats]=regress(Y,X,alpha)三種方式的主要區(qū)別在輸出項參數(shù)多少上，第3種方式可稱為全參數(shù)方式。以第3種為例來說明regeress命令的輸入與輸出參數(shù)的含義。輸入?yún)?shù)：輸入量Y表示模型（3.1.1）中因變量的觀測向量；X是一個的矩陣，其中第一列元全部是數(shù)“1”，第j列是自變量Xj的觀測向量，即對一元線性回歸，取p=1即可；alpha為顯著性水平輸出參數(shù)：輸出向量b為回歸系數(shù)估計值，bint為回歸系數(shù)的(1-alpha)置信區(qū)間；輸出向量r表示殘差列向量輸出rint為模型的殘差的(1-

)的置信區(qū)間；輸出stats是用于檢驗回歸模型的統(tǒng)計量，有4個分量值：第一個是R2，其中R是相關(guān)系數(shù)，第二個是F統(tǒng)計量值，第三個是與統(tǒng)計量F對應(yīng)的概率P，當(dāng)P<

時拒絕H0，即認(rèn)為線性回歸模型有意義，第四個是方差

2的無偏估計.例3.2.1某銷售公司將庫存占用資金情況、廣告投入的費用、員工薪酬以及銷售額等方面的數(shù)據(jù)作了匯總,該公司試圖根據(jù)這些數(shù)據(jù)找到銷售額與其他變量之間的關(guān)系，以便進行銷售額預(yù)測并為工作決策提供參考依據(jù)。(1)建立銷售額的回歸模型；(2)如果未來某月庫存資金額為150萬元，廣告投入預(yù)算為45萬元，員工薪酬總額為27萬元，試根據(jù)建立的回歸模型預(yù)測該月的銷售額。

表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解：為了確定銷售額與庫存占用資金、廣告投入、員工薪酬之間的關(guān)系，分別作出y與x1,x2,x3的散點圖，若散點圖顯示它們之間近似線性關(guān)系，則可設(shè)定y與x1,x2,x3的關(guān)系為三元線性回歸模型%輸入數(shù)據(jù)并作散點圖（圖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與庫存資金、廣告投入、員工薪酬具有線性關(guān)系，因此可以建立三元線性回歸模型.圖3.18銷售額與庫存、廣告、薪酬散點圖%調(diào)用命令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,%輸出結(jié)果程序運行結(jié)果b=162.06327.273913.9575-4.3996bint=-580.3603904.48674.373410.17437.164920.7501-46.779637.9805

stats=0.9574804050105.08665208910.000000000810077.9867891125輸出結(jié)果說明：b就是模型中的參數(shù)

0,1,2

，因此回歸模型為b就是模型中的參數(shù)

0,1,2

，因此回歸模型為bint的各行分別為參數(shù)

0,1,2的95%的置信區(qū)間。stats的第一列為模型可決系數(shù)，第二列為F統(tǒng)計量的觀測值，第三列得到概率p，最后一列為模型的殘差平方和2.多元回歸輔助圖形命令（1）殘差圖命令rcoplot，其調(diào)用格式rcoplot(r,rint)其中，輸入?yún)?shù)r,rint是多元回歸建模命令regress輸出的結(jié)果，運行該命令后展示了殘差與置信區(qū)間的圖形。該命令有助于對建立的模型進行分析，如果圖形中出現(xiàn)紅色的點，則可以認(rèn)作異常點，此時可刪除異常點，重新建模，最終得到改進的回歸模型。在上面的程序中加入

rcoplot(r,rint)得到如下圖形圖3.19殘差與置信區(qū)間圖

從圖形中可以看到第五個點為異常點，實際上從表3.7可以發(fā)現(xiàn)第5個月庫存占用資金、廣告投入、員工薪酬均比3月份少，為何銷售額反而增加？這就可以促使該公司的經(jīng)理找出原因，尋找對策。下面的例題介紹如何刪除異常點，對模型進行改進的方法。例3.2.2葛洲壩機組發(fā)電耗水率的主要影響因素為庫水位,出庫流量。數(shù)據(jù)如表3.8所示，利用多元線性回歸分析方法建立耗水率與出庫流量、庫水位的模型。表3.8某天耗水率與出庫流量、庫水位的數(shù)據(jù)

時間年-月-天-時

庫水位（米)出庫流量（立方米）

機組發(fā)電耗水率(立方米/萬千瓦)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解：%輸入原始數(shù)據(jù)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所示，從圖中可以看到無論是庫水位還是出庫流量都與機組發(fā)電耗水率具有線性關(guān)系，因此，可以建立機組發(fā)電耗水率與庫水位和出庫流量的二元線性回歸模型。圖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輸出回歸模型的系數(shù)、系數(shù)置信區(qū)間與統(tǒng)計量如表3.9所示

表3.9回歸模型的系數(shù)、系數(shù)置信區(qū)間與統(tǒng)計量回歸系數(shù)回歸系數(shù)估計值回歸系數(shù)置信區(qū)間

0373.8698[340.082,407.6577]

1-4.9759[-5.4642,-4.4875]

20.0007[0.0004,0.0009]R2=0.9863，F(xiàn)=468.4118，p<0.0001，s2=0.0278由此可得模型為：%模型改進

rcoplot(r,rint);得到圖形如圖3.21所示，發(fā)現(xiàn)有一個異常點，下面給出刪除異常點后，重新建模的程序。由此可得模型為：圖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所示，此時沒有異常點，改進回歸模型的系數(shù)、系數(shù)置信區(qū)間與統(tǒng)計量參見表3.10表3.10改進回歸模型的系數(shù)、系數(shù)置信區(qū)間與統(tǒng)計量回歸系數(shù)回歸系數(shù)估計值回歸系數(shù)置信區(qū)間

0328.4616[290.6145,366.3087]

1-4.3594[-4.8880,-3.8308]

20.0010[0.00073,0.0012]R2=0.9931，F(xiàn)=858.5846，p<0.0001，s2=0.0150我們將表3.9與表3.10加以比較，可以發(fā)現(xiàn)：可決系數(shù)從0.9863提高到0.9931，F(xiàn)統(tǒng)計量從468.4118提高到858.5846，由此可知改進后的模型顯著性提高。圖

3-22刪除異常點后殘差示意圖圖3.21殘差示意圖3.2.3多元線性回歸實例例3.2.3現(xiàn)代服務(wù)業(yè)是社會分工不斷深化的產(chǎn)物，隨著經(jīng)濟的發(fā)展，科學(xué)技術(shù)的進步，現(xiàn)代服務(wù)業(yè)的發(fā)展受到多種因素和條件的影響。不僅受到經(jīng)濟總體發(fā)展水平的影響，還受到第二產(chǎn)業(yè)、就業(yè)、投入等因素的影響，從這幾個主要方面出發(fā)，利用江蘇省統(tǒng)計年鑒的有關(guān)數(shù)據(jù)，通過建立多元線性回歸模型對1990-2008年各種因素對現(xiàn)代服務(wù)業(yè)的影響進行回歸分析。假如構(gòu)建如下江蘇省服務(wù)業(yè)增長模型：Y代表江蘇省服務(wù)業(yè)的增