回歸分析方法

1、第八章回歸分析方法當(dāng)人們對研究對象的內(nèi)在特性和各因素間的關(guān)系有比較充分的認(rèn)識時(shí)，一般用機(jī)理分析方法建立數(shù)學(xué)模型。 如果由于客觀事物內(nèi)部規(guī)律的復(fù)雜性及人們認(rèn)識程度的限制，無法分析實(shí)際對象內(nèi)在的因果關(guān)系，建立合乎 機(jī)理規(guī)律的數(shù)學(xué)模型，那么通常的辦法是搜集大量數(shù)據(jù)，基于對數(shù)據(jù)的統(tǒng)計(jì)分析去建立模型。本章討論其中 用途非常廣泛的一類模型一一統(tǒng)計(jì)回歸模型。回歸模型常用來解決預(yù)測、控制、生產(chǎn)工藝優(yōu)化等問題。變量之間的關(guān)系可以分為兩類：一類叫確定性關(guān)系，也叫函數(shù)關(guān)系，其特征是：一個變量隨著其它變量 的確定而確定。另一類關(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)選定因變量和自變量之間的模型，即一個數(shù)學(xué)式子，利用數(shù)據(jù)按照最小二乘準(zhǔn)則計(jì)算模型中的系 數(shù)；(3)利用統(tǒng)計(jì)分析方法對不同的模型進(jìn)行比較，找出與數(shù)據(jù)擬合得最好的模型；(4)判斷得到的模型是否適合于這組數(shù)據(jù)；(5)利用模型對因變量作出預(yù)測或解釋。應(yīng)用統(tǒng)計(jì)分析特別是多元統(tǒng)計(jì)分析方法一般都要處理大量數(shù)據(jù)，工作量非常大，所以在計(jì)算機(jī)普及以 前，這些方法大都是停留在理論研究上。運(yùn)用一般計(jì)算語言編程也要

3、占用大量時(shí)間，而對于經(jīng)濟(jì)管理及社會 學(xué)等對高級編程語言了解不深的人來說要應(yīng)用這些統(tǒng)計(jì)方法更是不可能。MATLAB 等軟件的開發(fā)和普及大大減少了對計(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ì)而繁雜的理論。在此基

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

5、函數(shù)形式只是要確定其中的參 數(shù)則是擬合問題，可以使用MATLAB 軟件的curvefit命令或nlinfit命令擬合得到參數(shù)的估計(jì)并進(jìn)行統(tǒng)計(jì)分析。本節(jié)主要考察線性回歸模型。8.1.1 一元線性回歸模型的建立及其MATLAB實(shí)現(xiàn)其中P0, P1是待定系數(shù)，又t于不同的x y是相互獨(dú)立的隨機(jī)變量。假設(shè)對于x的n個值Xi ,得到y(tǒng)的n個相應(yīng)的值 中，確定P0, P1的方法是根據(jù)最小二乘準(zhǔn)則，要使取最小值。利用極值必要條件令£Q = 0,£Q=0,求 % 3的估計(jì)值 % 因，從而得到回歸直線" 冷 0 1 0 1y = %+f?x。只不過這個過程可以由軟件通過直線擬合完成

6、，而無須進(jìn)行繁雜的運(yùn)算。(1)參數(shù)的區(qū)間估計(jì)由于我們所計(jì)算出的鼠 田仍然是隨機(jī)變量，因此要對 鼠 因取值的區(qū)間進(jìn)行估計(jì)，如果區(qū)間估計(jì)值是一個較短的區(qū)間表示模型精度較高。(2)對誤差方差的估計(jì)設(shè)?為回歸函數(shù)的值，X為測量值，殘差平方和剩余方差s2 = -Q-n -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 ,ri

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

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

9、t(x,y,'r*')畫出散點(diǎn)圖，從圖形可以看出，這些點(diǎn)大致分布在一條直線的左右，因此，可 以考慮一元線性回歸。可編制程序如下：y=143 145 146147 149 150 153 154 155 156 157 158 159 160 162 164;x=8885 8891929393959698 97969899 100 102;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.3196

10、51.22291.08461.4960s = 0.9282 180.95310.00003.12772R =0.9282,由 finv(0.95,1,14)= 4.6001 ,即 8心(1-2)= 4.6001<F=180.9531 , p<0.0001 , 可以通過殘差圖發(fā)現(xiàn)，第二個數(shù)據(jù)為奇異數(shù)據(jù)，去掉該數(shù)據(jù)后運(yùn)行后得到 b = 17.65491.4363bint = -0.598635.90831.24451.6281s = 0.9527 261.63890.00001.93132R =0.9527,由 finv(0.95,1,13)= 4.6672,即 F1q(1,n2)=

11、4.6672<F=261.6389 , p<0.0001 ,說明模型有效且有改進(jìn)，因此我們得到身高與腿長的關(guān)系y =17.6549 +1.4363x。當(dāng)然，也可以利用直線擬合得到同一方程。只不過不能得到參數(shù)置信區(qū)間和對模型進(jìn)行檢驗(yàn)。擬合程 序如下：y=143 145 146147149150153154155156 157 158 159 160 162 164;x=8885 8891929393959698 97 969899100 102;a=polyfit(x,y,1)temp=polyval(a,x);plot(x,y,'r*',x,temp)注意：函數(shù)相同

12、，但輸出一次函數(shù)參數(shù)順序與回歸分析(升募排列)中不同。另一個差別是擬合不能發(fā)現(xiàn)奇異 數(shù)據(jù)。8.2多元線性回歸分析8.2.1 多元線性回歸模型的建模步驟及其MATLAB實(shí)現(xiàn)如果根據(jù)經(jīng)驗(yàn)和有關(guān)知識認(rèn)為與因變量有關(guān)聯(lián)的自變量不止一個，那么就應(yīng)該考慮用最小二乘準(zhǔn)則建立 多元線性回歸模型。設(shè)影響因變量 y的主要因素(自變量)有 m個，記x = (x1,H|,xm),假設(shè)它們有如下的線性關(guān)系式：y = p0 + 01為+111 +Pm% +注， 名N(B2 )如果對變量 y與自變量x1,x2,H|,xm同時(shí)作n次觀察(n>m)得n組觀察值，采用最小二乘估計(jì)求得回 歸方程?= ?0?x "I

13、 Gm .建立回歸模型是一個相當(dāng)復(fù)雜的過程，概括起來主要有以下幾個方面工作(1)根據(jù)研究目的收集數(shù)據(jù)和預(yù)分析；(2)根據(jù)散點(diǎn)圖是否具有線性關(guān)系建立基本回歸模型；(3)模型的精細(xì)分析；(4)模型的確認(rèn)與應(yīng)用等。收集數(shù)據(jù)的一個經(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)。下面通過一個案例探討MATLAB軟件在回歸分析建模各個環(huán)節(jié)中如何應(yīng)用。 多

14、元線性回歸的 MATLAB 實(shí)現(xiàn)仍然用命令regress(y , X),只是要注意矩陣 X的形式，將通過如下例子說明其用法。8.2.2 某類研究學(xué)者的年薪1 .問題例2工薪階層關(guān)心年薪與哪些因素有關(guān)，以此可制定出它們自己的奮斗目標(biāo)。某科學(xué)基金會希望估計(jì)從事某研究的學(xué)者的年薪Y(jié)與他們的研究成果(論文、著作等)的質(zhì)量指標(biāo)Xi、從事研究工作的時(shí)間 X2、能成功獲得資助的指標(biāo) X3之間的關(guān)系，為此按一定的實(shí)驗(yàn)設(shè)計(jì)方法調(diào)查了24位研究學(xué)者，得到如下數(shù)據(jù)(i為學(xué)者序號)：表8-2從事某種研究的學(xué)者的相關(guān)指標(biāo)數(shù)據(jù)i1234567891011123.55.35.15.84.26.06.85.53.17.24.

15、54.992018333113253054725116.16.47.46.77.55.96.04.05.88.35.06.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試建立丫與Xi,X2,X3之間關(guān)系的數(shù)學(xué)模型

16、，并得出有關(guān)結(jié)論和作統(tǒng)計(jì)分析。2 .作出因變量 Y與各自變量的樣本散點(diǎn)圖作散點(diǎn)圖的目的主要是觀察因變量 丫與各自變量間是否有比較好的線性關(guān)系，以便選擇恰當(dāng)?shù)臄?shù)學(xué)模型形式。下圖分別為年薪丫與成果質(zhì)量指標(biāo)X1、研究工作時(shí)間 X2、獲得資助的指標(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)圖丫與x

17、2的散點(diǎn)圖丫與x3的散點(diǎn)圖圖8.1 因變量Y與各自變量的樣本散點(diǎn)圖3 .利用MATLAB統(tǒng)計(jì)工具箱得到初步的回歸方程設(shè)回歸方程為：/= ¥ 秋藐 ？3乂3.建立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

18、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.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.845.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

19、-3所示。表8-3對初步回歸模型的計(jì)算結(jié)果回歸系數(shù)回歸系數(shù)的估計(jì)值回歸系數(shù)的置信區(qū)間18.015713.905222.12621.08170.39001.77330.32120.2440 0.39841.28350.66911.897922R =0.9106F=67.9195p<0.0001 S = 3.0719計(jì)算結(jié)果包括回歸系數(shù)b=( P0 , P1, P2, P3)=(18.0157, 1.0817 , 0.3212 , 1.2835),且置信區(qū)間均不包含零點(diǎn)，；殘差及其置信區(qū)間；統(tǒng)計(jì)變量stats，它包含四個檢驗(yàn)統(tǒng)計(jì)量：相關(guān)系數(shù)的平方R2，假設(shè)檢驗(yàn)統(tǒng)計(jì)量22F,與F對應(yīng)的概率p,

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

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

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

23、縱坐標(biāo)所得到的散點(diǎn)圖稱為時(shí)序殘差圖，畫出時(shí)序殘差圖的MATLAB 語句為rcoplot(r,rint)(圖8.2)。可以清楚看到殘差大都分布在零的附近，因此還是比較好的，不過第4、12、19這三個樣本點(diǎn)的殘差偏離原點(diǎn)較遠(yuǎn)，如果作為奇異點(diǎn)看待，去掉后重新擬合，則得回歸模型 為：且回歸系數(shù)的置信區(qū)間更小均不包含原點(diǎn)，統(tǒng)計(jì)變量stats包含的三個檢驗(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) 變量間的交互作用討論變量間的交互作用包括：不同自變量之間的交互作用以

24、及同一變量的自相關(guān)性。不同自變量之間的交互作用：有時(shí)，在實(shí)驗(yàn)中不僅單因素對指標(biāo)有影響，而且因素間還會聯(lián)合起來對指標(biāo)產(chǎn)生影響，常稱這種聯(lián)合作用為交互作用。處理兩個因素間交互作用的一個簡單辦法是加入這兩個自變量 的乘積項(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)也不理想，因此，本例中 其交互作用并不顯著，該模型不如前面兩個模型好。自相關(guān)性的診斷和處理：若數(shù)據(jù)是以時(shí)間為序的，稱為時(shí)間

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

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 逐步回歸方法建模逐步回歸就是一種從眾多自變量中有效地選擇重要變量的方法。逐步回歸的基本思路是，先確定一個包 含若干自變量的初始集合，然后每次從集合外的變量中引入一個對因變量影響最大的，再對集合中的變量進(jìn) 行檢驗(yàn)，從變得不顯著的變量中移出一個影響最小的，依此進(jìn)行，直到不能引入和移出為止。引入和移出都 以給定的顯著性水平為標(biāo)準(zhǔn)。MATLAB 統(tǒng)計(jì)工具箱

27、中逐步回歸的命令是stepwise,它提供了一個人機(jī)交互式畫面，通過此工具可以自由地選擇變量進(jìn)行統(tǒng)計(jì)分析。該命令的用法是：stepwise(X , Y , inmodel , alpha)其中X是自變量數(shù)據(jù)，排成 n黑m矩陣(m為自變量個數(shù)，n為每個變量的數(shù)據(jù)量)，丫是因變量數(shù)據(jù)，排成門父1向量，inmodel是自變量初始集合的指標(biāo)，缺省時(shí)為全部自變量，alpha為顯著水平，缺省時(shí)為0.05。運(yùn)彳f stepwise命令時(shí)產(chǎn)生圖形窗口 ：Stepwise Plot , Stepwise Table , Stepwise History.當(dāng)鼠標(biāo)移到圖形某個區(qū)域時(shí)，鼠標(biāo)點(diǎn)擊后產(chǎn)生交互作用。Step

28、wise 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窗口中列出了一個統(tǒng)計(jì)表，包括回歸系數(shù)及其置信區(qū)間，以及模型的統(tǒng)計(jì)量剩余標(biāo)準(zhǔn)差RMSE、相關(guān)系數(shù) R-square、F值、與F對應(yīng)的概率。關(guān)于本節(jié)案例2,如果引入新的自變量x4= X1X2,X5=X1X3, X6= X2X

29、3,也可以采用逐步回歸法解決，源程序如下：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.3

30、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.845.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),則輸出一個交互式畫面，畫面顯示回歸曲線及 其置信區(qū)間，通過圖左下方的export下拉式菜單，還可以輸出回歸系數(shù)估計(jì)值及其置信區(qū)間、殘差等。下面通過一個用多元多項(xiàng)式回歸的實(shí)例說明什么時(shí)候用多項(xiàng)式回歸以及如何通過MATLAB軟件進(jìn)行處理。例3 為了了解人口平均預(yù)期壽命與人均國內(nèi)生產(chǎn)總值和體質(zhì)得分的關(guān)系，我

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

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

34、5100601565.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.2

35、9,67.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)總值 Xi的關(guān)系難以確定，我們建立二次函數(shù)的回歸模型。 一般的多元二項(xiàng)式回歸模型可表為MATLAB統(tǒng)計(jì)工具箱提供了一個很方便的多元二項(xiàng)式回歸命令：Rstool（x,y, 'model',alpha）輸入x為自變量（nxm矩陣），y為因變量（n

36、維向量），alpha為顯著水平，model從下列4個模型中 選擇一個：linear （只包含線性項(xiàng)）purequadratic （包含線性項(xiàng)和純二次項(xiàng)）interaction （包含線性項(xiàng)和純交互項(xiàng)）quadratic （包含線性項(xiàng)和完全二次項(xiàng)）輸出一個交互式畫面，對例 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 69.87 67.41 78.14 76.10

37、74.91 72.91 70.17 66.03 64.37;x1=12857 24495 24250 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.25 70.135 65.125 69.99 65.765 67.29 67.71 66.525 67.13,69.505 56.775 66.01

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

39、型中 的模型最佳，即間。移動鼠標(biāo)可改變xi, x2的值，同時(shí)圖左邊給出y的預(yù)測值及其置信區(qū)間。如輸入圖的左下方有兩個下拉式菜單，上面的菜單Export用于輸出數(shù)據(jù)(包括：回歸系數(shù)x2 =66.165,則 y =70.6948,其置信區(qū)間 70.6948 ± 1.1079 。residuals,剩余標(biāo)準(zhǔn)差 RMSE等)，在MATLAB 工作空間中得到有關(guān)數(shù)據(jù)。通過下面的菜單在上述 變更選擇，最后確定 RMSE值較小的模型。例 3則是包含線性項(xiàng)和完全二次項(xiàng)( quadratic) 剩余標(biāo)準(zhǔn)差為1.2622,因此，所得回歸模型為：利用此模型我們可以根據(jù)國內(nèi)生產(chǎn)總值及體質(zhì)得分，預(yù)測壽命。8.

40、3非線性回歸分析8.3.1 非線性最小二乘擬合線性最小二乘擬合與線性回歸中的“線性”并非指y與x的關(guān)系，而是指y是系數(shù)P0,P1或P =(Po, Pi ,111, Pm)的線性函數(shù)。擬合如y = Po + PlX2的函數(shù)仍然是最小二乘擬合；如果擬合如y = P°e做的曲線，y對P0,也是非線性的，但取對數(shù)后ln y對系數(shù)P0, Pi是線性的，屬于可化為線性回歸的類型。下面討論非線性擬合的情形。非線性最小二乘擬合問題的提法是：已知模型y = f (x, B), X = (Xi,|,Xm), B =(瓦,、1,用,九)，其中f對P是非線性的，為了估計(jì)參數(shù)P ,收集n個獨(dú)立觀測數(shù)據(jù)(Xi,

41、 y),Xi =(Xii,|xim) (i =1,|,n), n >m。記擬合誤差 剪件)=yi - f (x, F),求 P 使誤差的平方和最小。作為無約束非線性規(guī)劃的特例，解非線性最小二乘擬合可用MATLAB優(yōu)化工具箱命令I(lǐng)sqnonlin和Isqcurvefit。8.3.2 非線性回歸模型非線性回歸模型記作其中f對回歸系數(shù)P是非線性的，名N(0,仃2)。求得回歸系數(shù)P的最小二乘估計(jì)。MATLAB統(tǒng)計(jì)工具箱中非線性回歸的命令是：b,R,J=nlinfit(x,y, 'model',bo)輸入x是自變量數(shù)據(jù)矩陣，每列一個向量；y是因變量數(shù)據(jù)向量；model是模型的函數(shù)名

42、(M文件)，形式為y=f(b,x), b為待估系數(shù)P； b0是回歸系數(shù)P的初值。輸出b是P的估計(jì)值，R是殘差，J是用 于估計(jì)預(yù)測誤差的 Jacobi矩陣。這個命令是依據(jù)高斯一牛頓法求解的。將上面的輸出作為命令Bi=nlparci(b,R,J)的輸入，得到的 bi是回歸系數(shù) P的置信區(qū)間。用命令nlintool(x,y, 'model',b)可以得到一個交互式畫面，其內(nèi)容和用法與多項(xiàng)式回歸的Polytool類似。例4酶促反應(yīng)速度與底物濃度酶促反應(yīng)動力學(xué)簡稱酶動力學(xué)，主要研究酶促反應(yīng)速度與底物(即反應(yīng)物)濃度以及其它因素的關(guān)系。 在底物濃度很低時(shí)酶促反應(yīng)是一級反應(yīng)；當(dāng)?shù)孜餄舛忍幱谥?/p>

43、間范圍時(shí)，是混合級反應(yīng)；當(dāng)?shù)孜餄舛仍黾訒r(shí)， 向零級反應(yīng)過渡。某生化系學(xué)生為了研究喋吟霉素在某項(xiàng)酶促反應(yīng)中對反應(yīng)速度與底物濃度之間關(guān)系的影 響，設(shè)計(jì)了兩個實(shí)驗(yàn)，一個實(shí)驗(yàn)中所使用的酶是經(jīng)過喋吟霉素處理的，而另一個實(shí)驗(yàn)所用的酶是未經(jīng)喋吟霉 素處理的。所得實(shí)驗(yàn)數(shù)據(jù)見表8-5 o試根據(jù)問題的背景和這些數(shù)據(jù)建立一個合適的數(shù)學(xué)模型，來反映這項(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/速度處理7647971071231

44、39159152191201207200分析與假設(shè)記酶促反應(yīng)的速度為y，底物濃度為x ,二者之間的關(guān)系寫作y = f (x, P )，其中B為參數(shù)(B可為一向量)。由酶促反應(yīng)的基本性質(zhì)可知，當(dāng)?shù)孜餄舛群艿蜁r(shí)酶促反應(yīng)是一級反應(yīng)，此時(shí)反應(yīng)速度大致與底物濃度 成正比；而當(dāng)?shù)孜餄舛群艽螅瑵u近飽和時(shí)，反應(yīng)速度將趨于一個固定值(即零級反應(yīng))。下面的兩個簡單模 型具有這種性質(zhì)：Michaelis-Menten 模型指數(shù)增長模型非線性模型的求解首先作出給出的經(jīng)過喋吟霉素處理和未經(jīng)處理的反應(yīng)速度與底物濃度的散點(diǎn)圖，可以看出，上述兩個模 型與實(shí)際數(shù)據(jù)得到的散點(diǎn)圖是大致符合的。我們將主要對前一模型即Michaeli

45、s-Menten模型進(jìn)行詳細(xì)的分析。首先對經(jīng)過喋吟酶素處理的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行分析，在此基礎(chǔ)上，再來討論是否有更一般的模型來統(tǒng)一刻畫處理前后的數(shù)據(jù)，進(jìn)而揭示其中的聯(lián)系。我們用非線性回歸的方法直接估計(jì)模型的參數(shù)P1, P2,模型的求解可利用MATLAB統(tǒng)計(jì)工具箱中的命令進(jìn)行，使用格式為：beta,R,J=nlinfit(x,y,'model',beta0)其中輸入x為自變量數(shù)據(jù)矩陣，每列一個變量；y為因變量數(shù)據(jù)向量； model為模型的M文件名，M函數(shù)形式為y=f (beta,x),beta為待估計(jì)參數(shù)；beta。為給定的參數(shù)初值。輸出 beta為參數(shù)估計(jì)值，R為殘差，J為用于 估計(jì)預(yù)

46、測誤差的 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 200;beta0=195.8027 0.04841;beta,R,J=nlinfit(x,y,'huaxue',beta0);betaci=nlparci(beta,R,J);beta,betaciyy=b

47、eta(1)*x./(beta(2)+x);plot(x,y,'o',x,yy,'m+'),pausenlintool(x,y,'huaxue',beta)得到的數(shù)值2果見表 8-6 oNlintool用于給出一個交互式畫面，可以得到因變量y的預(yù)測值和預(yù)測區(qū)間，左下方的Export可向工作區(qū)傳送剩余標(biāo)準(zhǔn)差等數(shù)據(jù)。表8-6模型參數(shù)的估計(jì)結(jié)果參數(shù)參數(shù)估計(jì)值置信區(qū)間212,6818197.2028228.16080.06410.04570.0826從上面的結(jié)果可以知道，對經(jīng)過喋吟霉素處理的實(shí)驗(yàn)數(shù)據(jù)，在用Michaelis-Menten模型進(jìn)行回歸分析時(shí)，最 終反應(yīng)速度為 P 1=212