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

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

2、中可線性化問題則通過變量變換將其線性化，從而歸結(jié)為前面的多元線性回歸問題來解決。若實際問題的曲線類型不易確定時，由于任意曲線皆可由多項式來逼近，故?？捎枚囗検交貧w來擬合曲線。若變量間非線性關(guān)系式已知（多數(shù)未知） ，且難以用變量變換法將其線性化，則進行數(shù)值迭代的非線性回歸分析。一、可變換成線性的非線性回歸在實際問題中一些非線性回歸模型可通過變量變換的方法化為線性回歸問題。例如，對非線性回歸模型2yt0ai cosixt bi sin ix tt(34.1)i1即可作變換：x1tcosxt , x2tsin xt , x3tcos2xt , x4tsin 2 xt將其化為多元線性回歸模型。一般地，

3、若非線性模型的表達式為：yt b0 b1 g1 xtb2 g2xtbm gm xt(34.2)則可作變量變換：x1*tg1 xt , x2*tg 2 xt, xmt*gm xt(34.3)將其化為線性回歸模型的表達式，從而用前面線性模型的方法來解決，其中式 (34.3)中的 xt 也可為自變量構(gòu)成的向量。這種變量變換法也適用于因變量和待定參數(shù)bi 。如：yta exp b1 x1tb2 x2 tb3 x1t x2t1(34.4)時上式兩邊取對數(shù)得：ln ytln a b1 x1tb2 x2 tb3x1t x2t1(34.5)現(xiàn)作變換：yt*ln yt , b0 ln a, x3 tx1t x2

4、t1(34.6)則可得線性表達式：yt*b0b1 x1tb2 x2tb3 x3t(34.7)利用前面方法確定了?,0,1,2,3，并由?)abiiexp(b0a得到 ? 的值。變量變換的線性化方法可推廣到下列形式的非線性模型：h( yt ) c0 (b0 )c1 (b1 )g1 xtcm (bm ) gm xt(34.8)其中 x=（ x1,x2, ,xp），而 h（ yt）、 ci （ bi）、 gi（ xt）則分別化為新的因變量、線性回歸參數(shù)和自變量，即可歸結(jié)為線性回歸模型來解。表34.1 給出了一些常見的可線性化的非線性模型。表 34.1典型的函數(shù)及線性化方法函數(shù)名稱函數(shù)表達式線性化方法

5、雙曲線函數(shù)1abv1u1yxyx冪函數(shù)yaxbvln yuln xyaebxvln yux指數(shù)函數(shù)1yaeb / xvln yux對數(shù)函數(shù)yab ln xvyuln xS 型函數(shù)y1v1uexabexy當曲線的函數(shù)類型未確定時，我們常采用上述非線性模型作為其擬合曲線，即將自變量的各種初等函數(shù)的組合作為新自變量，用逐步回歸法（或正交篩選法等）對新變量進行篩選，以確定一個項數(shù)不多的線性函數(shù)表達式。該方法對表達式形式?jīng)]限制且精度要求不高的問題頗為有效。二、多項式回歸分析在式 (34.2) 中，若取 g i xxi ，則為多項式回歸模型。由數(shù)學分析知識可知，一般函數(shù)都可用多項式來逼近，故多項式回歸分析

6、可用來處理相當廣泛的非線性問題。對觀測數(shù)據(jù)（xt，yt）（ t= 1 ， N），多項式回歸模型為：ytb0b1 xtb2 xt2bm xtmt ， t=1,2,N令y11x1x12x1mb01Yy2， X1x2x22x2mb1，2， ByN1xNxN2xNmbmN則模型可表示為：Y XB當 X 列滿秩時，由前面的討論知，其最小二乘估計為：?X X1X YB由此即可求得其多項式回歸方程。但由于X X 1的計算既復雜又不穩(wěn)定，故我們一般采用正交多項式法來進行多項式回歸。三、不可變換成線性的非線性回歸分析假設(shè)因變量 y 與自變量（ x1， x2,， xp ）之間滿足非線性模型：yF x1, x2 ,

7、 xp ;(34.9)其中，1, 2,m為未知參數(shù)， F 為已知表達式，為誤差項?，F(xiàn)將觀察數(shù)據(jù)：yt , x1t , x2 t , xpt，t=1,2,N代人式 (34.9) 得非線性回歸模型：ytFx1t , x2 t , x pt ;t ，t=1,2,N常記為：YF ()E其中， Yy1 , y2 , , yN為 y 的觀察向量，1 , m為非線性回歸系數(shù)， E= 1,2, ,N為觀察誤差向量， F 為未知參數(shù)的函數(shù)向量。 非線性回歸分析就是利用最小二乘準則來估計回歸系數(shù)，即求? 使得殘差平方和：Q1 E E1 YFYF22在? 處達到最小。非線性回歸分析一般用數(shù)值迭代法來進行，其共同特點

8、是：由選定的初值0 出發(fā)，通過逐步迭代：0t(34.10)即選擇適當?shù)牟介Lt ( >0 ) 及確定搜索方向向量（1 ，2，m），使得：QQ0(34.11)再由取代0 ，重復上述迭代過程，直至Q( )可認為達到最小值為止，即可將所得的作為其最小二乘估計?，從而得到非線性回歸方程?, x2 , , x p ;?y F x11. 下降方向和步長的選擇首先考察 Q1 E E1 Y FY F的梯度向量（即導數(shù)） ：22QFG Y FY F其中， GFF ,F為 F 的梯度矩陣。1m為使0迭代收斂到? ，其迭代公式應(yīng)滿足下降性質(zhì)(34.11) ?，F(xiàn)考慮一元函數(shù)t Q0t，它從0 出發(fā)以為方向的射線上

9、取值。由復合求導公式得：dt t 0QY FG可以證明，當d<0 時，在以為方向向量的射線上可以找到0t ，使得QQ0 。我們將滿足d<0 的稱為下降方向 ，Bard 于 1974 年給出了為下降方向的充要條件為：PGYF其中， P 為對稱正定陣，由此我們可得下降算法的迭代公式為：0tPG Y F(34.12)其中， P 為任意正定陣， G 為 F 的梯度， t 為滿足 QQ0 的正實數(shù)，即步長。如何計算以便修改參數(shù)向量有五種常用的非線性回歸迭代方法：高斯牛頓法（ Gauss-Newton ）、最速下降法 （梯度法， Gradient）、牛頓法（ Newton ）、麥夸特法（ Ma

10、rquardt ）、正割法（ DUD ）。以下我們介紹高斯牛頓法。2. 高斯牛頓法首先選取的一切初始近似值0 ，令0 ，則只要確定的值即可確定。為此，考慮 F () 在0 處的 Taylor 展開式，并略去二次以上的項得：FF0F0F0F0G其中，GF0為 F 的梯度。此時其殘差平方和：Q1 YF0GYF0G由 Q20 ，得其 的正則方程為：G GG YF0(34.13)故G G1G YF0(34.14)由此即可用前面線性回歸法求，只需將 G 、YF (0 ) 視為前面（34.1）式中的 X、Y 即可。此時，對給定精度1 、 2，當 maxi1 或 Q02時，即得最小i二乘法估計?0；否則用所

11、得的?代替0 ，重復上述步驟，直至i或 Q()滿足精度要求為止。該法稱為高斯牛頓法，其一般迭代公式為：i1iti(34.15)其中，為 Gi GiGiYFi的解， ti為 tQit的最小值點。0高斯牛頓法在初值 選取適當，且 G G 可逆時非常有效，但在其他情形，其求解較為困難，對此， Marguardt 對 (34.14)中 的正則系數(shù)陣作適當修正，得到了改進算法。四、 nlin 非線性回歸過程在很多場合，可以對非線性模型進行線性化處理，尤其是關(guān)于變量非線性的模型，以運用 OLS 進行推斷。對線性化后的線性模型，可以應(yīng)用SAS 的 reg 過程進行計算。多項式模型可以直接應(yīng)用glm （廣義線

12、性模型）求解。對于不能線性化的非線性模型。其估計不能直接運用經(jīng)典的最小二乘法，而需要運用其他估計方法，如直接搜索法、直接最優(yōu)法與 Taylor 級數(shù)展開法進行線性逼近。此時，可以利用SAS/STAT 的 nlin 過程實現(xiàn)相應(yīng)的計算。1. proc nlin 過程proc nlin 采用最小誤差平方法（ Least Squares Method）及循環(huán)推測法 （ Iterative EstimationMethod ）來建立一個非線性模型。一般而言，用戶必須自 定參數(shù)的名字、 參數(shù)的啟動值 （ startingva1ue）、非線性的模型與循環(huán)推測法所用的準則。若用戶不指明，則nlin 程序自動

13、以高斯牛頓迭代法（ Gauss-Newton iterativeprocedure）為估計參數(shù)的方法。另外此程序也備有掃描（Grid search）的功能來幫助讀者選擇合適的參數(shù)啟動值。由于非線性回歸分析十分不易處理，nlin 程序不保證一定可以算出符合最小誤差平方法之標準的參數(shù)估計值。nlin 過程的功能，計算非線性模型參數(shù)的最小二乘估計LS 及加權(quán)最小二乘估計。與過程不同的是：模型的參數(shù)要命名、賦初值、求偏導數(shù)；model 語句與參數(shù)名、解釋變量的表達式有關(guān)；可以使用賦值語句及條件語句。regnlin過程一般由下列語句控制：proc nlindata=數(shù)據(jù)集</選項列表>par

14、ameters參數(shù)名 =數(shù)值；model因變量表達式</選項列表>；bounds表達式；der .參數(shù)名參數(shù)名= 表達式；id變量列表；outputout=數(shù)據(jù)集</選項列表>；by變量列表；run ;其中， parameters 語句和 model 語句是必需的，而其余語句供用戶根據(jù)需要選擇。2. proc nlin 語句中的主要選擇項outest數(shù)據(jù)集名 指定存放參數(shù)估計的每步迭代結(jié)果的數(shù)據(jù)集名。bestn 要求過程只輸出網(wǎng)格點初始值可能組合中最好的n 組殘差平方和。method gauss | marquardt | newton| gradient| dud |

15、 設(shè)定參數(shù)估計的迭代方法。缺省時為 gauss，除非沒有der.語句。eformat 要求所有數(shù)值以科學記數(shù)法輸出。nopoint 抑制打印輸出。noinpoint 抑制迭代結(jié)果的輸出。3. parameters（ parms）語句用于對所有參數(shù)賦初值，項目之間以空格分隔。例如， parms b0=0b1=1 to 10b2=1to 10 by 2b3=1,10,100；4. model 語句表達式可以是獲得數(shù)值結(jié)果的任意有效SAS 表達式。這個表達式包括參數(shù)名字、輸入數(shù)據(jù)集中的變量名以及在nlin過程中用程序設(shè)計語句創(chuàng)建的新變量。例如，model y=b0*( 1exp(-b1*x)；5.

16、bounds 語句用于設(shè)定參數(shù)的約束，主要是不等式約束，約束間用逗號分隔。例如，boundsa<=20，b>30， 1<=c<=10 ；6. der. 語句除非在proc nlin語句中指明所用的迭代法是dud，使用選擇項method dud，否則der.語句是必需的。der.語句用于計算模型關(guān)于各個參數(shù)的偏導數(shù)，相應(yīng)的格式為：一階偏導數(shù)der參數(shù)名表達式；二階偏導數(shù)der參數(shù)名參數(shù)名表達式；例 如 ， 對 于model y=b0*( 1 exp(-b1*x)； der.語句的書寫格式為der. b0=1 exp(-b1*x)； der.b1=b0*x*exp(-b1*

17、x)；對于多數(shù)算法，都必須對每個被估計的參數(shù)給出一階偏導數(shù)表達式。對于newton 法，必須 給 出 一 、 二 階 偏 導 數(shù) 表 達 式 。 例 如 ， 二 階 偏 導 數(shù) 表 達 式 為 ，； ；7. output 語句用于把一些計算結(jié)果輸出到指定的數(shù)據(jù)集中。有關(guān)的關(guān)鍵字及其意義如表34.2 所示。表 34.2nlin 過程中 output 語句的關(guān)鍵字關(guān)鍵字意義關(guān)鍵字意義關(guān)鍵字意義predicted|p預測值stdpclm 的標準差u9595cli 上限r(nóng)esidual |r殘差stdr殘差的標準差l9595cu 下限parms參數(shù)估計值l95m95 clm 下限student學生氏殘

18、差sse | ess殘差平方和u95m95 clm 上限h杠桿點統(tǒng)計量 hi關(guān)于nlin 過程的其他選擇項及意義，詳見SAS STAT 的用戶手冊。五、實例分析發(fā)現(xiàn)例 34.1Y 隨負指數(shù)增長曲線的非線性回歸。根據(jù)對已有數(shù)據(jù)的 XY 散點圖的觀察和分析， X 增長趨勢是減緩的，并且 Y 趨向一個極限值，我們認為用負指數(shù)增長曲線y b0 (1 e b1 x ) 來擬合模型較為合適。程序如下：data expd;input x y ;cards;020 0.57 030 0.72 040 0.81 050 0.87 060 0.91 070 0.94080 0.95 090 0.97 100 0.

19、98 110 0.99 120 1.00 130 0.99140 0.99 150 1.00 160 1.00 170 0.99 180 1.00 190 1.00200 0.99 210 1.00;proc nlin data=expd best=10 method=gauss;parms b0=0 to 2 by 0.5 b1=0.01 to 0.09 by 0.01;model y=b0*(1-exp(-b1*x);der.b0=1-exp(-b1*x);der.b1=b0*x*exp(-b1*x);output out=expout p=ygs ;run;goptions reset=

20、global gunit=pct cback=white borderhtitle=6 htext=3 ftext=swissb colors=(back);proc gplot data=expout;plot y*x ygs*x /haxis=axis1 vaxis=axis2 overlay;symbol1 i=none v=plus cv=red h=2.5 w=2;symbol2 i=join v=none l=1 h=2.5 w=2;axis1 order=20 to 210 by 10;axis2 order=0.5 to 1.1 by 0.05;title1 'y=b0

21、*(1-exp(-b1*x)'title2 'proc nlin method=gauss'run ;程序說明：由于在nlin 過程中使用選項method=gauss，即指定用高斯牛頓迭代算法來尋找 model 語句中非線性表達式y(tǒng) b0 (1 e b1 x ) 中參數(shù) b0 和 b1 的最小二乘估計。 我們知道，參數(shù)初始值選取好壞， 對迭代過程是否收斂影響很大。 parms 語句設(shè)置了初始值網(wǎng)格值為b0 取0，0.5，1，1.5，2 共 5 個值， b1 取 0.01，0.02， 0.09 共 9 個值， 所有可能組合為5× 9=45,選項best=10 要

22、求輸出殘差平方和最小的前10 種組合。高斯牛頓迭代算法要求給出模型y b0 (1 e b1 x ) 對參數(shù) b0 和 b1 的一階偏導數(shù)表達式，我們知道：d y1e b1xd b0d yb0 xe b1 xd b1der.語句用以表示上面兩個一階偏導數(shù)表達式。output 語句輸出一個新數(shù)據(jù)集expout，包括原數(shù)據(jù)集和非線性回歸模型的預測值ygs。gplot過程的主要作用是繪制輸出數(shù)據(jù)集expout中的原始數(shù)據(jù)的散點圖及回歸曲線的平滑線。程序的輸出結(jié)果如圖34. 1 和表 34.3 所示。圖 34. 1 XY 散點圖和非線性回歸曲線表 34.3負指數(shù)增長曲線 :Gauss-Newton 方法

23、的輸出結(jié)果y=b0*(1-exp(-b1*x)proc nlin method=gaussNon-Linear Least Squares Grid SearchDependent Variable YB0B1Sum of Squares1.0000000.0400000.0014041.0000000.0500000.0168111.0000000.0600000.0551551.0000000.0300000.0665711.0000000.0700000.0972841.0000000.0800000.1365361.0000000.0900000.1708391.0000000.020

24、0000.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 met.Non-Linear Least Squares Summary StatisticsDependent Variable YSourceDF Sum of SquaresMean SquareRegression217.6717231898.835861595Residua