matlab回歸分析方法

1、第八章回歸分析方法當(dāng)人們對(duì)研究對(duì)象的內(nèi)在特性和各因素間的關(guān)系有比較充分的認(rèn)識(shí)時(shí)，一般用機(jī)理分析方法建立數(shù)學(xué)模型。如果由于客觀事物內(nèi)部規(guī)律的復(fù)雜性及人們認(rèn)識(shí)程度的限制，無(wú)法分析實(shí)際對(duì)象內(nèi)在的因果關(guān)系，建立合乎機(jī)理規(guī)律的數(shù)學(xué)模型，那么通常的辦法是搜集大量數(shù)據(jù)，基于對(duì)數(shù)據(jù)的統(tǒng)計(jì)分析去建立模型。本章討論其中用途非常廣泛的一類模型一一統(tǒng)計(jì)回歸模型?；貧w模型常用來(lái)解決預(yù)測(cè)、控制、生產(chǎn)工藝優(yōu)化等問(wèn)題。變量之間的關(guān)系可以分為兩類：一類叫確定性關(guān)系，也叫函數(shù)關(guān)系，其特征是：一個(gè)變量隨著其它變量的確定而確定。另一類關(guān)系叫相關(guān)關(guān)系，變量之間的關(guān)系很難用一種精確的方法表示出來(lái)。例如，通常人的年齡越大血壓越高，但人的年

2、齡和血壓之間沒(méi)有確定的數(shù)量關(guān)系，人的年齡和血壓之間的關(guān)系就是相關(guān)關(guān)系?；貧w分析就是處理變量之間的相關(guān)關(guān)系的一種數(shù)學(xué)方法。其解決問(wèn)題的大致方法、步驟如下：(1)收集一組包含因變量和自變量的數(shù)據(jù)；(2)選定因變量和自變量之間的模型，即一個(gè)數(shù)學(xué)式子，利用數(shù)據(jù)按照最小二乘準(zhǔn)則計(jì)算模型中的系數(shù)；(3)利用統(tǒng)計(jì)分析方法對(duì)不同的模型進(jìn)行比較，找出與數(shù)據(jù)擬合得最好的模型；(4)判斷得到的模型是否適合于這組數(shù)據(jù)；(5)利用模型對(duì)因變量作出預(yù)測(cè)或解釋。應(yīng)用統(tǒng)計(jì)分析特別是多元統(tǒng)計(jì)分析方法一般都要處理大量數(shù)據(jù)，工作量非常大，所以在計(jì)算機(jī)普及以前，這些方法大都是停留在理論研究上。運(yùn)用一般計(jì)算語(yǔ)言編程也要占用大量時(shí)間，而

3、對(duì)于經(jīng)濟(jì)管理及社會(huì)學(xué)等對(duì)高級(jí)編程語(yǔ)言了解不深的人來(lái)說(shuō)要應(yīng)用這些統(tǒng)計(jì)方法更是不可能。MATLAB等軟件的開(kāi)發(fā)和普及大大減少了對(duì)計(jì)算機(jī)編程的要求，使數(shù)據(jù)分析方法的廣泛應(yīng)用成為可能。MATLAB統(tǒng)計(jì)工具箱幾乎包括了數(shù)理統(tǒng)計(jì)方面主要的概念、理論、方法和算法。運(yùn)用MATLAB統(tǒng)計(jì)工具箱，我們可以十分方便地在計(jì)算機(jī)上進(jìn)行計(jì)算，從而進(jìn)一步加深理解，同時(shí)，其強(qiáng)大的圖形功能使得概念、過(guò)程和結(jié)果可以直觀地展現(xiàn)在我們面前。本章內(nèi)容通常先介紹有關(guān)回歸分析的數(shù)學(xué)原理，主要說(shuō)明建模過(guò)程中要做的工作及理由，如模型的假設(shè)檢驗(yàn)、參數(shù)估計(jì)等，為了把主要精力集中在應(yīng)用上，我們略去詳細(xì)而繁雜的理論。在此基礎(chǔ)上再介紹在建模過(guò)程中如何有

4、效地使用MATLAB軟件。沒(méi)有學(xué)過(guò)這部分?jǐn)?shù)學(xué)知識(shí)的讀者可以不深究其數(shù)學(xué)原理，只要知道回歸分析的目的，按照相應(yīng)方法通過(guò)軟件顯示的圖形或計(jì)算所得結(jié)果表示什么意思，那么，仍然可以學(xué)到用回歸模型解決實(shí)際問(wèn)題的基本方法。包括：一元線性回歸、多元線性回歸、非線性回歸、逐步回歸等方法以及如何利用MATLAB軟件建立初步的數(shù)學(xué)模型，如何透過(guò)輸出結(jié)果對(duì)模型進(jìn)行分析和改進(jìn)，回歸模型的應(yīng)用等。8.1 一元線性回歸分析回歸模型可分為線性回歸模型和非線性回歸模型。非線性回歸模型是回歸函數(shù)關(guān)于未知參數(shù)具有非線性結(jié)構(gòu)的回歸模型。某些非線性回歸模型可以化為線性回歸模型處理；如果知道函數(shù)形式只是要確定其中的參數(shù)則是擬合問(wèn)題，可

5、以使用MATLAB軟件的curvefit命令或nlinfit命令擬合得到參數(shù)的估計(jì)并進(jìn)行統(tǒng)計(jì)分析。本節(jié)主要考察線性回歸模型。8.1.1 一元線性回歸模型的建立及其MATLAB實(shí)現(xiàn)y=，：1x；n(0,二2)其中吃,久是待定系數(shù)，對(duì)于不同的x,y是相互獨(dú)立的隨機(jī)變量。假設(shè)于x的n個(gè)值Xi,得到y(tǒng)的n個(gè)相應(yīng)的值y,確定P。，用的方法是根據(jù)最小二乘準(zhǔn)則，要使nnQ-o；i）:八x-（f-Xi）2i1i1取最小值。利用極值必要條件令JQ_=0,EQ=0,求Po,Pi的估計(jì)值照邑，從而得到：F：1回歸直線y二區(qū)+x。只不過(guò)這個(gè)過(guò)程可以由軟件通過(guò)直線擬合完成，而無(wú)須進(jìn)行繁雜的運(yùn)算。(1)參數(shù)的區(qū)間估計(jì)由于

6、我們所計(jì)算出的照因仍然是隨機(jī)變量，因此要對(duì)陽(yáng)，因取值的區(qū)間進(jìn)行估計(jì)，如果區(qū)間估計(jì)值是一個(gè)較短的區(qū)間表示模型精度較高。(2)對(duì)誤差方差的估計(jì)設(shè)？為回歸函數(shù)的值，yi為測(cè)量值，殘差平方和nQ=£(y？)2i1剩余方差s2=-Q-n-2(3)線性相關(guān)性的檢驗(yàn)由于我們采用的是一元線性回歸，因此，如果模型可用的話，應(yīng)該具有較好的線性關(guān)系。反映模型是否具有良好線性關(guān)系可通過(guò)相關(guān)系數(shù)R的值及F值觀察(后面的例子說(shuō)明)。(4)一元線性回歸的MATLAB實(shí)現(xiàn)MATLAB工具箱中用命令regress實(shí)現(xiàn)，其用法是：b=regress(y,x)b,bint,r,rint,s=regress(y,x,alp

7、ha)輸入y(因變量，列向量)、x(1與自變量組成的矩陣，見(jiàn)下例)，alpha是顯著性水平(缺省時(shí)默認(rèn)0.05)o輸出b=(£),因)，注意：b中元素順序與擬合命令polyfit的輸出不同，bint是P0,Pi的置信區(qū)間，r是殘差(列向量)，rint是殘差的置信區(qū)間，s包含4個(gè)統(tǒng)計(jì)量：決定系數(shù)R2(相2關(guān)系數(shù)為R);F值；F(1,n-2)分布大于F值的I率p；剩余萬(wàn)差s的值(MATLAB7.0以后版本)。s2也可由程序sum(12)/(n-2)計(jì)算。其意義和用法如下：R2的值越接近1,變量的線性相關(guān)性越強(qiáng)，說(shuō)明模型有效；如果滿足FiX1，n-2)MF,則認(rèn)為變量y與x顯著地有線性關(guān)系

8、，其中冗口(1刀-2)的值可查F分布表，或直接用MATLAB命令finv（1-a,1,n-2）計(jì)算得到；如果p<a表示線性模型可用。這三個(gè)值可以相互印證。s2的值主要用來(lái)比較模型是否有改進(jìn)，其值越小說(shuō)明模型精度越高。8.1.2身高與腿長(zhǎng)例1測(cè)得16名成年女子身高y與腿長(zhǎng)x所得數(shù)據(jù)如下：表8-116名女子身高（cm）腿長(zhǎng)（cm）數(shù)據(jù)x8885889192939395969897969899100102y143145146147149150153154155156157158159160162164首先利用命令plot(x,y,'r*')畫出散點(diǎn)圖，從圖形可以看出，這些點(diǎn)大致

9、分布在一條直線的左右，因此，可以考慮一元線性回歸?？删幹瞥绦蛉缦拢簓=143145146147149150153154155156157158159160162164;x=8885889192939395969897969899100102;n=16;X=ones(n,1),x'b,bint,r,rint,s=regress(y',X,0.05);b,bint,s,rcoplot(r,rint)運(yùn)行后得到b=31.77131.2903bint=12.319651.22291.08461.4960s=0.9282180.95310.00003.12772R=0.9282,由fin

10、v(0.95,1,14)=4.6001,即FMX(1,n-2)=4.6001<F=180.9531,p<0.0001,可以通過(guò)殘差圖發(fā)現(xiàn)，第二個(gè)數(shù)據(jù)為奇異數(shù)據(jù)，去掉該數(shù)據(jù)后運(yùn)行后得到b=17.65491.4363bint=-0.598635.90831.24451.6281s=0.9527261.63890.00001.9313R2=0.9527,由finv(0.95,1,13)=4.6672,即Fig(1,n2)=4.6672<F=261.6389,p<0.0001,說(shuō)明模型有效且有改進(jìn)，因此我彳門得到身高與腿長(zhǎng)的關(guān)系y=17.6549+1.4363x。當(dāng)然，也可以利

11、用直線擬合得到同一方程。只不過(guò)不能得到參數(shù)置信區(qū)間和對(duì)模型進(jìn)行檢驗(yàn)。擬合程序如下：y=143145146147149150153154155156157158159160162164;x=8885889192939395969897969899100102;a=polyfit(x,y,1)temp=polyval(a,x);plot(x,y,'r*',x,temp)注意：函數(shù)相同，但輸出一次函數(shù)參數(shù)順序與回歸分析（升哥排列）中不同。另一個(gè)差別是擬合不能發(fā)現(xiàn)奇異數(shù)據(jù)。8.2多元線性回歸分析8.2.1 多元線性回歸模型的建模步驟及其MATLAB實(shí)現(xiàn)如果根據(jù)經(jīng)驗(yàn)和有關(guān)知識(shí)認(rèn)為與因變量

12、有關(guān)聯(lián)的自變量不止一個(gè)，那么就應(yīng)該考慮用最小二乘準(zhǔn)則建立多元線性回歸模型。設(shè)影響因變量y的主要因素（自變量）有m個(gè)，記x=（Xi,|,Xm），假設(shè)它們有如下的線性關(guān)系式：y=P。+口1為十|+PmXm+W,WN（82）如果對(duì)變量y與自變量Xi,X2,|,Xm同時(shí)作n次觀察（n>m）得n組觀察值，采用最小二乘估計(jì)求得回歸方程y?=?0?x巾Z%.建立回歸模型是一個(gè)相當(dāng)復(fù)雜的過(guò)程，概括起來(lái)主要有以下幾個(gè)方面工作（1）根據(jù)研究目的收集數(shù)據(jù)和預(yù)分析；（2）根據(jù)散點(diǎn)圖是否具有線性關(guān)系建立基本回歸模型；（3）模型的精細(xì)分析；（4）模型的確認(rèn)與應(yīng)用等。收集數(shù)據(jù)的一個(gè)經(jīng)驗(yàn)準(zhǔn)則是收集的數(shù)據(jù)量（樣本容量）至

13、少應(yīng)為可能的自變量數(shù)目的610倍。在建模過(guò)程中首先要根據(jù)所研究問(wèn)題的目的設(shè)置因變量，然后再選取與該因變量區(qū)可以在得到初步的模型后利用MATLAB軟件進(jìn)行相關(guān)性檢驗(yàn)。下有統(tǒng)計(jì)關(guān)系的一些變量作為自變量。我們當(dāng)然希望選擇與問(wèn)題關(guān)系密切的變量面通過(guò)一個(gè)案例探討MATLAB軟件在回歸分析建模各個(gè)環(huán)節(jié)中如何應(yīng)用。多元線性回歸的MATLAB實(shí)現(xiàn)仍然用命令regress（y,X）,只是要注意矩陣X的形式，將通過(guò)如下例子說(shuō)明其用法。8.2.2 某類研究學(xué)者的年薪1 .問(wèn)題例2工薪階層關(guān)心年薪與哪些因素有關(guān)，以此可制定出它們自己的奮斗目標(biāo)。某科學(xué)基金會(huì)希望估計(jì)從事某研究的學(xué)者的年薪Y(jié)與他們的研究成果（論文、著作等

14、）的質(zhì)量指標(biāo)Xi、從事研究工作的時(shí)間X2、能成功獲得資助的指標(biāo)X3之間的關(guān)系，為此按一定的實(shí)驗(yàn)設(shè)計(jì)方法調(diào)查了24位研究學(xué)者，得到如下數(shù)據(jù)（i為學(xué)者序號(hào)）：表8-2從事某種研究的學(xué)者的相關(guān)指標(biāo)數(shù)據(jù)i 1234569101112Xii3.55.35.15.84.26.06.85.53.17.24.54.9*9201833311325305472511X36.16.47.46.77.55.96.04.05.88.35.06.4小33.240.338.746.841.437.539.040.730.152.938.231.8i1314151617181920212223248.06.56.63.76.

15、27.04.04.55.95.64.83.9xi223353921740352333273415x37.67.05.04.45.57.06.03.54.94.38.05.8y43.344.142.533.634.248.038.035.940.436.845.235.1試建立Y與XX2,X3之間關(guān)系的數(shù)學(xué)模型，并得出有關(guān)結(jié)論和作統(tǒng)計(jì)分析。2 .作出因變量Y與各自變量的樣本散點(diǎn)圖作散點(diǎn)圖的目的主要是觀察因變量Y與各自變量間是否有比較好的線性關(guān)系，以便選擇恰當(dāng)?shù)臄?shù)學(xué)模型形式。下圖分別為年薪丫與成果質(zhì)量指標(biāo)X1、研究工作時(shí)間X2、獲得資助的指標(biāo)X3之間的散點(diǎn)圖，1subplot(1,3,1),plo

16、t(x1,Y,'g*'),%subplot是MATLABH勺函數(shù)。使用方法：subplot(m,n,p)或者subplot(mnp)。subplot是將多個(gè)圖畫到一個(gè)平面上的工具。其中，m表示是圖排成m行，n表示圖排成n歹U,也就是整個(gè)figure中有n個(gè)圖是排成一行的，一共m行，如果m=2就是表示2行圖。p表示圖所在的位置，p=1表示從左到右從上到下的第一個(gè)位置。subplot(1,3,2),plot(x2,Y,'k+'),subplot(1,3,3),plot(x3,Y,'ro'),丫與x1的散點(diǎn)圖丫與x2的散點(diǎn)圖55305045403505

17、10丫與x3的散點(diǎn)圖從圖可以看出這些點(diǎn)大致分布在一條直線旁邊，因此，有比較好的線性關(guān)系，可以采用線性回歸。圖8.1因變量Y與各自變量的樣本散點(diǎn)圖3.利用MATLAB統(tǒng)計(jì)工具箱得到初步的回歸方程設(shè)回歸方程為：?=?01x1?2x33x3.建立m-文件輸入如下程序數(shù)據(jù)：x1=3.55.35.15.84.26.06.85.53.17.24.54.98.06.56.53.76.27.04.04.55.95.64.83.9;x2=920183331132530547251123353921740352333273415;x3=6.16.47.46.77.55.96.04.05.88.35.06.47.6

18、7.05.04.05.57.06.03.54.94.38.05.0;Y=33.240.338.746.841.437.539.040.730.152.938.231.843.344.142.533.634.248.038.035.940.436.845.235.1;n=24;m=3;X=ones(n,1),x1',x2',x3'b,bint,r,rint,s=regress(Y',X,0.05);b,bint,r,rint,s,運(yùn)行后即得到結(jié)果如表8-3所示。表8-3對(duì)初步回歸模型的計(jì)算結(jié)果回歸系數(shù)回歸系數(shù)的的值回歸系數(shù)的置信區(qū)間B。18.015713.9052

19、22.1262向1.08170.39001.77330.32120.24400.39841.28350.66911.897922R=0.9106F=67.9195p<0.0001S=3.0719計(jì)算結(jié)果包括回歸系數(shù)b=(P0,A,P2,3)=(18.0157,1.0817,0.3212,1.2835),且置信區(qū)間均不包含零點(diǎn)，；殘差及其置信區(qū)間；統(tǒng)計(jì)變量stats，它包含四個(gè)檢驗(yàn)統(tǒng)計(jì)量：相關(guān)系數(shù)的平方222.R，假設(shè)檢驗(yàn)統(tǒng)計(jì)量F,與F對(duì)應(yīng)的I率p,S的值(7.0以前版本S也可由程序sum(r.A2)/(n-m-1)計(jì)算)。因此我們得到初步的回歸方程為：y?=18.01571.0817x1

20、0.3212x21.2835x3由結(jié)果對(duì)模型的判斷：回歸系數(shù)置信區(qū)間不包含零點(diǎn)表示模型較好，殘差在零點(diǎn)附近也表示模型較好，接著就是利用檢驗(yàn)統(tǒng)計(jì)量R,F,p的值判斷該模型是否可用。(1)相關(guān)系數(shù)R的評(píng)價(jià)：一般地，相關(guān)系數(shù)絕對(duì)值在0.81范圍內(nèi)，可判斷回歸自變量與因變量具有較強(qiáng)的線性相關(guān)性。本例R的絕對(duì)值為0.9542,表明線性相關(guān)性較強(qiáng)。（2） F檢驗(yàn)法：當(dāng)f>F1g(m,nm1),即認(rèn)為因變量y與自變量x1,x2,lll,xm之間顯著地有線性相關(guān)關(guān)系；否則認(rèn)為因變量y與自變量為，*2,111,xm之間線性相關(guān)關(guān)系不顯著。本例F=67.919>F1_0.05(3,20)=3.10(查

21、F分布表或輸入命令finv(0.95,3,20)計(jì)算)。（3） p值檢驗(yàn)：若p<a（a為預(yù)定顯著水平），則說(shuō)明因變量y與自變量七,用,|?之間顯著地有線性相關(guān)關(guān)系。本例輸出結(jié)果，p<0.0001，顯然滿足P<o（=0.05。以上三種統(tǒng)計(jì)推斷方法推斷的結(jié)果是一致的，說(shuō)明因變量y與自變量之間顯著地有線性相關(guān)關(guān)系，所得線性回歸模型可用。s2當(dāng)然越小越好，這主要在模型改進(jìn)時(shí)作為參考。4.模型的精細(xì)分析和改進(jìn)（1）殘差分析殘差e=yi?（i=1,2，川，n）,是各觀測(cè)值yi與回歸方程所對(duì)應(yīng)得到的擬合值y?之差，實(shí)際上，它是線性回歸模型中誤差君的估計(jì)值。wN（0產(chǎn)2）即有零均值和常值方差

22、,利用殘差的這種特性反過(guò)來(lái)考察原模型的合理性就是殘差分析的基本思想。利用MATLAB進(jìn)行殘差分析則是通過(guò)殘差圖或時(shí)序殘差圖。殘差圖是指以殘差為縱坐標(biāo)，以其他指定的量為橫坐標(biāo)的散點(diǎn)圖。主要包括：（1）橫坐標(biāo)為觀測(cè)時(shí)間或觀測(cè)值序號(hào)；（2）橫坐標(biāo)為某個(gè)自變量的觀測(cè)值；（3）橫坐標(biāo)為因變量的擬合值。通過(guò)觀察殘差圖，可以對(duì)奇異點(diǎn)進(jìn)行分析，還可以對(duì)誤差的等方差性以及對(duì)回歸函數(shù)中是否包含其他自變量、自變量的高次項(xiàng)及交叉項(xiàng)等問(wèn)題給出直觀的檢驗(yàn)。以觀測(cè)值序號(hào)為橫坐標(biāo)，殘差為縱坐標(biāo)所得到的散點(diǎn)圖稱為時(shí)序殘差圖，畫出時(shí)序殘差圖的MATLAB語(yǔ)句為rcoplot（r,rint）（圖8.2）?？梢郧宄吹綒埐畲蠖挤植荚?/p>

23、零的附近，因此還是比較好的，不過(guò)第4、12、19這三個(gè)樣本點(diǎn)的殘差偏離原點(diǎn)較遠(yuǎn)，如果作為奇異點(diǎn)看待，去掉后重新擬合，則得回歸模型為：y?=19.08080.8616x10.3176x21.3463x3且回歸系數(shù)的置信區(qū)間更小均不包含原點(diǎn)，統(tǒng)計(jì)變量stats包含的三個(gè)檢驗(yàn)統(tǒng)計(jì)量：相關(guān)系數(shù)的平方R2,假設(shè)檢驗(yàn)統(tǒng)計(jì)量F,概率P,分別為：0.9533;115.5586;0.0000,比較可知R,F均增加模型得到改進(jìn)。圖8.2時(shí)序殘差圖(2)變量間的交互作用討論變量間的交互作用包括：不同自變量之間的交互作用以及同一變量的自相關(guān)性。不同自變量之間的交互作用：有時(shí)，在實(shí)驗(yàn)中不僅單因素對(duì)指標(biāo)有影響，而且因素間

24、還會(huì)聯(lián)合起來(lái)對(duì)指標(biāo)產(chǎn)生影響，常稱這種聯(lián)合作用為交互作用。處理兩個(gè)因素間交互作用的一個(gè)簡(jiǎn)單辦法是加入這兩個(gè)自變量的乘積項(xiàng)。本文案例如果加入交互項(xiàng)則為：?=?0?X1裝么3Zx“彳23'2X3用表8.2的數(shù)據(jù)，利用MATLAB統(tǒng)計(jì)工具箱得到回歸系數(shù)分別為：27.0727,1.1147,-0.0215,-0.1843,0.0033,-0.0054,0.0511。但它們的置信區(qū)間均包含原點(diǎn)，其他指標(biāo)也不理想，因此，本例中其交互作用并不顯著，該模型不如前面兩個(gè)模型好。自相關(guān)性的診斷和處理：若數(shù)據(jù)是以時(shí)間為序的，稱為時(shí)間序列數(shù)據(jù)。在時(shí)間序列數(shù)據(jù)中，同一變量的順序觀測(cè)值之間出現(xiàn)的相關(guān)現(xiàn)象稱為自相關(guān)。

25、一旦數(shù)據(jù)中存在這種自相關(guān)序列，如果仍采用普通的回歸模型直接處理，將產(chǎn)生不良后果，使預(yù)測(cè)失去意義。自相關(guān)的診斷主要有圖示檢驗(yàn)法、相關(guān)系數(shù)法和DW檢驗(yàn)法。圖示檢驗(yàn)法是通過(guò)繪制殘差et散點(diǎn)圖觀察，如果散布點(diǎn)(et.ej,t=2,3,|1大部分點(diǎn)落在第1,印象限，表明存在著正的序列DW檢驗(yàn)法可以利相關(guān)；如果大部分點(diǎn)落在第n,iv象限，表明存在著負(fù)的序列相關(guān)。對(duì)用MATLAB軟件編程計(jì)算統(tǒng)計(jì)量：netet1DW定2(1?),?=i廠，:e2.然后查閱DW檢驗(yàn)上下界表，以決定模型的自相關(guān)狀態(tài)。當(dāng)一個(gè)回歸模型存在序列相關(guān)性時(shí)，首先要查明序列相關(guān)產(chǎn)生的原因。如果是回歸模型選用不當(dāng)，則應(yīng)改用適當(dāng)?shù)幕貧w模型；如果

26、是缺少重要的自變量，則應(yīng)增加自變量；如果以上方法都不能消除序列相關(guān)性，則需要采用差分法、迭代法等處理，更詳細(xì)內(nèi)容參見(jiàn)相關(guān)概率統(tǒng)計(jì)參考文獻(xiàn)。8.2.3逐步回歸方法建模逐步回歸就是一種從眾多自變量中有效地選擇重要變量的方法。逐步回歸的基本思路是，先確定一個(gè)包含若干自變量的初始集合，然后每次從集合外的變量中引入一個(gè)對(duì)因變量影響最大的，再對(duì)集合中的變量進(jìn)行檢驗(yàn)，從變得不顯著的變量中移出一個(gè)影響最小的，依此進(jìn)行，直到不能引入和移出為止。引入和移出都以給定的顯著性水平為標(biāo)準(zhǔn)。MATLAB統(tǒng)計(jì)工具箱中逐步回歸的命令是stepwise,它提供了一個(gè)人機(jī)交互式畫面，通過(guò)此工具可以自由地選擇變量進(jìn)行統(tǒng)計(jì)分析。該命

27、令的用法是：stepwise(X,Y,inmodel,alpha)其中X是自變量數(shù)據(jù)，排成nMm矩陣(m為自變量個(gè)數(shù)，n為每個(gè)變量的數(shù)據(jù)量)，Y是因變量數(shù)據(jù)，排成口父1向量，inmodel是自變量初始集合的指標(biāo)，缺省時(shí)為全部自變量，alpha為顯著水平，缺省時(shí)為0.05。運(yùn)行stepwise命令時(shí)產(chǎn)生圖形窗口：StepwisePlot,StepwiseTable,StepwiseHistory.當(dāng)鼠標(biāo)移到圖形某個(gè)區(qū)域時(shí)，鼠標(biāo)點(diǎn)擊后產(chǎn)生交互作用。StepwisePlot窗口中的虛線表示回歸系數(shù)的置信區(qū)間包含零點(diǎn)，即該回歸系數(shù)與零無(wú)顯著差異，一般應(yīng)將該變量移去；實(shí)線則表明該回歸系數(shù)與零有顯著差異，

28、應(yīng)保留在模型中(藍(lán)色表示該變量已進(jìn)入模型，紅色表示該變量已移出模型)。引入和移出變量還可參考StepwiseHistory窗口中剩余標(biāo)準(zhǔn)差RMSE是否在下降，剩余標(biāo)準(zhǔn)差RMSE最小的就是最好的模型。StepwiseTable窗口中列出了一個(gè)統(tǒng)計(jì)表，包括回歸系數(shù)及其置信區(qū)間，以及模型的統(tǒng)計(jì)量剩余標(biāo)準(zhǔn)差RMSE、相關(guān)系數(shù)R-square、F值、與F對(duì)應(yīng)的概率。關(guān)于本節(jié)案例2,如果引入新的自變量x4=x1x2,x5=x1x3,x6=x2x3.也可以采用逐步回歸法解決，源程序如下：A=3.55.35.15.84.26.06.85.53.17.24.54.98.06.56.53.76.27.04.04.

29、55.95.64.83.9;920183331132530547251123353921740352333273415;6.16.47.46.77.55.96.04.05.88.35.06.47.67.05.04.05.57.06.03.54.94.38.05.0'Y=33.240.338.746.841.437.539.040.730.152.938.231.843.344.142.533.634.248.038.035.940.436.845.235.1'x1=A(:,1);x2=A(:,2);x3=A(:,3);x4=x1.*x2;x5=x1.*x3;x6=x2.*x3;

30、X=A,x4,x5,x6;stepwise(X,Y)運(yùn)行并按上述步驟操作后可以得到本文前面線性回歸相同的結(jié)論，即不含交互項(xiàng)的模型是最好的。在此只介紹操作過(guò)程，其交互界面，只要在MATLAB軟件上一試便知。8.2.4多項(xiàng)式回歸多項(xiàng)式回歸仍然屬于多元線性回歸，可以是一元多項(xiàng)式回歸或多元多項(xiàng)式回歸。一元多項(xiàng)式回歸模型的一般形式為¥=-0：1x-III-mxm用MATLAB求解一元多項(xiàng)式回歸，除了使用命令polyfit(x,y,m)外，還可以使用如下命令：Polytool(x,y,m,alpha)輸入x,y,m同命令polyfit,alpha是顯著性水平(默認(rèn)0.05),則輸出一個(gè)交互式畫面

31、，畫面顯示回歸曲線及其置信區(qū)間，通過(guò)圖左下方的export下拉式菜單，還可以輸出回歸系數(shù)估計(jì)值及其置信區(qū)間、殘差等。下面通過(guò)一個(gè)用多元多項(xiàng)式回歸的實(shí)例說(shuō)明什么時(shí)候用多項(xiàng)式回歸以及如何通過(guò)MATLAB軟件進(jìn)行處理。例3為了了解人口平均預(yù)期壽命與人均國(guó)內(nèi)生產(chǎn)總值和體質(zhì)得分的關(guān)系，我們查閱了國(guó)家統(tǒng)計(jì)局資料，北京體育大學(xué)出版社出版的2000國(guó)民體質(zhì)監(jiān)測(cè)報(bào)告，表8-4是我國(guó)大陸31個(gè)省市的有關(guān)數(shù)據(jù)。我們希望通過(guò)這幾組數(shù)據(jù)考察它們是否具有良好的相關(guān)關(guān)系，并通過(guò)它們的關(guān)系從人均國(guó)內(nèi)生產(chǎn)總值(可以看作反映生活水平的一個(gè)指標(biāo))、體質(zhì)得分預(yù)測(cè)其壽命可能的變化范圍。體質(zhì)是指人體的質(zhì)量，是遺傳性和獲得性的基礎(chǔ)上表現(xiàn)出

32、來(lái)的人體形態(tài)結(jié)構(gòu)，生理機(jī)能和心理因素綜合的、相對(duì)穩(wěn)定的特征。體質(zhì)是人的生命活動(dòng)和工作能力的物質(zhì)基礎(chǔ)。它在形成、發(fā)展和消亡過(guò)程中，具有明顯的個(gè)體差異和階段性。中國(guó)體育科學(xué)學(xué)會(huì)體質(zhì)研究會(huì)研究表明，體質(zhì)應(yīng)包括身體形態(tài)發(fā)育水平、生理功能水平、身體素質(zhì)和運(yùn)動(dòng)能力發(fā)展水平、心理發(fā)育水平和適應(yīng)能力等五個(gè)方面。目前，體質(zhì)的綜合評(píng)價(jià)主要是形態(tài)、機(jī)能和身體素質(zhì)三類指標(biāo)按一定的權(quán)重進(jìn)行換算而得。表8-431個(gè)省市人口預(yù)期壽命與人均國(guó)內(nèi)生產(chǎn)總值和體質(zhì)得分?jǐn)?shù)據(jù)序號(hào)預(yù)期壽命體質(zhì)得分人均產(chǎn)值序號(hào)預(yù)期壽命體質(zhì)得分人均產(chǎn)值序號(hào)預(yù)期壽命體質(zhì)得分人均產(chǎn)值171.5466.165128571265.4956.77587442369.

33、8764.30517717273.9271.25244951368.9566.01114942467.4160.48515205373.2770.135242501473.3467.97204612578.1470.2970622471.2065.125100601565.9662.953822676.1069.34547319573.9169.99299311672.3766.1190702774.9168.41540643672.5465.765182431770.0764.51109352872.9166.49511781770.6667.29107631872.5568.3852200

34、72970.1765.76510658871.8567.7199071971.6566.205135943066.0363.2811587971.0866.525132552071.73,65.77114743164.3762.8497251071.29,67.1390882173.1067.065143351174.7069.505337722267.4763.6057898模型的建立和求解作表8-4數(shù)據(jù)（x1,y）,（x2,y）的散點(diǎn)圖如圖8.3圖8.3預(yù)期壽命與人均國(guó)內(nèi)生產(chǎn)總值和體質(zhì)得分的散點(diǎn)圖從圖8.3可以看出人口預(yù)期壽命y與體質(zhì)得分X2有較好的線性關(guān)系，y與人均國(guó)內(nèi)生產(chǎn)總值Xi的關(guān)系

35、難以確定，我們建立二次函數(shù)的回歸模型。一般的多元二項(xiàng)式回歸模型可表為y=0：iX1|mXm-、:jkXj%.；1<j,kimMATLAB統(tǒng)計(jì)工具箱提供了一個(gè)很方便的多元二項(xiàng)式回歸命令:Rstool（X,y,'model',alpha）,y為因變量（n維向量），alpha為顯著水平，model從下輸入x為自變量（nxm矩陣）列4個(gè)模型中選擇一個(gè):linear（只包含線性項(xiàng)）purequadratic（包含線性項(xiàng)和純二次項(xiàng)）interaction（包含線性項(xiàng)和純交互項(xiàng)）quadratic（包含線性項(xiàng)和完全二次項(xiàng)）輸出一個(gè)交互式畫面，對(duì)例3,編程如下：y=71.5473.927

36、3.2771.2073.9172.5470.6671.8571.0871.29,74.7065.4968.9573.3465.9672.3770.0772.5571.6571.73,73.1067.4769.8767.4178.1476.1074.9172.9170.1766.0364.37;x1=12857244952425010060299311824310763990713255908833772874411494204615382190701093522007135941147414335789817717152057062247319406431178110658115879725;

37、x2=66.16571.2570.13565.12569.9965.76567.2967.7166.52567.13,69.50556.77566.0167.9762.966.164.5168.38566.20565.77,67.06563.60564.30560.48570.2969.34568.41566.49565.76563.2862.84;x=x1',x2'rstool（x,y','purequadratic'）得到一個(gè)如圖8.4的交互式畫面圖8.4預(yù)期壽命與人均國(guó)內(nèi)生產(chǎn)總值和體質(zhì)得分的一個(gè)交互式畫面左邊一幅圖形是x2固定時(shí)的曲線y(x1)及其

38、置信區(qū)間，右邊一幅圖形是用固定時(shí)的曲線y(x2)及其置信區(qū)間。移動(dòng)鼠標(biāo)可改變,x2的值，同時(shí)圖左邊給出y的預(yù)測(cè)值及其置信區(qū)間。如輸入Xi=128757,x2=66.165,則y=70.6948,其置信區(qū)間70.6948±1.1079。圖的左下方有兩個(gè)下拉式菜單，上面的菜單Export用于輸出數(shù)據(jù)(包括：回歸系數(shù)parameters殘差residuals,剩余標(biāo)準(zhǔn)差RMSE等)，在MATLAB工作空間中得到有關(guān)數(shù)據(jù)。通過(guò)下面的菜單在上述4個(gè)模型中變更選擇，最后確定RMSE值較小的模型。例3則是包含線性項(xiàng)和完全二次項(xiàng)(quadratic)的模型最佳，即2-：-yY-yy：Y21Y2y01

39、Xi2X23X1X24X15X2-剩余標(biāo)準(zhǔn)差為1.2622,因此，所得回歸模型為：一一一一一一一一一_5-922y=195.360.0045x1-5.5753x2-6.733810X1X23.352910Xi0.055556x2利用此模型我們可以根據(jù)國(guó)內(nèi)生產(chǎn)總值及體質(zhì)得分，預(yù)測(cè)壽命。8.3非線性回歸分析8.3.1 非線性最小二乘擬合線性最小二乘擬合與線性回歸中的“線性”并非指y與x的關(guān)系，而是指y是系數(shù)0°,01或P=(Po,凱川，以)的線性函數(shù)。擬合如y=Po+PiX2的函數(shù)仍然是最小二乘擬合；如果擬合如y=P°e跳的曲線，y對(duì)P0,P1是非線性的，但取對(duì)數(shù)后Iny對(duì)系數(shù)

40、P0,P1是線性的，屬于可化為線性回歸的類型。下面討論非線性擬合的情形。非線性最小二乘擬合問(wèn)題的提法是：已知模型y=f(X,P),X=(Xi,|,Xm),P=(Po,Pi|,Pk),其中f對(duì)P是非線性的，為了估計(jì)參數(shù)P,收集n個(gè)獨(dú)立觀測(cè)數(shù)據(jù)(Xi,y)Xi=(Xii,|Xim)(i=l,|,n),n>m。記擬合誤差身(口)=yi-f(x,口)，求口使誤差的平方和nnQL)=2L)="yi-f(X)2i1i1最小。作為無(wú)約束非線性規(guī)劃的特例，解非線性最小二乘擬合可用MATLAB優(yōu)化工具箱命令I(lǐng)sqnonlin和Isqcurvefit。8.3.2 非線性回歸模型非線性回歸模型記作y

41、=f(X；)；,X=(Xi,I兒Xm)；=(-0,-1,HI,-k)其中f對(duì)回歸系數(shù)P是非線性的，名N(0,。2)。求得回歸系數(shù)P的最小二乘估計(jì)。MATLAB統(tǒng)計(jì)工具箱中非線性回歸的命令是：b,R,J=nIinfit(X,y,'model',bo)輸入x是自變量數(shù)據(jù)矩陣，每列一個(gè)向量；y是因變量數(shù)據(jù)向量；model是模型的函數(shù)名(M文件)，形式為y=f(b,X),b為待估系數(shù)P；b0是回歸系數(shù)P的初值。輸出b是P的估計(jì)值，R是殘差，J是用于估計(jì)預(yù)測(cè)誤差的Jacobi矩陣。這個(gè)命令是依據(jù)高斯一牛頓法求解的。將上面的輸出作為命令Bi=nIparci(b,R,J)的輸入，得到的bi是

42、回歸系數(shù)P的置信區(qū)間。用命令nIintooI(X,y,'model',b)可以得到一個(gè)交互式畫面，其內(nèi)容和用法與多項(xiàng)式回歸的PoIytooI類似。例4酶促反應(yīng)速度與底物濃度酶促反應(yīng)動(dòng)力學(xué)簡(jiǎn)稱酶動(dòng)力學(xué)，主要研究酶促反應(yīng)速度與底物(即反應(yīng)物)濃度以及其它因素的關(guān)系。在底物濃度很低時(shí)酶促反應(yīng)是一級(jí)反應(yīng)；當(dāng)?shù)孜餄舛忍幱谥虚g范圍時(shí)，是混合級(jí)反應(yīng)；當(dāng)?shù)孜餄舛仍黾訒r(shí)，向零級(jí)反應(yīng)過(guò)渡。某生化系學(xué)生為了研究喋吟霉素在某項(xiàng)酶促反應(yīng)中對(duì)反應(yīng)速度與底物濃度之間關(guān)系的影響，設(shè)計(jì)了兩個(gè)實(shí)驗(yàn)，一個(gè)實(shí)驗(yàn)中所使用的酶是經(jīng)過(guò)喋吟霉素處理的，而另一個(gè)實(shí)驗(yàn)所用的酶是未經(jīng)喋吟霉素處理的。所得實(shí)驗(yàn)數(shù)據(jù)見(jiàn)表8-5。試根據(jù)問(wèn)

43、題的背景和這些數(shù)據(jù)建立一個(gè)合適的數(shù)學(xué)模型，來(lái)反映這項(xiàng)酶促反應(yīng)的速度與底物濃度以及喋吟霉素處理與否之間的關(guān)系。表8-5喋吟霉素實(shí)驗(yàn)中的反應(yīng)速度與底物濃度數(shù)據(jù)底物濃度(ppm)0.020.060.110.220.561.10反應(yīng)速度未處理6751848698115131124144158160/處理764797107123139159152191201207200分析與假設(shè)記酶促反應(yīng)的速度為y,底物濃度為x,二者之間的關(guān)系寫作y=f(x,P),其中3為參數(shù)(3可為一向量)。由酶促反應(yīng)的基本性質(zhì)可知，當(dāng)?shù)孜餄舛群艿蜁r(shí)酶促反應(yīng)是一級(jí)反應(yīng)，此時(shí)反應(yīng)速度大致與底物濃度成正比；而當(dāng)?shù)孜餄舛群艽?，漸近飽和時(shí)，

44、反應(yīng)速度將趨于一個(gè)固定值(即零級(jí)反應(yīng))。下面的兩個(gè)簡(jiǎn)單模型具有這種性質(zhì)：Michaelis-Menten模型y=f(x，？=1x指數(shù)增長(zhǎng)模型y=f(x,P)=P(1e-吃)非線性模型的求解首先作出給出的經(jīng)過(guò)喋吟霉素處理和未經(jīng)處理的反應(yīng)速度與底物濃度的散點(diǎn)圖，可以看出，上述兩個(gè)模型與實(shí)際數(shù)據(jù)得到的散點(diǎn)圖是大致符合的。我們將主要對(duì)前一模型即Michaelis-Menten模型進(jìn)行詳細(xì)的分析。首先對(duì)經(jīng)過(guò)喋吟酶素處理的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行分析，在此基礎(chǔ)上，再來(lái)討論是否有更一般的模型來(lái)統(tǒng)一刻畫處理前后的數(shù)據(jù)，進(jìn)而揭示其中的聯(lián)系。我們用非線性回歸的方法直接估計(jì)模型的參數(shù)P1,P2,模型的求解可利用MATLAB統(tǒng)計(jì)

45、工具箱中的命令進(jìn)行，使用格式為：beta,R,J=nlinfit(x,y,'model',beta0)其中輸入x為自變量數(shù)據(jù)矩陣，每列一個(gè)變量；y為因變量數(shù)據(jù)向量；model為模型的M文件名，M函數(shù)形式為y=f(beta,x),beta為待估計(jì)參數(shù)；beta。為給定的參數(shù)初值。輸出beta為參數(shù)估計(jì)值，R為殘差，J為用于估計(jì)預(yù)測(cè)誤差的Jacobi矩陣。參數(shù)beta的置信區(qū)間用命令nlparci(beta,R,J)得到。首先建立函數(shù)M文件huaxue.m,非線性模型參數(shù)估計(jì)的源程序如下：x=0.020.020.060.060.110.110.220.220.560.561.101

46、.10;y=764797107123139159152191201207200;beta0=195.80270.04841;beta,R,J=nlinfit(x,y,'huaxue',beta0);betaci=nlparci(beta,R,J);beta,betaciyy=beta(1)*x./(beta(2)+x);plot(x,y,'o',x,yy,'m+'),pausenlintool(x,y,'huaxue',beta)得到的數(shù)值結(jié)果見(jiàn)表8-6。Nlintool用于給出一個(gè)交互式畫面，可以得到因變量y的預(yù)測(cè)值和預(yù)測(cè)區(qū)間，左下方的Export可向工作區(qū)傳送剩余標(biāo)準(zhǔn)差等數(shù)據(jù)。表8-6模型參數(shù)的估計(jì)結(jié)果參數(shù)參數(shù)估計(jì)值置信區(qū)間p1212.6818197.2028228.1608P20.06410.04570.0826從上面的結(jié)果可以知道，對(duì)經(jīng)過(guò)喋吟霉素處理的實(shí)