matlab回歸分析方法

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

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

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

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

6、的值，為測量值，殘差平方和剩余方差（3）線性相關(guān)性的檢驗由于我們采用的是一元線性回歸，因此，如果模型可用的話，應(yīng)該具有較好的線性關(guān)系。反映模型是否具有良好線性關(guān)系可通過相關(guān)系數(shù)R的值及F值觀察（后面的例子說明）。（4）一元線性回歸的MATLAB實現(xiàn) MATLAB工具箱中用命令regress實現(xiàn)，其用法是： b=regress(y,x) b ,bint , r ,rint , s=regress(y , x , alpha)輸入y（因變量，列向量）、x（1與自變量組成的矩陣，見下例），alpha是顯著性水平（缺省時默認0.05）。輸出，注意：b中元素順序與擬合命令polyfit的輸出不同，bin

7、t是的置信區(qū)間，r是殘差（列向量），rint是殘差的置信區(qū)間，s包含4個統(tǒng)計量：決定系數(shù)（相關(guān)系數(shù)為R）；F值；F(1,n-2)分布大于F值的概率p；剩余方差的值（MATLAB7.0以后版本）。也可由程序sum(r.2)/(n-2)計算。其意義和用法如下：的值越接近1，變量的線性相關(guān)性越強，說明模型有效；如果滿足，則認為變量與顯著地有線性關(guān)系，其中的值可查F分布表，或直接用MATLAB命令finv(1-,1, n-2)計算得到；如果表示線性模型可用。這三個值可以相互印證。的值主要用來比較模型是否有改進，其值越小說明模型精度越高。8.1.2身高與腿長例1 測得16名成年女子身高與腿長所得數(shù)據(jù)如下

8、： 表8-1 16名女子身高(cm)腿長(cm)數(shù)據(jù)88 85 88 91 92 93 93 95 96 98 97 96 98 99 100 102143 145 146 147 149 150 153 154 155 156 157 158 159 160 162 164 首先利用命令plot(x,y,'r*')畫出散點圖，從圖形可以看出，這些點大致分布在一條直線的左右，因此，可以考慮一元線性回歸?？删幹瞥绦蛉缦拢簓=143 145 146 147 149 150 153 154 155 156 157 158 159 160 162 164;x=88 85 88 91 9

9、2 93 93 95 96 98 97 96 98 99 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)運行后得到b = 31.7713 1.2903bint = 12.3196 51.2229 1.0846 1.4960 s = 0.9282 180.9531 0.0000 3.1277=0.9282，由finv(0.95,1,14)= 4.6001，即= 4.6001<F=180.9531，p<0.0001，可以通過殘差圖發(fā)現(xiàn)，第二個數(shù)據(jù)

10、為奇異數(shù)據(jù)，去掉該數(shù)據(jù)后運行后得到b = 17.6549 1.4363bint = -0.5986 35.9083 1.2445 1.6281 s = 0.9527 261.6389 0.0000 1.9313=0.9527，由finv(0.95,1,13)= 4.6672，即= 4.6672<F=261.6389，p<0.0001，說明模型有效且有改進，因此我們得到身高與腿長的關(guān)系。當(dāng)然，也可以利用直線擬合得到同一方程。只不過不能得到參數(shù)置信區(qū)間和對模型進行檢驗。擬合程序如下：y=143 145 146 147 149 150 153 154 155 156 157 158 15

11、9 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,temp)注意：函數(shù)相同,但輸出一次函數(shù)參數(shù)順序與回歸分析（升冪排列）中不同。另一個差別是擬合不能發(fā)現(xiàn)奇異數(shù)據(jù)。8.2 多元線性回歸分析8.2.1 多元線性回歸模型的建模步驟及其MATLAB實現(xiàn) 如果根據(jù)經(jīng)驗和有關(guān)知識認為與因變量有關(guān)聯(lián)的自變量不止一個，那么就應(yīng)該考慮用最小二乘準則建立多元線性回歸模型。 設(shè)影響因變量的主要因素（自變量）有m個，記，假

12、設(shè)它們有如下的線性關(guān)系式： ， 如果對變量與自變量 同時作n次觀察（n>m）得n組觀察值，采用最小二乘估計求得回歸方程.建立回歸模型是一個相當(dāng)復(fù)雜的過程，概括起來主要有以下幾個方面工作（1）根據(jù)研究目的收集數(shù)據(jù)和預(yù)分析；（2）根據(jù)散點圖是否具有線性關(guān)系建立基本回歸模型；（3）模型的精細分析；（4）模型的確認與應(yīng)用等。收集數(shù)據(jù)的一個經(jīng)驗準則是收集的數(shù)據(jù)量（樣本容量）至少應(yīng)為可能的自變量數(shù)目的610倍。在建模過程中首先要根據(jù)所研究問題的目的設(shè)置因變量，然后再選取與該因變量有統(tǒng)計關(guān)系的一些變量作為自變量。我們當(dāng)然希望選擇與問題關(guān)系密切的變量，同時這些變量之間相關(guān)性不太強，這可以在得到初步的模型

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

14、關(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.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.63

15、4.248.038.035.940.436.845.235.1試建立Y與之間關(guān)系的數(shù)學(xué)模型，并得出有關(guān)結(jié)論和作統(tǒng)計分析。2. 作出因變量Y與各自變量的樣本散點圖作散點圖的目的主要是觀察因變量Y與各自變量間是否有比較好的線性關(guān)系，以便選擇恰當(dāng)?shù)臄?shù)學(xué)模型形式。下圖分別為年薪Y(jié)與成果質(zhì)量指標(biāo)、研究工作時間、獲得資助的指標(biāo)之間的散點圖，subplot(1,3,1),plot(x1,Y,'g*'),% subplot是MATLAB中的函數(shù)。1 使用方法：subplot（m,n,p）或者subplot（m n p）。subplot是將多個圖畫到一個平面上的工具。其中，m表示是圖排成m行，n

16、表示圖排成n列，也就是整個figure中有n個圖是排成一行的，一共m行，如果m=2就是表示2行圖。p表示圖所在的位置，p=1表示從左到右從上到下的第一個位置。subplot(1,3,2),plot(x2,Y,'k+'),subplot(1,3,3),plot(x3,Y,'ro'),從圖可以看出這些點大致分布在一條直線旁邊，因此，有比較好的線性關(guān)系，可以采用線性回歸。 Y與x1的散點圖 Y與x2的散點圖 Y與x3的散點圖圖8.1 因變量Y與各自變量的樣本散點圖3. 利用MATLAB統(tǒng)計工具箱得到初步的回歸方程設(shè)回歸方程為：.建立m-文件輸入如下程序數(shù)據(jù)：x1=3.

17、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.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.

18、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,運行后即得到結(jié)果如表8-3所示。表8-3 對初步回歸模型的計算結(jié)果回歸系數(shù)回歸系數(shù)的估計值回歸系數(shù)的置信區(qū)間18.015713.9052 22.12621.08170.3900 1.77330.32120.2440 0.398

19、41.28350.6691 1.8979=0.9106 F=67.9195 p<0.0001 = 3.0719計算結(jié)果包括回歸系數(shù)b=()=(18.0157, 1.0817 , 0.3212 , 1.2835)，且置信區(qū)間均不包含零點，；殘差及其置信區(qū)間；統(tǒng)計變量stats ,它包含四個檢驗統(tǒng)計量：相關(guān)系數(shù)的平方，假設(shè)檢驗統(tǒng)計量，與F對應(yīng)的概率p，的值（7.0以前版本也可由程序sum(r.2)/(n-m-1)計算）。因此我們得到初步的回歸方程為：由結(jié)果對模型的判斷：回歸系數(shù)置信區(qū)間不包含零點表示模型較好，殘差在零點附近也表示模型較好，接著就是利用檢驗統(tǒng)計量，p的值判斷該模型是否可用。（）

20、相關(guān)系數(shù)的評價：一般地，相關(guān)系數(shù)絕對值在0.81范圍內(nèi)，可判斷回歸自變量與因變量具有較強的線性相關(guān)性。本例的絕對值為0.9542，表明線性相關(guān)性較強。（）F檢驗法：當(dāng)，即認為因變量與自變量之間顯著地有線性相關(guān)關(guān)系；否則認為因變量與自變量之間線性相關(guān)關(guān)系不顯著。本例67.919>= 3.10 （查F分布表或輸入命令finv(0.95,3,20)計算）。（）p值檢驗：若（為預(yù)定顯著水平），則說明因變量與自變量之間顯著地有線性相關(guān)關(guān)系。本例輸出結(jié)果，p<0.0001,顯然滿足P<=0.05。以上三種統(tǒng)計推斷方法推斷的結(jié)果是一致的，說明因變量與自變量之間顯著地有線性相關(guān)關(guān)系，所得線性

21、回歸模型可用。當(dāng)然越小越好，這主要在模型改進時作為參考。4. 模型的精細分析和改進（1） 殘差分析殘差，是各觀測值與回歸方程所對應(yīng)得到的擬合值之差，實際上，它是線性回歸模型中誤差的估計值。即有零均值和常值方差，利用殘差的這種特性反過來考察原模型的合理性就是殘差分析的基本思想。利用MATLAB進行殘差分析則是通過殘差圖或時序殘差圖。殘差圖是指以殘差為縱坐標(biāo)，以其他指定的量為橫坐標(biāo)的散點圖。主要包括：（1）橫坐標(biāo)為觀測時間或觀測值序號；（2）橫坐標(biāo)為某個自變量的觀測值；（3）橫坐標(biāo)為因變量的擬合值。通過觀察殘差圖，可以對奇異點進行分析，還可以對誤差的等方差性以及對回歸函數(shù)中是否包含其他自變量、自變

22、量的高次項及交叉項等問題給出直觀的檢驗。以觀測值序號為橫坐標(biāo)，殘差為縱坐標(biāo)所得到的散點圖稱為時序殘差圖，畫出時序殘差圖的MATLAB語句為rcoplot(r,rint)（圖8.2）?？梢郧宄吹綒埐畲蠖挤植荚诹愕母浇?，因此還是比較好的 ，不過第4、12、19這三個樣本點的殘差偏離原點較遠，如果作為奇異點看待，去掉后重新擬合，則得回歸模型為：且回歸系數(shù)的置信區(qū)間更小均不包含原點，統(tǒng)計變量stats包含的三個檢驗統(tǒng)計量：相關(guān)系數(shù)的平方，假設(shè)檢驗統(tǒng)計量，概率，分別為：0.9533 ； 115.5586 ； 0.0000 ，比較可知R，F(xiàn)均增加模型得到改進。 圖8.2 時序殘差圖（2） 變量間的交互作

23、用討論變量間的交互作用包括：不同自變量之間的交互作用以及同一變量的自相關(guān)性。不同自變量之間的交互作用：有時，在實驗中不僅單因素對指標(biāo)有影響，而且因素間還會聯(lián)合起來對指標(biāo)產(chǎn)生影響，常稱這種聯(lián)合作用為交互作用。處理兩個因素間交互作用的一個簡單辦法是加入這兩個自變量的乘積項。本文案例如果加入交互項則為：用表8.2的數(shù)據(jù)，利用MATLAB統(tǒng)計工具箱得到回歸系數(shù)分別為：27.0727 ，1.1147，-0.0215 ，-0.1843 ，0.0033 ，-0.0054 ，0.0511 。但它們的置信區(qū)間均包含原點，其他指標(biāo)也不理想，因此，本例中其交互作用并不顯著，該模型不如前面兩個模型好。自相關(guān)性的診斷和

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

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

26、地選擇變量進行統(tǒng)計分析。該命令的用法是：stepwise(X , Y , inmodel , alpha)其中X是自變量數(shù)據(jù)，排成矩陣（m為自變量個數(shù)，n為每個變量的數(shù)據(jù)量），Y是因變量數(shù)據(jù)，排成向量，inmodel 是自變量初始集合的指標(biāo)，缺省時為全部自變量，alpha為顯著水平，缺省時為0.05。運行stepwise命令時產(chǎn)生圖形窗口：Stepwise Plot , Stepwise Table , Stepwise History.當(dāng)鼠標(biāo)移到圖形某個區(qū)域時，鼠標(biāo)點擊后產(chǎn)生交互作用。Stepwise Plot窗口中的虛線表示回歸系數(shù)的置信區(qū)間包含零點，即該回歸系數(shù)與零無顯著差異，一般應(yīng)將該

27、變量移去；實線則表明該回歸系數(shù)與零有顯著差異，應(yīng)保留在模型中(藍色表示該變量已進入模型，紅色表示該變量已移出模型)。引入和移出變量還可參考Stepwise History窗口中剩余標(biāo)準差RMSE是否在下降，剩余標(biāo)準差RMSE最小的就是最好的模型。Stepwise Table窗口中列出了一個統(tǒng)計表，包括回歸系數(shù)及其置信區(qū)間，以及模型的統(tǒng)計量剩余標(biāo)準差RMSE、相關(guān)系數(shù)R-square、F值、與F對應(yīng)的概率。關(guān)于本節(jié)案例2，如果引入新的自變量 . 也可以采用逐步回歸法解決,源程序如下：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

28、 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 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.

29、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)運行并按上述步驟操作后可以得到本文前面線性回歸相同的結(jié)論，即不含交互項的模型是最好的。在此只介紹操作過程，其交互界面，只要在MATLAB軟件上一試便知。8.2.4 多項式回歸多項式回歸仍然屬于多元線性回歸，可以是一元多項式回歸或多元多項式回歸。一元多項式回歸模型的一般形式為用MATLAB求解一元多項式回歸，除了使用命令polyfit(x,y,m)外，還可

30、以使用如下命令： Polytool(x,y,m,alpha)輸入x,y,m同命令polyfit，alpha是顯著性水平（默認0.05），則輸出一個交互式畫面，畫面顯示回歸曲線及其置信區(qū)間，通過圖左下方的export下拉式菜單，還可以輸出回歸系數(shù)估計值及其置信區(qū)間、殘差等。下面通過一個用多元多項式回歸的實例說明什么時候用多項式回歸以及如何通過MATLAB軟件進行處理。例3 為了了解人口平均預(yù)期壽命與人均國內(nèi)生產(chǎn)總值和體質(zhì)得分的關(guān)系，我們查閱了國家統(tǒng)計局資料，北京體育大學(xué)出版社出版的2000國民體質(zhì)監(jiān)測報告，表8-4是我國大陸31個省市的有關(guān)數(shù)據(jù)。我們希望通過這幾組數(shù)據(jù)考察它們是否具有良好的相關(guān)關(guān)

31、系，并通過它們的關(guān)系從人均國內(nèi)生產(chǎn)總值（可以看作反映生活水平的一個指標(biāo)）、體質(zhì)得分預(yù)測其壽命可能的變化范圍。體質(zhì)是指人體的質(zhì)量，是遺傳性和獲得性的基礎(chǔ)上表現(xiàn)出來的人體形態(tài)結(jié)構(gòu)，生理機能和心理因素綜合的、相對穩(wěn)定的特征。體質(zhì)是人的生命活動和工作能力的物質(zhì)基礎(chǔ)。它在形成、發(fā)展和消亡過程中，具有明顯的個體差異和階段性。中國體育科學(xué)學(xué)會體質(zhì)研究會研究表明，體質(zhì)應(yīng)包括身體形態(tài)發(fā)育水平、生理功能水平、身體素質(zhì)和運動能力發(fā)展水平、心理發(fā)育水平和適應(yīng)能力等五個方面。目前，體質(zhì)的綜合評價主要是形態(tài)、機能和身體素質(zhì)三類指標(biāo)按一定的權(quán)重進行換算而得。 表8-4 31個省市人口預(yù)期壽命與人均國內(nèi)生產(chǎn)總值和體質(zhì)得分數(shù)據(jù)

32、序號預(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.125100601565.9662.953822676.1069.34547319573.9169.99299311672.3766.1190702774.9168.41540643672.

33、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.1390882173.1067.065143351174.7069 .505337722267.4763.6057898模型的建立和求解 作表8-4數(shù)據(jù)的散點圖如圖8.3圖8.

34、3 預(yù)期壽命與人均國內(nèi)生產(chǎn)總值和體質(zhì)得分的散點圖 從圖8.3可以看出人口預(yù)期壽命與體質(zhì)得分有較好的線性關(guān)系，與人均國內(nèi)生產(chǎn)總值的關(guān)系難以確定，我們建立二次函數(shù)的回歸模型。一般的多元二項式回歸模型可表為 MATLAB統(tǒng)計工具箱提供了一個很方便的多元二項式回歸命令：Rstool(x,y, 'model',alpha)輸入x為自變量（n×m矩陣），y為因變量（n維向量），alpha為顯著水平，model從下列4個模型中選擇一個：linear（只包含線性項）purequadratic（包含線性項和純二次項）interaction（包含線性項和純交互項）quadratic（包含

35、線性項和完全二次項）輸出一個交互式畫面，對例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.55 71.65 71.73,73.10 67.47 69.87 67.41 78.14 76.10 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

36、 5382 19070 10935 22007 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 67.97 62.9 66.1 64.51 68.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

37、.84;x=x1',x2'rstool(x,y','purequadratic')得到一個如圖8.4的交互式畫面 圖8.4 預(yù)期壽命與人均國內(nèi)生產(chǎn)總值和體質(zhì)得分的一個交互式畫面左邊一幅圖形是固定時的曲線及其置信區(qū)間，右邊一幅圖形是固定時的曲線及其置信區(qū)間。移動鼠標(biāo)可改變，的值，同時圖左邊給出的預(yù)測值及其置信區(qū)間。如輸入=128757，=66.165，則=70.6948，其置信區(qū)間70.6948±1.1079。圖的左下方有兩個下拉式菜單，上面的菜單Export用于輸出數(shù)據(jù)（包括：回歸系數(shù)parameters,殘差residuals,剩余標(biāo)準差RM

38、SE等）, 在MATLAB工作空間中得到有關(guān)數(shù)據(jù)。通過下面的菜單在上述4個模型中變更選擇，最后確定RMSE值較小的模型。例3則是包含線性項和完全二次項（quadratic）的模型最佳，即剩余標(biāo)準差為1.2622，因此，所得回歸模型為：利用此模型我們可以根據(jù)國內(nèi)生產(chǎn)總值及體質(zhì)得分，預(yù)測壽命。8.3 非線性回歸分析8.3.1 非線性最小二乘擬合線性最小二乘擬合與線性回歸中的“線性”并非指與的關(guān)系，而是指是系數(shù)或的線性函數(shù)。擬合如的函數(shù)仍然是最小二乘擬合；如果擬合如的曲線，對是非線性的，但取對數(shù)后對系數(shù)是線性的，屬于可化為線性回歸的類型。下面討論非線性擬合的情形。非線性最小二乘擬合問題的提法是：已知

39、模型，其中對是非線性的，為了估計參數(shù)，收集n個獨立觀測數(shù)據(jù)。記擬合誤差，求使誤差的平方和最小。作為無約束非線性規(guī)劃的特例，解非線性最小二乘擬合可用MATLAB優(yōu)化工具箱命令lsqnonlin和lsqcurvefit。8.3.2 非線性回歸模型非線性回歸模型記作其中對回歸系數(shù)是非線性的，。求得回歸系數(shù)的最小二乘估計。MATLAB統(tǒng)計工具箱中非線性回歸的命令是：b,R,J=nlinfit(x,y, 'model',bo)輸入x是自變量數(shù)據(jù)矩陣，每列一個向量；y是因變量數(shù)據(jù)向量；model是模型的函數(shù)名（M文件），形式為，b為待估系數(shù)；b0是回歸系數(shù)的初值。輸出b是的估計值，R是殘差

40、，J是用于估計預(yù)測誤差的Jacobi矩陣。這個命令是依據(jù)高斯牛頓法求解的。將上面的輸出作為命令Bi=nlparci(b,R,J)的輸入，得到的bi是回歸系數(shù)的置信區(qū)間。用命令nlintool(x,y, 'model',b)可以得到一個交互式畫面，其內(nèi)容和用法與多項式回歸的Polytool類似。例4 酶促反應(yīng)速度與底物濃度 酶促反應(yīng)動力學(xué)簡稱酶動力學(xué)，主要研究酶促反應(yīng)速度與底物（即反應(yīng)物）濃度以及其它因素的關(guān)系。在底物濃度很低時酶促反應(yīng)是一級反應(yīng)；當(dāng)?shù)孜餄舛忍幱谥虚g范圍時，是混合級反應(yīng)；當(dāng)?shù)孜餄舛仍黾訒r，向零級反應(yīng)過渡。某生化系學(xué)生為了研究嘌呤霉素在某項酶促反應(yīng)中對反應(yīng)速度與底物

41、濃度之間關(guān)系的影響，設(shè)計了兩個實驗，一個實驗中所使用的酶是經(jīng)過嘌呤霉素處理的，而另一個實驗所用的酶是未經(jīng)嘌呤霉素處理的。所得實驗數(shù)據(jù)見表8-5。試根據(jù)問題的背景和這些數(shù)據(jù)建立一個合適的數(shù)學(xué)模型，來反映這項酶促反應(yīng)的速度與底物濃度以及嘌呤霉素處理與否之間的關(guān)系。表8-5 嘌呤霉素實驗中的反應(yīng)速度與底物濃度數(shù)據(jù)底物濃度（ppm）0.020.060.110.220.561.10反應(yīng)速度未處理6751848698115131124144158160/處理764797107123139159152191201207200分析與假設(shè)記酶促反應(yīng)的速度為,底物濃度為，二者之間的關(guān)系寫作,其中為參數(shù)(可為一向量

42、)。由酶促反應(yīng)的基本性質(zhì)可知，當(dāng)?shù)孜餄舛群艿蜁r酶促反應(yīng)是一級反應(yīng)，此時反應(yīng)速度大致與底物濃度成正比；而當(dāng)?shù)孜餄舛群艽螅瑵u近飽和時，反應(yīng)速度將趨于一個固定值（即零級反應(yīng)）。下面的兩個簡單模型具有這種性質(zhì)：Michaelis-Menten 模型指數(shù)增長模型非線性模型的求解首先作出給出的經(jīng)過嘌呤霉素處理和未經(jīng)處理的反應(yīng)速度與底物濃度的散點圖，可以看出，上述兩個模型與實際數(shù)據(jù)得到的散點圖是大致符合的。我們將主要對前一模型即Michaelis-Menten模型進行詳細的分析。首先對經(jīng)過嘌呤酶素處理的實驗數(shù)據(jù)進行分析，在此基礎(chǔ)上，再來討論是否有更一般的模型來統(tǒng)一刻畫處理前后的數(shù)據(jù)，進而揭示其中的聯(lián)系。我們

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

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