回歸分析MATLAB實(shí)現(xiàn)

1、普通高等院校計(jì)算機(jī)課程規(guī)劃教材MATLAB數(shù)據(jù)分析方法 李柏年 吳禮斌 主編 張孔生 丁 華 參編 回歸分析MATLAB實(shí)現(xiàn)第1頁(yè) 回歸分析是最慣用數(shù)據(jù)分析方法之一。它是依據(jù)已得試驗(yàn)結(jié)果以及以往經(jīng)驗(yàn)來(lái)建立統(tǒng)計(jì)模型，并研究變量間相關(guān)關(guān)系，建立起變量之間關(guān)系近似表示式即經(jīng)驗(yàn)公式，并由此對(duì)對(duì)應(yīng)變量進(jìn)行預(yù)測(cè)和控制等.3.1一元回歸模型 3.1.1一元線性回歸模型1.一元線性回歸基本概念回歸分析MATLAB實(shí)現(xiàn)第2頁(yè)通常，我們對(duì)總體(x,Y)進(jìn)行n次獨(dú)立觀察，取得n組數(shù)據(jù)（稱為樣本觀察值） (x1,y1),(x2,y2),(xn,yn)利用最小二乘法能夠得到回歸模型參數(shù)0，1最小二乘預(yù)計(jì) 設(shè)Y是一個(gè)可觀

2、察隨機(jī)變量，它受到一個(gè)非隨機(jī)變量原因x和隨機(jī)誤差影響。若Y與x有以下線性關(guān)系：(3.1.1)且E=0,D=2,則稱(3.1.1)為一元線性回歸模型.其中0，1為回歸系數(shù)，x為自變量，Y為因變量.回歸分析MATLAB實(shí)現(xiàn)第3頁(yè) (3.1.2) 其中 于是建立經(jīng)驗(yàn)公式模型：(3.1.3)一元線性回歸分析主要任務(wù)：一是利用樣本觀察值對(duì)回歸系數(shù)0，1和作點(diǎn)預(yù)計(jì)；二是對(duì)方程線性關(guān)系即1作顯著性檢驗(yàn)；三是在x=x0處對(duì)Y作預(yù)測(cè)等.以下舉例說(shuō)明建立經(jīng)驗(yàn)公式（3.1.3）方法?；貧w分析MATLAB實(shí)現(xiàn)第4頁(yè)例3.1.1 近來(lái)，某市社會(huì)商品零售總額與職員工資總額（單位：億元）數(shù)據(jù)以下表3.1。表3.1 商品零售

3、總額與職員工資表 （單位：億元）建立社會(huì)商品零售總額與職員工資總額數(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;回歸分析MATLAB實(shí)現(xiàn)第5頁(yè) % 然后作散點(diǎn)圖plot(x,y,*) %

4、作散點(diǎn)圖xlabel(x(職員工資總額) %橫坐標(biāo)名ylabel(y(商品零售總額) %縱坐標(biāo)名圖3.1商品零售總額與職員工資總額數(shù)據(jù)散點(diǎn)圖回歸分析MATLAB實(shí)現(xiàn)第6頁(yè)% 計(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);運(yùn)行后得到：b1 = 2.7991，b0 = -23.5493所以，回歸模型為 問(wèn)題1：當(dāng)x=0，得到y(tǒng)=-23.5493億元怎樣了解？問(wèn)題2：怎樣檢驗(yàn)E=0? D=2?回歸分析MATLAB實(shí)現(xiàn)第7頁(yè)2. 一元多項(xiàng)式回歸模型在一元回歸模型中，

5、假如變量y與x關(guān)系是n次多項(xiàng)式，即其中是隨機(jī)誤差，服從正態(tài)分布N(0,2)a0,a1,an為回歸系數(shù)，則稱(3.1.4)為多項(xiàng)式回歸模型. (3.1.4)(1)多項(xiàng)式曲線擬合在MATLAB7統(tǒng)計(jì)工具箱中，有多項(xiàng)式曲線擬合命令polyfit，其調(diào)用格式有以下三種：回歸分析MATLAB實(shí)現(xiàn)第8頁(yè)p=polyfit(x,y,n) p,S=polyfit(x,y,n) p,S,mu=polyfit(x,y,n)其中，輸入x,y分別為自變量與因變量樣本觀察數(shù)據(jù)向量；n是多項(xiàng)式階數(shù)，對(duì)于一元線性回歸則取n=1;輸出p是按照降冪排列多項(xiàng)式系數(shù)向量，S是一個(gè)矩陣，用于預(yù)計(jì)預(yù)測(cè)誤差或供MATLAB其它函數(shù)調(diào)用

6、。例3.1.2 某種合金中主要成份為A,B兩種金屬，經(jīng)過(guò)試驗(yàn)發(fā)覺(jué)：這兩種金屬成份之和x與合金膨脹系數(shù)y有以下關(guān)系，建立描述這種關(guān)系數(shù)學(xué)表示式.回歸分析MATLAB實(shí)現(xiàn)第9頁(yè)表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; %其次做散點(diǎn)圖plot(x,y,*)xlabel(x(兩種合金之和) %橫坐標(biāo)名ylabel(y(合金膨脹系數(shù)) %縱坐標(biāo)名 %然后依據(jù)散點(diǎn)圖猜測(cè)曲線類別（2.1.7） x3737.53838.53939.54040.54141.54242.543

7、y3.4332.272.11.831.52.352.542.9回歸分析MATLAB實(shí)現(xiàn)第10頁(yè)因?yàn)樯Ⅻc(diǎn)圖呈拋物線，故選擇二次函數(shù)曲線進(jìn)行擬合.p = polyfit(x,y,2) %注意取n=2運(yùn)行得到回歸系數(shù)：p=0.1660 -13.3866 271.6231即二次回歸模型為：回歸分析MATLAB實(shí)現(xiàn)第11頁(yè)多項(xiàng)式曲線擬合預(yù)測(cè)命令polyval，其調(diào)用格式有以下兩種：Y=polyval(p,x0)Y,Delta=polyconf(p,x0,S,alpha)其中，輸入p,S是由多項(xiàng)式擬合命p,S=polyfit(x,y,n)輸出，x0是要預(yù)測(cè)自變量值.輸出Y是polyfi

8、t所得回歸多項(xiàng)式在x處預(yù)測(cè)值。 (2) 多項(xiàng)式回歸預(yù)測(cè)與置信區(qū)間假如輸入數(shù)據(jù)誤差相互獨(dú)立，且方差為常數(shù)，則YDelta最少包含95%預(yù)測(cè)值；alpha缺省時(shí)為0.05。(Y-Delta, Y+Delta)即95%置信區(qū)間回歸分析MATLAB實(shí)現(xiàn)第12頁(yè) (3) 多項(xiàng)式回歸GUI界面命令多項(xiàng)式回歸GUI界面命令polytool，其經(jīng)典調(diào)用格式 polytool(x,y,n,alpha)其中，輸入x,y分別為自變量與因變量樣本觀察數(shù)據(jù)向量；n是多項(xiàng)式階數(shù)；置信度為(1-alpha)%，alpha缺省時(shí)為0.05。 該命令能夠繪出總體擬合圖形以及(1-alpha)上、下置信區(qū)間直線（屏幕上顯示為紅色

9、）.另外，用鼠標(biāo)拖動(dòng)圖中縱向虛線，就能夠顯示出對(duì)于不一樣自變量數(shù)值所對(duì)應(yīng)預(yù)測(cè)情況，與此同時(shí)圖形左端數(shù)值框中會(huì)伴隨自變量改變而得到預(yù)報(bào)數(shù)值以及(1-alpha) 置信區(qū)間長(zhǎng)度二分之一數(shù)值?；貧w分析MATLAB實(shí)現(xiàn)第13頁(yè)回歸分析MATLAB實(shí)現(xiàn)第14頁(yè)例3.1.3為了分析X射線殺菌作用，用200千伏X射線來(lái)照射細(xì)菌，每次照射6分鐘用平板計(jì)數(shù)法預(yù)計(jì)尚存活細(xì)菌數(shù)，照射次數(shù)記為t，照射后細(xì)菌數(shù)y如表3.3所表示。t123456789101112131415y3522111971601421061046056383632211915表3.3 X射線照射次數(shù)與殘留細(xì)菌數(shù)試求： 給出y與t二次函數(shù)回歸模型；

10、 在同一坐標(biāo)系內(nèi)做出原始數(shù)據(jù)與擬合結(jié)果散點(diǎn)圖 預(yù)測(cè)t=16時(shí)殘留細(xì)菌數(shù)； 依據(jù)問(wèn)題實(shí)際意義選擇多項(xiàng)式函數(shù)是否適當(dāng)？數(shù)據(jù)起源：http/hadi/RABE回歸分析MATLAB實(shí)現(xiàn)第15頁(yè)解：% 輸入原始數(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); % 作二次多項(xiàng)式回歸y1= polyval(p,t); % 模型預(yù)計(jì)與作圖plot(t,y,-*,t,y1,-o); legend(原始數(shù)據(jù),二次函數(shù)) xlabel(t(照射次數(shù)) ylabel(y(殘留細(xì)菌數(shù))t0=16;yc1= p

11、olyconf(p,t0) % 預(yù)測(cè)t0=16時(shí)殘留細(xì)菌數(shù)運(yùn)行結(jié)果為p =1.9897 -51.1394 347.8967，yc1 =39.0396即二次回歸模型為回歸分析MATLAB實(shí)現(xiàn)第16頁(yè)yc1 =39.0396,表明照射16次后，用二次函數(shù)計(jì)算出細(xì)菌殘留數(shù)為39.0396，顯然與實(shí)際不相符合。調(diào)用多項(xiàng)式回歸GUI界面命令polytool,如圖3.4原始數(shù)據(jù)與擬合結(jié)果散點(diǎn)圖如圖3.3所表示，從圖形可知擬合效果很好.圖 3.3 原始數(shù)據(jù)與擬合結(jié)果散點(diǎn)圖回歸分析MATLAB實(shí)現(xiàn)第17頁(yè)依據(jù)實(shí)際問(wèn)題意義可知：盡管二次多項(xiàng)式擬合效果很好，不過(guò)用于預(yù)測(cè)并不理想。所以怎樣依據(jù)原始數(shù)據(jù)散點(diǎn)圖規(guī)律，選

12、擇適當(dāng)回歸曲線是非常主要，所以有必要研究非線性回歸分析. 圖 3.4 二次函數(shù)預(yù)測(cè)交互圖回歸分析MATLAB實(shí)現(xiàn)第18頁(yè) 3.1.2一元非線性回歸模型 為了便于正確地選擇適當(dāng)函數(shù)進(jìn)行回歸分析建模，我們給出通常選擇六類曲線以下所表示：1. 非線性曲線選擇（1）雙曲線1/y=a+b/x(見(jiàn)圖3.5)。圖3.5雙曲線圖3.5雙曲線回歸分析MATLAB實(shí)現(xiàn)第19頁(yè)(2) 冪函數(shù)曲線y=axb, 其中x0,a0(圖3.6)。圖3.6 冪函數(shù)曲線（3）指數(shù)曲線y=aebx，其中參數(shù)a0(見(jiàn)圖3.7)。圖3.7 指數(shù)曲線回歸分析MATLAB實(shí)現(xiàn)第20頁(yè)（4）倒指數(shù)曲線 ,其中a0(圖3.8)。圖3.8 倒指

13、數(shù)曲線（5）y=a+blnx (見(jiàn)圖3.9)。圖3.9 對(duì)數(shù)曲線回歸分析MATLAB實(shí)現(xiàn)第21頁(yè)（6）S型曲線 (見(jiàn)圖3.10)。圖3.10 S型曲線 對(duì)于非線性回歸建模通常有兩種方法：一是經(jīng)過(guò)適當(dāng)變換轉(zhuǎn)化為線性回歸模型，比如雙曲線模型(圖3.5)。假如無(wú)法實(shí)現(xiàn)線性化，能夠利用最小二乘法直接建立非線性回歸模型，求解最正確參數(shù)。回歸分析MATLAB實(shí)現(xiàn)第22頁(yè)2.非線性回歸MATLAB命令MATLAB統(tǒng)計(jì)工具箱中實(shí)現(xiàn)非線性回歸命令有nlinfit、nlparci、lpredci和nlintool。下面逐一介紹調(diào)用格式。 非線性擬合命令nlinfit，調(diào)用格式：beta,r,J = nlinfit

14、(x,y,model,beta0)其中，輸人數(shù)據(jù)x，y分別為nm矩陣和n維列向量，對(duì)一元非線性回歸，x為n維列向量，model是事先用M文件定義非線性函數(shù)，beta0是回歸系數(shù)初值(需要經(jīng)過(guò)解方程組得到)，beta是預(yù)計(jì)出最正確回歸系數(shù)，r是殘差，J是Jacobian矩陣，它們是預(yù)計(jì)預(yù)測(cè)誤差需要數(shù)據(jù)。回歸分析MATLAB實(shí)現(xiàn)第23頁(yè)非線性回歸預(yù)測(cè)命令nlpredci，調(diào)用格式： ypred = nlpredci(FUN,inputs,beta,r,J)其中，輸入?yún)?shù)beta,r,J是非線性回歸命令nlinfit輸出結(jié)果, FUN 是擬合函數(shù)，inputs是需要預(yù)測(cè)自變量；輸出量ypred是in

15、puts預(yù)測(cè)值。非線性回歸置信區(qū)間命令nlparci，調(diào)用格式： ci = nlparci(beta,r,J,alpha)輸入?yún)?shù)beta,r,J就是非線性回歸命令nlinfit輸出結(jié)果，輸出ci是一個(gè)矩陣，每一行分別為每個(gè)參數(shù)(1-alpha)% 置信區(qū)間，alpha缺省時(shí)默認(rèn)為0.05.回歸分析MATLAB實(shí)現(xiàn)第24頁(yè)非線性回歸GUI界面命令nlintool，經(jīng)典調(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=in

16、line(b(1)*(1-b(2)*exp(-b(3)*x),b,x);此處，將b看成參變量，b(1),b(2),b(3)為其分量.例3.1.5 煉鋼廠出鋼時(shí)所用盛鋼水鋼包，因?yàn)殇撍畬?duì)耐火材料侵蝕，容積不停增大，我們希望找出使用次數(shù)與增大容積之間函數(shù)關(guān)系.試驗(yàn)數(shù)據(jù)如表3.4。 回歸分析MATLAB實(shí)現(xiàn)第25頁(yè)使用次數(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

17、；（2）預(yù)測(cè)鋼包使用x0=17次后增大容積y0；（3）計(jì)算回歸模型參數(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;回歸分析MATLAB實(shí)現(xiàn)第26頁(yè)%建立非線性雙曲線回歸模型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); % 擬合曲線pl

18、ot(x,y,*,x,y1,-or)legend(原始數(shù)據(jù),擬合曲線)注意：初始值要先計(jì)算后，才能得到上面程序中b0，因?yàn)榇_定兩個(gè)參數(shù)值，所以我們選擇已知數(shù)據(jù)中兩點(diǎn)（2,6.42）和（16,10.76）代入設(shè)定方程，得到方程組回歸分析MATLAB實(shí)現(xiàn)第27頁(yè)上述方程組有兩種解法：手工方法與Matlab方法。下面用Matlab方法解方程組：a,b=solve(6.42*(2*a+b)=2,10.76*(16*a+b)=16)輸出為a =.83961597702347450462657355615004e-1b =.14360328434608391527406223581049回歸分析MATLA

19、B實(shí)現(xiàn)第28頁(yè)圖3.11鋼包使用次數(shù)與增大容積非線性擬合圖在例3.1.5中，預(yù)測(cè)鋼包使用17次后增大容積，可在執(zhí)行上面程序中，繼續(xù)輸入命令ypred=nlpredci(fun,17,beta,r,J)得到：ypred =10.9599即鋼包使用17次后增大容積10.9599。回歸分析MATLAB實(shí)現(xiàn)第29頁(yè)求回歸模型參數(shù)95%置信區(qū)間，只要繼續(xù)添加程序ci = nlparci(beta,r,J)運(yùn)行后得到ci =0.0814 0.0876 0.0934 0.1370即回歸模型 中參數(shù)a,b95%置信區(qū)間分別為(0.0814 ,0.0876) 與（0.0934，0.1370）.我們求出最正確參數(shù)

20、分別為 a=0.0845,b=0.1152均屬于上述置信區(qū)間?；貧w分析MATLAB實(shí)現(xiàn)第30頁(yè)圖3.12給出鋼包使用次數(shù)與增大容積非線性擬合交互圖形，圖中圓圈是試驗(yàn)原始數(shù)據(jù)點(diǎn)，兩條虛線為95%上、下置信區(qū)間曲線（屏幕上顯示為紅色），中間實(shí)線（屏幕上顯示為綠色）是回歸模型曲線，縱向藍(lán)色虛線顯示了在自變量為8.9502，橫向虛線給出了對(duì)應(yīng)預(yù)測(cè)值為10.2734.圖3.12 鋼包使用次數(shù)與增大容積非線性擬合交互圖回歸分析MATLAB實(shí)現(xiàn)第31頁(yè)例3.1.6 對(duì)例題3.1.3進(jìn)行非線性回歸，并預(yù)測(cè)照射16次后細(xì)菌殘留數(shù)目，給出模型參數(shù)95%置信區(qū)間，繪出模型交互圖形.解：我們選取函數(shù)y=aebt進(jìn)行非

21、線性回歸，該方程兩個(gè)參數(shù)含有簡(jiǎn)單物了解釋，a表示試驗(yàn)開(kāi)始時(shí)細(xì)菌數(shù)目，b表示細(xì)菌死亡（或衰變）速率。MATLAB腳本程序以下：t=1:15;y=352 211 197 160 142 106 104 60 56 38 36 32 21 19 15;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ù)值計(jì)算回歸分析MATLAB實(shí)現(xiàn)第32頁(yè)plot(t,

22、y,*,t,y1,-or),legend(原始數(shù)據(jù),非線性回歸)xlabel(t(照射次數(shù)) ylabel(y(殘留細(xì)菌數(shù))ypred = nlpredci(fun,16,beta,r,J) % 預(yù)測(cè)殘留細(xì)菌數(shù)ci = nlparci(beta,r,J) % 參數(shù)95%區(qū)間預(yù)計(jì)nlintool(t,y,fun,beta0) % 作出交互圖形運(yùn)行后結(jié)果以下：beta = 400.0904 -0.2240即，最正確參數(shù)為：a=400.0904，b=-0.2240故非線性回歸模型為回歸分析MATLAB實(shí)現(xiàn)第33頁(yè) 預(yù)測(cè)為：ypred =11.1014即，照射16次后細(xì)菌殘留數(shù)目為11.1014，該預(yù)

23、測(cè)符合實(shí)際，顯然百分比3.1.3中多項(xiàng)式回歸結(jié)果合理。ci =355.2481 444.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ū)間。回歸分析MATLAB實(shí)現(xiàn)第34頁(yè)圖3.13原始數(shù)據(jù)與非線性回歸圖形圖3.14 原始數(shù)據(jù)與非線性回歸GUI圖形回歸分析MATLAB實(shí)現(xiàn)第35頁(yè)從交互圖形3.14能夠看出：圓圈為原始數(shù)據(jù)，兩條虛線（屏幕上顯示紅色）是

24、置信區(qū)間曲線；兩條虛線內(nèi)實(shí)線（屏幕上顯示綠色）是回歸模型曲線；縱向虛線指示照射8次，此時(shí)對(duì)應(yīng)水平虛線表示模型得到殘留細(xì)菌數(shù)為：66.6451。圖3.14 原始數(shù)據(jù)與非線性回歸GUI圖形回歸分析MATLAB實(shí)現(xiàn)第36頁(yè)3.1.3一元回歸建模實(shí)例例3.1.7在四川白鵝生產(chǎn)性能研究中，得到以下一組關(guān)于雛鵝重（g）與70日齡重(g)數(shù)據(jù)，試建立70日齡重(y)與雛鵝重(x)直線回歸方程，計(jì)算模型誤差平方和以及可決系數(shù)，當(dāng)雛鵝重分別為：85,95,115時(shí)預(yù)測(cè)其70日齡重，以及置信區(qū)間。 表3.5 四川白鵝重與70日齡重測(cè)定結(jié)果 （單位：g）編號(hào)123456789101112雛鵝重(x)80869890

25、120102958311310511010070日齡重(Y)235024002720250031502680263024003080292029602860回歸分析MATLAB實(shí)現(xiàn)第37頁(yè)解：（1）作散點(diǎn)圖。以雛鵝重（x）為橫坐標(biāo)，70日齡重（y）為縱坐標(biāo)作散點(diǎn)圖，如圖2-14。在MATLAB命令窗口中輸入:x=80 86 98 90 120 102 95 83 113 105 110 100; % 雛鵝重y=2350 2400 2720 2500 3150 2680 2630 2400 3080 2920 2960 2860; %70日齡重plot(x,y,*) %作散點(diǎn)圖xlabel(x(

26、雛鵝重) %橫坐標(biāo)名ylabel(y(70日齡重) %縱坐標(biāo)名圖3.15 四川白鵝雛鵝重與70日齡重散點(diǎn)圖和回歸直線圖回歸分析MATLAB實(shí)現(xiàn)第38頁(yè)由圖形3.15可見(jiàn)白鵝70日齡重與雛鵝重間存在直線關(guān)系，且70日齡重隨雛鵝重增大而增大。所以，可認(rèn)為y與x符合一元線性回歸模型。（2）建立直線回歸方程。在MATLAB中調(diào)用命令polyfit，從而求出參數(shù)0，1最小二乘預(yù)計(jì). 在MATLAB命令窗口中繼續(xù)輸入:n= size(x,1) % 計(jì)算樣本容量p,s=polyfit(x,y,1); % 調(diào)用命令polyfit計(jì)算回歸參數(shù)y1=polyval(p,x); % 計(jì)算回歸模型函數(shù)值hold on

27、plot(x,y1) % 作回歸方程圖形，結(jié)果如圖3.15p % 顯示參數(shù)最小二乘預(yù)計(jì)結(jié)果p=582.1850 21.7122回歸分析MATLAB實(shí)現(xiàn)第39頁(yè)即參數(shù) 最小二乘預(yù)計(jì)為所以70日齡重(y)與雛鵝重(x)直線回歸經(jīng)驗(yàn)方程為（3）誤差預(yù)計(jì)與決定系數(shù)。在MATLAB命令窗口中繼續(xù)輸入:TSS=sum(y-mean(y).2) %計(jì)算總離差平方和RSS=sum(y1-mean(y).2) %計(jì)算回歸平方和ESS=sum(y-y1).2) %計(jì)算殘差平方和R2=RSS/TSS; %計(jì)算樣本決定系數(shù)R2.回歸分析MATLAB實(shí)現(xiàn)第40頁(yè)輸出:TSS =8.314917e+005RSS =7.9

28、43396e+005ESS =3.715217e+004R2= 0.9553TSS =8.314917e+005RSS =7.943396e+005ESS =3.715217e+004R2=0.9553因?yàn)闃颖緵Q定系數(shù)R2=0.9553靠近于1，所以模型擬合效果很好。回歸分析MATLAB實(shí)現(xiàn)第41頁(yè)（4）回歸方程關(guān)系顯著性F檢驗(yàn)。在MATLAB命令窗口中繼續(xù)輸入:F=(n-2)*RSS/ESS %計(jì)算F統(tǒng)計(jì)量F1=finv(0.95,1,n-2) %查F統(tǒng)計(jì)量0.05分位數(shù)F2=finv(0.99,1,n-2) %查F統(tǒng)計(jì)量0.01分位數(shù)輸出結(jié)果:F=2.138e+002，F(xiàn)1 =4.9646

29、，F(xiàn)2 =10.0442為了方便，將以上計(jì)算結(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回歸分析MATLAB實(shí)現(xiàn)第42頁(yè)因?yàn)?表明四川白鵝70日齡重與雛鵝重間存在顯著線性關(guān)系。（5）回歸關(guān)系顯著性t檢驗(yàn)。在MATLAB命令窗口中繼續(xù)輸入:T=p(2)/sqrt(ESS/(n-2)*sqrt(sum(x-mean(x).2)%計(jì)算T統(tǒng)計(jì)量T1=tinv(0.975,

30、n-2) %t統(tǒng)計(jì)量0.05分位數(shù)T2=tinv(0.995,n-2) %t統(tǒng)計(jì)量0.01分位數(shù)輸出:T =14.622, T1 =2.228, T2 =3.169因?yàn)門(mén)=14.62t0.01(10)，否定H0，接收H1即四川白鵝70日齡重(y)與雛鵝重(x)線性回歸系數(shù)是顯著,可用所建立回歸方程進(jìn)行預(yù)測(cè)和控制。回歸分析MATLAB實(shí)現(xiàn)第43頁(yè) (6)預(yù)測(cè)x1=85,95,115; % 輸入自變量yc=polyval(p,x1) % 計(jì)算預(yù)測(cè)值Y,Delta=polyconf(p,x1,s); I1=Y-Delta,Y+Delta % 置信區(qū)間輸出:yc = 2427.72 2644.84 3

31、079.08I1 = 2279.47 2575.96 2503.01 2786.67 2927.55 3230.62所以當(dāng)雛鵝重分別為85,95,115時(shí)，白鵝70日齡重分別為2427.72, 2644.84, 3079.08;且95%置信區(qū)間分別為：2279.47 ,2575.96，2503.01,2786.67, 2927.55,3230.62.回歸分析MATLAB實(shí)現(xiàn)第44頁(yè)在程序中加入：polytool(x,y) % 交互功效bar(x,y-y1), % 殘差圖legend(殘差)h=lillietest(y-y1) % 殘差正態(tài)性檢驗(yàn)輸出h = 0得到交互圖形如圖3.16所表示，能夠

32、看出當(dāng)雛鵝重為100時(shí)，模型給出70日齡鵝重為2753.4016.圖3.16 四川白鵝70日齡重與雛鵝重線性模型交互圖回歸分析MATLAB實(shí)現(xiàn)第45頁(yè)3.2多元線性回歸模型3.2.1多元線性回歸模型及其表示回歸分析MATLAB實(shí)現(xiàn)第46頁(yè)對(duì)于總體n組觀察值它應(yīng)滿足式（3.2.1），即其中i (i=1,2,n)相互獨(dú)立，且設(shè) 記, , , 回歸分析MATLAB實(shí)現(xiàn)第47頁(yè)則模型(3.2.2)可用矩陣形式表示為 Y=X+ (3.2.3)其中Y稱為觀察向量，X稱為設(shè)計(jì)矩陣，稱為待預(yù)計(jì)向量，是不可觀察n維隨機(jī)向量，它分量相互獨(dú)立，假定 .2. 多元線性回歸建?；静襟E(1)對(duì)問(wèn)題進(jìn)行直觀分析，選擇因變

33、量與解釋變量，作出與因變量與各解釋變量散點(diǎn)圖，初步設(shè)定多元線性回歸模型參數(shù)個(gè)數(shù)；(2)輸入因變量與自變量觀察數(shù)據(jù)(y,X)調(diào)用命令 b,bint,r,rint,s=regress(y,X,alpha)，計(jì)算參數(shù)預(yù)計(jì)。(3)調(diào)用命令rcoplot(r,rint)，分析數(shù)據(jù)異常點(diǎn)情況?；貧w分析MATLAB實(shí)現(xiàn)第48頁(yè)(4)作顯著性檢驗(yàn)，若檢驗(yàn)經(jīng)過(guò)，則用模型作預(yù)測(cè)。(5)對(duì)模型深入研究：如殘差正態(tài)性檢驗(yàn)，殘差異方差檢驗(yàn)，殘差進(jìn)行自相關(guān)性檢驗(yàn)等。3.2.2 MATLAB回歸分析命令在MATLAB7.0統(tǒng)計(jì)工具箱中，與多元回歸模型相關(guān)命令有多個(gè)，下面逐一介紹。1.多元回歸建模命令regeress，其調(diào)用

34、格式有以下三種：（1）b = regress(y,X)（2）b,bint,r,rint,stats = regress(Y,X)（3）b,bint,r,rint,stats = regress(Y,X,alpha) 回歸分析MATLAB實(shí)現(xiàn)第49頁(yè)三種方式主要區(qū)分在輸出項(xiàng)參數(shù)多少上，第3種方式可稱為全參數(shù)方式。以第3種為例來(lái)說(shuō)明regeress命令輸入與輸出參數(shù)含義。輸入?yún)?shù)：輸入量Y表示模型（3.1.1）中因變量觀察向量；X是一個(gè)矩陣，其中第一列元全部是數(shù)“1”，第j列是自變量Xj觀察向量，即對(duì)一元線性回歸，取p=1即可；alpha為顯著性水平輸出參數(shù)：輸出向量b為回歸系數(shù)預(yù)計(jì)值，bint為

35、回歸系數(shù)(1-alpha)置信區(qū)間；輸出向量r表示殘差列向量回歸分析MATLAB實(shí)現(xiàn)第50頁(yè)輸出rint為模型殘差 (1- )置信區(qū)間；輸出stats是用于檢驗(yàn)回歸模型統(tǒng)計(jì)量，有4個(gè)分量值：第一個(gè)是R2，其中R是相關(guān)系數(shù)，第二個(gè)是F統(tǒng)計(jì)量值，第三個(gè)是與統(tǒng)計(jì)量F對(duì)應(yīng)概率P，當(dāng)P時(shí)拒絕H0，即認(rèn)為線性回歸模型有意義，第四個(gè)是方差2無(wú)偏預(yù)計(jì).例3.2.1某銷售企業(yè)將庫(kù)存占用資金情況、廣告投入費(fèi)用、員工薪酬以及銷售額等方面數(shù)據(jù)作了匯總,該企業(yè)試圖依據(jù)這些數(shù)據(jù)找到銷售額與其它變量之間關(guān)系，方便進(jìn)行銷售額預(yù)測(cè)并為工作決議提供參考依據(jù)。(1)建立銷售額回歸模型；(2)假如未來(lái)某月庫(kù)存資金額為150萬(wàn)元，廣告

36、投入預(yù)算為45萬(wàn)元，員工薪酬總額為27萬(wàn)元，試依據(jù)建立回歸模型預(yù)測(cè)該月銷售額。 回歸分析MATLAB實(shí)現(xiàn)第51頁(yè)表3.7 占用資金、廣告投入、員工薪酬、銷售額（單位：萬(wàn)元）月份庫(kù)存資金額(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

37、.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回歸分析MATLAB實(shí)現(xiàn)第52頁(yè)解：為了確定銷售額與庫(kù)存占用資金、廣告投入、員工薪酬之間關(guān)系，分別作出y與x1,x2,x3散點(diǎn)圖，若散點(diǎn)圖顯示它們之間近似線性關(guān)系，則可設(shè)定y與x1,x2,x3關(guān)系為三元線性回歸模型%輸入數(shù)據(jù)并作散點(diǎn)圖（圖3.18）A=75.2 30.6 21.

38、1 1090.4;77.6 31.3 21.4 113380.7 33.9 22.9 1242.1;76 29.6 21.4 1003.279.5 32.5 21.5 1283.2;81.8 27.9 21.7 1012.298.3 24.8 21.5 1098.8;67.7 23.6 21 826.374 33.9 22.4 1003.3;151 27.7 24.7 1554.690.8 45.5 23.2 1199;102.3 42.6 24.3 1483.1115.6 40 23.1 1407.1;125 45.8 29.1 1551.3137.8 51.7 24.6 1601.2;17

39、5.6 67.2 27.5 2311.7155.2 65 26.5 2126.7;174.3 65.4 26.8 2256.5;回歸分析MATLAB實(shí)現(xiàn)第53頁(yè)m,n=size(A);subplot(3,1,1),plot(A(:,1),A(:,4),+),xlabel(x1(庫(kù)存資金額) 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.1

40、8所表示，可見(jiàn)銷售額y與庫(kù)存資金、廣告投入、員工薪酬含有線性關(guān)系，所以能夠建立三元線性回歸模型.回歸分析MATLAB實(shí)現(xiàn)第54頁(yè)圖3.18銷售額與庫(kù)存、廣告、薪酬散點(diǎn)圖% 調(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é)果 回歸分析MATLAB實(shí)現(xiàn)第55頁(yè)程序運(yùn)行結(jié)果b =162.0632 7.2739 13.9575 -4.3996bint =-580.3603 904.4867 4.3734 10.1743

41、 7.1649 20.7501 -46.7796 37.9805 stats =0.9574804050 105.0866520891 0.0000000008 10077.9867891125輸出結(jié)果說(shuō)明：b就是模型中參數(shù)0 ,1 ,2 ，所以回歸模型為b就是模型中參數(shù)0 ,1 ,2 ，所以回歸模型為bint各行分別為參數(shù)0 ,1 ,295%置信區(qū)間。stats第一列為模型可決系數(shù)，第二列為F統(tǒng)計(jì)量觀察值，第三列得到概率p，最終一列為模型殘差平方和回歸分析MATLAB實(shí)現(xiàn)第56頁(yè)2.多元回歸輔助圖形命令（1）殘差圖命令rcoplot，其調(diào)用格式rcoplot(r,rint)其中，輸入?yún)?shù)r,

42、rint 是多元回歸建模命令regress 輸出結(jié)果，運(yùn)行該命令后展示了殘差與置信區(qū)間圖形。該命令有利于對(duì)建立模型進(jìn)行分析，假如圖形中出現(xiàn)紅色點(diǎn)，則能夠認(rèn)作異常點(diǎn)，此時(shí)可刪除異常點(diǎn)，重新建模，最終得到改進(jìn)回歸模型。在上面程序中加入 rcoplot(r,rint)得到以下列圖形回歸分析MATLAB實(shí)現(xiàn)第57頁(yè)圖3.19 殘差與置信區(qū)間圖 從圖形中能夠看到第五個(gè)點(diǎn)為異常點(diǎn)，實(shí)際上從表3.7能夠發(fā)覺(jué)第5個(gè)月庫(kù)存占用資金、廣告投入、員工薪酬均比3月份少，為何銷售額反而增加？這就能夠促使該企業(yè)經(jīng)理找出原因，尋找對(duì)策。下面例題介紹怎樣刪除異常點(diǎn)，對(duì)模型進(jìn)行改進(jìn)方法?；貧w分析MATLAB實(shí)現(xiàn)第58頁(yè)例3.2

43、.2 葛洲壩機(jī)組發(fā)電耗水率主要影響原因?yàn)閹?kù)水位,出庫(kù)流量。數(shù)據(jù)如表3.8所表示，利用多元線性回歸分析方法建立耗水率與出庫(kù)流量、庫(kù)水位模型。表3.8 某天耗水率與出庫(kù)流量、庫(kù)水位數(shù)據(jù) 時(shí)間年-月-天-時(shí) 庫(kù)水位（米)出庫(kù)流量（立方米） 機(jī)組發(fā)電耗水率 (立方米/萬(wàn)千瓦) -10-15:0065.081560760.46-10-15:0265.101556560.28-10-15:0465.121554060.10-10-15:0665.171550759.78-10-15:0865.211543259.44-10-15:1065.371561959.25-10-15:1265.381553658

44、.91-10-15:1465.391551458.76-10-15:1665.401551958.73-10-15:1865.431551058.63-10-15:2065.471548958.48-10-15:2265.531543758.31-10-16:0065.621635557.96-10-16:0265.581470857.06-10-16:0465.701439356.43-10-16:0665.841429655.83回歸分析MATLAB實(shí)現(xiàn)第59頁(yè)解：% 輸入原始數(shù)據(jù)A=65.08 15607 60.4665.10 15565 60.2865.12 15540 60.1065

45、.17 15507 59.7865.21 15432 59.4465.37 15619 59.2565.38 15536 58.9165.39 15514 58.7665.40 15519 58.7365.43 15510 58.6365.47 15489 58.4865.53 15437 58.3165.62 16355 57.9665.58 14708 57.0665.70 14393 56.4365.84 14296 55.83;% 做散點(diǎn)圖subplot(1,2,1),plot(A(:,1),A(:,3),+)xlabel(x1(庫(kù)水位) ylabel(y(耗水率)subplot(1,

46、2,2),plot(A(:,2),A(:,3),o)xlabel(x2(出庫(kù)流量) ylabel(y(耗水率)運(yùn)行后得到圖形如圖3.20所表示，從圖中能夠看到不論是庫(kù)水位還是出庫(kù)流量都與機(jī)組發(fā)電耗水率含有線性關(guān)系，所以，能夠建立機(jī)組發(fā)電耗水率與庫(kù)水位和出庫(kù)流量二元線性回歸模型?；貧w分析MATLAB實(shí)現(xiàn)第60頁(yè)圖3.20 庫(kù)水位、出庫(kù)流量與耗水率散點(diǎn)圖% 建立模型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)計(jì)量如表3.9所表示回歸分析

47、MATLAB實(shí)現(xiàn)第61頁(yè) 表3.9回歸模型系數(shù)、系數(shù)置信區(qū)間與統(tǒng)計(jì)量回歸系數(shù)回歸系數(shù)預(yù)計(jì)值回歸系數(shù)置信區(qū)間0373.8698340.082 ,407.65771-4.9759-5.4642 ,-4.4875 20.00070.0004,0.0009R2= 0.9863，F(xiàn)= 468.4118，p0.0001，s2 =0.0278由此可得模型為：% 模型改進(jìn) rcoplot(r,rint);得到圖形如圖3.21所表示，發(fā)覺(jué)有一個(gè)異常點(diǎn)，下面給出刪除異常點(diǎn)后，重新建模程序。由此可得模型為：回歸分析MATLAB實(shí)現(xiàn)第62頁(yè)圖3.21 殘差示意圖% 刪除異常點(diǎn)程序并建模b1,bint1,r1,rint

48、1,stats1=regress(y(1:12);y(14:m),ones(m-1,1),x(1:12,:);x(14:m,:)rcoplot(r1,rint1);刪除異常點(diǎn)后，殘差示意圖如圖2-21所表示，此時(shí)沒(méi)有異常點(diǎn)，改進(jìn)回歸模型系數(shù)、系數(shù)置信區(qū)間與統(tǒng)計(jì)量參見(jiàn)表3.10回歸分析MATLAB實(shí)現(xiàn)第63頁(yè)表3.10改進(jìn)回歸模型系數(shù)、系數(shù)置信區(qū)間與統(tǒng)計(jì)量回歸系數(shù)回歸系數(shù)預(yù)計(jì)值回歸系數(shù)置信區(qū)間0328.4616290.6145 ,366.30871-4.3594-4.8880 ,-3.830820.00100.00073,0.0012R2= 0.9931，F(xiàn)= 858.5846，p F,或者由F

49、0p值為0小于0.05，所以拒絕H0，接收備擇假設(shè)H1，說(shuō)明總體回歸系數(shù)i不全為零，即表明模型線性關(guān)系在95%置信水平下顯著成立。回歸分析MATLAB實(shí)現(xiàn)第73頁(yè)t檢驗(yàn)：對(duì)于=0.05，從表3.13 最終一列概率能夠看出均小于0.05，所以拒絕H0，接收備擇假設(shè)，即回歸系數(shù)i（i=0,1,2,3,4），顯著不為零。 回歸分析MATLAB實(shí)現(xiàn)第74頁(yè)方法2：從包含全部變量回歸方程中逐次剔除不顯著因子。首先建立包含全部變量回歸方程，然后對(duì)每一個(gè)因子作顯著性檢驗(yàn)，剔除不顯著因子中偏回歸平方和最小一個(gè)因子，重新建立包含全部變量（剔除除外）回歸方程。然后重復(fù)上面過(guò)程，對(duì)新建立回歸方程每一個(gè)因子作顯著性檢

50、驗(yàn)，剔除不顯著因子中偏回歸平方和最小因子，再重新建立回歸方程。如此，當(dāng)新建立回歸方程中全部因子都顯著時(shí)，回歸方程就是“最優(yōu)”了。這種方法在因子，尤其是不顯著因子不多時(shí)，能夠采取。但計(jì)算工作量依然可能較大?；貧w分析MATLAB實(shí)現(xiàn)第75頁(yè)3.3逐步回歸 在建立經(jīng)濟(jì)預(yù)測(cè)問(wèn)題數(shù)學(xué)模型時(shí)，經(jīng)常從可能影響預(yù)測(cè)量Y許多原因中挑選一批原因作為自變量，應(yīng)用回歸分析方法建立回歸方程作預(yù)報(bào)或控制用。問(wèn)題是怎樣在為數(shù)眾多原因中挑選變量，以建立我們稱為是這批觀察數(shù)據(jù)“最優(yōu)”回歸方程。3.3.1最優(yōu)回歸方程選擇回歸分析MATLAB實(shí)現(xiàn)第76頁(yè)選擇“最優(yōu)”回歸方程有以下幾個(gè)方法：方法1：從全部可能變量組合回歸方程中挑選最

51、優(yōu)者，即把全部包含1個(gè)，2個(gè)，直至全部變量線性回歸方程全部計(jì)算，對(duì)每個(gè)方程及自變量作顯著性檢驗(yàn)，然后從中挑選一個(gè)方程，要求全部變量全部顯著，且剩下均方和MSE較小。這種方法只在變量較少時(shí)可用?；貧w分析MATLAB實(shí)現(xiàn)第77頁(yè)方法3：從一個(gè)變量開(kāi)始，把變量逐一引入回歸方程。這一方法首先計(jì)算各因子與Y相關(guān)系數(shù)，將絕對(duì)值最大一個(gè)因子引入方程，并對(duì)回歸平方和進(jìn)行檢驗(yàn)，若顯著則引入。然后找出余下因子中與Y偏相關(guān)系數(shù)最大那個(gè)因子，將其引入方程，檢驗(yàn)顯著性，等等，當(dāng)引入因子建立方程檢驗(yàn)不顯著時(shí)，該因子就不再引入。這種方法盡管工作量較小，但并不確保最終所得到方程是“最優(yōu)”，還得深入作檢驗(yàn)，剔除不顯著因子。同時(shí)

52、這種方法每一步要計(jì)算偏相關(guān)系數(shù)，也較麻煩?；貧w分析MATLAB實(shí)現(xiàn)第78頁(yè)結(jié)合方法與，產(chǎn)生了一個(gè)建立“最優(yōu)”回歸方程方法逐步回歸分析。逐步回歸基本思想是，將變量一個(gè)一個(gè)引入，引入變量條件是偏回歸平方和經(jīng)檢驗(yàn)是顯著，同時(shí)每引入一個(gè)新變量后，對(duì)已選入變量要進(jìn)行逐一檢驗(yàn)，將不顯著變量剔除， 詳細(xì)做法是將變量一個(gè)一個(gè)引入，當(dāng)每引入一個(gè)自變量后，對(duì)已選變量要進(jìn)行逐一檢驗(yàn)，當(dāng)原引入變量因?yàn)楹竺孀兞恳攵兊貌辉亠@著時(shí)，要將其剔除。引入一個(gè)變量或從回歸方程中剔除一個(gè)變量，為逐步回歸一步，每一步都要進(jìn)行F檢驗(yàn)，以確保每次引入新變量之前回歸方程中只包含顯著變量?；貧w分析MATLAB實(shí)現(xiàn)第79頁(yè)這個(gè)過(guò)程重復(fù)進(jìn)行，直到既無(wú)顯著自變量選入回歸方程，也無(wú)不顯著自變量從回歸方程中剔除為止。這么確保最終所得變量子集中全部變量都是顯著。這么經(jīng)若干步以后便得“最優(yōu)”變量子集?；貧w分析MATLAB實(shí)現(xiàn)第80頁(yè)3.3.2逐步回歸MAT