應(yīng)用MATLAB進(jìn)行非線性回歸分析

1、應(yīng)用MATLAB進(jìn)行非線性回歸分析 摘 要早在十九世紀(jì)，英國生物學(xué)家兼統(tǒng)計學(xué)家高爾頓在研究父與子身高的遺傳問題時，發(fā)現(xiàn)子代的平均高度又向中心回歸大的意思，使得一段時間內(nèi)人的身高相對穩(wěn)定。之后回歸分析的思想滲透到了數(shù)理統(tǒng)計的其他分支中。隨著計算機的發(fā)展，各種統(tǒng)計軟件包的出現(xiàn)，回歸分析的應(yīng)用就越來越廣泛?；貧w分析處理的是變量與變量間的關(guān)系。有時，回歸函數(shù)不是自變量的線性函數(shù)，但通過變換可以將之化為線性函數(shù)，從而利用一元線性回歸對其進(jìn)行分析，這樣的問題是非線性回歸問題。下面的第一題：煉鋼廠出鋼水時用的鋼包，在使用過程中由于鋼水及爐渣對耐火材料的侵蝕，使其容積不斷增大。要找出鋼包的容積用盛滿鋼水時的質(zhì)

2、量與相應(yīng)的實驗次數(shù)的定量關(guān)系表達(dá)式，就要用到一元非線性回歸分析方法。首先我們要對數(shù)據(jù)進(jìn)行分析，描出數(shù)據(jù)的散點圖，判斷兩個變量之間可能的函數(shù)關(guān)系，對題中的非線性函數(shù)，參數(shù)估計是最常用的“線性化方法”，即通過某種變換，將方程化為一元線性方程的形式，接著我們就要對得到的一些曲線回歸方程進(jìn)行選擇，找出到底哪一個才是更好一點的。此時我們通?？刹捎脙蓚€指標(biāo)進(jìn)行選擇，第一個是決定系數(shù)，第二個是剩余標(biāo)準(zhǔn)差。進(jìn)而就得到了我們想要的定量關(guān)系表達(dá)式。第二題：給出了某地區(qū)19712000年的人口數(shù)據(jù)，對該地區(qū)的人口變化進(jìn)行曲線擬合。也用到了一元非線性回歸的方法。首先我們也要對數(shù)據(jù)進(jìn)行分析，描出數(shù)據(jù)的散點圖，然后用MA

3、TLAB編程進(jìn)行回歸分析擬合計算輸出利用 Logistic模型擬合曲線。 關(guān)鍵詞：參數(shù)估計， Logistic模型，MATLAB正文一、一元非線性回歸分析的求解思路：求解函數(shù)類型并檢驗。求解未知參數(shù)?？苫€回歸為直線回歸，用最小二乘法求解；可化曲線回歸為多項式回歸。二、回歸曲線函數(shù)類型的選取和檢驗1、直接判斷法2、作圖觀察法，與典型曲線比較，確定其屬于何種類型，然后檢驗。3、直接檢驗法（適應(yīng)于待求參數(shù)不多的情況）4、表差法（適應(yīng)于多想式回歸，含有常數(shù)項多于兩個的情況）三、化曲線回歸為直線回歸問題 用直線檢驗法或表差法檢驗的曲線回歸方程都可以通過變量代換轉(zhuǎn)化為直線回歸方程，

4、利用線性回歸分析方法可求得相應(yīng)的參數(shù)估計值。題目： 例 煉鋼廠出鋼水時用的鋼包，在使用過程中由于鋼水及爐渣對耐火材料的浸蝕，其容積不斷增大。現(xiàn)在鋼包的容積用盛滿鋼水時的重量y (kg)表示，相應(yīng)的試驗次數(shù)用x表示。數(shù)據(jù)見表，要找出y 與x的定量關(guān)系表達(dá)式。 表8.5.1 鋼包的重量y與試驗次數(shù)x數(shù)據(jù) 序號xy序號xy12106.42811110.5923108.20914110.6034109.581015110.9045109.501116110.7657110.001218111.0068109.931319111.20710110.49   1) 1/y=a

5、+b/x y=a+blnxformat longx=2 3 4 5 7 8 10 11 14 15 16 18 19;y=106.42 108.20 109.58 109.5 110 109.93 110.49 110.59 110.60 110.9 110.76 111 111.20;plot(x,y,k+);% 數(shù)據(jù)的散點圖x1=1./x;y1=1./y;plot(x1,y1,k+); %變換后數(shù)據(jù)的散點圖x2=ones(13,1) x1'b,bint,rint,stats=regress(y1',x2);z=b(1)+b(2)*x1;yc=1./z;plot(x1,y1,

6、k+,x1,z,r)%變換后數(shù)據(jù)的散點圖和回歸直線圖變換后數(shù)據(jù)的散點圖及回歸直線圖R2=1-sum(y-yc).2)/lyy;%模型的擬合優(yōu)度系數(shù)plot(x,y,k+,x,yc,r)%數(shù)據(jù)的散點圖和回歸曲線圖legend('散點圖','回歸函數(shù)')b = 0.00896662968057 0.00082917436336R2 =0.97292374957556第一種方法的程序：format longx=2 3 4 5 7 8 10 11 14 15 16 18 19;y=106.42 108.20 109.58 109.5 110 109.93 110.49

7、110.59 110.60 110.9 110.76 111 111.20;plot(x,y,k+);%數(shù)據(jù)的散點圖x1=1./x;y1=1./y;plot(x1,y1,k+); %變換后數(shù)據(jù)的散點圖x2=ones(13,1) x1'b,bint,rint,stats=regress(y1',x2);z=b(1)+b(2)*x1;yc=1./z;plot(x1,y1,k+,x1,z,r)%變換后數(shù)據(jù)的散點圖和回歸直線圖n=length(x);lyy=sum(y.2)-n*(mean(y)2;R2=1-sum(y-yc).2)/lyy;%模型的擬合優(yōu)度系數(shù)b = 0.008966

8、62968057 0.00082917436336R2 =0.97292374957556用類似的方法可以得出其它三個曲線回歸方程，它們分別是： 第二種方法的程序：format longx=2 3 4 5 7 8 10 11 14 15 16 18 19;y=106.42 108.20 109.58 109.5 110 109.93 110.49 110.59 110.60 110.9 110.76 111 111.20;x1=log(x);y1=y;x2=ones(13,1) x1'b,bint,rint,stats=regress(y1',x2);bz=b(1)+b(2)*

9、x1;yc=z;n=length(x);lyy=sum(y.2)-n*(mean(y)2;R2=1-sum(y-yc).2)/lyy;plot(x,y,'k+',x,yc,c');legend('散點圖','回歸函數(shù)')b = 1.0e+002 *R2 =0.87731500489620第三種方法的程序：format longx=2 3 4 5 7 8 10 11 14 15 16 18 19;y=106.42 108.20 109.58 109.5 110 109.93 110.49 110.59 110.60 110.9 110.76

10、 111 111.20;x1=sqrt(x);y1=y;x2=ones(13,1) x1'b,bint,rint,stats=regress(y1',x2);bz=b(1)+b(2)*x1;yc=z;n=length(x);lyy=sum(y.2)-n*(mean(y)2;R2=1-sum(y-yc).2)/lyy;plot(x,y,'k+',x,yc,'k');legend('散點圖','回歸函數(shù)')b = 1.0e+002 * 1.06301275014382三種方法的擬合效果比較：R2 =0.97292374

11、957556R2 =0.877315004896201.原始數(shù)據(jù)下表給出了某地區(qū)19712000年的人口數(shù)據(jù)（表1）。試分別用Matlab和SPSS軟件，對該地區(qū)的人口變化進(jìn)行曲線擬合。表1 某地區(qū)人口變化數(shù)據(jù)年份時間變量t=年份-1970人口y/人1971133 8151972233 9811973334 0041974434 1651975534 2121976634 3271977734 3441978834 4581979934 49819801034 47619811134 48319821234 48819831334 51319841434 49719851534 5111986

12、1634 52019871734 50719881834 50919891934 52119902034 51319912134 51519922234 51719932334 51919942434 51919952534 52119962634 52119972734 52319982834 52519992934 52520003034 527根據(jù)上表中的數(shù)據(jù)，做出散點圖，見圖1。圖1 某地區(qū)人口隨時間變化的散點圖從圖1可以看出，人口隨時間的變化呈非線性過程，而且存在一個與橫坐標(biāo)軸平行的漸近線，故可以用Logistic曲線模型進(jìn)行擬合。因為Logistic曲線模型的基本形式為：所以，只要

13、令：，就可以將其轉(zhuǎn)化為直線模型：下面，我們分別用Matlab和SPSS軟件進(jìn)行回歸分析擬合計算。2用Matlab編程進(jìn)行回歸分析擬合計算源程序（Nonlinear-Regression-Model.m），如下：clearclc% 讀入人口數(shù)據(jù)（19712000年）y = 33815 33981 34004 34165 34212 34327 34344 34458 34498 34476 34483 34488 34513 34497 34511 34520 34507 34509 34521 34513 34515 34517 34519 34519 34521 34521 34523 34

14、525 34525 34527;% 讀入時間變量數(shù)據(jù)（t年份1970）t=1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30; % 線性化處理for t = 1:30, x(t)=exp(-t); y(t)=1/y(t,1);end% 計算，并輸出回歸系數(shù)Bc=zeros(30,1)+1;X=c,x'B=inv(X'*X)*X'*yfor i=1:30,% 計算回歸擬合值 z(i)=B(1,1)+B(2,1)*x(i);% 計算離差 s(i)=y(i)-sum(y)/30;% 計算誤差 w(i)=z(i)-y(i);end% 計算離差平方和SS=s*s'% 回歸誤差平方和QQ=w*w'% 計算回歸平方和UU=S-Q;% 計算，并輸出F檢驗值F=28*U/Q% 計算非線性回歸模型的擬合值for j=1:30, Y(j)=1/(B(1,1)+B(2,1)*exp(-j);end% 輸出非線性回歸模型的擬合曲線（Logisic曲線）plot(T,Y)上述程序運行后，輸出（1）輸出回歸系數(shù)B