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

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

2、化問題則通過變量變換將其線性化，從而歸結(jié)為前面的多元線性回歸問題來解決。若實際問題的曲線類型不易確定時，由于任意曲線皆可由多項式來逼近，故?？捎枚囗検交貧w來擬合曲線。若變量間非線性關(guān)系式已知(多數(shù)未知)，且難以用變量變換法將其線性化，則進行數(shù)值迭代的非線性回歸分析。一、可變換成線性的非線性回歸在實際問題中一些非線性回歸模型可通過變量變換的方法化為線性回歸問題。例如，對非線性回歸模型2yt0_i.aicosiXtb,sinixt”(34.1)11即可作變換x1t=cosxt,x2t=sinxt,x3t=cos2xt,x4t=sin2xt將其化為多元線性回歸模型。一般地，若非線性模型的表達式為yt

3、=bbg1b2g2/-bmgmxt(34.2)則可作變量變換*xit=gixt,x2t=g2X,xmt=gmxt(34.3)將其化為線性回歸模型的表達式，從而用前面線性模型的方法來解決，其中(34.3)中的 xt也可為自變量構(gòu)成的向量。這種變量變換法也適用于因變量和待定參數(shù) bi。如yt=aexphxtb?x2tb3x1x2t-11(34.4)時上式兩邊取對數(shù)得Inyt=Inab2x2tb3xx2t一1(34.5)現(xiàn)作變換*yt=lnyt,b=Ina,x3t=xx2t-1(34.6)則可得線性表達式y(tǒng)t=bobix1tb2X2tb3x3t利用前面方法確定了?,i=0,1,2,3,并由j?=ex

4、p(b0)得到a的值。變量變換的線性化方法可推廣到下列形式的非線性模型h(yt)=co(bo)G(bi)gxcm(bm)gmXt(34.8)其中 x=(xi,x2，,xp),而 h(yt)、G(bi)、gi(為)則分別化為新的因變量、線性回歸參數(shù)和自變量，即可歸結(jié)為線性回歸模型來解。見表 34.1 所示給出了一些常見的可線性化的非線性模型。表 34.1 典型的函數(shù)及線性化方法函數(shù)名稱函數(shù)表達式線性化方法雙曲線函數(shù)1b一=a+一yxvuyx募函數(shù)by=axv=lnyu=Inx指數(shù)函數(shù)bxy=aev=Inyu=xb/xy=aein1v=Inyu=一x對數(shù)函數(shù)y=a+blnxv=yu=InxS 型函

5、數(shù)1y-，a+be1-v=u=ey當曲線的函數(shù)類型未確定時，我們常采用上述非線性模型作為其擬合曲線，即將自變量的各種初等函數(shù)的組合作為新自變量，用逐步回歸法(或正交篩選法等)對新變量進行篩選，以確定一個項數(shù)不多的線性函數(shù)表達式。該方法對表達式形式?jīng)]限制且精度要求不高的問題頗為有效。多項式回歸分析在式(34.2)中，若取 gdxx，則為多項式回歸模型。由數(shù)學分析知識可知，一般函數(shù)都可用多項式來逼近，故多項式回歸分析可用來處理相當廣泛的非線性問題。對觀測數(shù)據(jù)(xt,yj(t=1,N),多項式回歸模型為yt=bobxtb?x2bmxt，t=i,2,N(34.7)則模型可表示為Y=XB；當 X 列滿秩

6、時，由前面的討論知，其最小二乘估計為E?=XXXYi由此即可求得其多項式回歸萬程。但由于（XX）的計算既復雜又不穩(wěn)定，故我們一般采用正交多項式法來進行多項式回歸。三、不可變換成線性的非線性回歸分析假設因變量 y 與自變量（xi,X2,xp）之間滿足非線性模型y=FXI,X2,Xp；:；（34.9）其中P=#1,22,，Pm）為未知參數(shù)，F(xiàn) 為已知表達式，僅P未知的非線性函數(shù)，8為誤差項?，F(xiàn)將觀察數(shù)據(jù)（yt,Xit,X2t,，Xpt）,t=1,2,N代人上式（34.9）得非線性回歸模型ytF（Xit,X2t,Xpt；0）+%,t=i,2，N常記為Y=F）E其中Y=3。2,，yN）為y的觀察向量，

7、P=（叫,，心為非線性回歸系數(shù)，E=g,%，,心）為觀察誤差向量，F(xiàn) 為未知參數(shù)B的函數(shù)向量。非線性回歸分析就是利用最小二乘準則來估計回歸系數(shù)P,即求伊使得殘差平方和Q-gEE=JY-F-二Y-F122在B=胃處達到最小。非線性回歸分析一般來用數(shù)值迭代法來進行，其共同特點是：由選定P的初值P0出發(fā)，-yily2JNm-X1mX2一blbimXNJbm-2通過逐步迭代=-t(34.10)即選擇適當?shù)牟介L t(0)及確定搜索方向向量=(&,&,，,Am),使得Q-Q-(5.4.11)再 B 由取代 P,重復上述迭代過程，直至 Q(p)可認為達到最小值為止，即可將所得的 P作為其最小

8、二乘估計從而得到非線性回歸方程?=FXi，x2,xp;?1.下降方向和步長的選擇首先考察Q便)=3EE(Y-FfP)(Y-F(3)的梯度向量(即導數(shù))22FY-F：YF一：cPI心P)F-F5.F；F其中 G=-1=7,，一后一|為 F 的梯度矩陣。他,相:m為使 B迭代收斂到骨，其迭代公式應滿足下降性質(zhì)(5.4.11)?，F(xiàn)考慮一元函數(shù)邛(t)=Q(B+t它從 P出發(fā)以為方向的射線上取值。由復合求導公式得rd=t=坐:=YFG可以證明，當 d時，在以為方向向量的射線上可以找到 B=P+t,使得Q(P)Q(B1 我們將滿足 d的稱為下降方向，Bard 于 1974 年給出了 A 為下降方向的充要

9、條件為,PGY-F:其中 P 為對稱正定陣，由此我們可得下降算法的迭代公式為:=tPGY-Fif-ii(34.12)其中 P 為任意正定陣，G 為 F 的梯度，t 為滿足Q伊)Q倬)的正實數(shù)，即步長。如何計算以便修改參數(shù)向量 P 有五種常用的非線性回歸迭代方法：高斯-牛頓法(Gauss-Newton)、 最速下降法(梯度法， Gradient)、 牛頓法(Newton)、 麥夸特法(Marquardt)、 正割法(DUD)。以下我們介紹其中高斯-牛頓法。2.GaussNewton 法首先選取 P 的一切初始近似值00,令=P-P0,則只要確定的值即可確定 P。為此，考慮 F(P)在 P0處的

10、Taylor 展開式，并略去二次以上的項得F(P)=F(P+A)=F(B*方呼=F 伊0)+GA 萬 F 其中 G 邛為 F 的梯度。此時其殘差平方和Q=1YF0-G.=Y-F0-G2 由 0=0,得其的正則方程為必GG:.=GY-F0(34.13)故.：=GG，GY-F0(34.14)由此即可用前面線性回歸法求&只需將 G、Y-F(B0)視為前面(5.2.1)式中的 X、丫即可。此時，對給定精度叫、與，當 maxf*)5 或 Q(P0+)君2時，即得的 P 最小i二乘法估計陰=P0+；否則用所得的印代替 P,重復上述步驟，直至或 Q(B)滿足精度要求為止。該法稱為Gauss-Newt

11、on 法，其一般迭代公式為-i1-i-ti：(34.15)1procnlin 過程其中，parameters 語句和 model 語句是必需的，而其余語句供用戶根據(jù)需要選擇。2.procnlin 語句中的主要選擇項。outes 仁數(shù)據(jù)集名一一指定存放參數(shù)估計的每步迭代結(jié)果的數(shù)據(jù)集名。best=n要求過程只輸出網(wǎng)格點初始值可能組合中最好的 n 組殘差平方和。method=gauss|marquardt|newton|gradient|dud|設定參數(shù)估計的迭代方法。缺省時為 gauss,除非沒有 der.語句。eformat 要求所有數(shù)值以科學記數(shù)法輸出。nopoint 抑制打印輸出。noinp

12、oint 抑制迭代結(jié)果的輸出。3.parameters(parms)語句。用于對所有參數(shù)賦初值，項目之間以空格分隔。例如，parmsb0=0b1=1to10b2=1to10by2b3=1,10,100;4.model 語句。表達式可以是獲得數(shù)值結(jié)果的任意有效 SAS 表達式。這個表達式包括參數(shù)名字、輸入數(shù)據(jù)集中的變量名以及在 nlin 過程中用程序設計語句創(chuàng)建的新變量。例如，modely=b0*(1exp(-b1*x)；其中：為 GG3=Ge 卜F(卬)的解，ti 為甲(t)nQ(Bi+t 的最小值點。Gauss-Newton 法在初值口。選取適當，且 GG 可逆時非常有效，但在其他情形，其求

13、解較為困難，對此，Marguardt(34.14)中 A 的正則系數(shù)陣作適當修正，得到了改進算法。四、nlin非線性回歸過程在很多場合，可以對非線性模型進行線性化處理，尤其是關(guān)于變量非線性的模型，以運用 OLS 進行推斷。對線性化后的線性模型，可以應用 SAS 的 reg 過程進行計算。多項式模型可以直接應用 glm(廣義線性模型)求解。對于不能線性化的非線性模型。其估計不能直接運用經(jīng)典的最小二乘法，而需要運用其他估計方法，如直接搜索法、直接最優(yōu)法與 Taylor級數(shù)展開法進行線性逼近。此時，可以利用 SAS/STAT 的 nlin 過程實現(xiàn)相應的計算。5.bounds 語句。用于設定參數(shù)的約

14、束，主要是不等式約束，約束間用逗號分隔。例如，boundsa30,1=c=10;6.der.語句。除非在 procnlin 語句中指明所用的迭代法是 dud,使用選擇項 method=dud,否則 der 語句是必需的。der.語句用于計算模型關(guān)于各個參數(shù)的偏導數(shù)，相應的格式為：一階偏導數(shù) der.參數(shù)名=表達式；dataexpd;inputxy;cards;0200.570300.720400.810500.870600.910700.940800.950900.971000.981100.991201.001300.991400.991501.001601.001700.991801.00

15、1901.002000.992101.00；procnlindata=expdbest=10method=gauss;parmsb0=0to2by0.5b1=0.01to0.09by0.01;modely=b0*(1-exp(-b1*x);der.b0=1-exp(-b1*x);der.b1=b0*x*exp(-b1*x);outputout=expoutp=ygs;run;goptionsreset=globalgunit=pctcback=whiteborderhtitle=6htext=3ftext=swissbcolors=(back);procgplotdata=expout;plo

16、ty*xygs*x/haxis=axis1vaxis=axis2overlay;symbol1i=nonev=pluscv=redh=2.5w=2;procnlin 采用最小誤差平方法(LeastSquaresMethod)及循環(huán)推測法(IterativeEstimationMethod)來建立一個非線性模型。一般而言，用戶必須自訂參數(shù)的名字、參數(shù)的啟動值(startingva1ue)、非線性的模型與循環(huán)推測法所用的準則。若用戶不指明，則 nlin 程序自動以高斯-牛頓迭代法(Gauss-Newtoniterativeprocedure)為估計參數(shù)的方法。另外此程序也備有掃描(Gridsear

17、ch)的功能來幫助讀者選擇合適的參數(shù)啟動值。由于非線性回歸分析十分不易處理，nlin 程序不保證一定可以算出符合最小誤差平方法之標準的參數(shù)估計值。nlin 過程的功能，計算非線性模型參數(shù)的最小二乘估計 LS 及加權(quán)最小二乘估計。與 reg過程不同的是：模型的參數(shù)要命名、賦初值、求偏導數(shù)；model 語句與參數(shù)名、解釋變量的表達式有關(guān)；可以使用賦值語句及條件語句。nlin 過程一般由下列語句控制：procnlindata=#集;parameters 參數(shù)名=數(shù)值；model 因變量=表達式;bounds 表達式；der.參數(shù)名.參數(shù)名=表達式；id 變量列表；outputout=數(shù)據(jù)集;by 變

18、量列表；run;二階偏導數(shù) der.參數(shù)名.參數(shù)名=表達式；例如，對于 modely=b0*(1exp(-b1*x);der.語句的書寫格式為 der.b0=1exp(-b1*x);der.b1=b0*x*exp(-b1*x);對于多數(shù)算法，都必須對每個被估計的參數(shù)給出一階偏導數(shù)表達式。對于 newton 法，必須給出一、二階偏導數(shù)表達式。例如，二階偏導數(shù)表達式為，der.b0.b0.=0;der.b0.b1=x*exp(-b1*x)；der.b1.b1=-der.b1*x；7.output 語句用于把一些計算結(jié)果輸出到指定的數(shù)據(jù)集中。有關(guān)的關(guān)鍵字及其意義如表 34.2 所示：表 34.2nl

19、in 過程中 output 語句的關(guān)鍵字關(guān)鍵字意義關(guān)鍵字意義關(guān)鍵字意義predicted|p預測值stdpclm 的標準差u9595%cli 上限r(nóng)esidual|r殘差stdr殘差的標準差l9595%cu 下限parms參數(shù)估計值l95m95%clm 下限student學生氏殘差sse|ess殘差平方和u95m95%clm 上限h杠桿點統(tǒng)計量 hi關(guān)于 nlin 過程的其他選擇項及意義，詳見 SAS/STAT 的用戶手冊。五、實例分析例 34.1 負指數(shù)增長曲線的非線性回歸。根據(jù)對已有數(shù)據(jù)的 XY 散點圖的觀察和分析,發(fā)現(xiàn) Y 隨 X 增長趨勢是減緩的，并且 Y 趨向一個極限值，我們認為用負

20、指數(shù)增長曲線y=b0(1-ex)來擬合模型較為合適。程序如下：symbol2i=joinv=nonel=1h=2.5w=2;axislorder=20to210by10;axis2order=0.5to1.1by0.05;title1y=b0*(1-exp(-b1*x);title2procnlinmethod=gauss;run;程序說明：由于在 nlin 過程中使用選項 method=gauss,即指定用高斯-牛頓迭代算法來尋b/x、找 model 語句中非線性表達式y(tǒng)=b0(1-e1)中參數(shù)b0和b1的最小二乘估計。我們知道參數(shù)初始值選取好壞，對迭代過程是否收斂影響很大。parms 語句

21、設置了初始值網(wǎng)格值為b0取 0,0.5,1,1.5,2 共 5 個值，”取 0.01,0.02,0.09 共 9 個值，所有可能組合為 5X9=45,選項 best=10 要求輸出殘差平方和最小的前 10 種組合。高斯-牛頓迭代算法要求給出模型y=b0(1-4)對參數(shù)b0和b1的一階偏導數(shù)表達式，我們知道der 語句用以表示上面兩個一階偏導數(shù)表達式。output 語句輸出一個新數(shù)據(jù)集 expout,包括原數(shù)據(jù)集和非線性回歸模型的預測值 ygs。gplot 過程的主要作用是繪制輸出數(shù)據(jù)集 expout中的原始數(shù)據(jù)的散點圖及回歸曲線的平滑線。程序的輸出結(jié)果見圖 34-1 和見表 34.3 所示。y

22、=b0*(1-exp(b1*x)procnlin圖 34-1XY 散點圖和非線性回歸曲線dydb0=1-bp-edydh二b0 xe4x表 34.3 負指數(shù)增長曲線：Gauss-Newton 方法的輸出結(jié)果y=b0*(1-exp(-b1*x)procnlinmethod=gaussNon-LinearLeastSquaresGridSearchDependentVariable丫B0B1SumofSquares1.0000000.0400000.0014041.0000000.0500000.0168111.0000000.0600000.0551551.0000000.0300000.066

23、5711.0000000.0700000.0972841.0000000.0800000.1365361.0000000.0900000.1708391.0000000.0200000.4192851.5000000.0100000.9757241.0000000.0100002.165290Non-LinearLeastSquaresIterativePhaseDependentVariableYMethod:Gauss-NewtonIterB0B1SumofSquares01.0000000.0400000.00140410.9961390.0418570.00058020.9961920.0419520.00057730.9961890.0419540.00057740.9961890.0419540.000577NOTE:Convergencecriterionmet.Non-LinearLeastSquaresSummaryStatisticsDependentVariableYSourceDFSumofSquaresMeanSquareRegression217.6717231898.835861595Residual180