SAS講義 第三十四課非線性回歸分析

1、第三十四課 非線性回歸分析現(xiàn)實(shí)世界中嚴(yán)格的線性模型并不多見，它們或多或少都帶有某種程度的近似；在不少情 況下，非線性模型可能更加符合實(shí)際。由于人們?cè)趥鹘y(tǒng)上常把“非線性”視為畏途，非線性 回歸的應(yīng)用在國(guó)內(nèi)還不夠普及。事實(shí)上，在計(jì)算機(jī)與統(tǒng)計(jì)軟件十分發(fā)達(dá)的令天，非線性回歸 的基本統(tǒng)計(jì)分析已經(jīng)與線性回歸一樣切實(shí)可行。在常見的軟件包中(諸如SAS、SPSS等等), 人們已經(jīng)可以像線性回歸一樣，方便的對(duì)非線性回歸進(jìn)行統(tǒng)計(jì)分析。因此，在國(guó)內(nèi)回歸分析 方法的應(yīng)用中，已經(jīng)到了“更上一層樓”，線性回歸與非線性回歸同時(shí)并重的時(shí)候。對(duì)變量間非線性相關(guān)問題的曲線擬合，處理的方法主要有：首先決定非線性模型的函數(shù)類型，對(duì)于

2、其中可線性化問題則通過變量變換將其線 性化，從而歸結(jié)為前面的多元線性回歸問題來解決。若實(shí)際問題的曲線類型不易確定時(shí)，由于任意曲線皆可由多項(xiàng)式來逼近，故?？?用多項(xiàng)式回歸來擬合曲線。若變量間非線性關(guān)系式已知(多數(shù)未知)，且難以用變量變換法將其線性化，則進(jìn) 行數(shù)值迭代的非線性回歸分析。一、可變換成線性的非線性回歸在實(shí)際問題中一些非線性回歸模型可通過變量變換的方法化為線性回歸問題。例如，對(duì) 非線性回歸模型(34.1)y =a + (a cosix + b sinix )+s(34.1)0it itti=1即可作變換二 cos x , x 二 sin x , x 二 cos 2x , x 二 sin

3、2xt 二 cos x , x 二 sin x , x 二 cos 2x , x 二 sin 2xt 2tt3tt4tt將其化為多元線性回歸模型。一般地，若非線性模型的表達(dá)式為y 二 b + b g (jc )+ b g (jc )F b g (jc )t則可作變量變換x1t011 t(34.2)x *1t二 g (x ) x *1 t2t*mt二 g (x )m t(34.3)將其化為線性回歸模型的表達(dá)式，從而用前面線性模型的方法來解決，其中(34.3)中的xt也 可為自變量構(gòu)成的向量。這種變量變換法也適用于因變量和待定參數(shù)b.。如iyt這種變量變換法也適用于因變量和待定參數(shù)b.。如iyt+

4、 b x22t-1)31t 2t(34.4)時(shí)上式兩邊取對(duì)數(shù)得ln yt=ln a ln yt=ln a + b x + b x1 1t2 2t31t 2t-1)(34.5)現(xiàn)作變換(34.6)y* = ln y , b = ln a, x = x x -1 tt 03t 1t 2t(34.6)則可得線性表達(dá)式y(tǒng) * 二 b + b x + b x + b TOC o 1-5 h z t 01 1t2 2t3 3t利用前面方法確定了b ,i = 0,1,2,3，并由a = exp(b )得到a的值。i0變量變換的線性化方法可推廣到下列形式的非線性模型h(y )二 c (b ) + c (b )

5、g (x )+ + c (b )g (x )(34.8)t00111 tmmmt其中X=(X,x2，,x )，而h (y)、c. (b.)、g. (x)則分別化為新的因變量、線性回歸參數(shù)12ptiiit和自變量，即可歸結(jié)為線性回歸模型來解。見表34.1所示給出了一些常見的可線性化的非 線性模型。表341典型的函數(shù)及線性化方法函數(shù)名稱函數(shù)表達(dá)式線性化方法雙曲線函數(shù)1b=a + yx11v = u yx冪函數(shù)y 二 axbv = In y u = ln x指數(shù)函數(shù)y 二 aebxv = ln y u = xy 二 aeb/xv = ln y u =x對(duì)數(shù)函數(shù)y = a + b In xv = yu

6、 = ln xS型函數(shù)1y 二 a + be - x1v = u = e - xy當(dāng)曲線的函數(shù)類型未確定時(shí)，我們常采用上述非線性模型作為其擬合曲線，即將自變量 的各種初等函數(shù)的組合作為新自變量，用逐步回歸法(或正交篩選法等)對(duì)新變量進(jìn)行篩選, 以確定一個(gè)項(xiàng)數(shù)不多的線性函數(shù)表達(dá)式。該方法對(duì)表達(dá)式形式?jīng)]限制且精度要求不高的問題 頗為有效。二、多項(xiàng)式回歸分析在式(34.2 )中，若取印匸”，則為多項(xiàng)式回歸模型。由數(shù)學(xué)分析知識(shí)可知，一般函數(shù)Ii1都可用多項(xiàng)式來逼近，故多項(xiàng)式回歸分析可用來處理相當(dāng)廣泛的非線性問題。 對(duì)觀測(cè)數(shù)據(jù)(x，y) (t= 1,，N)，多項(xiàng)式回歸模型為t J ty 二 b + b

7、x + b x2 + b xm +e ，仁 1,2，,Nt0 1_t 2_tm_tJ， ，Y =-I y2-J11 X Y =-I y2-J11 X =.1x1x2XNX 21X 22Xm1Xm2XmN則模型可表示為Y 二 XB + s當(dāng)X列滿秩時(shí)，由前面的討論知，其最小二乘估計(jì)為B =(X X )-1 X Y由此即可求得其多項(xiàng)式回歸方程。但由于(XX人的計(jì)算既復(fù)雜又不穩(wěn)定，故我們一般采 用正交多項(xiàng)式法來進(jìn)行多項(xiàng)式回歸。不可變換成線性的非線性回歸分析假設(shè)因變量y與自變量(x1，x2,，xp)之間滿足非線性模型(34.9):,X，,X ； 0 )+ (34.9)12p 其中卩=(0,0，0 )為

8、未知參數(shù)，F(xiàn)為已知表達(dá)式，僅0未知的非線性函數(shù)，為誤差項(xiàng)?，F(xiàn)將觀察數(shù)據(jù)(誤差項(xiàng)?，F(xiàn)將觀察數(shù)據(jù)(y , x , x，, xt It2tptt=l,2,N代人上式(34.9)得非線性回歸模型代人上式(34.9)得非線性回歸模型y 二Fk ,X，,Xt1t 2tpt；0)+匕t=l,2,N常記為其中Y = (y ,y ,y )為y的觀察向量，0二，,0 )為非線性回歸系數(shù)，E12N1m=(1，。，)為觀察誤差向量，F(xiàn)為未知參數(shù)0的函數(shù)向量。非線性回歸分析就是利用 1 2最小二乘準(zhǔn)則來估計(jì)回歸系數(shù)0，即求同 使得殘差平方和Q(0 )= - EE = - (Y - F(0)(Y F(0 )在0 = 0處

9、達(dá)到最小。非線性回歸分析一般來用數(shù)值迭代法來進(jìn)行，其共同特點(diǎn)是：由選定0的初值冋出發(fā),通過逐步迭代即選擇適當(dāng)?shù)牟介L(zhǎng)t ( 0 )及確定搜索方向向量A=(A1，A2，A )，使得12m(5.4.11)QW qC 0 (5.4.11)再P由取代因，重復(fù)上述迭代過程，直至Q(P)可認(rèn)為達(dá)到最小值為止，即可將所得的P 作為其最小二乘估計(jì)同，從而得到非線性回歸方程1( )y = fX , x ,x ; p 12p1.下降方向和步長(zhǎng)的選擇首先考察Q )= 2 E E=(y -F (p) (y -1.下降方向和步長(zhǎng)的選擇首先考察Q )= 2 E EaF其中G =araQaaF其中G =araQapf aF、

10、ap(Y - F(P)=-Gr(Y - F(p )為F的梯度矩陣。ap (舌那 1迭代收斂到同，其迭代公式應(yīng)滿足下降性質(zhì)(5.4.11)?，F(xiàn)考慮一元函數(shù) )=qG 0 +1 - A j，它從因出發(fā)以A為方向的射線上取值。由復(fù)合求導(dǎo)公式得t=0Vq、鄰丿-A = -(Y - F (卩 G-AAA為方向向量的射線上可以找到p = p0 +1-A，使得可以證明，當(dāng)d0時(shí)，在以Q(p)v qG R 我們將滿足d0的A稱為下降方向，Bard于1974年給出了 A為下降方向的充要條件為A = PGf(Y - F(p )其中P為對(duì)稱正定陣，由此我們可得下降算法的迭代公式為(34.12)p = p o + t

11、PG - F (p)(34.12)其中P為任意正定陣，G為F的梯度，t為滿足Q(p)v 其中P為任意正定陣，G為F的梯度，t如何計(jì)算A以便修改參數(shù)向量P有五種常用的非線性回歸迭代方法：高斯-牛頓法(Gauss-Newton)、最速下降法(梯度法，Gradient)、牛頓法(Newton)、麥夸特法(Marquardt)、 正割法(DUD)。以下我們介紹其中高斯-牛頓法。GaussNewton 法首先選取0的一切初始近似值阿，令A(yù)=0-0 0，則只要確定A的值即可確定0。為A = F G o)+ G A此，考慮F(0 A = F G o)+ G AF(0)= F C +A)= F(00)+30

12、0=0其中G = 其中G = 0F00 I 0=00為F的梯度。此時(shí)其殘差平方和(34.13)(34.14)由此即可用前面線性回歸法求A，只需將G(34.13)(34.14)由此即可用前面線性回歸法求A，只需將G、Y - F(0 0)視為前面(5.2.1)式中的X、Y=或qG0 +a)2時(shí)，即得的0最小，當(dāng) max a . ii代替網(wǎng)，重復(fù)上述步驟，直至或Q(0)滿i足精度要求為止。該法稱為GaussNewton法，其一般迭代公式為即可。此時(shí)，對(duì)給定精度 1、* 2A入二乘法估計(jì)卩=卩0 + A ；否則用所得的0Q =丄 C - F G 0)- G A)Y - F G 0)- G A 20A

13、= 0，得其A的正則方程為(GG)A 二 G&- F4 A =(GG)-1 G&-FC (34.15)0 i+1 = 0 i + t A(34.15)i的最小值點(diǎn)。其中：a為 G G i L G i h=g G i - f G i的解,t.為啊)=qG+t 胡 i的最小值點(diǎn)。Gauss-Newton 法在初值網(wǎng)選取適當(dāng)，且G G可逆時(shí)非常有效，但在其他情形，其求解 較為困難，對(duì)此，Marguardt對(duì)(34.14)中4的正則系數(shù)陣作適當(dāng)修正，得到了改進(jìn)算法。四、nlin非線性回歸過程在很多場(chǎng)合，可以對(duì)非線性模型進(jìn)行線性化處理，尤其是關(guān)于變量非線性的模型，以運(yùn) 用OLS進(jìn)行推斷。對(duì)線性化后的線性

14、模型，可以應(yīng)用SAS的reg過程進(jìn)行計(jì)算。多項(xiàng)式模型可以直接應(yīng)用glm (廣義線性模型)求解。對(duì)于不能線性化的非線性模型。 其估計(jì)不能直接運(yùn)用經(jīng)典的最小二乘法，而需要運(yùn)用其他估計(jì)方法，如直接搜索法、直接最 優(yōu)法與Taylor級(jí)數(shù)展開法進(jìn)行線性逼近。此時(shí)，可以利用SAS/STAT的nlin過程實(shí)現(xiàn)相應(yīng)的 計(jì)算。1. proc nlin 過程proc nlin 采用最小誤差平方法(Least Squares Method)及循環(huán)推測(cè)法(Iterative Estimation Method)來建立一個(gè)非線性模型。一般而言，用戶必須自訂參數(shù)的名字、參數(shù)的啟動(dòng)值(starting value)、非線性

15、的模型與循環(huán)推測(cè)法所用的準(zhǔn)則。若用戶不指明，則nlin程序自動(dòng)以高斯- 牛頓迭代法(Gauss-Newton iterative procedure)為估計(jì)參數(shù)的方法。另外此程序也備有掃描 (Grid search)的功能來幫助讀者選擇合適的參數(shù)啟動(dòng)值。由于非線性回歸分析十分不易處理, nlin程序不保證一定可以算出符合最小誤差平方法之標(biāo)準(zhǔn)的參數(shù)估計(jì)值。nlin過程的功能，計(jì)算非線性模型參數(shù)的最小二乘估計(jì)LS及加權(quán)最小二乘估計(jì)。與reg 過程不同的是：模型的參數(shù)要命名、賦初值、求偏導(dǎo)數(shù)；model語句與參數(shù)名、解釋變量的 表達(dá)式有關(guān)；可以使用賦值語句及條件語句。nliproc nlindata

16、=數(shù)據(jù)集/選項(xiàng)列表；parameters 參數(shù)名=數(shù)值；model因變量=表達(dá)式/選項(xiàng)列表;bounds表達(dá)式；der.參數(shù)名.參數(shù)名=表達(dá)式；id變量列表；outputout=數(shù)據(jù)集/選項(xiàng)列表；by變量列表；其中，parameters語句和model語句是必需的，而其余語句供用戶根據(jù)需要選擇。proc nlin語句中的主要選擇項(xiàng)。0%他時(shí)=數(shù)據(jù)集名定存放參數(shù)估計(jì)的每步迭代結(jié)果的數(shù)據(jù)集名。best=n一 求過程只輸出網(wǎng)格點(diǎn)初始值可能組合中最好的n組殘差平方和。method= gauss I marquardt I newtonl gradientl dud I設(shè)定參數(shù)估計(jì)的迭代方法。缺省時(shí)為g

17、auss,除非沒有der語句。eformat 求所有數(shù)值以科學(xué)記數(shù)法輸出。nopoint 制打印輸出。noinpoint 制迭代結(jié)果的輸出。parameters (parms) 語句。用于對(duì)所有參數(shù)賦初值，項(xiàng)目之間以空格分隔。例如，parms b0=0b1=1 to 10b2=1to 10 by 2b3=1,10,100；model 語句。表達(dá)式可以是獲得數(shù)值結(jié)果的任意有效SAS表達(dá)式。這個(gè)表達(dá)式包括參數(shù)名字、輸入數(shù) 據(jù)集中的變量名以及在nlin過程中用程序設(shè)計(jì)語句創(chuàng)建的新變量。例如，model y=b0*(1 exp(-bl*x)；bounds 語句。用于設(shè)定參數(shù)的約束，主要是不等式約束，約

18、束間用逗號(hào)分隔。例如，bounds a30， 1=c=10；der.語句。除非在proc nlin語句中指明所用的迭代法是dud,使用選擇項(xiàng)method=dud，否則der語 句是必需的。der.語句用于計(jì)算模型關(guān)于各個(gè)參數(shù)的偏導(dǎo)數(shù)，相應(yīng)的格式為：一階偏導(dǎo)數(shù)der.參數(shù)名=表達(dá)式；二階偏導(dǎo)數(shù)der.參數(shù)名.參數(shù)名=表達(dá)式；例如，對(duì)于 model y=bO*(l exp(-bl*x) ； der.語句 的書寫格式為 der. bO=l exp(-bl*x) ； der. bl=bO*x*exp(-bl*x)；對(duì)于多數(shù)算法，都必須對(duì)每個(gè)被估計(jì)的參數(shù)給出一階偏導(dǎo)數(shù)表達(dá)式。對(duì)于new ton法，必 須

19、給出一、二階偏導(dǎo)數(shù)表達(dá)式。例如，二階偏導(dǎo)數(shù)表達(dá)式為，der.b0.b0.=0 ； der.bO.bl二x*exp(bl*x) ； der.bl.bl二der.bl*x；7. output 語句用于把一些計(jì)算結(jié)果輸出到指定的數(shù)據(jù)集中。有關(guān)的關(guān)鍵字及其意義如表34.2所示: 表 342 nlin 過程中output語句的關(guān)鍵字關(guān)鍵字意義關(guān)鍵字意義關(guān)鍵字意義predicted | p預(yù)測(cè)值stdpclm的標(biāo)準(zhǔn)差u9595%cli上限r(nóng)esidual | r殘差stdr殘差的標(biāo)準(zhǔn)差l9595%cu下限parms參數(shù)估計(jì)值l95m95%clm下限student學(xué)生氏殘差sse | ess殘差平方和u95

20、m95%clm上限h杠桿點(diǎn)統(tǒng)計(jì)量關(guān)于nlin過程的其他選擇項(xiàng)及意義，詳見SAS/STAT的用戶手冊(cè)。五、實(shí)例分析例341負(fù)指數(shù)增長(zhǎng)曲線的非線性回歸。根據(jù)對(duì)已有數(shù)據(jù)的XY散點(diǎn)圖的觀察和分析, 發(fā)現(xiàn)Y隨X增長(zhǎng)趨勢(shì)是減緩的，并且Y趨向一個(gè)極限值，我們認(rèn)為用負(fù)指數(shù)增長(zhǎng)曲線 y = b (1 - e-b1x)來擬合模型較為合適。程序如下：0data expd ;input x y ;cards;0200.570300.720400.810500.870600.910700.940800.950900.971000.981100.991201.001300.991400.991501.001601.001

21、700.991801.001901.002000.992101.00proc nlin data二expd best=10 method二gauss;parms b0=0 to 2 by 0.5 b1=001 to 0.09 by 0.01;model y=b0*(1exp(b1*x);derb0=1exp(b1*x);der b1二b0 *x* exp(b1 *x);output out二expout p=ygs ;run;goptions reset二global gunit=pct cback二white borderhtitle=6 htext=3 ftext=swissb color

22、s=(back);proc gplot data=expout;plot y*x ygs*x /haxis=axis1 vaxis=axis2 overlay;symboll i=none v=plus cv=red h=2.5 w=2;symbol2 i=join v=none l=1 h=2.5 w=2;axisl order=20 to 210 by 10;axis2 order=0.5 to 1.1 by 0.05;titlel y=b0*（1-exp（-b1*x）；title2 proc nlin method=gauss;run ;程序說明：由于在nlin過程中使用選項(xiàng)method

23、=gauss,即指定用高斯-牛頓迭代算法來尋 找model語句中非線性表達(dá)式y(tǒng)二b （1 -e-卡）中參數(shù)b和b的最小二乘估計(jì)。我們知道參 0 0 1數(shù)初始值選取好壞，對(duì)迭代過程是否收斂影響很大oparms語句設(shè)置了初始值網(wǎng)格值為b取0, 0.5, 1, 1.5, 2共5個(gè)值，耳取0.01, 0.02,，0.09共9個(gè)值，所有可能組合為5X9=45, 選項(xiàng)best=10要求輸出殘差平方和最小的前10種組合。高斯-牛頓迭代算法要求給出模型 y = b （1 - e-x）對(duì)參數(shù)b和b的一階偏導(dǎo)數(shù)表達(dá)式，我們知道d yd b0d yd b0d y d b1=1 e -bx=b xe -qx0der.

24、語句用以表示上面兩個(gè)一階偏導(dǎo)數(shù)表達(dá)式。output語句輸出一個(gè)新數(shù)據(jù)集expout，包 括原數(shù)據(jù)集和非線性回歸模型的預(yù)測(cè)值ygso gplot過程的主要作用是繪制輸出數(shù)據(jù)集expout 中的原始數(shù)據(jù)的散點(diǎn)圖及回歸曲線的平滑線。程序的輸出結(jié)果見圖34-1和見表34.3所示。y= b0*(1exp(b1*X)prac hlin1,1Qr1,W0.551,1Qr1,W0.55：O.KV曲 3Q 朋閃胡 R3 口舞 1QQ 110 100 13Q 匕口 150 怡口 170 19 口 190 20口 21Q圖34-1 XY散點(diǎn)圖和非線性回歸曲線表34.3負(fù)指數(shù)增長(zhǎng)曲線:Gauss-Newton方法的輸

25、出結(jié)果y=bO*(l-exp(-bl*x)proc nlin method二gaussNon-Linear LeastSquares Grid SearchDependent VariableYB0B1Sum of Squares1.0000000.0400000.0014041.0000000.0500000.0168111.0000000.0600000.0551551.0000000.0300000.0665711.0000000.0700000.0972841.0000000.0800000.1365361.0000000.0900000.1708391.0000000.0200000

26、.4192851.5000000.0100000.9757241.0000000.0100002.165290Non-Linear Least Squares Iterative PhaseDependent Variable YMethod:Gauss-NewtonIterB0B1Sum of Squares01.0000000.0400000.00140410.9961390.0418570.00058020.9961920.0419520.00057730.9961890.0419540.00057740.9961890.0419540.000577NOTE: Convergence criterion metNon-Linear Least Squares Summary StatisticsDepend