回歸分析方法

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

2、血壓之間沒有確定的數(shù)量關(guān)系，人的年齡和血壓之間的關(guān)系就是相關(guān)關(guān)系?；貧w分析就是處理變量之間的相關(guān)關(guān)系的一種數(shù)學(xué)方法。其解決問題的大致方法、步驟如下：（ 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ì)算語言編程也要占用大量時(shí)

3、間，而對(duì)于經(jīng)濟(jì)管理及社會(huì)學(xué)等對(duì)高級(jí)編程語言了解不深的人來說要應(yīng)用這些統(tǒng)計(jì)方法更是不可能。 MATLAB 等軟件的開發(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)大的圖形功能使得概念、過程和結(jié)果可以直觀地展現(xiàn)在我們面前。本章內(nèi)容通常先介紹有關(guān)回歸分析的數(shù)學(xué)原理，主要說明建模過程中要做的工作及理由，如模型的假設(shè)檢驗(yàn)、參數(shù)估計(jì)等，為了把主要精力集中在應(yīng)用上，我們略去詳細(xì)而繁雜的理論。在此基礎(chǔ)上再介紹

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

5、中的參數(shù)則是擬合問題，可以使用MATLAB軟件的 curvefit命令或 nlinfit命令擬合得到參數(shù)的估計(jì)并進(jìn)行統(tǒng)計(jì)分析。本節(jié)主要考察線性回歸模型。8.1.1 一元線性回歸模型的建立及其MATLAB實(shí)現(xiàn)其中0，1 是待定系數(shù)，對(duì)于不同的x, y是相互獨(dú)立的隨機(jī)變量。假設(shè)對(duì)于 x 的 n 個(gè)值 xi，得到y(tǒng) 的 n 個(gè)相應(yīng)的值yi ，確定0，1 的方法是根據(jù)最小二乘準(zhǔn)則，要使取最小值。利用極值必要條件令Q0, Q0，求0，1 的估計(jì)值?， ? ，從而得到回歸直線0101y?01x 。只不過這個(gè)過程可以由軟件通過直線擬合完成，而無須進(jìn)行繁雜的運(yùn)算。（ 1）參數(shù)的區(qū)間估計(jì)由于我們所計(jì)算出的?0，

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

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

8、，則認(rèn)為變量 y 與 x 顯著地有線性關(guān)系，其中F1 (1,n 2) 的值可查 F 分布表，或直接用 MATLAB命令 finv(1-,1, n-2) 計(jì)算得到；如果 p表示線性模型可用。這三個(gè)值可以相互印證。s2的值主要用來比較模型是否有改進(jìn)，其值越小說明模型精度越高。例 1 測(cè)得 16 名成年女子身高 y 與腿長(zhǎng) x 所得數(shù)據(jù)如下：表 8-1 16 名女子身高 (cm) 腿長(zhǎng) (cm) 數(shù)據(jù)8885889192939395969897969899100102143145146147 149150153154155156157158159160162164首先利用命令 plot(x,y,&#

9、39;r*') 畫出散點(diǎn)圖，從圖形可以看出，這些點(diǎn)大致分布在一條直線的左右，因此，可以考慮一元線性回歸?？删幹瞥绦蛉缦拢簓=143 145 146147149150153154155156157158159160162164;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

10、.9282180.95310.00003.1277R2 =0.9282 ，由 finv(0.95,1,14)= 4.6001，即 F1(1,n2) = 4.6001<F=180.9531 ， p<0.0001 ，可以通過殘差圖發(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 ，即 F1(1,n2) = 4.6672<F=261.6389 ， p&l

11、t;0.0001 ，說明模型有效且有改進(jìn)，因此我們得到身高與腿長(zhǎng)的關(guān)系y 17.65491.4363 x 。當(dāng)然，也可以利用直線擬合得到同一方程。只不過不能得到參數(shù)置信區(qū)間和對(duì)模型進(jìn)行檢驗(yàn)。擬合程序如下：y=143 145 146147149150153154155156157158159160162164;x=8885889192939395969897969899100102;a=polyfit(x,y,1)temp=polyval(a,x);plot(x,y,'r*',x,temp)注意：函數(shù)相同,但輸出一次函數(shù)參數(shù)順序與回歸分析（升冪排列）中不同。另一個(gè)差別是擬合不能發(fā)

12、現(xiàn)奇異數(shù)據(jù)。8.2多元線性回歸分析8.2.1多元線性回歸模型的建模步驟及其MATLAB實(shí)現(xiàn)如果根據(jù)經(jīng)驗(yàn)和有關(guān)知識(shí)認(rèn)為與因變量有關(guān)聯(lián)的自變量不止一個(gè)，那么就應(yīng)該考慮用最小二乘準(zhǔn)則建立多元線性回歸模型。設(shè)影響因變量y 的主要因素（自變量）有m 個(gè)，記x( x1 , xm )，假設(shè)它們有如下的線性關(guān)系式：y01 x1m xm，N(0,2)如果對(duì)變量y 與自變量x1 , x2 , xm同時(shí)作n 次觀察（ n>m）得n 組觀察值，采用最小二乘估計(jì)求得回歸方程y?0? x11? xkm.建立回歸模型是一個(gè)相當(dāng)復(fù)雜的過程，概括起來主要有以下幾個(gè)方面工作（1）根據(jù)研究目的收集數(shù)據(jù)和預(yù)分析；（2）根據(jù)散點(diǎn)

13、圖是否具有線性關(guān)系建立基本回歸模型；（3）模型的精細(xì)分析；（4）模型的確認(rèn)與應(yīng)用等。收集數(shù)據(jù)的一個(gè)經(jīng)驗(yàn)準(zhǔn)則是收集的數(shù)據(jù)量（樣本容量）至少應(yīng)為可能的自變量數(shù)目的610倍。在建模過程中首先要根據(jù)所研究問題的目的設(shè)置因變量，然后再選取與該因變量有統(tǒng)計(jì)關(guān)系的一些變量作為自變量。我們當(dāng)然希望選擇與問題關(guān)系密切的變量，同時(shí)這些變量之間相關(guān)性不太強(qiáng)，這可以在得到初步的模型后利用MATLAB軟件進(jìn)行相關(guān)性檢驗(yàn)。下面通過一個(gè)案例探討MATLAB軟件在回歸分析建模各個(gè)環(huán)節(jié)中如何應(yīng)用。多元線性回歸的MATLAB實(shí)現(xiàn)仍然用命令regress(y , X)，只是要注意矩陣X 的形式，將通過如下例子說明其用法。8.2.2

14、某類研究學(xué)者的年薪1. 問題例 2 工薪階層關(guān)心年薪與哪些因素有關(guān)，以此可制定出它們自己的奮斗目標(biāo)。某科學(xué)基金會(huì)希望估計(jì)從事某研究的學(xué)者的年薪Y(jié) 與他們的研究成果 (論文、著作等 )的質(zhì)量指標(biāo) X1、從事研究工作的時(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ù)i1234567891011123.55.35.15.84.26.06.85.53.17.24.54.992018333113253054725116.16.47.46.77.55.96.04.05.88.35.06

15、.433.240.338.746.841.437.539.040.730.152.938.231.8i1314151617181920212223248.06.56.63.76.27.04.04.55.95.64.83.9233539217403523332734157.67.05.04.45.57.06.03.54.94.38.05.843.344.142.533.634.248.038.035.940.436.845.235.1試建立Y 與 X1 , X 2 , X 3 之間關(guān)系的數(shù)學(xué)模型，并得出有關(guān)結(jié)論和作統(tǒng)計(jì)分析。2. 作出因變量 Y 與各自變量的樣本散點(diǎn)圖作散點(diǎn)圖的目的主要是觀察因變

16、量Y 與各自變量間是否有比較好的線性關(guān)系，以便選擇恰當(dāng)?shù)臄?shù)學(xué)模型形式。下圖分別為年薪 Y 與成果質(zhì)量指標(biāo)X1 、研究工作時(shí)間 X 2、獲得資助的指標(biāo)X3 之間的散點(diǎn)圖，subplot(1,3,1),plot(x1,Y,'g*'),subplot(1,3,2),plot(x2,Y,'k+'),subplot(1,3,3),plot(x3,Y,'ro'),從圖可以看出這些點(diǎn)大致分布在一條直線旁邊，因此，有比較好的線性關(guān)系，可以采用線性回歸。Y 與 x1 的散點(diǎn)圖Y 與 x2 的散點(diǎn)圖Y 與 x3的散點(diǎn)圖圖 8.1 因變量 Y 與各自變量的樣本散點(diǎn)圖3

17、. 利用 MATLAB 統(tǒng)計(jì)工具箱得到初步的回歸方程設(shè)回歸方程為：?.01 x12 x33x3y建立 m-文件輸入如下程序數(shù)據(jù)：x1=3.5 5.3 5.1 5.8 4.2 6.0 6.8 5.5 3.1 7.2 4.5 4.9 8.0 6.5 6.5 3.7 6.2 7.0 4.0 4.5 5.9 5.6 4.8 3.9; x2=9 20 18 33 31 13 25 30 5 47 25 11 23 35 39 21 7 40 35 23 33 27 34 15;x3=6.1 6.4 7.4 6.7 7.5 5.9 6.0 4.0 5.8 8.3 5.0 6.4 7.6 7.0 5.0 4

18、.0 5.5 7.0 6.0 3.5 4.9 4.3 8.0 5.0;Y=33.2 40.3 38.7 46.8 41.4 37.5 39.0 40.7 30.1 52.9 38.2 31.8 43.3 44.1 42.5 33.6 34.2 48.0 38.0 35.9 40.4 36.8 45.2 35.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ù)的

19、估計(jì)值回歸系數(shù)的置信區(qū)間18.015713.905222.12621.08170.39001.77330.32120.24400.39841.28350.66911.8979R2 =0.9106F=67.9195p<0.0001s2 = 3.0719計(jì)算結(jié)果包括回歸系數(shù)b=(0 ,1 ,2 , 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ù)的平方R2 ，假設(shè)檢驗(yàn)統(tǒng)計(jì)量，與 F 對(duì)應(yīng)的概率p， s2的值（ 7.0以前版本 s2也可由程序sum(r.2)/(n-m-

20、1) 計(jì)算） 。因此我們得到初步的回歸方程為：由結(jié)果對(duì)模型的判斷：回歸系數(shù)置信區(qū)間不包含零點(diǎn)表示模型較好，殘差在零點(diǎn)附近也表示模型較好，接著就是利用檢驗(yàn)統(tǒng)計(jì)量， p 的值判斷該模型是否可用。（） 相關(guān)系數(shù)的評(píng)價(jià)：一般地，相關(guān)系數(shù)絕對(duì)值在0.81范圍內(nèi)，可判斷回歸自變量與因變量具有較強(qiáng)的線性相關(guān)性。本例的絕對(duì)值為0.9542，表明線性相關(guān)性較強(qiáng)。（） F 檢驗(yàn)法：當(dāng)FF1(m, nm 1) ，即認(rèn)為因變量 y 與自變量 x1 , x2 , , xm 之間顯著地有線性相關(guān)關(guān) 系 ； 否 則 認(rèn) 為 因 變 量 y 與 自 變 量 x1 , x2 , xm 之 間 線 性 相 關(guān) 關(guān) 系 不 顯 著

21、 。 本 例 67.919> F1 0.05 (3,20) = 3.10（查 F 分布表或輸入命令finv(0.95,3,20) 計(jì)算）。（） p 值檢驗(yàn)：若p（為預(yù)定顯著水平），則說明因變量y 與自變量 x1 , x2 , xm 之間顯著地有線性相關(guān)關(guān)系。本例輸出結(jié)果，p<0.0001, 顯然滿足 P<=0.05 。以上三種統(tǒng)計(jì)推斷方法推斷的結(jié)果是一致的，說明因變量y 與自變量之間顯著地有線性相關(guān)關(guān)系，所得線性回歸模型可用。s2 當(dāng)然越小越好，這主要在模型改進(jìn)時(shí)作為參考。4. 模型的精細(xì)分析和改進(jìn)（ 1） 殘差分析殘差eiyi yi(i1,2, ,n) ，是各觀測(cè)值yi 與

22、回歸方程所對(duì)應(yīng)得到的擬合值y?i 之差，實(shí)際上，它是線性回歸模型中誤差的估計(jì)值。 N(0,2 ) 即有零均值和常值方差，利用殘差的這種特性反過來考察原模型的合理性就是殘差分析的基本思想。利用MATLAB 進(jìn)行殘差分析則是通過殘差圖或時(shí)序殘差圖。殘差圖是指以殘差為縱坐標(biāo)，以其他指定的量為橫坐標(biāo)的散點(diǎn)圖。主要包括：（ 1）橫坐標(biāo)為觀測(cè)時(shí)間或觀測(cè)值序號(hào)；（2）橫坐標(biāo)為某個(gè)自變量的觀測(cè)值；（3）橫坐標(biāo)為因變量的擬合值。通過觀察殘差圖，可以對(duì)奇異點(diǎn)進(jìn)行分析，還可以對(duì)誤差的等方差性以及對(duì)回歸函數(shù)中是否包含其他自變量、自變量的高次項(xiàng)及交叉項(xiàng)等問題給出直觀的檢驗(yàn)。以觀測(cè)值序號(hào)為橫坐標(biāo)，殘差為縱坐標(biāo)所得到的散點(diǎn)

23、圖稱為時(shí)序殘差圖，畫出時(shí)序殘差圖的MATLAB語句為 rcoplot(r,rint) （圖 8.2）。可以清楚看到殘差大都分布在零的附近，因此還是比較好的，不過第 4、 12 、 19 這三個(gè)樣本點(diǎn)的殘差偏離原點(diǎn)較遠(yuǎn)，如果作為奇異點(diǎn)看待，去掉后重新擬合，則得回歸模型為：且回歸系數(shù)的置信區(qū)間更小均不包含原點(diǎn)，統(tǒng)計(jì)變量stats 包含的三個(gè)檢驗(yàn)統(tǒng)計(jì)量：相關(guān)系數(shù)的平方R2，假設(shè)檢驗(yàn)統(tǒng)計(jì)量，概率，分別為：0.9533；115.5586；0.0000，比較可知R， F 均增加模型得到改進(jìn)。圖 8.2時(shí)序殘差圖（ 2） 變量間的交互作用討論變量間的交互作用包括：不同自變量之間的交互作用以及同一變量的自相關(guān)

24、性。不同自變量之間的交互作用 ：有時(shí)，在實(shí)驗(yàn)中不僅單因素對(duì)指標(biāo)有影響，而且因素間還會(huì)聯(lián)合起來對(duì)指標(biāo)產(chǎn)生影響，常稱這種聯(lián)合作用為交互作用。處理兩個(gè)因素間交互作用的一個(gè)簡(jiǎn)單辦法是加入這兩個(gè)自變量的乘積項(xiàng)。本文案例如果加入交互項(xiàng)則為：用表 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ù)中

25、，同一變量的順序觀測(cè)值之間出現(xiàn)的相關(guān)現(xiàn)象稱為自相關(guān)。一旦數(shù)據(jù)中存在這種自相關(guān)序列，如果仍采用普通的回歸模型直接處理，將產(chǎn)生不良后果，使預(yù)測(cè)失去意義。自相關(guān)的診斷主要有圖示檢驗(yàn)法、相關(guān)系數(shù)法和DW 檢驗(yàn)法。圖示檢驗(yàn)法是通過繪制殘差et 散點(diǎn)圖觀察，如果散布點(diǎn)(et 1, et ) , t2,3, ,n 大部分點(diǎn)落在第，象限，表明存在著正的序列相關(guān)；如果大部分點(diǎn)落在第，象限，表明存在著負(fù)的序列相關(guān)。對(duì)DW 檢驗(yàn)法可以利用 MATLAB 軟件編程計(jì)算統(tǒng)計(jì)量：nee?tt 1DWt 22(1 ),nn，et2et21t2t 2然后查閱DW 檢驗(yàn)上下界表，以決定模型的自相關(guān)狀態(tài)。當(dāng)一個(gè)回歸模型存在序列相

26、關(guān)性時(shí)，首先要查明序列相關(guān)產(chǎn)生的原因。如果是回歸模型選用不當(dāng)，則應(yīng)改用適當(dāng)?shù)幕貧w模型；如果是缺少重要的自變量，則應(yīng)增加自變量；如果以上方法都不能消除序列相關(guān)性，則需要采用差分法、迭代法等處理，更詳細(xì)內(nèi)容參見相關(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ì)工具箱中逐步回歸的命令是s

27、tepwise，它提供了一個(gè)人機(jī)交互式畫面，通過此工具可以自由地選擇變量進(jìn)行統(tǒng)計(jì)分析。該命令的用法是：stepwise(X , Y , inmodel , alpha)其中成 nX 是自變量數(shù)據(jù)，排成nm 矩陣（ 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)生圖形窗口：Stepwise Plot , Stepwise Table , Stepwise History.當(dāng)鼠標(biāo)移到圖形某個(gè)區(qū)域時(shí)，鼠標(biāo)點(diǎn)擊后產(chǎn)生交互作用。Stepwise

28、 Plot 窗口中的虛線表示回歸系數(shù)的置信區(qū)間包含零點(diǎn)，即該回歸系數(shù)與零無顯著差異，一般應(yīng)將該變量移去；實(shí)線則表明該回歸系數(shù)與零有顯著差異，應(yīng)保留在模型中(藍(lán)色表示該變量已進(jìn)入模型，紅色表示該變量已移出模型)。引入和移出變量還可參考Stepwise History 窗口中剩余標(biāo)準(zhǔn)差RMSE 是否在下降，剩余標(biāo)準(zhǔn)差RMSE 最小的就是最好的模型。Stepwise Table 窗口中列出了一個(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，如果引入新的自變量x4x1x2 , x5x1 x3 , x6x2

29、x3 . 也可以采用逐步回歸法解決,源程序如下：A=3.5 5.3 5.1 5.8 4.2 6.0 6.8 5.5 3.1 7.2 4.5 4.9 8.0 6.5 6.5 3.7 6.2 7.0 4.0 4.5 5.9 5.6 4.8 3.9;9 20 18 33 31 13 25 30 5 47 25 11 23 35 39 21 7 40 35 23 33 27 34 15;6.1 6.4 7.4 6.7 7.5 5.9 6.0 4.0 5.8 8.3 5.0 6.4 7.6 7.0 5.0 4.0 5.5 7.0 6.0 3.5 4.9 4.3 8.0 5.0'Y=33.2 40

30、.3 38.7 46.8 41.4 37.5 39.0 40.7 30.1 52.9 38.2 31.8 43.3 44.1 42.5 33.6 34.2 48.0 38.0 35.9 40.4 36.8 45.2 35.1'x1=A(:,1);x2=A(:,2);x3=A(:,3);x4=x1.*x2;x5=x1.*x3;x6=x2.*x3;X=A,x4,x5,x6;stepwise(X,Y)運(yùn)行并按上述步驟操作后可以得到本文前面線性回歸相同的結(jié)論，即不含交互項(xiàng)的模型是最好的。在此只介紹操作過程，其交互界面，只要在MATLAB軟件上一試便知。8.2.4多項(xiàng)式回歸多項(xiàng)式回歸仍然屬于多元

31、線性回歸，可以是一元多項(xiàng)式回歸或多元多項(xiàng)式回歸。一元多項(xiàng)式回歸模型的一般形式為用 MATLAB求解一元多項(xiàng)式回歸，除了使用命令polyfit(x,y,m)外，還可以使用如下命令：Polytool(x,y,m,alpha)輸入 x,y,m 同命令 polyfit ，alpha 是顯著性水平（默認(rèn)0.05 ），則輸出一個(gè)交互式畫面，畫面顯示回歸曲線及其置信區(qū)間，通過圖左下方的export 下拉式菜單，還可以輸出回歸系數(shù)估計(jì)值及其置信區(qū)間、殘差等。下面通過一個(gè)用多元多項(xiàng)式回歸的實(shí)例說明什么時(shí)候用多項(xiàng)式回歸以及如何通過MATLAB軟件進(jìn)行處理。例3為了了解人口平均預(yù)期壽命與人均國內(nèi)生產(chǎn)總值和體質(zhì)得分的

32、關(guān)系，我們查閱了國家統(tǒng)計(jì)局資料，北京體育大學(xué)出版社出版的2000 國民體質(zhì)監(jiān)測(cè)報(bào)告，表8-4 是我國大陸31 個(gè)省市的有關(guān)數(shù)據(jù)。我們希望通過這幾組數(shù)據(jù)考察它們是否具有良好的相關(guān)關(guān)系，并通過它們的關(guān)系從人均國內(nèi)生產(chǎn)總值（可以看作反映生活水平的一個(gè)指標(biāo)）、體質(zhì)得分預(yù)測(cè)其壽命可能的變化范圍。體質(zhì)是指人體的質(zhì)量，是遺傳性和獲得性的基礎(chǔ)上表現(xiàn)出來的人體形態(tài)結(jié)構(gòu)，生理機(jī)能和心理因素綜合的、相對(duì)穩(wěn)定的特征。體質(zhì)是人的生命活動(dòng)和工作能力的物質(zhì)基礎(chǔ)。它在形成、發(fā)展和消亡過程中，具有明顯的個(gè)體差異和階段性。中國體育科學(xué)學(xué)會(huì)體質(zhì)研究會(huì)研究表明，體質(zhì)應(yīng)包括身體形態(tài)發(fā)育水平、生理功能水平、身體素質(zhì)和運(yùn)動(dòng)能力發(fā)展水平、心

33、理發(fā)育水平和適應(yīng)能力等五個(gè)方面。目前，體質(zhì)的綜合評(píng)價(jià)主要是形態(tài)、機(jī)能和身體素質(zhì)三類指標(biāo)按一定的權(quán)重進(jìn)行換算而得。表 8-4 31個(gè)省市人口預(yù)期壽命與人均國內(nèi)生產(chǎn)總值和體質(zhì)得分?jǐn)?shù)據(jù)序預(yù)期體質(zhì)得分人均序預(yù)期體質(zhì)得人均序預(yù)期體質(zhì)得人均號(hào)壽命產(chǎn)值號(hào)壽命分產(chǎn)值號(hào)壽命分產(chǎn)值171.5466.165128571265.4956.77587442369.8764.30517717273.9271.25244951368.9566.01114942467.4160.48515205373.2770.135242501473.3467.97204612578.1470.2970622471.2065.125100

34、601565.9662.953822676.1069.34547319573.9169.99299311672.3766.1190702774.9168.41540643672.5465.765182431770.0764.51109352872.9166.49511781770.6667.29107631872.5568.385220072970.1765.76510658871.8567.7199071971.6566.205135943066.0363.2811587971.0866.525132552071.73,65.77114743164.3762.8497251071.29,67

35、.1390882173.1067.065143351174.7069 .505337722267.4763.6057898模型的建立和求解作表 8-4 數(shù)據(jù) (x1, y),( x2 , y) 的散點(diǎn)圖如圖 8.3圖 8.3預(yù)期壽命與人均國內(nèi)生產(chǎn)總值和體質(zhì)得分的散點(diǎn)圖從圖 8.3 可以看出人口預(yù)期壽命y 與體質(zhì)得分 x2 有較好的線性關(guān)系，y 與人均國內(nèi)生產(chǎn)總值x1 的關(guān)系難以確定，我們建立二次函數(shù)的回歸模型。一般的多元二項(xiàng)式回歸模型可表為MATLAB統(tǒng)計(jì)工具箱提供了一個(gè)很方便的多元二項(xiàng)式回歸命令：Rstool(x,y, 'model',alpha)輸入 x 為自變量（ n&#

36、215; m矩陣）， y 為因變量（ n 維向量），alpha 為顯著水平，model 從下列 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.54 73.92 73.27 71.20 73.91 72.54 70.66 71.85 71.08 71.29,74.70 65.49 68.95 73.34 65.96 72.37 70.07 72.5571.65 71.73,73.10 67.47 6

37、9.87 67.41 78.1476.10 74.91 72.91 70.17 66.03 64.37;x1=12857 2449524250 10060 29931 18243 10763 9907 13255 9088 33772 8744 11494 20461 5382 19070 1093522007 13594 11474 14335 7898 17717 15205 70622 47319 40643 11781 10658 11587 9725;x2=66.165 71.2570.135 65.125 69.99 65.765 67.29 67.71 66.525 67.13,

38、69.505 56.775 66.01 67.97 62.9 66.1 64.5168.385 66.205 65.77,67.065 63.605 64.30560.485 70.29 69.345 68.415 66.495 65.765 63.28 62.84;x=x1',x2'rstool(x,y','purequadratic')得到一個(gè)如圖8.4 的交互式畫面圖 8.4預(yù)期壽命與人均國內(nèi)生產(chǎn)總值和體質(zhì)得分的一個(gè)交互式畫面左邊一幅圖形是x2 固定時(shí)的曲線y( x1 ) 及其置信區(qū)間，右邊一幅圖形是x1 固定時(shí)的曲線y(x2 ) 及其置信區(qū)間。移

39、動(dòng)鼠標(biāo)可改變x1 ， x2 的值，同時(shí)圖左邊給出y 的預(yù)測(cè)值及其置信區(qū)間。如輸入x1 =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ù)。通過下面的菜單在上述4 個(gè)模型中變更選擇，最后確定RMSE 值較小的模型。例3 則是包含線性項(xiàng)和完全二次項(xiàng)（quadratic ）的模型最佳，即剩余標(biāo)準(zhǔn)差為1.2622，因此，所得回歸模型為：利用此模

40、型我們可以根據(jù)國內(nèi)生產(chǎn)總值及體質(zhì)得分，預(yù)測(cè)壽命。8.3非線性回歸分析8.3.1非線性最小二乘擬合線性 最 小二 乘 擬合 與 線 性回 歸 中的 “ 線性 ” 并 非指 y 與 x 的關(guān) 系 ，而 是 指 y 是系 數(shù) 0 , 1或( 0 , 1, , m ) 的 線 性 函 數(shù) 。 擬 合 如 y01x2 的 函 數(shù) 仍 然 是 最 小 二 乘 擬 合 ； 如 果 擬 合 如y0 e 1x 的曲線，y 對(duì) 0 ,1 是非線性的，但取對(duì)數(shù)后ln y 對(duì)系數(shù)0 ,1 是線性的，屬于可化為線性回歸的類型。下面討論非線性擬合的情形。非線性最小二乘擬合問題的提法是：已知模型yf ( x, ), x(x

41、1, , xm ) ,(0,1, k ) ，其中 f 對(duì)是非線性的，為了估計(jì)參數(shù)，收集 n 個(gè)獨(dú)立觀測(cè)數(shù)據(jù)( xi , yi ), xi(xi1,xim ) (i1, n), nm 。記擬合誤差i ( )yif ( xi , ) ，求使誤差的平方和最小。作為無約束非線性規(guī)劃的特例，解非線性最小二乘擬合可用MATLAB優(yōu)化工具箱命令lsqnonlin和lsqcurvefit 。8.3.2非線性回歸模型非線性回歸模型記作其中 f 對(duì)回歸系數(shù)是非線性的， N(0,2 ) 。求得回歸系數(shù)的最小二乘估計(jì)。MATLAB 統(tǒng)計(jì)工具箱中非線性回歸的命令是：b,R,J=nlinfit(x,y, 'mod

42、el',bo)輸入 x 是自變量數(shù)據(jù)矩陣，每列一個(gè)向量；y 是因變量數(shù)據(jù)向量；model 是模型的函數(shù)名（M 文件），形式為 yf (b, x) ， b 為待估系數(shù)； b0 是回歸系數(shù)的初值。輸出b 是的估計(jì)值， R 是殘差， J 是用于估計(jì)預(yù)測(cè)誤差的Jacobi 矩陣。這個(gè)命令是依據(jù)高斯牛頓法求解的。將上面的輸出作為命令Bi=nlparci(b,R,J)的輸入，得到的bi是回歸系數(shù)的置信區(qū)間。用命令nlintool(x,y, 'model',b)可以得到一個(gè)交互式畫面，其內(nèi)容和用法與多項(xiàng)式回歸的Polytool類似。例4酶促反應(yīng)速度與底物濃度酶促反應(yīng)動(dòng)力學(xué)簡(jiǎn)稱酶動(dòng)力學(xué)

43、，主要研究酶促反應(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)過渡。某生化系學(xué)生為了研究嘌呤霉素在某項(xiàng)酶促反應(yīng)中對(duì)反應(yīng)速度與底物濃度之間關(guān)系的影響，設(shè)計(jì)了兩個(gè)實(shí)驗(yàn)，一個(gè)實(shí)驗(yàn)中所使用的酶是經(jīng)過嘌呤霉素處理的，而另一個(gè)實(shí)驗(yàn)所用的酶是未經(jīng)嘌呤霉素處理的。所得實(shí)驗(yàn)數(shù)據(jù)見表8-5 。試根據(jù)問題的背景和這些數(shù)據(jù)建立一個(gè)合適的數(shù)學(xué)模型，來反映這項(xiàng)酶促反應(yīng)的速度與底物濃度以及嘌呤霉素處理與否之間的關(guān)系。表 8-5嘌呤霉素實(shí)驗(yàn)中的反應(yīng)速度與底物濃度數(shù)據(jù)底物濃度0.020.060.110.220.561.10

44、（ ppm ）反應(yīng)未處理6751848698115131124144158160/速度處理76 47 97107 123 139 159152 191201 207 200分析與假設(shè)記酶促反應(yīng)的速度為y ,底物濃度為x ，二者之間的關(guān)系寫作yf ( x,) ,其中 為參數(shù) ( 可為一向量 ) 。由酶促反應(yīng)的基本性質(zhì)可知，當(dāng)?shù)孜餄舛群艿蜁r(shí)酶促反應(yīng)是一級(jí)反應(yīng)，此時(shí)反應(yīng)速度大致與底物濃度成正比；而當(dāng)?shù)孜餄舛群艽?，漸近飽和時(shí)，反應(yīng)速度將趨于一個(gè)固定值（即零級(jí)反應(yīng)）。下面的兩個(gè)簡(jiǎn)單模型具有這種性質(zhì)：Michaelis-Menten 模型指數(shù)增長(zhǎng)模型非線性模型的求解首先作出給出的經(jīng)過嘌呤霉素處理和未經(jīng)處理

45、的反應(yīng)速度與底物濃度的散點(diǎn)圖，可以看出，上述兩個(gè)模型與實(shí)際數(shù)據(jù)得到的散點(diǎn)圖是大致符合的。我們將主要對(duì)前一模型即Michaelis-Menten模型進(jìn)行詳細(xì)的分析。首先對(duì)經(jīng)過嘌呤酶素處理的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行分析，在此基礎(chǔ)上，再來討論是否有更一般的模型來統(tǒng)一刻畫處理前后的數(shù)據(jù)，進(jìn)而揭示其中的聯(lián)系。我們用非線性回歸的方法直接估計(jì)模型的參數(shù)1， 2 ，模型的求解可利用MATLAB 統(tǒng)計(jì)工具箱中的命令進(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ù)；beta0 為給定的參數(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.02 0.02 0.06 0.06 0.11 0.11 0.22 0.22 0.56 0.56 1.10 1.10;y=76 47 97 107 123 139 159 152 191 201 207 20