SAS讲义第三十四课非线性回归分析

1、第三十四课非线性回归分析现实世界中严格的线性模型并不多见，它们或多或少都带有某种程度的近似；在不少情 况下，非线性模型可能更加符合实际。由于人们在传统上常把“非线性”视为畏途，非线性 回归的应用在国内还不够普及。事实上，在计算机与统计软件十分发达的令天，非线性回归 的基本统计分析已经与线性回归一样切实可行。在常见的软件包中（诸如SAS、SPSS等等），人们已经可以像线性回归一样，方便的对非线性回归进行统计分析。因此，在国内回归分析 方法的应用中，已经到了“更上一层楼”，线性回归与非线性回归同时并重的时候。对变量间非线性相关问题的曲线拟合，处理的方法主要有：首先决定非线性模型的函数类型，对于其中

2、可线性化问题则通过变量变换将其线性化，从而归结为前面的多元线性回归问题来解决。若实际问题的曲线类型不易确定时，由于任意曲线皆可由多项式来逼近，故常可 用多项式回归来拟合曲线。若变量间非线性关系式已知（多数未知），且难以用变量变换法将其线性化，则进行数值迭代的非线性回归分析。一、可变换成线性的非线性回归在实际问题中一些非线性回归模型可通过变量变换的方法化为线性回归问题。例如，对 非线性回归模型2yt = ： o 亠 向 cosixt bi sin ixt i 亠和（34.1）=i即可作变换：x1t = cosxt, x2t 二 sin 人，x3t = cos2xt, x4t = si门2人将其化

3、为多元线性回归模型。一般地，若非线性模型的表达式为：yt =bo bigi xtb2g2 Xbmgm xt（34.2）则可作变量变换：x；t =gi 人,x2t 二 g2 人，x：二 gm 焉（34.3）将其化为线性回归模型的表达式，从而用前面线性模型的方法来解决，其中式（34.3）中的xt也可为自变量构成的向量。这种变量变换法也适用于因变量和待定参数bi。如：yt =aexpbx1t +b2X?t + 匕3仗朮乂玄 一1 9（34.4）时上式两边取对数得：ln yt = In abzXzt b?-1（34.5）现作变换:(34.6)(34.7)(34.8)yt = In yt,b°

4、 二 In a,X3t =-1则可得线性表达式：yt - b0b1X1t b2 X2tb3X3t利用前面方法确定了 ！?，i =0,1,2,3，并由?= exp(?)得到?的值。变量变换的线性化方法可推广到下列形式的非线性模型：h(yt) = Co(b。)G(bi)giXtCm(bm)gmXt其中x= (Xi,X2,Xp),而h (yj、Ci (bi)、gi (xj则分别化为新的因变量、线性回归 参数和自变量，即可归结为线性回归模型来解。表34.1给出了一些常见的可线性化的非线性模型。表34.1典型的函数及线性化方法函数名称函数表达式线性化方法双曲线函数1+b=a + yx11v =u =yx

5、幂函数by =axv = I n y u = I n x指数函数bxy =aev = I n y u = xb/ xy =ae. 1v = I n y u =X对数函数y = a + bl n xv = yu = I n xS型函数1y _ - a + be1 一 v _u _ ey当曲线的函数类型未确定时，我们常采用上述非线性模型作为其拟合曲线，即将自变量 的各种初等函数的组合作为新自变量，用逐步回归法(或正交筛选法等)对新变量进行筛选， 以确定一个项数不多的线性函数表达式。该方法对表达式形式没限制且精度要求不高的问题 颇为有效。二、多项式回归分析在式(34.2)中，若取gj x = Xi，

6、则为多项式回归模型。由数学分析知识可知，一般函数 都可用多项式来逼近，故多项式回归分析可用来处理相当广泛的非线性问题。对观测数据(xt, yt) (t= 1，N),多项式回归模型为：yt = b0 Rxt b2x2 亠 亠 bmxtmt, t=1,2, ,Nyily2>N则模型可表示为:m -mxbi?m 一；2&NY 二 XB ；当X列满秩时，由前面的讨论知，其最小二乘估计为:B = XX JXY由此即可求得其多项式回归方程。但由于XX '的计算既复杂又不稳定，故我们一般采用正交多项式法来进行多项式回归。三、不可变换成线性的非线性回归分析假设因变量y与自变量(xi, X

7、2,，xp )之间满足非线性模型：y 二 F Xi,X2, ,Xp； ：；(34.9)其中=：1, ：2,，：m 为未知参数，F为已知表达式，:为误差项。现将观察数据：yt,Xit,X2t,Xpt ,t=i,2,，N代人式(34.9)得非线性回归模型：% 二 F公戌；：；t,t=1,2, ,N常记为：Y = F ( J E其中，丫二y1, y2/ , yN为y的观察向量，一：二m为非线性回归系数，E=；i,；2,，；N为观察误差向量，F为未知参数一：的函数向量。非线性回归分析就是利用最小二乘准则来估计回归系数1 ,即求?使得残差平方和：Q ： = - E 丄丫 FY F :2 2在-_ ?处达

8、到最小。非线性回归分析一般用数值迭代法来进行，其共同特点是：由选定1的初值1 0出发,通过逐步迭代:t即选择适当的步长t ( >0 )及确疋搜索方向向量 也=(A1，也2,，釧),使得:Q： Q -0(34.10)(34.11)再由取代一：0，重复上述迭代过程，直至qC)可认为达到最小值为止，即可将所得的1作为其最小二乘估计？，从而得到非线性回归方程?=F 花公2, ,Xp； ?1. 下降方向和步长的选择11f首先考察Q 二;EE=?Y_F ： Y-F: 的梯度向量(即导数)：i - Y -F 川川 G Y - F -千；:F为F的梯度矩阵。为使1°迭代收敛到?，其迭代公式应满

9、足下降性质(34.11)。现考虑一元函数t =Q '0：，它从一：0出发以 厶为方向的射线上取值。由复合求导公式得:d 八 t y 二=：Q八 Y F ： G :可以证明，当d<0时，在以厶为方向向量的射线上可以找到一：=' 0 t ，使得j Q 0。我们将满足d<0的厶称为下降方向，Bard于1974年给出了厶为下降方向的 充要条件为：二 PG Y - F 1其中，P为对称正定阵，由此我们可得下降算法的迭代公式为：tPG Y -F- (34.12)其中，P为任意正定阵，G为F的梯度，t为满足Q 1 : Q - 0的正实数，即步长。如何计算厶以便修改参数向量1有五种

10、常用的非线性回归迭代方法：高斯一牛顿法(Gauss-Newton)、最速下降法(梯度法，Gradient)、牛顿法(Newton )、麦夸特法(Marquardt )、 正割法(DUD、。以下我们介绍高斯牛顿法。2. 咼斯一牛顿法首先选取1的一切初始近似值1 0，令.-： - - - -0，则只要确定厶的值即可确定。为此，考虑F(J在一：0处的Taylor展开式，并略去二次以上的项得：F(P)=F(P°+A)=F(P0)J 手 p申也= F(P0)+G迪'0_ j其中，G =务| p上 为F的梯度。此时其残差平方和：Q =丄 Y 一 F " G 厶 Y 一 F -

11、0 -G .=2由 2=o，得其.泊勺正则方程为：gAGG 厶二G Y - F - 0(34.13)故 ：二 GG °G Y - F - 0(34.14)由此即可用前面线性回归法求上只需将G、YF0)视为前面(34.1)式中的X、Y即可。此时，对给定精度 r、 ；2，当=可或Q * 锐2时，即得：最小 二乘法估计閉=P° +;否则用所得的黒代替P 0，重复上述步骤，直至|纠|或Q(0)满足 精度要求为止。该法称为高斯-牛顿法，其一般迭代公式为：卩宀=P +皿(34.15)其中，厶为G ' i G 1=G 1 i Y - F “的解，ti为、t二Q - i的最小值点。

12、高斯-牛顿法在初值1 0选取适当，且G G可逆时非常有效，但在其他情形，其求解较为困难，对此，Marguardt对(34.14)中二的正则系数阵作适当修正，得到了改进算法。四、nlin非线性回归过程在很多场合，可以对非线性模型进行线性化处理，尤其是关于变量非线性的模型，以运用OLS进行推断。对线性化后的线性模型，可以应用SAS的reg过程进行计算。多项式模型可以直接应用glm (广义线性模型)求解。对于不能线性化的非线性模型。其估计不能直接运用经典的最小二乘法，而需要运用其他估计方法，如直接搜索法、直接最 优法与Taylor级数展开法进行线性逼近。此时，可以利用SAS/STAT的nlin过程实

13、现相应的计算。1. proc nlin 过程proc nlin 采用最小误差平方法(Least Squares Method)及循环推测法(Iterative EstimationMethod)来建立一个非线性模型。一般而言，用户必须自 定参数的名字、参数的启动值(startingva1ue)、非线性的模型与循环推测法所用的准则。若用户不指明，则nlin程序自动以高斯一牛顿迭代法(Gauss-Newt on iterative procedure)为估计参数的方法。另外此程序也备有扫描(Grid search)的功能来帮助读者选择合适的参数启动值。由于非线性回归分析十分不易处理，nlin程序不

14、保证一定可以算出符合最小误差平方法之标准的参数估计值。nlin过程的功能，计算非线性模型参数的最小二乘估计LS及加权最小二乘估计。与reg过程不同的是：模型的参数要命名、赋初值、求偏导数；model语句与参数名、解释变量的表达式有关；可以使用赋值语句及条件语句。nlin过程一般由下列语句控制:proc nlindata据集 /选项列表;parameters参数名=数值；model因变量=表达式/选项列表；bounds表达式 ；der.参数名参数名=表达式；id变量列表；output out=数据集/选项列表;by变量列表；run ;其中，parameters语句和model语句是必需的，而其余

15、语句供用户根据需要选择。2. proc nlin语句中的主要选择项outesA数据集名一一指定存放参数估计的每步迭代结果的数据集名。best= n要求过程只输出网格点初始值可能组合中最好的n组残差平方和。method= gauss | marquardt | newton| gradient| dud |设定参数估计的迭代方法。缺省时为gauss,除非没有der.语句。eformat要求所有数值以科学记数法输出。nopoint抑制打印输出。noinpoint抑制迭代结果的输出。3. parameters(parms)语句用于对所有参数赋初值，项目之间以空格分隔。 例如，parms b0=0b1

16、=1 to 10b2=1to 10 by 2b3=1,10,100;4. model 语句表达式可以是获得数值结果的任意有效SAS表达式。这个表达式包括参数名字、输入数据集中的变量名以及在nlin过程中用程序设计语句创建的新变量。例如， model y=b0*(1 exp(-b1*x)；bounds a<=20.5. bounds 语句用于设定参数的约束，主要是不等式约束，约束间用逗号分隔。例如,b>30， l<=c<=10 ；6. der.语句除非在proc nlin语句中指明所用的迭代法是dud，使用选择项 method = dud，否则der语句是必需的。der.

17、语句用于计算模型关于各个参数的偏导数，相应的格式为：一阶偏导数der .参数名=表达式 ；二阶偏导数der .参数名.参数名=表达式；例如，对于 model y=bO*( 1 exp(-b1*x)； der.语句 的书写格式为 der. b0=1 exp(-b1*x) ; der. b1=b0*x*exp(-b1*x);对于多数算法，都必须对每个被估计的参数给出一阶偏导数表达式。对于newton法，必须给出一、二阶偏导数表达式。例如，二阶偏导数表达式为，der.bO.bO.=O ;der.bO.b仁x*exp(-b1*x)； der.b1.b仁-der.b1*x；7. output 语句用于把

18、一些计算结果输出到指定的数据集中。有关的关键字及其意义如表34.2所示。表34.2nlin过程中output语句的关键字关键字意义关键字意义关键字意义predicted | p预测值stdpclm的标准差u9595 % cli上限residual | r残差stdr残差的标准差l9595%cu下限parms参数估计值l95m95% clm下限stude nt学生氏残差sse | ess残差平方和u95m95% clm上限h杠杆点统计量hi关于nlin过程的其他选择项及意义，详见SAS/ STAT的用户手册。五、实例分析例34.1负指数增长曲线的非线性回归。根据对已有数据的XY散点图的观察和分析

19、,发现Y随X增长趋势是减缓的，并且Y趋向一个极限值，我们认为用负指数增长曲线y二b°(1 - ex)来拟合模型较为合适。程序如下：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.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=

20、10 method=gauss;parms bO=O 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=global gunit=pct cback=white border htitle=6 htext=3 ftext=swissb colors=(back);proc gplot data=expout;plot y*x ygs*x /haxis

21、=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*(1-exp(-b1*x)'title2 'proc nlin method=gauss'run ;程序说明：由于在 nlin过程中使用选项 method=gauss,即指定用高斯一牛顿迭代算法来寻找model语句中非

22、线性表达式 y二b0 (1 - e 中参数b0和b1的最小二乘估计。我们知道， 参数初始值选取好坏，对迭代过程是否收敛影响很大。parms语句设置了初始值网格值为b0取0, 0.5, 1, 1.5, 2共5个值，bi取0.01, 0.02,，0.09共9个值，所有可能组合为 5X 9=45,选项best=10要求输出残差平方和最小的前10种组合。高斯牛顿迭代算法要求给出模型y = b° (1 - e_b,x)对参数b°和d的一阶偏导数表达式，我们知道：db°der.语句用以表示上面两个一阶偏导数表达式。output语句输出一个新数据集 expout，包括原数据集和

23、非线性回归模型的预测值ygs。gplot过程的主要作用是绘制输出数据集 expout中的原始数据的散点图及回归曲线的平滑线。程序的输出结果如图34. 1和表34.3所示。y= bo*(1exp(b1*x)pro© Ein nw恤皿卯谣$1.00 J05S：0,5Q;I ' I ''' I7 !|T1-p1-. p .-.( I1=!=1> I I-I -|-I|. s .|. -rI1-. J ,-,|30 3Q -109060 TO 8Q 汕 1Q0 TIQ 120 130 刚 1S0 KO 170 WO 190 2Q0 210图34. 1 X

24、Y散点图和非线性回归曲线表34.3负指数增长曲线:Gauss-Newton方法的输出结果y=bO*(1-exp(-b1*x)proc nlin method=gaussNon-Linear Least Squares Grid Search Dependent Variable YB0B1Sum of Squares1.0000000.0400000.0014041.0000000.0500000.0168111.0000000.0600000.0551551.0000000.0300000.0665711.0000000.0700000.0972841.0000000.0800000.136

25、5361.0000000.0900000.1708391.0000000.0200000.4192851.5000000.0100000.9757241.0000000.0100002.165290Non-Linear Least Squares Iterative Phase Dependent Variable Y Method: 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 Statistics Dependent Variable Y SourceDF Sum of Squares