非线性拟合技巧:分离线性参数与数据标准化
例题
对非线性模型
y = p 1 + p 2 x 1 + p 3 x 2 + p 4 x 3 + p 5 x 4 1 + p 6 x 1 + p 7 x 2 + p 8 x 3 + p 9 x 4 y = \frac{p_1 + p_2 x_1 + p_3 x_2 + p_4 x_3 + p_5 x_4}{1 + p_6 x_1 + p_7 x_2 + p_8 x_3 + p_9 x_4} y=1+p6x1+p7x2+p8x3+p9x4p1+p2x1+p3x2+p4x3+p5x4
的一次观测数据的 csv
格式文本如下所示,各列分别为
x
1
,
x
2
,
x
3
,
x
4
,
y
x_1, x_2, x_3, x_4, y
x1,x2,x3,x4,y,以最小化残差平方和为目标拟合参数
p
1
p_1
p1 至
p
9
p_9
p9。一个较好结果的残差平方和为 1.650522
。
15100,29000,508,180,3.4
20500,43350,453.7,141,3
80000,92610,487.9,132,2.7
91500,142775,572.3,182,3.37
82500,2123160,455.7,113,6.89
20000,227800,481.3,170,5.03
17800,140000,541.3,179,3.55
3900,15980,538.6,186,2.72
17300,223200,460.6,100,4.05
25700,229400,393.1,133,3.22
49400,424500,373.9,106,2.65
40700,561700,482.8,107,1.91
77000,563600,482.1,140,3
92900,557600,415.1,121,1.31
63300,528300,536.7,144,2.33
51600,488940,385.1,154,3.55
60000,480500,412.2,111,3.37
70000,530500,567.2,139,2.55
解
使用 R
语言进行求解,使用 nloptr
包中的 neldermead
函数,以参数全为 1
为起始点尝试求解,具体代码如下所示。
# 代码不能直接运行,详情请参阅各处注释
library(nloptr)
# 在此处编写代码将数据读入至 data 中
loss = function(p) {
sum(((p[1] + p[2] * data[, 1] + p[3] * data[, 2] + p[4] * data[, 3] + p[5] * data[, 4])
/ (1 + p[6] * data[, 1] + p[7] * data[, 2] + p[8] * data[, 3] + p[9] * data[, 4]) - data[, 5]
) ^ 2)
}
s = neldermead(rep(1, 9), loss, control = list(maxeval = 1e4))
然后反复执行 s = neldermead(s$par, loss, control = list(maxeval = 1e4))
直到结果不再进一步优化为止,得到一组结果的残差平方和为 3.493562
。
分离线性参数
设非线性模型
y = f ( x ⃗ ; p ⃗ ) y = f(\vec x; \vec p) y=f(x;p)
其中 x ⃗ ∈ R n \vec x \in \mathbb R^n x∈Rn 是自变量, y ∈ R y \in \mathbb R y∈R 是因变量, f : R n → R f:\mathbb R^n \rightarrow \mathbb R f:Rn→R 是一映射, p ⃗ ∈ R m \vec p \in \mathbb R^m p∈Rm 是 f f f 的参数向量。若 y y y 对 p ⃗ \vec p p 的某分量 p i p_i pi 可视作一个线性函数,则称 p i p_i pi 是线性参数。若固定 p ⃗ \vec p p 的所有非线性参数,则 p ⃗ \vec p p 余下的参数均为线性参数,这些参数可以使用线性拟合的方法求解。根据线性代数的知识,当解存在时,线性拟合方法总是求得最优解。
在该例题中,注意到
p
1
p_1
p1 至
p
5
p_5
p5 是线性参数,则只需拟合
p
6
p_6
p6 至
p
9
p_9
p9,而
p
1
p_1
p1 至
p
5
p_5
p5 则交由 lm
方法拟合即可。
数据标准化
在该例题中,观察到各变量之间的数量级相差较大。使用 x 1 x_1 x1 至 x 4 x_4 x4 的标准分进行拟合,具体代码如下所示。
# 代码不能直接运行,详情请参阅各处注释
# 在此处编写代码将数据读入至 data 中
sd = scale(data[, 1:4])
model = function(p) {
sd[, 1:4] = sd[, 1:4] / (1 + p[1] * sd[, 1] + p[2] * sd[, 2] + p[3] * sd[, 3] + p[4] * sd[, 4])
lm(data[, 5] ~ sd[, 1] + sd[, 2] + sd[, 3] + sd[, 4])
}
loss = function(p) {
sum(model(p)$residuals ^ 2)
}
s = neldermead(rep(1, 4), loss)
e = s$par / attr(sd, "scaled:scale")
p = e / (1 - c(e %*% attr(sd, "scaled:center")))
无需反复使用 neldermead
求解,可以得到符合题目要求的结果。实际上
p
6
p_6
p6 至
p
9
p_9
p9 的一个较好的结果是 -2.235640774844147663555e-06, 7.302384225577161910251e-08, -2.128725631770886621980e-03, 3.161796617995193568393e-04
。可以看到它们的数量级都很小。事实上对于大部分的局部优化算法来说,起始点都必须选择在这个解的附近才能得到该结果。以下代码将起始点的每个分量加上一个均匀分布于
[
−
3
×
1
0
−
7
4
,
3
×
1
0
−
7
4
]
\left[-\dfrac{3 \times 10^{-7}}{4}, \dfrac{3 \times 10^{-7}}{4}\right]
[−43×10−7,43×10−7] 的随机数,成功求解的概率已经低于 95%
,而随机数的区间若选择
[
−
1
0
−
5
,
1
0
−
5
]
[-10^{-5},10^{-5}]
[−10−5,10−5],则成功概率不大于 5%
。
# 代码不能直接运行,详情请参阅各处注释
# 在此处编写代码将数据读入至 data 中
model = function(p) {
data[, 1:4] = data[, 1:4] / (1 + p[1] * data[, 1] + p[2] * data[, 2] + p[3] * data[, 3] + p[4] * data[, 4])
lm(data[, 5] ~ data[, 1] + data[, 2] + data[, 3] + data[, 4])
}
loss = function(p) {
sum(model(p)$residuals ^ 2)
}
cl = makeCluster(max(detectCores() - 1), 1)
clusterEvalQ(cl, library(nloptr))
clusterExport(cl, c('p', 'loss', 'model', 'data'))
sum(parSapply(cl, 1:100, function(x)
neldermead(p + runif(4,-3e-7 / 4, 3e-7 / 4), loss)$value < 1.66))
stopCluster(cl)
若数据中各变量之间的数量级差别较大,此时可以在保持模型等价的情况下使用一个或几个变量的标准分,增加求解成功的概率。
在一组实数据中,实数据点 x x x 的标准分定义为
z x = x − x ‾ σ x z_x = \frac{x - \overline x}{\sigma_x} zx=σxx−x
其中 x ‾ \overline x x 是该组数据的平均值, σ x \sigma_x σx 是该组数据的标准差。
使用标准分后,模型也相应地发生变化。例如使用自变量与因变量的标准分,则模型变为
z y = f ( z ⃗ ; p ′ ⃗ ) z_y = f(\vec z;\vec{p'}) zy=f(z;p′)
其中 z ⃗ \vec z z 的各分量是 x ⃗ \vec x x 各分量的标准分, z y z_y zy 是 y y y 的标准分,拟合的参数向量 p ′ ⃗ \vec{p'} p′ 维度与 p ⃗ \vec p p 相同,但 p ′ ⃗ ≠ p ⃗ \vec{p'}\neq\vec p p′=p。即将原来模型的自变量与因变量都变为对应的标准分,参数及映射表达式不变。要注意此时拟合所得到的参数向量 p ′ ⃗ \vec{p'} p′ 不是原来模型所需要的参数向量 p ⃗ \vec p p,需要在拟合出 p ′ ⃗ \vec{p'} p′ 后再求 p ⃗ \vec p p。
根据实际情况可以不对所有变量都使用标准分;根据模型的表达式,也不是任何变量都可以使用标准分。例如对于模型 y = exp ( a x 1 ) + b x 2 + c y = \exp (a x_1) + b x_2 + c y=exp(ax1)+bx2+c,如果使用 y y y 的标准分,那么拟合所得的残差平方和将与拟合原来模型所得的不同,因此可以考虑不使用 y y y 的标准分。若使用 x 1 x_1 x1 的标准分,则模型变为 y = exp ( a ′ ( x 1 − x 1 ‾ ) σ x 1 ) + b x 2 + c = exp ( a ′ σ x 1 x 1 − a ′ x 1 ‾ σ x 1 ) + b x 2 + c y = \exp \left(\dfrac{a'(x_1 - \overline{x_1})}{\sigma_{x_1}}\right) + b x_2 + c = \exp \left(\dfrac{a'}{\sigma_{x_1}}x_1 - \dfrac{a'\overline{x_1}}{\sigma_{x_1}}\right) + b x_2 + c y=exp(σx1a′(x1−x1))+bx2+c=exp(σx1a′x1−σx1a′x1)+bx2+c,对于任意的 a ′ ≠ 0 a'\neq 0 a′=0,不存在 a a a 能够使两个模型等价,因此不能使用 x 1 x_1 x1 的标准分。可以使用 x 2 x_2 x2 的标准分,则模型变为 y = exp ( a x 1 ) + b ′ x 2 − x 2 ‾ σ x 2 + c ′ = exp ( a x 1 ) + b ′ σ x 2 x 2 − b ′ x 2 ‾ σ x 2 + c ′ y = \exp(a x_1) + b' \dfrac{x_2 - \overline{x_2}}{\sigma_{x_2}} + c' = \exp(a x_1) + \dfrac{b'}{\sigma_{x_2}}x_2 - \dfrac{b' \overline{x_2}}{\sigma_{x_2}} + c' y=exp(ax1)+b′σx2x2−x2+c′=exp(ax1)+σx2b′x2−σx2b′x2+c′,从而有 b = b ′ σ x 2 , c = c ′ − b ′ x 2 ‾ σ x 2 b = \dfrac{b'}{\sigma_{x_2}}, c = c' - \dfrac{b' \overline{x_2}}{\sigma_{x_2}} b=σx2b′,c=c′−σx2b′x2,则拟合出 b ′ , c ′ b', c' b′,c′ 后即可代入上述等式求出 b , c b, c b,c。
上述标准分其实是统计学中的一种标准化方法。由于标准分对 x x x 进行平移,而模型中的参数 a a a 只能对 x x x 进行缩放,因此导致对 x 1 x_1 x1 使用标准分后出现模型不等价的问题。如果使用一种变体
z x 1 ′ = x σ x z'_{x_1} = \dfrac{x}{\sigma_x} zx1′=σxx
则模型变为 y = exp ( a ′ x 1 σ x 1 ) + b x 2 + c y = \exp \left(\dfrac{a' x_1}{\sigma_{x_1}}\right) + b x_2 + c y=exp(σx1a′x1)+bx2+c,则 a = a ′ σ x 1 a = \dfrac{a'}{\sigma_{x_1}} a=σx1a′,拟合出 a ′ a' a′ 后就能求出 a a a。或者在标准分模型上修改:
y = exp ( a ′ σ x 1 x 1 − a ′ x 1 ‾ σ x 1 + d ) + b x 2 + c y = \exp \left(\dfrac{a'}{\sigma_{x_1}}x_1 - \dfrac{a'\overline{x_1}}{\sigma_{x_1}} + d\right) + b x_2 + c y=exp(σx1a′x1−σx1a′x1+d)+bx2+c
在指数内添加一个 d d d。为使两个模型等价,应令 d − a ′ x 1 ‾ σ x 1 = 0 d - \dfrac{a'\overline{x_1}}{\sigma_{x_1}} = 0 d−σx1a′x1=0,则模型变为上述所示形式。
在例题中,选择对 x 1 x_1 x1 至 x 4 x_4 x4 使用标准分,标准分模型为
y = p 1 ′ + p 2 ′ z x 1 + p 3 ′ z x 2 + p 4 ′ z x 3 + p 5 ′ z x 4 1 + p 6 ′ z x 1 + p 7 ′ z x 2 + p 8 ′ z x 3 + p 9 ′ z x 4 y = \frac{p_1' + p_2' z_{x_1} + p_3' z_{x_2} + p_4' z_{x_3} + p_5' z_{x_4}}{1 + p_6' z_{x_1} + p_7' z_{x_2} + p_8' z_{x_3} + p_9' z_{x_4}} y=1+p6′zx1+p7′zx2+p8′zx3+p9′zx4p1′+p2′zx1+p3′zx2+p4′zx3+p5′zx4
p 1 ′ p_1' p1′ 至 p 5 ′ p_5' p5′ 可线性拟合,将分母整理,使其与原模型的分母形式相似
( 1 − p 6 ′ σ x 1 x 1 ‾ − p 7 ′ σ x 2 x 2 ‾ − p 8 ′ σ x 3 x 3 ‾ − p 9 ′ σ x 4 x 4 ‾ ) + p 6 ′ σ x 1 x 1 + p 7 ′ σ x 2 x 2 + p 8 ′ σ x 3 x 3 + p 9 ′ σ x 4 x 4 \left(1 - \dfrac{p_6'}{\sigma_{x_1}}\overline{x_1} - \dfrac{p_7'}{\sigma_{x_2}}\overline{x_2} - \dfrac{p_8'}{\sigma_{x_3}}\overline{x_3} - \dfrac{p_9'}{\sigma_{x_4}}\overline{x_4}\right) + \dfrac{p_6'}{\sigma_{x_1}}x_1 + \dfrac{p_7'}{\sigma_{x_2}}x_2 + \dfrac{p_8'}{\sigma_{x_3}}x_3 + \dfrac{p_9'}{\sigma_{x_4}}x_4 (1−σx1p6′x1−σx2p7′x2−σx3p8′x3−σx4p9′x4)+σx1p6′x1+σx2p7′x2+σx3p8′x3+σx4p9′x4
为了使该分母相对于 x 1 x_1 x1 至 x 4 x_4 x4 的常数项也化为 1 1 1,记 e 1 = 1 − p 6 ′ σ x 1 x 1 ‾ − p 7 ′ σ x 2 x 2 ‾ − p 8 ′ σ x 3 x 3 ‾ − p 9 ′ σ x 4 x 4 ‾ e_1 = 1 - \dfrac{p_6'}{\sigma_{x_1}}\overline{x_1} - \dfrac{p_7'}{\sigma_{x_2}}\overline{x_2} - \dfrac{p_8'}{\sigma_{x_3}}\overline{x_3} - \dfrac{p_9'}{\sigma_{x_4}}\overline{x_4} e1=1−σx1p6′x1−σx2p7′x2−σx3p8′x3−σx4p9′x4,标准分模型的分子分母同除 e 1 e_1 e1,再比较两模型分母得
p i = p i ′ e 1 σ x i − 5 , i = 6 , 7 , 8 , 9 p_i = \dfrac{p_i'}{e_1 \sigma_{x_i-5}}, i = 6, 7, 8, 9 pi=e1σxi−5pi′,i=6,7,8,9
可以使用向量符号来简化记号。设标准分模型参数向量 s ⃗ = ( p 6 ′ , ⋯ , p 9 ′ ) \vec s = (p_6', \cdots ,p_9') s=(p6′,⋯,p9′),平均值向量 x ‾ ⃗ = ( x 1 ‾ , ⋯ , x 4 ‾ ) \vec{\overline x} = (\overline{x_1}, \cdots, \overline{x_4}) x=(x1,⋯,x4),标准差向量 σ ⃗ = ( σ x 1 , ⋯ , σ x 4 ) \vec\sigma = (\sigma_{x_1}, \cdots, \sigma_{x_4}) σ=(σx1,⋯,σx4),记 e 2 ⃗ = s ⃗ σ ⃗ \vec{e_2} = \dfrac{\vec s}{\vec\sigma} e2=σs,则 e 1 = 1 − e 2 ⃗ ⋅ x ‾ ⃗ e_1 = 1 - \vec{e_2}\cdot\vec{\overline x} e1=1−e2⋅x,设原模型参数向量 p ⃗ = ( p 6 , ⋯ , p 9 ) \vec p = (p_6, \cdots, p_9) p=(p6,⋯,p9),则 p ⃗ = e 2 ⃗ e 1 = e 2 ⃗ 1 − e 2 ⃗ ⋅ x ‾ ⃗ \vec p = \dfrac{\vec{e_2}}{e_1} = \dfrac{\vec{e_2}}{1 - \vec{e_2}\cdot\vec{\overline x}} p=e1e2=1−e2⋅xe2。