SAS学习系列25. 非线性回归

25.非线性回归现实世界中严格的线性模型并不多见，它们或多或少都带有某种程度的近似；在不少情况下，非线性模型可能更加符合实际。对变量间非线性相关问题的曲线拟合，处理的方法主要有：（1）首先确定非线性模型的函数类型，对于其中可线性化问题则通过变量变换将其线性化，从而归结为前面的多元线性回归问题来解决；（2）若实际问题的曲线类型不易确定时，由于任意曲线皆可由多项式来逼近，故常可用多项式回归来拟合曲线；（3）若变量间非线性关系式已知（多数未知），且难以用变量变换法将其线性化，则进行数值迭代的非线性回归分析。一）可变换为线性的非线性回归表1典型的函数及线性化方法"窗数窑称.函数盍达式中线性化方法.取曲线函数J1b—=d——.■y工V——u——V工V=111JU-111X指数函数.y=ae^©v=Inyu=Xry=◎严-X对数函数.y=戊+&Ln兀"v=yu=InXh，W型函数十v=—u=e~x4

在很多场合，可以对非线性模型进行线性化处理，尤其是可变换为线性的非线性回归，运用最小二乘法进行推断，对线性化后的线性模型，可以应用REG过程步进行计算。例1有实验数据如下：X1.11.21.31.41.51.61.71.81.922.12.22.32.4Y109.9540.4520.0924.5311.027.394.952.721.821.490.820.30.20.22试分别采用指数回归（y=aebx）方法进行回归分析。代码：dataexam25_1;inputxy;cards•;1.1109.951.240.451.320.091.424.531.511.021.67.391.74.951.82.721.91.8221.492.10.822.20.32.30.22.4•0.22;run;procsgplotdata=exam251;scatterx=xy=y;run;proccorrdata=exam251;varxy;run;datanewl;setexam25_1;v=log(y);run;procsgplotdata=newl;scatterx=xy=v;title'变量代换后数据'；run;procregdata=newl;varxv;modelv=x;printcli;title'残差图';plotresidual.*predicted.;run;datanew2;setexam25_1;y1=14530・28*exp(-4.73895*x);run;procgplotdata=new2;ploty*x=1y1*x=2/overlay;symbolv=doti=nonecv=red;symbol2i=smcolor=blue;title'指数回归图'；run;运行结果：

CORR过程2变量：kV筒单统计量N均值标准差总和最小值最大值XI41.750000.4133324500001100002.40000y16.13929294SD3722595DOOo.?noon109.95000Pearson相关系数，N_14

Prob>|r|underHO:Hho=0Xy1.00000■0.711040.0044-[)./11IMrooooo0.0()14變咼忙膜后歆据读取的观测数14使用的观测数E方差分析源自由度平方和均方F值Pr>F模型151.0910651.09106750.31<.0001误差1?0.810630.血麻校正合计1351.90170均方根逞差0.25991R方09844国塑量均值1.29084调整R方0.9631变异系数20,13498参数怙计值自由度参数估计值准差标误L值Pr>|t|Intercept19.583990.309-1530.97<.oaoiX1-473895().1m?-儿b()<.()001输出统计星观测因变量标准误差均值预测95%GL预测残差1-1.70001S7I20.I31S3.了砂5.00810.32892U./U013.89730.11753.27584.5137-0.1^/233.000230.KM12.0133-1一（）畑】-0,4231d3.19992.9+950.03202.34S73.55020.250452.39972.-1/bG0.08171.88193.0692-0.U/bLI62.0001?.00170.07411.4l?fi2.5305-0.00154971.59941.527B0.U/000.㈣132.11430.0/1681.00061.05390.(17000.-IfiM1.fi40d-0.0533甘0.59300.58000.0741-0.0088761.ICS90.OIOS100.3988[].10N10.0811-0.48750.6997(1.292711-0.1995-0.367(50.OJ1PO-096S50.23290.IC9312-1.2040-0.841f0.1041-1.4517-0.2317-0.362313-1.6094-1.3156o.urn-1.-o.cfi-n-0.加豹1d-1.5141-I.73950.1310-2.4244-1.15450.2754□.Hi-aa.H44-□.H31TOC\o"1-5"\h\z测U境掘2唉芒口由屢12N5EDDB75F?方0呛44出車R冉DWJ11.82.12.2程序说明：（1）调整后的R2=0.9831，说明拟合程度很好；F检验的P值=0.0001<a=0.05，拒绝原假设，故直线回归的斜率不为0；（2）将线性回归系数代入，得到原回归方程y=14530.28*e-4.73895x（3）残差图趋势，符合残差随机正态分布的假设（不带其它明显趋势）。二、多项式回归一般函数都可用多项式来逼近，故多项式回归分析可用来处理相当广泛的非线性问题。对观测数据（xt,yt）,t=l,...,N.多项式回归模型为：yt—+力］兀+玄X；+…++str匸12…N"令则模型可写为：Y=XB+s当X列满秩时，用最小二乘估计B=（XX）-1XY可求得其多项式回归方程。但由于（XX）-1的计算既复杂又不稳定，故一般采用正交多项式法来进行多项式回归。多项式模型可以直接应用GLM（广义线性模型）求解。例2重庆市种畜场奶牛群1—12月份(xl)，产犊母牛平均产奶量(y)的资料如程序数据步中，试对该资料配置一个合适的回归方程。代码：dataexam25_2;inputxly@@;x2=x1*x1;datalines;3833.433476.763811.583466.223769.473395.423565.743807.083481.993817.033372.823884.52•run;procsgplotdata=exam25_2;scatterx=x1y=y;title'原始数据散点图'；run;procregdata=exam25_2;modely=x1x2;run;运行结果：

原蜡妁埠敬点图3WD3TOD2.57.510.0原蜡妁埠敬点图3WD3TOD2.57.510.0REfl过程

核型:MODEI1

因娈量：/

读取的观测数吃

使用的观测裁壮方差分析源自由奩平方和均方F值Pr>F模型233259416629/1G.890.0009误差988S019844.56254校正合计II42II95均方根溟差99.21977R方0【伽因銮量均值3640.17167调整R方0.7429蛮异系数2.72509参数估计值翌量自由度参数怙计値标准溟差土值Pr>|t|Intercept14117.20136102.54G4940.15<.0001xl1-2049366836.2S837-5.650.0003115.785702,715895.810.om

GLM过程源自由度平方和均方1值Pr>1模型332594.3259166297,1630n.nans*1=1主恢差88601.06299844,5625楼正合计11421195.3388R方变异系数根MSt丫均值0.7896442.72550999.219773040.172源自由度I型SS均方F值Pr>FX1I11.005011.00500.000.9741x1*x1I332583.3210332583.321033.780.0003源自由庚I11型均方卜值Pr>1-xl1314325.5414314325.541431.930.0003x1*x11332583.3210332593.32100.0003参数怙计值标准误差tfiPr>|t|IntercEipt4117.201364102.546490540.15<.0001xl-204.93607836.2683725-A650.CKXJ3xl+x115.7856992.7158977L).81U.ZUS程序说明：（1）观察数据的散点图，更适合二次多项式拟合，也可以测试几种不同次数的多项式拟合选择其中最优的；（2）将回归系数代入多项式方程得到：y=4117.20136-204.93668x1+15.78570x12三、不能变换为线性的非线性回归该类非线性回归分析就是利用最小二乘准则来估计回归系数卩，使得残差平方和最小。一般来用数值迭代法来进行，先选定回归系数的初值卩0,按照给定的步长和搜索方向逐步迭代，直到残差平方和达到最小。有5种常用的非线性回归迭代方法:高斯-牛顿法（Gauss-Newton）、最速下降法（梯度法）、牛顿法（Newton）、麦夸特法（Marquardt）、正割法（DUD）。高斯-牛顿法在初值选取适当，且可逆时非常有效，但在其他情形，其求解较为困难，对此，Marguardt对其中的正则系数阵作适当修正，得到了改进算法。(二)PROCNLIN过程步对于不能线性化的非线性模型。其估计不能直接运用经典的最小二乘法，而需要运用其他估计方法，如加权最小二乘法、直接搜索法、直接最优法与Taylor级数展开法进行线性逼近。此时，可以利用NLIN过程步实现相应的计算，它是采用最小误差平方法及迭代推测法来建立一个非线性模型，估计参数默认采用高斯-牛顿迭代法。NLIN过程不保证一定可以算出符合最小误差平方法之标准的参数估计值。基本语法：PROCNLINdata=数据集</可选项>;PARMS参数名=数值;MODEL因变量=表达式</可选项>;BOUNDS不等式;>CONTROL变量;>der.参数名=表达式;><OUTPUTout=输出数据集</可选统计量>;>说明：NLIN的可选项包括outest=输出数据集输出每步迭代的结果；best=n只输出最好的n组残差平方和；method=gauss|marquardt|newton|gradient]dud|设定参数估计的迭代方法，默认为gauss(没有der.语句)；(2)PARMS语句指定参数并赋值，一般包括参数名、初始值(GridSearch可以帮助选择合适的初始值)、迭代准则；例如：parmsb0=0b1=1to10b2=1to10by2b3=1,10,100;bounds语句用于设定参数的约束，主要是不等式约束，约束间用逗号分隔。例如，boundsa<=20,b>30,1<=c<=10;der.语句用于计算模型关于各参数的偏导数，相应格式为:一阶偏导数：der.参数名=表达式；二阶偏导数：der.参数名.参数名=表达式；例如，对于模型modely=b0*(1—exp(-b1*x));二阶偏导数表达式：der・b0・b1=x*exp(-b1*x);例3根据对已有数据的XY散点图的观察和分析，发现Y随X增长趋势是减缓的，并且Y趋向一个极限值，我们认为用负指数增长曲线来拟合模型较为合适。代码：dataexpd;inputxy@@;datalines;0200.570300.720400.810500.870600.910700.940800.950900.971000.981100.991201.001300.991400.991501.001601.001700.991801.001901.002000.992101.00procnlindata=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;ploty*xygs*x/haxis=axis1vaxis=axis2overlay;symbol1i=nonev=pluscv=redh=2.5w=2;symbol2i=joinv=nonel=1h=2.5w=2;axis1order=20to210by10;axis2order=0.5to1・1by0.05;title1'y=b0*(1-exp(-b1*x)';title2'procnlinmethod=gauss';run;运行结果：NLIN过程因姿量网格搜索hOhi平方和1.【)1灿00400I).[)01d()1.0000o.oroo0.1.0000n.0.05521.00000.03000.06661.00000.07000.09731.00000.08000.13651.00000.09000.17081.£)0(10t).矽E】o.lias1.MXN)0.01000.1.()0(J()001002.If曲

NLIN过程因变量7

方法：Gauss-Newton迭代阶段迭代bOb1平方和0I00000.04000.0014010.99(510.0-1190.000^8020.S9620.00057730.S9620.04200.000577斗0,9962NOIL;Convergencecriterionmet.估计汇总方法Uauss-NciMton進代1R2.E52E1PPG(b1)1.O2ftE-0RPC(b1)394E7Object2.b(5F-1()□标0.000577读取的观测2014用的观测2()缺失的观测UNoteNote;AninterueplwtisriuLspeci1iedlotthismodeI.源自由度平方和均方FVa1ue近似Pr>F模型217.67178.8359275733<.0001误差180.0005770.000032未校正合计2017.6723参数怙计值近似标准淒差近似95%査信限bO0.驭为0.001610.99280