回歸分析方法

1、第八章 回歸分析方法當(dāng)人們對研究對象的內(nèi)在特性和各因素間的關(guān)系有比較充分的認(rèn)識(shí)時(shí)，一般用機(jī)理分析方法建立數(shù)學(xué)模型。如果由于客觀事物內(nèi)部規(guī)律的復(fù)雜性及人們認(rèn)識(shí)程度的限制，無法分析實(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)系，其特征是：一個(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ì)分析方法對不同的模型進(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ì)管理及社會(huì)學(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ì)而繁雜的理論。在此基礎(chǔ)上再介紹在建模過

4、程中如何有效地使用MATLAB軟件。沒有學(xué)過這部分?jǐn)?shù)學(xué)知識(shí)的讀者可以不深究其數(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ù)形式只是要確定其中的參數(shù)則是擬

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

6、值，為測量值，殘差平方和剩余方差（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 , alpha)輸入y（因變量，列向量）、x（1與自變量組成的矩陣，見下例），alpha是顯著性水平（缺省時(shí)默認(rèn)0.05）。輸出，注意：b中元素順序與擬合命令polyfit的輸出不同，bint

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

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

9、 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)運(yùn)行后得到b = 31.7713 1.2903bint = 12.3196 51.2229 1.0846 1.4960s = 0.9282 180.9531 0.0000 3.1277=0.9282，由finv(0.95,1,14)= 4.6001，即= 4.6001F=180.9531，p0.0001，可以通過殘差圖發(fā)現(xiàn)，第二個(gè)數(shù)據(jù)為奇異數(shù)據(jù)，去掉該數(shù)據(jù)后運(yùn)行后得到b = 17.6549 1.4363bint

10、= -0.5986 35.9083 1.2445 1.6281s = 0.9527 261.6389 0.0000 1.9313=0.9527，由finv(0.95,1,13)= 4.6672，即= 4.6672F=261.6389，pm）得n組觀察值，采用最小二乘估計(jì)求得回歸方程.建立回歸模型是一個(gè)相當(dāng)復(fù)雜的過程，概括起來主要有以下幾個(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ù)量（樣本容量）至少應(yīng)為可能的自變量數(shù)目的610倍。在建模過程中首先要根據(jù)所研究問題的

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

12、作的時(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.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.8

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

14、),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. 利用MATLAB統(tǒng)計(jì)工具箱得到初步的回歸方程設(shè)回歸方程為：.建立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 2

15、5 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.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

16、,rint,s=regress(Y,X,0.05);b,bint,r,rint,s,運(yùn)行后即得到結(jié)果如表8-3所示。表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.8979=0.9106 F=67.9195 p= 3.10 （查F分布表或輸入命令finv(0.95,3,20)計(jì)算）。（）p值檢驗(yàn)：若（為預(yù)定顯著水平），則說明因變量與自變量之間顯著地有線性相關(guān)關(guān)系。本例輸出結(jié)果，p0.0001,顯然滿足P=0.05。以上

17、三種統(tǒng)計(jì)推斷方法推斷的結(jié)果是一致的，說明因變量與自變量之間顯著地有線性相關(guān)關(guān)系，所得線性回歸模型可用。當(dāng)然越小越好，這主要在模型改進(jìn)時(shí)作為參考。4. 模型的精細(xì)分析和改進(jìn)（1） 殘差分析殘差，是各觀測值與回歸方程所對應(yīng)得到的擬合值之差，實(shí)際上，它是線性回歸模型中誤差的估計(jì)值。即有零均值和常值方差，利用殘差的這種特性反過來考察原模型的合理性就是殘差分析的基本思想。利用MATLAB進(jìn)行殘差分析則是通過殘差圖或時(shí)序殘差圖。殘差圖是指以殘差為縱坐標(biāo)，以其他指定的量為橫坐標(biāo)的散點(diǎn)圖。主要包括：（1）橫坐標(biāo)為觀測時(shí)間或觀測值序號；（2）橫坐標(biāo)為某個(gè)自變量的觀測值；（3）橫坐標(biāo)為因變量的擬合值。通過觀察殘差

18、圖，可以對奇異點(diǎn)進(jìn)行分析，還可以對誤差的等方差性以及對回歸函數(shù)中是否包含其他自變量、自變量的高次項(xiàng)及交叉項(xiàng)等問題給出直觀的檢驗(yàn)。以觀測值序號為橫坐標(biāo)，殘差為縱坐標(biāo)所得到的散點(diǎn)圖稱為時(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ù)的平方，假設(shè)檢驗(yàn)統(tǒng)計(jì)量，概率，分別為：0.9533 ； 115.5586 ； 0.0

19、000 ，比較可知R，F(xiàn)均增加模型得到改進(jìn)。 圖8.2 時(shí)序殘差圖（2） 變量間的交互作用討論變量間的交互作用包括：不同自變量之間的交互作用以及同一變量的自相關(guān)性。不同自變量之間的交互作用：有時(shí)，在實(shí)驗(yàn)中不僅單因素對指標(biāo)有影響，而且因素間還會(huì)聯(lián)合起來對指標(biāo)產(chǎn)生影響，常稱這種聯(lián)合作用為交互作用。處理兩個(gè)因素間交互作用的一個(gè)簡單辦法是加入這兩個(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)，其他

20、指標(biāo)也不理想，因此，本例中其交互作用并不顯著，該模型不如前面兩個(gè)模型好。自相關(guān)性的診斷和處理：若數(shù)據(jù)是以時(shí)間為序的，稱為時(shí)間序列數(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)法是通過繪制殘差散點(diǎn)圖觀察，如果散布點(diǎn)大部分點(diǎn)落在第，象限，表明存在著正的序列相關(guān)；如果大部分點(diǎn)落在第，象限，表明存在著負(fù)的序列相關(guān)。對DW檢驗(yàn)法可以利用MATLAB軟件編程計(jì)算統(tǒng)計(jì)量：，然后查閱DW檢驗(yàn)上下界表，以決定模型的自相關(guān)狀態(tài)。當(dāng)

21、一個(gè)回歸模型存在序列相關(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è)對因變量影響最大的，再對集合中的變量進(jìn)行檢驗(yàn)，從變得不顯著的變量中移出一個(gè)影響最小的，依此進(jìn)行，直到不能引入和移出為止。引入和移出都以給定的顯著性水平為標(biāo)準(zhǔn)。MATLAB統(tǒng)計(jì)工

22、具箱中逐步回歸的命令是stepwise，它提供了一個(gè)人機(jī)交互式畫面，通過此工具可以自由地選擇變量進(jìn)行統(tǒng)計(jì)分析。該命令的用法是：stepwise(X , Y , inmodel , alpha)其中X是自變量數(shù)據(jù)，排成矩陣（m為自變量個(gè)數(shù)，n為每個(gè)變量的數(shù)據(jù)量），Y是因變量數(shù)據(jù)，排成向量，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 Pl

23、ot窗口中的虛線表示回歸系數(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對應(yīng)的概率。關(guān)于本節(jié)案例2，如果引入新的自變量 . 也可以采用逐步回歸法解決,源程序如下：A=3.5 5.3 5.1

24、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 38.7 46.8 41.4 37.5 39.0 40.7 30.1 52.9

25、 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)式回歸仍然屬于多元線性回歸，可以是一元多項(xiàng)式回歸或多元多項(xiàng)式回歸。一元多項(xiàng)式回歸模型的一般形式為用MATLA

26、B求解一元多項(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ì)得分的關(guān)系，我們查閱了國家統(tǒng)計(jì)局資料，北京體育大學(xué)出版社出版的2000國民體質(zhì)監(jiān)測報(bào)告，表8-4是我國大陸3

27、1個(gè)省市的有關(guān)數(shù)據(jù)。我們希望通過這幾組數(shù)據(jù)考察它們是否具有良好的相關(guān)關(guān)系，并通過它們的關(guān)系從人均國內(nèi)生產(chǎn)總值（可以看作反映生活水平的一個(gè)指標(biāo)）、體質(zhì)得分預(yù)測其壽命可能的變化范圍。體質(zhì)是指人體的質(zhì)量，是遺傳性和獲得性的基礎(chǔ)上表現(xiàn)出來的人體形態(tài)結(jié)構(gòu)，生理機(jī)能和心理因素綜合的、相對穩(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ā)展水平、心理發(fā)育水平和適應(yīng)能力等五個(gè)方面。目前，體質(zhì)的綜合評價(jià)主要是形態(tài)、機(jī)能和身體素質(zhì)三類指標(biāo)按一定的權(quán)重進(jìn)行換算而

28、得。 表8-4 31個(gè)省市人口預(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.125100601565.9662.953822676.1069.34547319573.9169.9929931167

29、2.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.1390882173.1067.065143351174.7069 .505337722267.4763.

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

31、和純交互項(xiàng)）quadratic（包含線性項(xiàng)和完全二次項(xiàng)）輸出一個(gè)交互式畫面，對例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 3377

32、2 8744 11494 20461 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.4

33、95 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è)交互式畫面左邊一幅圖形是固定時(shí)的曲線及其置信區(qū)間，右邊一幅圖形是固定時(shí)的曲線及其置信區(qū)間。移動(dòng)鼠標(biāo)可改變，的值，同時(shí)圖左邊給出的預(yù)測值及其置信區(qū)間。如輸入=，=66.165，則=70.6948，其置信區(qū)間70.69481.1079。圖的左下方有兩個(gè)下拉式菜單，上面的菜單Export用于輸出數(shù)據(jù)（包括：回歸系數(shù)parameters,殘差residuals,剩余標(biāo)準(zhǔn)差RMSE等）, 在MATLAB工作空間中

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

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

36、高斯牛頓法求解的。將上面的輸出作為命令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é)簡稱酶動(dòng)力學(xué)，主要研究酶促反應(yīng)速度與底物（即反應(yīng)物）濃度以及其它因素的關(guān)系。在底物濃度很低時(shí)酶促反應(yīng)是一級反應(yīng)；當(dāng)?shù)孜餄舛忍幱谥虚g范圍時(shí)，是混合級反應(yīng)；當(dāng)?shù)孜餄舛仍黾訒r(shí)，向零級反應(yīng)過渡。某生化系學(xué)生為了研究嘌呤霉素在某項(xiàng)酶促反應(yīng)中對反應(yīng)速度與底物濃度之間關(guān)系的影響，設(shè)計(jì)了兩個(gè)實(shí)驗(yàn)，一個(gè)實(shí)驗(yàn)中所使用的酶是經(jīng)過嘌呤霉素處理的

37、，而另一個(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ù)底物濃度（ppm）0.020.060.110.220.561.10反應(yīng)速度未處理6751848698115131124144158160/處理764797107123139159152191201207200分析與假設(shè)記酶促反應(yīng)的速度為,底物濃度為，二者之間的關(guān)系寫作,其中為參數(shù)(可為一向量)。由酶促反應(yīng)的基本性質(zhì)可知，當(dāng)?shù)孜餄舛群艿蜁r(shí)酶促反應(yīng)是一級反應(yīng)，此時(shí)反應(yīng)速

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

39、中的命令進(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ù)測誤差的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.

40、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用于給出一個(gè)交互式畫面，可以得到因變量y的預(yù)測值和預(yù)測區(qū)間，左下方的Export可向工作區(qū)傳送剩余標(biāo)準(zhǔn)差等數(shù)據(jù)。表8-6 模型參數(shù)的估計(jì)結(jié)果參 數(shù)參 數(shù) 估 計(jì) 值置 信 區(qū) 間2126818197.2028 228.1608006410.0457 0.0826從上面的結(jié)果可以知道,對經(jīng)過嘌呤霉素處理的實(shí)驗(yàn)數(shù)據(jù),在用Michaelis-Menten模型進(jìn)行回歸分析時(shí),最終反應(yīng)速度為=212.6818,反應(yīng)的半速度點(diǎn)(達(dá)到最終反應(yīng)速度的一半時(shí)的底物濃度x值)恰為=0.06412。混合反應(yīng)模型由酶動(dòng)力學(xué)知識(shí)我們知道,酶促反應(yīng)的濃度依賴于底物濃度,并且可以假定,嘌呤霉素的處理會(huì)影響最