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

【精品文档】如有侵权，请联系网站删除，仅供学习与交流第三十四课第三十五课第三十六课第三十七课 SAS讲义第三十四课非线性回归分析.精品文档.第三十八课 非线性回归分析现实世界中严格的线性模型并不多见，它们或多或少都带有某种程度的近似；在不少情况下，非线性模型可能更加符合实际。由于人们在传统上常把“非线性”视为畏途，非线性回归的应用在国内还不够普及。事实上，在计算机与统计软件十分发达的令天，非线性回归的基本统计分析已经与线性回归一样切实可行。在常见的软件包中（诸如SAS、SPSS等等），人们已经可以像线性回归一样，方便的对非线性回归进行统计分析。因此，在国内回归分析方法的应用中，已经到了“更上一层楼”，线性回归与非线性回归同时并重的时候。对变量间非线性相关问题的曲线拟合，处理的方法主要有：l 首先决定非线性模型的函数类型，对于其中可线性化问题则通过变量变换将其线性化，从而归结为前面的多元线性回归问题来解决。l 若实际问题的曲线类型不易确定时，由于任意曲线皆可由多项式来逼近，故常可用多项式回归来拟合曲线。l 若变量间非线性关系式已知（多数未知），且难以用变量变换法将其线性化，则进行数值迭代的非线性回归分析。一、 可变换成线性的非线性回归在实际问题中一些非线性回归模型可通过变量变换的方法化为线性回归问题。例如，对非线性回归模型(34.1)即可作变换将其化为多元线性回归模型。一般地，若非线性模型的表达式为(34.2)则可作变量变换(34.3)将其化为线性回归模型的表达式，从而用前面线性模型的方法来解决，其中(34.3)中的xt 也可为自变量构成的向量。这种变量变换法也适用于因变量和待定参数 bi 。如(34.4)时上式两边取对数得(34.5)现作变换(34.6)则可得线性表达式(34.7)利用前面方法确定了，并由得到 的值。变量变换的线性化方法可推广到下列形式的非线性模型(34.8)其中x=（x1,x2,L,xp），而h（yt）、ci（bi）、gi（xt）则分别化为新的因变量、线性回归参数和自变量，即可归结为线性回归模型来解。见表34.1所示给出了一些常见的可线性化的非线性模型。表34.1 典型的函数及线性化方法函数名称函数表达式线性化方法双曲线函数 幂函数 指数函数 对数函数 S型函数 当曲线的函数类型未确定时，我们常采用上述非线性模型作为其拟合曲线，即将自变量的各种初等函数的组合作为新自变量，用逐步回归法（或正交筛选法等）对新变量进行筛选，以确定一个项数不多的线性函数表达式。 该方法对表达式形式没限制且精度要求不高的问题颇为有效。二、 多项式回归分析在式(34.2)中，若取，则为多项式回归模型。由数学分析知识可知，一般函数都可用多项式来逼近，故多项式回归分析可用来处理相当广泛的非线性问题。对观测数据（xt，yt）（t= 1，N），多项式回归模型为 ，t=1,2,L,N令则模型可表示为当X列满秩时，由前面的讨论知，其最小二乘估计为由此即可求得其多项式回归方程。但由于的计算既复杂又不稳定，故我们一般采用正交多项式法来进行多项式回归。三、 不可变换成线性的非线性回归分析假设因变量y与自变量（x1，x2,，xp）之间满足非线性模型(34.9)其中 为未知参数，F为已知表达式，仅未知的非线性函数，e 为误差项。现将观察数据 ， t=1,2,L,N代人上式(34.9)得非线性回归模型 ， t=1,2,L,N常记为其中为y的观察向量，为非线性回归系数，E =为观察误差向量，F为未知参数的函数向量。非线性回归分析就是利用最小二乘准则来估计回归系数，即求 使得残差平方和在 处达到最小。非线性回归分析一般来用数值迭代法来进行，其共同特点是：由选定的初值出发，通过逐步迭代(34.10)即选择适当的步长t ( >0 ) 及确定搜索方向向量D（D1，D2，Dm），使得(5.4.11)再由取代，重复上述迭代过程，直至 Q(b)可认为达到最小值为止，即可将所得的作为其最小二乘估计，从而得到非线性回归方程1. 下降方向和步长的选择首先考察的梯度向量（即导数）其中为F的梯度矩阵。为使迭代收敛到，其迭代公式应满足下降性质(5.4.11)。现考虑一元函数，它从出发以 D为方向的射线上取值。由复合求导公式得可以证明，当 d<0 时，在以 D为方向向量的射线上可以找到，使得。我们将满足 d<0 的D称为下降方向，Bard于1974年给出了D为下降方向的充要条件为其中P为对称正定阵，由此我们可得下降算法的迭代公式为(34.12)其中P为任意正定阵，G为F的梯度，t为满足的正实数，即步长。如何计算D以便修改参数向量有五种常用的非线性回归迭代方法：高斯-牛顿法（Gauss-Newton）、最速下降法（梯度法，Gradient）、牛顿法（Newton）、麦夸特法（Marquardt）、正割法（DUD）。以下我们介绍其中高斯-牛顿法。2. GaussNewton法首先选取的一切初始近似值，令，则只要确定D的值即可确定。为此，考虑在处的Taylor展开式，并略去二次以上的项得其中为F的梯度。此时其残差平方和由，得其D的正则方程为(34.13)故(34.14)由此即可用前面线性回归法求D，只需将、视为前面（5.2.1）式中的X、Y即可。此时，对给定精度、 ，当或时，即得的最小二乘法估计；否则用所得的代替，重复上述步骤，直至或Q(b)满足精度要求为止。该法称为GaussNewton法，其一般迭代公式为(34.15)其中：D为的解，ti为的最小值点。Gauss-Newton法在初值选取适当，且可逆时非常有效，但在其他情形，其求解 较为困难，对此，Marguardt对(34.14)中D的正则系数阵作适当修正，得到了改进算法。四、 nlin非线性回归过程在很多场合，可以对非线性模型进行线性化处理，尤其是关于变量非线性的模型，以运用OLS进行推断。对线性化后的线性模型，可以应用SAS的reg过程进行计算。多项式模型可以直接应用glm（广义线性模型）求解。对于不能线性化的非线性模型。其估计不能直接运用经典的最小二乘法，而需要运用其他估计方法，如直接搜索法、直接最优法与Taylor级数展开法进行线性逼近。此时，可以利用SAS/STAT的nlin过程实现相应的计算。1. proc nlin过程proc nlin采用最小误差平方法（Least Squares Method）及循环推测法（Iterative Estimation Method）来建立一个非线性模型。一般而言，用户必须自订参数的名字、参数的启动值（starting va1ue）、非线性的模型与循环推测法所用的准则。若用户不指明，则nlin 程序自动以高斯-牛顿迭代法（Gauss-Newton iterative procedure）为估计参数的方法。另外此程序也备有扫描（Grid search）的功能来帮助读者选择合适的参数启动值。由于非线性回归分析十分不易处理，nlin程序不保证一定可以算出符合最小误差平方法之标准的参数估计值。nlin过程的功能，计算非线性模型参数的最小二乘估计LS及加权最小二乘估计。与reg过程不同的是：模型的参数要命名、赋初值、求偏导数；model语句与参数名、解释变量的表达式有关；可以使用赋值语句及条件语句。nlin过程一般由下列语句控制：proc nlin data=数据集 </选项列表>parameters 参数名=数值 ；model 因变量表达式 </选项列表>；bounds 表达式 ；der.参数名参数名= 表达式；id 变量列表；output out=数据集 </选项列表>；by 变量列表；run ;其中，parameters语句和model语句是必需的，而其余语句供用户根据需要选择。2. proc nlin语句中的主要选择项。l outest数据集名指定存放参数估计的每步迭代结果的数据集名。l bestn要求过程只输出网格点初始值可能组合中最好的n组残差平方和。l methodgauss | marquardt | newton| gradient| dud |设定参数估计的迭代方法。缺省时为gauss，除非没有der.语句。l eformat要求所有数值以科学记数法输出。l nopoint抑制打印输出。l noinpoint抑制迭代结果的输出。3. parameters（parms）语句。用于对所有参数赋初值，项目之间以空格分隔。例如，parms b0=0 b1=1 to 10 b2=1 to 10 by 2 b3=1,10,100；4. model语句。表达式可以是获得数值结果的任意有效SAS表达式。这个表达式包括参数名字、输入数据集中的变量名以及在nlin过程中用程序设计语句创建的新变量。例如，model y=b0*(1exp(-b1*x)；5. bounds语句。用于设定参数的约束，主要是不等式约束，约束间用逗号分隔。例如，bounds a<=20，b>30，1<=c<=10；6. der.语句。除非在proc nlin语句中指明所用的迭代法是dud，使用选择项methoddud，否则der.语句是必需的。der.语句用于计算模型关于各个参数的偏导数，相应的格式为： 一阶偏导数 der参数名表达式 ； 二阶偏导数 der参数名参数名表达式；例如，对于model y=b0*(1exp(-b1*x)；der.语句的书写格式为der.b0=1exp(-b1*x)；der.b1=b0*x*exp(-b1*x)；对于多数算法，都必须对每个被估计的参数给出一阶偏导数表达式。对于newton法，必须给出一、二阶偏导数表达式。例如，二阶偏导数表达式为，der.b0.b0.=0；der.b0.b1=x*exp(-b1*x)；der.b1.b1=-der.b1*x；7. output语句用于把一些计算结果输出到指定的数据集中。有关的关键字及其意义如表34.2所示：表34.2 nlin过程中output语句的关键字关键字意义关键字意义关键字意义predicted | p预测值stdpclm的标准差u9595cli上限residual | r残差stdr残差的标准差l9595cu下限parms参数估计值l95m95clm下限student学生氏残差sse | ess残差平方和u95m95clm上限h杠杆点统计量hi关于 nlin过程的其他选择项及意义，详见SASSTAT的用户手册。五、 实例分析例34.1 负指数增长曲线的非线性回归。根据对已有数据的XY散点图的观察和分析，发现Y随X增长趋势是减缓的，并且Y趋向一个极限值，我们认为用负指数增长曲线来拟合模型较为合适。程序如下：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=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=global gunit=pct cback=white border htitle=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*(1-exp(-b1*x)'title2 'proc nlin method=gauss'run ;程序说明：由于在nlin过程中使用选项method=gauss，即指定用高斯-牛顿迭代算法来寻找model语句中非线性表达式中参数和的最小二乘估计。我们知道参数初始值选取好坏，对迭代过程是否收敛影响很大。parms语句设置了初始值网格值为取0，0.5，1，1.5，2共5个值，取0.01，0.02，0.09共9个值，所有可能组合为5×9=45,选项best=10要求输出残差平方和最小的前10种组合。高斯-牛顿迭代算法要求给出模型对参数和的一阶偏导数表达式，我们知道der.语句用以表示上面两个一阶偏导数表达式。output语句输出一个新数据集expout，包括原数据集和非线性回归模型的预测值ygs。gplot过程的主要作用是绘制输出数据集expout中的原始数据的散点图及回归曲线的平滑线。程序的输出结果见图341和见表34.3所示。图341 XY散点图和非线性回归曲线 y=b0*(1-exp(-b1*x) proc nlin method=gauss Non-Linear Least Squares Grid Search Dependent Variable Y B0 B1 Sum of Squares 1.000000 0.040000 0.001404 1.000000 0.050000 0.016811 1.000000 0.060000 0.055155 1.000000 0.030000 0.066571 1.000000 0.070000 0.097284 1.000000 0.080000 0.136536 1.000000 0.090000 0.170839 1.000000 0.020000 0.419285 1.500000 0.010000 0.975724 1.000000 0.010000 2.165290Non-Linear Least Squares Iterative Phase Dependent Variable Y Method: Gauss-Newton Iter B0 B1 Sum of Squares 0 1.000000 0.040000 0.001404 1 0.996139 0.041857 0.000580 2 0.996192 0.041952 0.000577 3 0.996189 0.041954 0.000577 4 0.996189 0.041954 0.000577NOTE: Convergence criterion met. Non-Linear Least Squares Summary Statistics Dependent Variable Y Source DF Sum of Squares Mean Square Regression 2 17.671723189 8.835861595 Residual 18 0.000576811 0.000032045 Uncorrected Total 20 17.672300000 (Corrected Total) 19 0.243855000 Parameter Estimate Asymptotic Asymptotic 95 % Std. Error Confidence Interval Lower Upper B0 0.9961885657 0.00161380015 0.99279811976 0.99957901159 B1 0.0419538868 0.00039822900 0.04111724424 0.04279052935 Asymptotic Correlation Matrix Corr B0 B1 - B0 1 -0.555895746 B1 -0.555895746 1表34.3 负指数增长曲线:Gauss-Newton方法的输出结果输出结果分析：首先列表输出了10个最好的迭代初始值组合，按回归模型的残差平方和从小到大排列。因此，最好的迭代初始值为B0=1.000000，B1=0.040000，此时回归模型的ESS=0.001404，与其他初始值组合的ESS相比是最小的。紧接着输出迭代过程分析表，从这个最好迭代初始值开始经过4次迭代误差平方和的变化就满足收敛准则，停止迭代，例如，经过第1次迭代，B0值从1修正为0.996139，B1值从0.040000修正为0.041857，而ESS从0.001404减小到0.000580，而从第2次迭代开始，B0 的值只有在小数点四位以后才有变化，B1的值只有在小数点五位以后才有变化，而ESS值几乎不变了，因此是满足了收敛性标准而停止迭代的。方差分析表给出，回归平方和为17.671723189，残差平方和为0.000576811，总平方和为17.672300000，总偏平方和为0.243855000。参数估计表给出了，B0和B1的渐近估计值，因此得到的非线性回归模型为同时还给出B0和B1参数估计的渐近有效标准差（Asymptotic Std. Error）和渐近95%置信区间（Asymptotic 95 % Confidence Interval）。最后还给出了参数的渐近相关阵（Asymptotic Correlation Matrix）。