回歸分析方法

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

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

4、程中如何有效地使用MATLAB軟件。沒有學過這部分數(shù)學知識的讀者可以不深究其數(shù)學原理，只要知道回歸分析的目的，按照相應方法通過軟件顯示的圖形或計算所得結(jié)果表示什么意思，那么，仍然可以學到用回歸模型解決實際問題的基本方法。包括：一元線性回歸、多元線性回歸、非線性回歸、逐步回歸等方法以及如何利用MATLAB軟件建立初步的數(shù)學模型，如何透過輸出結(jié)果對模型進行分析和改進，回歸模型的應用等。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個相應的值，確定的方法是根據(jù)最小二乘準則，要使取最小值。利用極值必要條件令，求的估計值，從而得到回歸直線。只不過這個過程可以由軟件通過直線擬合完成，而無須進行繁雜的運算。（1）參數(shù)的區(qū)間估計由于我們所計算出的仍然是隨機變量，因此要對取值的區(qū)間進行估計，如果區(qū)間估計值是一個較短的區(qū)間表示模型精度較高。（2）對誤差方差的估計設(shè)為回歸函數(shù)

6、的值，為測量值，殘差平方和剩余方差（3）線性相關(guān)性的檢驗由于我們采用的是一元線性回歸，因此，如果模型可用的話，應該具有較好的線性關(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 92 93 93 95

9、 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.6001F=180.9531，p0.0001，可以通過殘差圖發(fā)現(xiàn)，第二個數(shù)據(jù)為奇異數(shù)據(jù)，去掉該數(shù)據(jù)后運行后得到b = 17.654

10、9 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.6672F=261.6389，pm）得n組觀察值，采用最小二乘估計求得回歸方程.建立回歸模型是一個相當復雜的過程，概括起來主要有以下幾個方面工作（1）根據(jù)研究目的收集數(shù)據(jù)和預分析；（2）根據(jù)散點圖是否具有線性關(guān)系建立基本回歸模型；（3）模型的精細分析；（4）模型的確認與應用等。收集數(shù)據(jù)的一個經(jīng)驗準則是收集的數(shù)據(jù)量（樣本容量）至少應為可能的自變量數(shù)目的610倍。在建模

11、過程中首先要根據(jù)所研究問題的目的設(shè)置因變量，然后再選取與該因變量有統(tǒng)計關(guān)系的一些變量作為自變量。我們當然希望選擇與問題關(guān)系密切的變量，同時這些變量之間相關(guān)性不太強，這可以在得到初步的模型后利用MATLAB軟件進行相關(guān)性檢驗。下面通過一個案例探討MATLAB軟件在回歸分析建模各個環(huán)節(jié)中如何應用。多元線性回歸的MATLAB實現(xiàn) 仍然用命令regress(y , X)，只是要注意矩陣X的形式，將通過如下例子說明其用法。8.2.2 某類研究學者的年薪1. 問題例2 工薪階層關(guān)心年薪與哪些因素有關(guān)，以此可制定出它們自己的奮斗目標。某科學基金會希望估計從事某研究的學者的年薪Y(jié)與他們的研究成果(論文、著作等

12、)的質(zhì)量指標X1、從事研究工作的時間X2、能成功獲得資助的指標X3之間的關(guān)系，為此按一定的實驗設(shè)計方法調(diào)查了24位研究學者，得到如下數(shù)據(jù)（i為學者序號）：表8-2 從事某種研究的學者的相關(guān)指標數(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

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

14、),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.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 1

15、3 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.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,

16、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.39841.28350.6691 1.8979=0.9106 F=67.9195 p= 3.10 （查F分布表或輸入命令finv(0.95,3,20)計算）。（）p值檢驗：若（為預定顯著水平），則說明因變量與自變量之間顯著地有線性相關(guān)關(guān)系。本例輸出結(jié)果，p0.0001

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

18、因變量的擬合值。通過觀察殘差圖，可以對奇異點進行分析，還可以對誤差的等方差性以及對回歸函數(shù)中是否包含其他自變量、自變量的高次項及交叉項等問題給出直觀的檢驗。以觀測值序號為橫坐標，殘差為縱坐標所得到的散點圖稱為時序殘差圖，畫出時序殘差圖的MATLAB語句為rcoplot(r,rint)（圖8.2）。可以清楚看到殘差大都分布在零的附近，因此還是比較好的 ，不過第4、12、19這三個樣本點的殘差偏離原點較遠，如果作為奇異點看待，去掉后重新擬合，則得回歸模型為：且回歸系數(shù)的置信區(qū)間更小均不包含原點，統(tǒng)計變量stats包含的三個檢驗統(tǒng)計量：相關(guān)系數(shù)的平方，假設(shè)檢驗統(tǒng)計量，概率，分別為：0.9533 ；

19、115.5586 ； 0.0000 ，比較可知R，F(xiàn)均增加模型得到改進。 圖8.2 時序殘差圖（2） 變量間的交互作用討論變量間的交互作用包括：不同自變量之間的交互作用以及同一變量的自相關(guān)性。不同自變量之間的交互作用：有時，在實驗中不僅單因素對指標有影響，而且因素間還會聯(lián)合起來對指標產(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 。但它

20、們的置信區(qū)間均包含原點，其他指標也不理想，因此，本例中其交互作用并不顯著，該模型不如前面兩個模型好。自相關(guān)性的診斷和處理：若數(shù)據(jù)是以時間為序的，稱為時間序列數(shù)據(jù)。在時間序列數(shù)據(jù)中，同一變量的順序觀測值之間出現(xiàn)的相關(guān)現(xiàn)象稱為自相關(guān)。一旦數(shù)據(jù)中存在這種自相關(guān)序列，如果仍采用普通的回歸模型直接處理，將產(chǎn)生不良后果，使預測失去意義。自相關(guān)的診斷主要有圖示檢驗法、相關(guān)系數(shù)法和DW檢驗法。圖示檢驗法是通過繪制殘差散點圖觀察，如果散布點大部分點落在第，象限，表明存在著正的序列相關(guān)；如果大部分點落在第，象限，表明存在著負的序列相關(guān)。對DW檢驗法可以利用MATLAB軟件編程計算統(tǒng)計量：，然后查閱DW檢驗上下界表

21、，以決定模型的自相關(guān)狀態(tài)。當一個回歸模型存在序列相關(guān)性時，首先要查明序列相關(guān)產(chǎn)生的原因。如果是回歸模型選用不當，則應改用適當?shù)幕貧w模型；如果是缺少重要的自變量，則應增加自變量；如果以上方法都不能消除序列相關(guān)性，則需要采用差分法、迭代法等處理，更詳細內(nèi)容參見相關(guān)概率統(tǒng)計參考文獻。8.2.3 逐步回歸方法建模 逐步回歸就是一種從眾多自變量中有效地選擇重要變量的方法。逐步回歸的基本思路是，先確定一個包含若干自變量的初始集合，然后每次從集合外的變量中引入一個對因變量影響最大的，再對集合中的變量進行檢驗，從變得不顯著的變量中移出一個影響最小的，依此進行，直到不能引入和移出為止。引入和移出都以給定的顯著性

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

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

24、：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 38.7 46.8 41.4 37.5 39.0

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

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

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

28、三類指標按一定的權(quán)重進行換算而得。 表8-4 31個省市人口預期壽命與人均國內(nèi)生產(chǎn)總值和體質(zhì)得分數(shù)據(jù)序號預期壽命體質(zhì)得分人均產(chǎn)值序號預期壽命體質(zhì)得分人均產(chǎn)值序號預期壽命體質(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.

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

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

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

32、 13255 9088 33772 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

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

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

35、，其中對是非線性的，為了估計參數(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是殘差，J是用于估計預測誤差的

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

37、中所使用的酶是經(jīng)過嘌呤霉素處理的，而另一個實驗所用的酶是未經(jīng)嘌呤霉素處理的。所得實驗數(shù)據(jù)見表8-5。試根據(jù)問題的背景和這些數(shù)據(jù)建立一個合適的數(shù)學模型，來反映這項酶促反應的速度與底物濃度以及嘌呤霉素處理與否之間的關(guān)系。表8-5 嘌呤霉素實驗中的反應速度與底物濃度數(shù)據(jù)底物濃度（ppm）0.020.060.110.220.561.10反應速度未處理6751848698115131124144158160/處理764797107123139159152191201207200分析與假設(shè)記酶促反應的速度為,底物濃度為，二者之間的關(guān)系寫作,其中為參數(shù)(可為一向量)。由酶促反應的基本性質(zhì)可知，當?shù)孜餄舛群艿?/p>

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

39、求解可利用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為用于估計預測誤差的Jacobi矩陣。參數(shù)beta的置信區(qū)間用命令nlparci(beta,R,J)得到。首先建立函數(shù)M文件huaxue.m，非線性模型參數(shù)估計的源程序如下：x=0.02 0.02 0.06 0.06 0.11 0.11 0.22 0.22 0

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