回歸分析方法

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ì)回歸模型?；貧w模型常用來解決預(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計(jì)工具箱幾乎包括了數(shù)理統(tǒng)計(jì)方面主要的概念、理 論、方法和算法。運(yùn)用 MATLAB計(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)的回歸模型。某些非線性回歸模型可以化為線性回歸模型處理；如 果知道函數(shù)

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

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

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

8、系，其中F1 (1,n 2)的值可查F分布表， 或直接用MATLA命令finv(1-,1, n-2)計(jì)算得到；如果p表示線性模型可用。這三個值可以相互印證。s2的值主要用來比較模型是否有改進(jìn)，其值越小說明模型精度越高。例1測得16名成年女子身高y與腿長x所得數(shù)據(jù)如下：表8-1 16名女子身高(cm)腿長(cm)數(shù)據(jù)88 85 88 91 92 93 93 95 96 98 97 96 9899 100 102143 145 146 147 149 150 153 154 155 156 157 158 159160 162 164首先利用命令plot(x,y,'r*')畫出散點(diǎn)

9、圖，從圖形可以看出，這些點(diǎn)大致分布在一條直線的左右，因此，可以考慮一元線性回歸?？删幹瞥绦蛉缦拢簓二143 145 146 147 149 150 153 154 155 156 157 158 159 160 162164;x=88 85 88 9192 93 93 95 96 98 97 96 98 99 100102;n=16;X=o nes( n,1),x'b,b in t,r,ri nt,s=regress(y',X,0.05);b,bi nt,s,rcoplot(r,r int)運(yùn)行后得到b = 31.77131.2903bint = 12.3196 51.2229

10、1.08461.4960s = 0.9282 180.95310.00003.1277R2 =0.9282，由 finv(0.95,1,14)=4.6001，即 F1 (1,n 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.5986 35.90831.24451.6281s = 0.9527 261.63890.0000 1.9313R2 =0.9527，由 finv(0.95,1,13)=4.6672，即 R (1,n 2)= 4.6672<

11、F=261.6389 , p<0.0001，y 17.6549 1.4363 x 。說明模型有效且有改進(jìn)，因此我們得到身高與腿長的關(guān)系當(dāng)然，也可以利用直線擬合得到同一方程。 只不過不能得到參數(shù)置信區(qū)間和對模型進(jìn)行 檢驗(yàn)。擬合程序如下： y=143 145 146 147 149 150 153 154 155 156 157 158 159 160 162 164;x=88 85 88 91 92 93 93 95 96 98 97 96 98 99 100 102;a=polyfit（x,y,1）temp=polyval（a,x）;plot（x,y,'r*',x,tem

12、p）注意：函數(shù)相同 , 但輸出一次函數(shù)參數(shù)順序與回歸分析（升冪排列）中不同。另一個差別是 擬合不能發(fā)現(xiàn)奇異數(shù)據(jù)。8.2 多元線性回歸分析821多元線性回歸模型的建模步驟及其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 （xiL n），假設(shè)它們有如下的線 性關(guān)系式：2y 0 1x1 Lmxm， N（0, ）如果對變量y與自變量x,x2丄，xm同時(shí)作n次觀察（n>m得n組觀察值，采用最小二 乘估計(jì)求得回歸方程y? ?0 ?1x1 L ?k xm .建立回歸模型是一

13、個相當(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)用。多元線性回歸的 MA

14、TLA或現(xiàn)仍然用命令regress（y , X），只是要注意矩陣X的形式，將通過如下例子說明其用法。822 某類研究學(xué)者的年薪1. 問題例2工薪階層關(guān)心年薪與哪些因素有關(guān)，以此可制定出它們自己的奮斗目標(biāo)。某科學(xué)基金會希望估計(jì)從事某研究的學(xué)者的年薪Y(jié)與他們的研究成果（論文、著作等）的質(zhì)量指標(biāo)X、從事研究工作的時(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.54.99201833

15、3113253054725116.16.47.46.77.55.96.04.05.88.35.06.433.40.38.46.41.37.39.40.30.52.38.31.237845071928i1314151617181920212223248.06.56.63.76.27.04.04.55.95.64.83.9233539217403523332734157.67.05.04.45.57.06.03.54.94.38.05.843.44.42.33.34.48.38.35.40.36.45.35.315620094821試建立Y與X!,X2,X3之間關(guān)系的數(shù)學(xué)模型，并得出有關(guān)結(jié)論和作統(tǒng)

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

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

18、 5.0 6.4 7.6 7.0 5.0 4.0 5.5 7.0 6.03.5 4.9 4.3 8.0 5.0;Y二33.2 40.3 38.746.8 41.4 37.5 39.0 40.7 30.1 52.9 38.2 31.8 43.3 44.1 42.5 33.634.2 48.0 38.0 35.9 40.4 36.8 45.2 35.1;n=24; m=3;X=o nes( n,1),x1',x2',x3'b,bi nt,r,ri nt,s=regress( Y',X,0.05);b,bi nt,r,r in t,s,運(yùn)行后即得到結(jié)果如表 8-3所示

19、。表8-3對初步回歸模型的計(jì)算結(jié)果回歸系數(shù)回歸系數(shù)的估計(jì)值回歸系數(shù)的置信區(qū)間18.015713.9052 22.12621.08170.3900 1.77330.32120.2440 0.39841.28350.6691 1.8979R2 =0.9106F=67.9195p<0.0001s2 = 3.0719y 與自變量之間顯著地有線以上三種統(tǒng)計(jì)推斷方法推斷的結(jié)果是一致的，說明因變量性相關(guān)關(guān)系，所得線性回歸模型可用。s2 當(dāng)然越小越好，這主要在模型改進(jìn)時(shí)作為參考。4. 模型的精細(xì)分析和改進(jìn)（ 1） 殘差分析殘差e yi ? （i 1,2丄,n）,是各觀測值y與回歸方程所對應(yīng)得到的擬合值？

20、之差，實(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)，殘差為縱坐標(biāo)所得到的散點(diǎn)圖稱為時(shí)序殘差圖，

21、畫出時(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） 變量間的交互作用討論 變量間的交互作用包括：不同自變量之間的交互作用以及同一變量的自相關(guān)性

22、。不同自變量之間的交互作用 ：有時(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í)間序列數(shù)據(jù)。在時(shí)間序列數(shù)據(jù)中，同一變量的順

23、序觀測值之間出現(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)(eM,el) ,t 2,3,L ,n大部分點(diǎn)落在第I,川象限，表明存在著正的序列相 關(guān)；如果大部分點(diǎn)落在第U，W象限，表明存在著負(fù)的序列相關(guān)。對DW檢驗(yàn)法可以利用MATLAB軟件編程計(jì)算統(tǒng)計(jì)量：etet 1DW 2(1然后查閱DW僉驗(yàn)上下界表，以決定模型的自相關(guān)狀態(tài)。當(dāng)一個回歸模型存在序列相關(guān)性時(shí)，首先要查明序列相關(guān)產(chǎn)生的原因。如果是回歸模型選用不當(dāng)，則應(yīng)

24、改用適當(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)行，直到不能引入和移出為止。弓I入和移出都以給定的顯著性水平為標(biāo)準(zhǔn)。MATLAB統(tǒng)計(jì)工具箱中逐步回歸的命令是stepwise，它提供了一個人機(jī)交互式畫面，通過此工

25、具可以自由地選擇變量進(jìn)行統(tǒng)計(jì)分析。該命令的用法是：stepwise(X , Y , inmodel , alpha)其中X是自變量數(shù)據(jù)，排成n m矩陣(m為自變量個數(shù)，n為每個變量的數(shù)據(jù)量)，Y是因 變量數(shù)據(jù)， 排成 n 1向量， inmodel 是自變量初始集合的指標(biāo)， 缺省時(shí)為全部自變量， alpha 為顯著水平，缺省時(shí)為 0.05 。運(yùn)行 stepwise 命令時(shí)產(chǎn)生圖形窗口： Stepwise Plot , Stepwise Table , StepwiseHistory. 當(dāng)鼠標(biāo)移到圖形某個區(qū)域時(shí)，鼠標(biāo)點(diǎn)擊后產(chǎn)生交互作用。 Stepwise Plot 窗口中的 虛線表示回歸系數(shù)的置信

26、區(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)差RMS、E相關(guān)系數(shù)R-square、F值、與F對應(yīng)的概率。關(guān)于本節(jié)案例 2，如果引入新的自變量 x4 x1x2, x5 x1x3,x6 x2x3 . 也可以采用逐步回 歸法

27、解決 , 源程序如下：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.55.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 3415;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.54.9 4.3 8.0 5.0'Y=33.2 40.3 38.7 46.8 41.

28、4 37.5 39.0 40.7 30.1 52.9 38.2 31.8 43.3 44.1 42.5 33.634.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)式回歸仍然屬于多元線性回歸，可以是一元多項(xiàng)式回歸或多

29、元多項(xiàng)式回歸。 一元多項(xiàng)式回歸模型的一般形式為用MATLA求解一元多項(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)系，我們查閱 了國

30、家統(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ā)展水平

31、、心理發(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.566.1612851265.456.787442369.64.3177145797587057273.971.2524491368.966.011492467.60.415202551441855373.270.1324251473.367.920462578.70.270627504711492471.265.12100

32、61565.962.953822676.69.34731050610459573.969.9929931672.366.119072774.68.44064117091153672.565.7618241770.064.510932872.66.4117845371591951770.667.2910761872.568.322002970.65.7106563585717658871.867.7199071971.666.213593066.63.21158550540387971.066.5213252071.765.711473164.62.897258553,743741071.267

33、.1390882173.167.014339,06551174.769 .533772267.463.678980052705模型的建立和求解作表8-4數(shù)據(jù)(xi, 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計(jì)工具箱提供了一個很方便的多元二項(xiàng)式回歸命令：Rstool（x,y, 'model',alpha）輸入x為自變量（nxm矩陣），y為因變量（n維向量

34、），alpha為顯著水平，model從 下列4個模型中選擇一個：lin ear（只包含線性項(xiàng)）purequadratic（包含線性項(xiàng)和純二次項(xiàng)）in teracti on（包含線性項(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.9573.34 65.96 72.37 70.07 72.55 71.65 71.73,73.1067.47 69.87 67.41 78.14 76.10 74.917

35、2.91 70.17 66.03 64.37;x1=12857 24495 24250 10060 29931 18243 10763 9907 13255 9088 33772 8744 11494 204615382 19070 10935 22007 13594 11474 14335 7898 17717 15205 70622 47319 40643 1178110658 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.77566.01 67.97 62

36、.9 66.1 64.51 68.385 66.205 65.77,67.065 63.605 64.305 60.485 70.2969.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ì)得分的一個交互式畫面左邊一幅圖形是X2固定時(shí)的曲線y(xj及其置信區(qū)間，右邊一幅圖形是xi固定時(shí)的曲線y(X2) 及其置信區(qū)間。移動鼠標(biāo)可改變xi，X2的值，同時(shí)圖左邊給出 y的預(yù)測值及其置信區(qū)

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

38、分析8.3.1 非線性最小二乘擬合線性最小二乘擬合與線性回歸中的“線性”并非指y與x的關(guān)系，而是指y是系數(shù)0, 1或 (0, 1丄，m)的線性函數(shù)。擬合如y 0低2的函數(shù)仍然是最小二乘擬合；如果擬合 如y°e1x的曲線，y對°, 1是非線性的，但取對數(shù)后Iny對系數(shù)°, 1是線性的，屬于可化為線性回歸的類型。下面討論非線性擬合的情形。非線性最小二乘擬合問題的提法是：已知模型y f (x, ), x (x1,L ,xm),( 0, 1,L , k) ，其中 f 對 是非線性的，為了估計(jì)參數(shù) ，收集 n 個獨(dú)立觀測數(shù)據(jù) (xi,yi),xi (xi1,L xim) (

39、i 1,L ,n),n m 。記擬合誤差 i( ) yi f (xi, )，求 使誤差的平方和 最小。作為無約束非線性規(guī)劃的特例，解非線性最小二乘擬合可用MATLAB優(yōu)化工具箱命令lsqnonlin 和 lsqcurvefit 。8.3.2 非線性回歸模型非線性回歸模型記作其中 f 對回歸系數(shù) 是非線性的， N(0, 2 ) 。求得回歸系數(shù)的最小二乘估計(jì)。MATLAB計(jì)工具箱中非線性回歸的命令是：b,R,J=nlinfit(x,y, 'model',bo)輸入 x 是自變量數(shù)據(jù)矩陣，每列一個向量； y 是因變量數(shù)據(jù)向量； model 是模型的函數(shù) 名(M文件)，形式為y f(b

40、,x)，b為待估系數(shù)；bO是回歸系數(shù)的初值。輸出b是 的估計(jì)值， R 是殘差， J 是用于估計(jì)預(yù)測誤差的 Jacobi 矩陣。這個命令是依據(jù)高斯牛頓法 求解的。將上面的輸出作為命令Bi=nlparci(b,R,J)的輸入，得到的 bi 是回歸系數(shù) 的置信區(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ù)孜餄舛忍幱谥虚g范圍時(shí)， 是混合級

41、反應(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í)來反映這項(xiàng)酶促驗(yàn)數(shù)據(jù)見表8-5。試根據(jù)問題的背景和這些數(shù)據(jù)建立一個合適的數(shù)學(xué)模型, 反應(yīng)的速度與底物濃度以及嘌呤霉素處理與否之間的關(guān)系表8-5嘌呤霉素實(shí)驗(yàn)中的反應(yīng)速度與底物濃度數(shù)據(jù)底物濃度0.020.060.110.220.561.10(ppm反未處6584869811131121441516/應(yīng)理715480速處理74971012313159151912020200

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

43、分析。首先對經(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二n li nfit(x,y,'model',betaO)其中輸入x為自變量數(shù)據(jù)矩陣，每列一個變量；y為因變量數(shù)據(jù)向量；model為模型的M文件名，M函數(shù)形式為y=f (beta,x),beta 為待估計(jì)參數(shù)；betaO為給定的參數(shù)初值。輸出beta 為參數(shù)估計(jì)值，R為殘差，J為用于估計(jì)預(yù)測誤差的Jacobi矩陣。參

44、數(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二n li nfit(x,y,'huaxue',beta0);betaci 二n lparci(beta,R,J);beta,betaciyy=beta(1)*x.

45、/(beta(2)+x);plot(x,y,'o',x,yy,'m+'),pausenli ntool(x,y,'huaxue',beta)得到的數(shù)值結(jié)果見表8-6。Nlintool用于給出一個交互式畫面，可以得到因變量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)速度

46、為 =212.6818,反應(yīng)的半速度點(diǎn)（達(dá)到最終反應(yīng)速度的一半 時(shí)的底物濃度x值）恰為 =0.06412 o混合反應(yīng)模型由酶動力學(xué)知識我們知道，酶促反應(yīng)的濃度依賴于底物濃度，并且可以假定，嘌呤霉素 的處理會影響最終反應(yīng)速度參數(shù)，而基本上不影響半速度參數(shù).表8-5的數(shù)據(jù)也印證了這種看法。Michaelis-Me nte n模型的形式可以分別描述經(jīng)過嘌呤霉素處理和未處理的反應(yīng)速度與底物濃度的關(guān)系（兩個模型的參數(shù) 會不同）,為了在同一個模型中考慮嘌呤霉素處 理的影響，我們采用對未經(jīng)嘌呤霉素處理的模型附加增量的方法，考察如下的混合反應(yīng)模型：其中自變量X為底物濃度，X2為一示性變量（0-1變量），用來表示是否經(jīng)嘌呤霉素處理， X2 =1表示經(jīng)過處理，X2=0表示未經(jīng)處理；參數(shù)是未處理的反應(yīng)的最終反應(yīng)速度，是經(jīng)處理后最終反應(yīng)速度的增長值，是未經(jīng)處理的反應(yīng)的半速度點(diǎn)，是經(jīng)處理后反應(yīng)的半 速度點(diǎn)的增長值。混合模型的求解和分析為了給出初始迭代值，從實(shí)驗(yàn)數(shù)據(jù)我們注意到，未經(jīng)處理的反應(yīng)速度的最大實(shí)驗(yàn)值為160,經(jīng)過處理的最大實(shí)驗(yàn)值為207,于是可取參數(shù)初值10 170, ； 60 ;又從數(shù)據(jù)可大致估計(jì)未經(jīng)處理的半速度點(diǎn)約為0.05 ,經(jīng)過處理的半速度點(diǎn)約為0.06 ,我們?nèi)?0.05,