第三十四课 非线性回归分析
现实世界中严格的线性模型并不多见，它们或多或少都带有某种程度的近似；在不少情况下，非线性模型可能更加符合实际。由于人们在传统上常把“非线性”视为畏途，非线性回归的应用在国内还不够普及。事实上，在计算机与统计软件十分发达的令天，非线性回归的基本统计分析已经与线性回归一样切实可行。在常见的软件包中（诸如SAS 、SPSS 等等），人们已经可以像线性回归一样，方便的对非线性回归进行统计分析。因此，在国内回归分析方法的应用中，已经到了“更上一层楼”，线性回归与非线性回归同时并重的时候。
对变量间非线性相关问题的曲线拟合，处理的方法主要有：
● 首先决定非线性模型的函数类型，对于其中可线性化问题则通过变量变换将其线
性化，从而归结为前面的多元线性回归问题来解决。
● 若实际问题的曲线类型不易确定时，由于任意曲线皆可由多项式来逼近，故常可
用多项式回归来拟合曲线。
● 若变量间非线性关系式已知（多数未知），且难以用变量变换法将其线性化，则进
行数值迭代的非线性回归分析。
一、 可变换成线性的非线性回归
在实际问题中一些非线性回归模型可通过变量变换的方法化为线性回归问题。例如，对非线性回归模型
()t i t i t i t ix b ix a y εα+++=∑=2
1
0sin cos
(34.1)
即可作变换
t t t t t t t t x x x x x x x x 2sin ,2cos ,sin ,cos 4321====
将其化为多元线性回归模型。一般地，若非线性模型的表达式为
()()()t m m t t t x g b x g b x g b b y ++++= 22110
(34.2)
则可作变量变换
()()()
t m m t t t t t x g x x g x x g x ===*
2*21*1,,, (34.3)
将其化为线性回归模型的表达式，从而用前面线性模型的方法来解决，其中(34.3)中的x t 也
可为自变量构成的向量。
这种变量变换法也适用于因变量和待定参数 b i 。如
()[]1exp 2132211-++=t t t t t x x b x b x b a y
(34.4)
时上式两边取对数得
()1ln ln 2132211-+++=t t t t t x x b x b x b a y
(34.5)
现作变换
1,ln ,ln 2130*-===t t t t t x x x a b y y
(34.6)
则可得线性表达式
t t t t x b x b x b b y 3322110*+++=
(34.7)
利用前面方法确定了3,2,1,0,?=i b i ，并由)?exp(?0
b a =得到a ? 的值。 变量变换的线性化方法可推广到下列形式的非线性模型
()()t m m m t t x g b c x g b c b c y h )()()()(11100+++=
(34.8)
其中x =（x 1,x 2, ,x p ），而h （y t ）、c i （b i ）、g i （x t ）则分别化为新的因变量、线性回归参数
和自变量，即可归结为线性回归模型来解。见表34.1所示给出了一些常见的可线性化的非线性模型。
表34.1 典型的函数及线性化方法
函数名称
函数表达式
线性化方法
双曲线函数
x
b a y +=1 y
v 1
=
x
u 1=
幂函数
b ax y =
y v ln = x u ln = 指数函数
bx ae y =
y v ln = x u = x b ae y /=
y v ln = x
u 1=
对数函数
x b a y ln += y v = x u ln =
S 型函数 x
be a y -+=
1
y
v 1=
x e u
-=
当曲线的函数类型未确定时，我们常采用上述非线性模型作为其拟合曲线，即将自变量的各种初等函数的组合作为新自变量，用逐步回归法（或正交筛选法等）对新变量进行筛选，以确定一个项数不多的线性函数表达式。 该方法对表达式形式没限制且精度要求不高的问题颇为有效。
二、 多项式回归分析
在式(34.2)中，若取()i i x x g =，则为多项式回归模型。由数学分析知识可知，一般函数都可用多项式来逼近，故多项式回归分析可用来处理相当广泛的非线性问题。
对观测数据（x t ，y t ）（t = 1，…，N ），多项式回归模型为
t m t m t t t x b x b x b b y ε+++++= 2210，t =1,2, ,N
令
?
?
??
??
??????=N y y y Y 2
1，???
???
?????
???=m N N
N
m m x x x x x x x x x X 2222212
11111，????????????=m b b b B 10，??????
??????=N εεεε 21 则模型可表示为
ε+=XB Y
当X 列满秩时，由前面的讨论知，其最小二乘估计为
()Y X X X B ''=-1?
由此即可求得其多项式回归方程。但由于()1
-'X X 的计算既复杂又不稳定，故我们一般采
用正交多项式法来进行多项式回归。
三、 不可变换成线性的非线性回归分析
假设因变量y 与自变量（x 1，x 2,…，x p ）之间满足非线性模型
()εβ+=;,,,21p x x x F y
(34.9)
其中()'
=m ββββ,,,21 为未知参数，F 为已知表达式，仅β未知的非线性函数，ε 为误差项。
现将观察数据
()
pt t t t x x x y ,,,,21 ， t =1,2, ,N
代人上式(34.9)得非线性回归模型
()
t pt t t t x x x F y εβ+=;,,,21 ， t =1,2, ,N
常记为
E F Y +=)(β
其中()'=N y y y Y ,,,21 为y 的观察向量，()'
=m βββ,,1 为非线性回归系数，E =()'
N εεε,,,21 为观察误差向量，F 为未知参数β的函数向量。非线性回归分析就是利用
最小二乘准则来估计回归系数β，即求β
? 使得残差平方和
()()()()()βββF Y F Y E E Q -'
-='=
2
121 在
β
β?= 处达到最小。 非线性回归分析一般来用数值迭代法来进行，其共同特点是：由选定β的初值0
β出发，
通过逐步迭代
??+=t 0ββ
(34.10)
即选择适当的步长t ( >0 ) 及确定搜索方向向量?＝（?1，?2，…，?m ），使得
()()
0ββQ Q <
(5.4.11)
再β由取代0β，重复上述迭代过程，直至 Q (β)可认为达到最小值为止，即可将所得的β
作为其最小二乘估计β
?，从而得到非线性回归方程
()
β?;,,,?21p x x x F y
= 1. 下降方向和步长的选择
首先考察()()()()()βββF Y F Y E E Q -'
-='=
2
121的梯度向量（即导数）
()()()()ββββF Y G F Y F Q
-'-=-'???
? ????-=?? 其中'
???
?
??????=??=m F F F G βββ,,1 为F 的梯度矩阵。 为使
0β迭代收敛到β
?，其迭代公式应满足下降性质(5.4.11)。现考虑一元函数()()??+=t Q t 0β?，它从0β出发以 ?为方向的射线上取值。由复合求导公式得
()()()??'
--=??'
???
? ????='==G F Y Q t d t ββ?0
可以证明，当 d <0 时，在以 ?为方向向量的射线上可以找到??+=t 0
ββ，使得
()()
0ββQ Q <。我们将满足 d <0 的?称为下降方向，Bard 于1974年给出了?为下降方向的
充要条件为
()()βF Y G P -'=?
其中P 为对称正定阵，由此我们可得下降算法的迭代公式为
()()βββF Y G tP -'+=0
(34.12)
其中P 为任意正定阵，G 为F 的梯度，t 为满足()()0
ββQ
Q <的正实数，即步长。
如何计算?以便修改参数向量β有五种常用的非线性回归迭代方法：高斯-牛顿法（Gauss-Newton ）、最速下降法（梯度法，Gradient ）、牛顿法（Newton ）、麦夸特法（Marquardt ）、
正割法（DUD ）。以下我们介绍其中高斯-牛顿法。
2. Gauss －Newton 法
首先选取β的一切初始近似值0β，令0ββ-=?，则只要确定?的值即可确定β。为此，考虑)(βF 在0β处的Taylor 展开式，并略去二次以上的项得
()()()
()
??'+=??'?
??
? ?
???+=?+==G F F F F F 0
0βββ
ββββ 其中0
βββ
=??=
F G 为F 的梯度。此时其残差平方和
()()()()
?'--'?'--=
G F Y G F Y Q 002
1
ββ 由
0=?
??Q
，得其?的正则方程为 ()()()0βF Y G G G -'=?'
(34.13)
故
()()()
01
βF Y G G G -''=?-
(34.14)
由此即可用前面线性回归法求?，只需将G 、)(0
βF Y -视为前面（5.2.1）式中的X 、Y
即可。此时，对给定精度1ε、2ε ，当{}
1max ε
或()
20
εβ
二乘法估计?+=0?ββ；否则用所得的β?代替0β，重复上述步骤，直至i
?或Q (β)满足精度要求为止。该法称为Gauss －Newton 法，其一般迭代公式为
?+=+i i i t ββ1
(34.15)
其中：?为()()()()()i
i
i
i
F Y
G G G ββββ-'=?'
的解，t i
为()()
??+=t Q t i
β
?的最小值点。
Gauss-Newton 法在初值0
β选取适当，且G G '可逆时非常有效，但在其他情形，其求解 较为困难，对此，Marguardt 对(34.14)中?的正则系数阵作适当修正，得到了改进算法。
四、 nlin 非线性回归过程
在很多场合，可以对非线性模型进行线性化处理，尤其是关于变量非线性的模型，以运
用OLS 进行推断。对线性化后的线性模型，可以应用SAS 的reg 过程进行计算。
多项式模型可以直接应用glm （广义线性模型）求解。对于不能线性化的非线性模型。其估计不能直接运用经典的最小二乘法，而需要运用其他估计方法，如直接搜索法、直接最优法与Taylor 级数展开法进行线性逼近。此时，可以利用SAS/STA T 的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语句中的主要选择项。
●outest＝数据集名——指定存放参数估计的每步迭代结果的数据集名。
●best＝n——要求过程只输出网格点初始值可能组合中最好的n组残差平方和。
●method＝gauss | marquardt | newton| gradient| dud |——设定参数估计的迭代方法。缺省时为gauss，除非没有der.语句。
●eformat——要求所有数值以科学记数法输出。
●nopoint——抑制打印输出。
●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*(1－exp(-b1*x))；
5. bounds语句。
用于设定参数的约束，主要是不等式约束，约束间用逗号分隔。例如，bounds a<=20，b>30，1<=c<=10；
6. der.语句。
除非在proc nlin语句中指明所用的迭代法是dud，使用选择项method＝dud，否则der.语句是必需的。der.语句用于计算模型关于各个参数的偏导数，相应的格式为：
一阶偏导数der．参数名＝表达式；
二阶偏导数 der ．参数名．参数名＝表达式；
例如，对于model y=b0*(1－exp(-b1*x))；der.语句的书写格式为der.b0=1－exp(-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 预测值 stdp clm 的标准差 u 95 95％cli 上限 residual | r 残差 stdr 残差的标准差 l95 95％cu 下限 parms 参数估计值 l95m 95％clm 下限 student 学生氏残差 sse | ess
残差平方和
u95m
95％clm 上限
h
杠杆点统计量h i
关于 nlin 过程的其他选择项及意义，详见SAS ／STAT 的用户手册。
五、 实例分析
例34.1 负指数增长曲线的非线性回归。根据对已有数据的XY 散点图的观察和分析，发现Y 随X 增长趋势是减缓的，并且Y 趋向一个极限值，我们认为用负指数增长曲线
)1(10x b e b y --=来拟合模型较为合适。程序如下：
data expd ;
input x y @@;; cards;
020 0.57 030 0.72 040 0.81 050 0.87 060 0.91 070 0.94 080 0.95 090 0.97 100 0.98 110 0.99 120 1.00 130 0.99 140 0.99 150 1.00 160 1.00 170 0.99 180 1.00 190 1.00 200 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 语句中非线性表达式)1(10x
b e
b y --=中参数0b 和1b 的最小二乘估计。我们知道参
数初始值选取好坏，对迭代过程是否收敛影响很大。parms 语句设置了初始值网格值为0b 取0，0.5，1，1.5，2共5个值，1b 取0.01，0.02，…，0.09共9个值，所有可能组合为5×9=45,选项best=10要求输出残差平方和最小的前10种组合。高斯-牛顿迭代算法要求给出模型
)1(10x b e b y --=对参数0b 和1b 的一阶偏导数表达式，我们知道
x
b x
b xe b b y
e b y
1101
d d 1d d --=-=
der.语句用以表示上面两个一阶偏导数表达式。output 语句输出一个新数据集expout ，包括原数据集和非线性回归模型的预测值ygs 。gplot 过程的主要作用是绘制输出数据集expout 中的原始数据的散点图及回归曲线的平滑线。程序的输出结果见图34-1和见表34.3所示。
图34-1 XY 散点图和非线性回归曲线
表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的渐近估计值，因此得到的非线性回归模型为
)1(9961885657.00419538868.0x e y --=
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.165290
Non-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.000577 NOTE: 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
同时还给出B0和B1参数估计的渐近有效标准差（Asymptotic Std. Error）和渐近95%置信区间（Asymptotic 95 % Confidence Interval）。最后还给出了参数的渐近相关阵（Asymptotic Correlation Matrix）。