第七章 非线性回归模型
前面六章讲的各种回归模型,不管变量是一元还是多元、连续还是离散,不管误差是什
么形态,不管有没有什么约束,都是线性回归模型,即变量是以一次函数形态出现。这一章
我们将研究非线性回归模型。
经济模型本来就存在许多非线性形式,我们在第一章与第二章就曾经处理过“可以线性
化的非线性模型”,即经过简单函数变换后可以化为一元或多元线性回归模型的非线性回归
模型。但是在一般情况下,非线性模型难以精确线性化,这就需要予以特别的考虑。
一般的非线性回归模型可以表示为
Y
=
(
,Xf
) εβ +
(7.0.1)
这里 X 是可观察的独立随机变量,β是待估的参数向量,Y 是独立观察变量,它的平均数
依赖于 X 与β,ε是随机误差。函数形式 (·)是已知的。
f
Cobb-Douglas 生产函数是非线性回归模型的典型例子:
Q
=
α
ββ
1 KL
2
+
ε
(7.0.2)
这里Q 是经济部门的产出, L 是劳动力投入, K 是资本投入,待估参数是α、 1β 与 2β 。
定义
Y Q=
,
X
′ =
L K
( ,
)
,
)
β αβ β ′
=
(
,
,
1
2
,以及 (
f X
,
β α=
L Kβ β
1
1
1
2
)
,则 Cobb-Douglas
生产函数就可以写为(7.0.1)的形式。
另一个例子是消费函数
C
=
ββ
2
+
1
Y
β +
3
ε
(7.0.3)
这里Y 是居民收入, 是居民消费。其中参数
C
3β 的估计问题就很有必要。如果贸然假定
3 1β = ,那就是线性函数了,可是实际数据也许会否定 3 1β = 。
有些经济模型到底能不能线性化,取决于误差项的假定。例如 Cobb-Douglas 生产函数,
如果将误差假定为与函数部分相乘,即
Q
α=
L K e
β β ε
1
2
则取对数后可以线性化:
ln
Q
=
ln
βα
1
+
ln
L
+
β
2
ln
K
+
ε
272
(7.0.4)
(7.0.5)
另一方面,有些线性回归模型也可以视为非线性问题,例如广义最小二乘问题
XY
=
E
,
εβ
+
( )
ε
=
,0
Var
( )
2
Ψ
σε
=
(7.0.6)
的极大似然估计就可以看作非线性问题。
本章就讨论这些非线性回归模型的性质与计算问题,涉及到一些大样本理论,介绍了非
线性强度度量的几何意义。作为特别的非线性回归模型,重点介绍了增长曲线模型与失效率
模型。
非线性回归模型的解法主要有两个。一个是最小二乘法,即求向量Y 与集合
( βXf
)
,
的
最短距离:
||
Y
−
ββXf
(
⎯→⎯
||)
,
min
一个是极大似然法,即假设误差的分布密度函数
合密度函数),再求其最大值:
2σβXg
(
,
,
)
已知,作似然函数(样本的联
L
2
(
σβ
,
)
n
= ∏
i
1
=
iXg
(
,
,
σβσβ
⎯⎯ →⎯
)
2
,
2
max
可见解法都归结到求极值问题。本章第一节介绍了两种求极值的计算方法,这两种方法都需
y
Δ
要计算目标函数的导函数 x
Δ
。如果是等距的 xΔ ,自然没有什么问题,可是数理统计里的
数据是随机的, 未必等距,甚至有时为 0,这将由于分母为 0 而使计算溢出。我们在实
xΔ
际编程计算时求极值是采用 Sargent 改进的 Powell 算法,它不需要计算导函数或差分,计
算稳定。
因此似乎可以说,本章前三节更多的是其理论上的意义。本章比较有实际意义的是第四
节的增长曲线模型和第五节的失效率模型。
增长曲线模型可以用来拟合增长曲线,特点是增长速率开始较慢,中期较快,后期又较
慢且存在极限。上一章的 Logit 模型就属于这一类模型,这一章我们将了解更多的增长曲线
模型。
失效率模型用于生存分析,关键是理解寿命数据
ti
,{
i
,2,1
L=
n
},
、它的密度函数
f
)(t
和分布函数
)(tF
,以及失效率
tr
)(
=
f
−
t
)(
tF
)(
1
。本书作者利用自己对密度的核估计的研究
成果,提出了失效率的图估计方法,再利用
)(tr
在分段 Weibull 分布上的特点,拟合浴盆曲
线,在生存分析领域有较好的实际意义。
273
第一节 非线性回归模型最小二乘估计的计算
为了引进非线性回归的最小二乘方法(Gauss-Newton 算法),我们先考虑一个简单的单
参数模型:
=
Y
i
(
Xf
,
i
定义残差平方和
)
X
βεβ
=
+
i
2
β
X
+
i
1
+
ε
i
,
2
i
E
(
ε
i
)
=
,0
Var
(
2
σε
=
)
i
(7.1.1)
S
(
)
β
=
=
n
∑
i
1
=
n
∑
i
1
=
2
ε
i
=
n
∑
i
1
=
[
Y
i
−
(
Xf
i
]
)
2
,β
(
Y
i
−
2
X
1 ββ
−
i
2
)
X
2
i
(7.1.2)
回归的原则还是要使残差平方和最小,于是对 (
S β 求导得:
dS
d
β
=
2
n
∑
i
1
=
[
Y
i
−
(
Xf
,
i
⎛
]
)
β
⎜⎜
⎝
−
(
Xdf
,
i
d
β
)
)
β
⎞
⎟⎟
⎠
n
= ∑
2
i
1
=
Y
[
i
−
整理得:
2
X
ββ
−
i
1
X
2
i
[ ]
−
X
i
1
−
X
2
β
0]
=
2
i
(7.1.3)
3
2
β
n
∑
i
1
=
X
2
i
2
2
3
β
+
n
∑
i
1
=
XX
i
1
2
i
+
⎛
⎜
β
⎜
⎝
n
∑
i
1
=
X
2
i
1
−
2
n
∑
i
1
=
YX
i
i
2
⎞
⎟
−⎟
⎠
n
∑
i
1
=
YX
i
i
1
=
0
(7.1.4)
这是关于β的三次方程,β有 3 个可能的解,将这三个解分别代入 (
S β ,取使 (
S β 最小
)
)
的那个β为回归模型的最终解。
这 一 节 就 专 门 介 绍 非 线 性 回 归 模 型 参 数 最 小 二 乘 估 计 的 计 算 方 法 , 一 个 是
Gauss-Newton 算法,一个是 Newton-Raphson 算法。为了清楚易懂,我们都是先介绍单参数
的,再介绍多参数的。
一、非线性模型 LSE 的 Gauss-Newton 算法
先考虑单参数的非线性回归模型( (·)已知):
f
Y
i
=
(
Xf
,
)
εβ +
i
i
(7.1.5)
其残差平方和函数为
274
S
(
)
β
=
n
∑
i
1
=
[
Y
i
−
(
Xf
i
]
)
2
,β
(7.1.6)
要使 (
S β 取极小值,其一阶条件为
)
dS
d
β
n
= ∑
2
i
1
=
[
Y
i
−
(
Xf
,
i
⎛
]
)
β
⎜⎜
⎝
−
)
β
(
Xdf
,
i
d
β
⎞
=⎟⎟
⎠
0
(7.1.7)
现在的问题是要求出上述方程的解β,并且判断出整体最小值解β。
一个近似办法是用 (
f X β 的一阶 Taylor 近似展开去代替 (
f X β 。设β的初值为
)
)
,
,
i
i
1β ,则在 1β 点附近函数 (
(
f X
,
i
于是求得导数值为
i
,
)
f X β 有近似 Taylor 展式
)
β
+
(
df X
d
,
i
β
)
β
≈
(
f X
i
,
β
1
)
(
)1
β β−
β
1
)
β
(
df X
d
,
i
β
≈
β
1
(
f X
i
,
)
(
f X
β
β β
1
−
−
,i
β
1
)
(7.1.8)
(7.1.9)
简记
则
Z
i
(
β
1
)
=
)
β
(
df X
d
,i
β
β
1
(7.1.10)
S
(
)
β
=
n
∑
i
1
=
[
Y
i
−
(
Xf
,
β
1
i
)
−
Z
i
(
βββ
1
)(
−
1
]
)
2
=
这里
n
[∑
~
Y
i
i
1
=
(
β
1
)
−
Z
i
]
(
ββ
)
⋅
1
2
(7.1.11)
)
~
(
Y
β
i
1
(~
1βiY
残差平方和正是线性回归
对于给定的初值 1β ,
=
Y
i
−
(
Xf
,
β
1
i
)
+
Z
i
(
ββ
1
)
1
⋅
(7.1.12)
)
以及
iZ β 都是确定的,可计算的。于是(7.1.11)所表达的
1(
)
~
Y
i
(
β
1
)
=
Z
i
(
εββ
i
+
)
1
(7.1.13)
的残差平方和。Malinvaud(1980)将上式称作拟线性模型,其最小二乘估计是
275
~
Y
i
n
∑
i
1
=
n
∑
i
1
=
(
β
1
)
⋅
Z
i
(
β
1
)
(
Z
i
(
β
1
)
)
2
β
2
=
这里
(
Z
=
(
ββ
1
1
Z
(
′
)
1
−
)
)
Z
(
)1
ββ
1
(
′
)
~
Y
(7.1.14)
Z
(
β
1
)
=
Z
(
β
1
)
1
M
(
β
1
)
nZ
⎞
⎟
⎟
⎟
⎠
⎛
⎜
⎜
⎜
⎝
(~
Y
β
1
)
=
,
(~
Y
β
1
1
M
(~
nY
β
1
)
)
⎛
⎜
⎜
⎜⎜
⎝
⎞
⎟
⎟
⎟⎟
⎠
(7.1.15)
于是我们看到,如果我们有待估参数β的一个初值 1β ,就可以得到β的一个新值 2β 。
重复使用这个方法,又有一个拟线性模型
) εββ
(
~
Y
Z
+
=
)
(
β
2
2
其解为
β
3
⎡
= ⎢
⎣
Z
(
β
2
′
)
Z
(
β
2
)
1
−
⎤
⎥
⎦
Z
(
β
2
′
)
Y
%
(
)2
β
(7.1.16)
(7.1.17)
这样得到一个序列 1
β β βL L (注意它们都是向量)。我们可以写出一般迭代表达式
n
2
,
,
,
这里
f X
(
,
)
β
=
[
f X
(
1
,
),
β
f X
(
,
),
β
2
,
f X
(
n
L
,
由于 (
)
S β 取极小值的一阶条件(7.1.7)可写作
′
(
)
β
[
Y
Z
−
(
Xf
]
)
β
,
=
0
故若在迭代过程中有 1n
β β+ = ,则由(7.1.18)知必有(7.1.19)成立,即
n
(7.1.19)
0=
,此时 (
S β
)
dS
βd
取得一个极值。
下面要考虑的问题是,这样的迭代解是使 (
S β 取极小值还是极大值?如果是极小值,
)
它是不是整体最小值?我们可以用不同的迭代初值去计算,如果不同的迭代初值导致不同的
极小值 ,则其中最小的那个 是整体最小值。至于为什么这些迭代过程只会导致极小而
S
S
不是极大,可以分析
dS
βd
的表达式
276
β
n
1
+
=
Z
(
ββ
n
Z
(
n
(
(
)n
ββ
(
n
′
)
~
Y
′
)
′
)
1
−
Z
1
−
)
)
)
⎤
⎥
⎦
=
Z
⎡
⎢
⎣
β
n
Z
(
β
n
=
β
n
+
Z
(
ββ
n
Z
(
n
′
)
1
−
)
)
Z
(
(
Z
(
β
n
′
)
Y
⎡
⎣
−
(
f X
,
β
n
)
+
Z
(
β
n
)
⋅
β
n
⎤
⎦
[
Y
−
(
Xf
,
β
n
] (7.1.18)
)
′
)
(
β
n
]
)
β 。
dS
d
β
2
−=
Z
′
)
(
β
[
Y
−
(
Xf
]
)
β
,
于是迭代公式(7.1.18)可写作
ββ
n
1
+
=
n
(
Z
1
2
−
(
ββ
n
Z
(
n
′
)
1
−
)
)
(7.1.20)
(7.1.21)
为正时, 1n
β β+ < ,于
n
n
dS
d
ββ
dS
d
ββ
n
由于 Z Z′ 是一些平方和,它总是正的,
(
)Z Z −
1
′ 也是正的。当
是函数值 将走向极小;当
S
dS
d
ββ
n
为负时, 1n
β β+ > ,于是函数值还是走向极小。
n
图 7.1.1.1
为了避免迭代时间过长或迭代来回反复,可以引进步长控制函数 ,
nt
ββ
n
1
+
=
n
(
Zt
n
−
(
ββ
n
Z
(
n
′
)
1
−
)
)
dS
d
ββ
n
=
β
n
+
⎡
t Z
2
⎢
n
⎣
(
β
n
′
)
Z
(
β
n
)
Z
(
β
n
′
)
−
1
⎤
⎥
⎦
Y
⎡
⎣
−
(
f X
,
β
n
)
⎤
⎦
(7.1.22)
nt 由计算程序根据误差自动调整。上述算法一般称为 Gauss-Newton 算法。
由于非线性模型函数形式复杂,一般难以建立有限样本的统计性质。然而我们可以考虑
它的渐近性质。一般来说 , 是一致估计,
βˆ
2
(
*
ZZσ
'
,方差为
0
分布去作近似:
) 1
−n
*
/
,其中
Z
* =
lim
(
)ββ ˆ−
n
)βZ
(
的极限分布为正态分布,平均数为
。于是在作假设检验时,可以用渐近正态
~ˆ
β
N
⎛
⎜⎜
⎝
2ˆ,
⎤
)
ββσβ
⎥⎦
⎡
⎢⎣
Z
Z
(
(
′
)
ˆ 2
σ
=
(
)
S β
n
1
−
−1
⎞
⎟⎟
⎠
(7.1.23)
(7.1.24)
Z
具体检验过程与线性回归的假设检验是一样的,不过以 (
Zβ β′
)
)
(
代替 X X′ 。
277
再考虑多参数非线性回归模型( (·)已知):
f
Y
i
=
(
Xf
,
)
εβ
i
+
i
,
E
(
ε
i
)
=
,0
Var
(
2
σε
=
)
i
其中β是
1m×
未知参数。采用矩阵记法,模型是
Y
=
(
,Xf
) εβ +
E
( )
ε
=
,0
Var
( )
2
σε
=
nI
残差平方和为
(
)
=′=
εεβ
它取极小值的一阶条件为
S
[
Y
−
(
Xf
′
]
)
β
[
Y
,
−
(
Xf
]β
)
,
S
∂
∂
β
−=
2
)
β
(
Xf
,
∂
∂
β
⎡
⎢
⎣
′
⎤
⎥
⎦
[
Y
−
(
Xf
]
)
β
,
=
0
这里
′
⎛
⎜⎜
⎝
f
⎞
∂
⎟⎟
∂
β
⎠
是一个 (
m n×
)
的矩阵。记
(n m× 矩阵
)
(7.1.25)
(7.1.26)
(7.1.27)
(7.1.28)
Z
(
)
β
=
)
β
Xf
(
,
∂
∂
β
=
⎛
⎜
⎜
⎜
⎜
⎜⎜
⎝
)
β
)
β
Xf
(
,
∂
1
∂
β
1
L
Xf
,
(
∂
n
∂
β
1
L
O
L
)
β
)
β
Xf
(
,
∂
1
∂
β
m
L
Xf
(
,
∂
n
∂
β
m
(7.1.29)
⎞
⎟
⎟
⎟
⎟
⎟⎟
⎠
注意这里带下标的 iβ 是分量,则残差平方和取极小值的一阶条件为
[
Z
′
]
(
)
β
[
Y
−
(
Xf
]
)
β
,
=
0
(7.1.30)
为了使用 Gauss-Newton 算法,需要对多元函数 (
f X β 在初值 1β 附近作多元 Taylor
)
,
展开:
(
Xf
i
)
β
,
≈
(
Xf
,
β
1
i
)
+
)
β
(
Xf
,
∂
i
∂
β
⎡
⎢
⎢
⎣
β
1
′
⎤
⎥
⎥
⎦
(
ββ
1
−
)
(7.1.31)
将这些展式合并在一起有
(
Xf
)
β
,
≈
(
Xf
,
β
1
)
+
Z
(
)1
βββ
)(
−
1
这样非线性回归模型(7.1.26)就成为
Y
≈
(
Xf
,
β
1
)
+
Z
(
) εβββ
)(
+
−
1
1
(7.1.32)
(7.1.33)
278
如果记
~
Y
(
β
1
)
=
Y
−
(
Xf
,
β
1
)
+
Z
(
) 1
ββ
1
则非线性模型就线性化了
(
Y
%
n
1
×
)
Zβ
1
n k
×
=
(
)
β β ε
n
1
×
+
1
×
1
k
(7.1.34)
(7.1.35)
我们可以写出这个线性回归模型的最小二乘解,也就是原非线性回归模型的第一次迭代解
β
2
⎡
= ⎢
⎣
Z
(
β
1
′
)
Z
(
β
1
)
−
1
⎤
⎥
⎦
Z
(
β
1
′
)
Y
%
(
β
1
)
=
β
1
+
Z
(
β
1
′
)
Z
(
β
1
)
⎡
⎢
⎣
Z
(
β
1
′
)
−
1
⎤
⎥
⎦
Y
−⎡
⎣
(
f X 1
β
,
)
⎤
⎦
(7.1.36)
继续这个过程,我们得到一般迭代形式
β β
n
=
1
+
n
Z
(
β
n
′
)
Z
(
β
n
)
+
⎡
⎢
⎣
Z
(
β
n
′
)
−
1
⎤
⎥
⎦
Y
−⎡
⎣
(
f X
,
β
n
)
⎤
⎦
如果出现了 1n
β β+ = ,则意味着
n
Z
(
β
n
′
)
[
Y
−
(
Xf
,
β
n
]
)
=
0
(7.1.37)
(7.1.38)
此即意味着
S
∂
β
∂
0=
。同样可以分析出,极值必定是极小值。虽然我们不能保证整体最小值
已经被发现,但总是可以通过改变迭代初值而获得不同的极小值并加以比较。这样整体最小
值被遗漏的机会就大大减少了,计算程序对多元情形也有迭代步长的设置选择。
在合适条件下,β的最小二乘估计 ˆβ的渐近分布是正态分布,平均数为β,方差阵为
ˆ
=Σ
)ˆ(
ˆ
ββσ
)ˆ(
′
Z
2
(
Z
−
) 1
(7.1.39)
这里
ˆ 2
σ
=
( )
ˆ
S
β
mn
−
(7.1.40)
这样我们可以作出假设检验。
βˆ
关于 的渐近分布为正态的合适条件,可以从三个方面考虑。
首先是残差假定。我们已经假定 iε是 i.i.d.样本。平均数为 0,方差为 2σ ,这对于保证
βˆ
的渐近分布为正态已经够了。
其次是函数 (
f X β 的假定。从分析过程可以看到,我们需要假定 (
)
,
f X β 关于 iX 是
)
,
i
i
连续的,关于β有二阶连续导数。
279