非线性拟合技巧：分离线性参数与数据标准化

非线性拟合技巧：分离线性参数与数据标准化

例题

对非线性模型

$$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}$$

的一次观测数据的 csv 格式文本如下所示，各列分别为 $x_1,x_2,x_3,x_4,y$，以最小化残差平方和为目标拟合参数 $p_1$ 至 $p_9$。一个较好结果的残差平方和为 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(\vec{x};\vec{p})$$

其中 $\vec{x}\in {\mathbb{R}}^n$ 是自变量，$y\in \mathbb{R}$ 是因变量，$f:{\mathbb{R}}^n\to \mathbb{R}$ 是一映射，$\vec{p}\in {\mathbb{R}}^m$ 是 $f$ 的参数向量。若 $y$ 对 $\vec{p}$ 的某分量 $p_i$ 可视作一个线性函数，则称 $p_i$ 是线性参数。若固定 $\vec{p}$ 的所有非线性参数，则 $\vec{p}$ 余下的参数均为线性参数，这些参数可以使用线性拟合的方法求解。根据线性代数的知识，当解存在时，线性拟合方法总是求得最优解。

在该例题中，注意到 $p_1$ 至 $p_5$ 是线性参数，则只需拟合 $p_6$ 至 $p_9$，而 $p_1$ 至 $p_5$ 则交由 lm 方法拟合即可。

数据标准化

在该例题中，观察到各变量之间的数量级相差较大。使用 $x_1$ 至 $x_4$ 的标准分进行拟合，具体代码如下所示。

# 代码不能直接运行，详情请参阅各处注释
# 在此处编写代码将数据读入至 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_9$ 的一个较好的结果是 -2.235640774844147663555e-06, 7.302384225577161910251e-08, -2.128725631770886621980e-03, 3.161796617995193568393e-04。可以看到它们的数量级都很小。事实上对于大部分的局部优化算法来说，起始点都必须选择在这个解的附近才能得到该结果。以下代码将起始点的每个分量加上一个均匀分布于 $\left[-\frac{3\times 10^{-7}}{4},\frac{3\times 10^{-7}}{4}\right]$ 的随机数，成功求解的概率已经低于 95%，而随机数的区间若选择 $[-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$ 的标准分定义为

$$z_x=\frac{x-\overline{x}}{{\sigma }_x}$$

其中 $\overline{x}$ 是该组数据的平均值，${\sigma }_x$ 是该组数据的标准差。

使用标准分后，模型也相应地发生变化。例如使用自变量与因变量的标准分，则模型变为

$$z_y=f(\vec{z};\overrightarrow{p^{\prime }})$$

其中 $\vec{z}$ 的各分量是 $\vec{x}$ 各分量的标准分，$z_y$ 是 $y$ 的标准分，拟合的参数向量 $\overrightarrow{p^{\prime }}$ 维度与 $\vec{p}$ 相同，但 $\overrightarrow{p^{\prime }}\ne \vec{p}$。即将原来模型的自变量与因变量都变为对应的标准分，参数及映射表达式不变。要注意此时拟合所得到的参数向量 $\overrightarrow{p^{\prime }}$ 不是原来模型所需要的参数向量 $\vec{p}$，需要在拟合出 $\overrightarrow{p^{\prime }}$ 后再求 $\vec{p}$。

根据实际情况可以不对所有变量都使用标准分；根据模型的表达式，也不是任何变量都可以使用标准分。例如对于模型 $y=\exp (a{x}_1)+b{x}_2+c$，如果使用 $y$ 的标准分，那么拟合所得的残差平方和将与拟合原来模型所得的不同，因此可以考虑不使用 $y$ 的标准分。若使用 $x_1$ 的标准分，则模型变为 $y=\exp \left(\frac{a^{\prime }(x_1-\overline{x_1})}{{\sigma }_{x_1}}\right)+b{x}_2+c=\exp \left(\frac{a^{\prime }}{{\sigma }_{x_1}}{x}_1-\frac{a^{\prime }\overline{x_1}}{{\sigma }_{x_1}}\right)+b{x}_2+c$，对于任意的 $a^{\prime }\ne 0$，不存在 $a$ 能够使两个模型等价，因此不能使用 $x_1$ 的标准分。可以使用 $x_2$ 的标准分，则模型变为 $y=\exp (a{x}_1)+b^{\prime }\frac{x_2-\overline{x_2}}{{\sigma }_{x_2}}+c^{\prime }=\exp (a{x}_1)+\frac{b^{\prime }}{{\sigma }_{x_2}}{x}_2-\frac{b^{\prime }\overline{x_2}}{{\sigma }_{x_2}}+c^{\prime }$，从而有 $b=\frac{b^{\prime }}{{\sigma }_{x_2}},c=c^{\prime }-\frac{b^{\prime }\overline{x_2}}{{\sigma }_{x_2}}$，则拟合出 $b^{\prime },c^{\prime }$ 后即可代入上述等式求出 $b,c$。

上述标准分其实是统计学中的一种标准化方法。由于标准分对 $x$ 进行平移，而模型中的参数 $a$ 只能对 $x$ 进行缩放，因此导致对 $x_1$ 使用标准分后出现模型不等价的问题。如果使用一种变体

$$z_{x_1}^{\prime }=\frac{x}{{\sigma }_x}$$

则模型变为 $y=\exp \left(\frac{a^{\prime }{x}_1}{{\sigma }_{x_1}}\right)+b{x}_2+c$，则 $a=\frac{a^{\prime }}{{\sigma }_{x_1}}$，拟合出 $a^{\prime }$ 后就能求出 $a$。或者在标准分模型上修改：

$$y=\exp \left(\frac{a^{\prime }}{{\sigma }_{x_1}}{x}_1-\frac{a^{\prime }\overline{x_1}}{{\sigma }_{x_1}}+d\right)+b{x}_2+c$$

在指数内添加一个 $d$。为使两个模型等价，应令 $d-\frac{a^{\prime }\overline{x_1}}{{\sigma }_{x_1}}=0$，则模型变为上述所示形式。

在例题中，选择对 $x_1$ 至 $x_4$ 使用标准分，标准分模型为

$$y=\frac{p_1^{\prime }+p_2^{\prime }{z}_{x_1}+p_3^{\prime }{z}_{x_2}+p_4^{\prime }{z}_{x_3}+p_5^{\prime }{z}_{x_4}}{1+p_6^{\prime }{z}_{x_1}+p_7^{\prime }{z}_{x_2}+p_8^{\prime }{z}_{x_3}+p_9^{\prime }{z}_{x_4}}$$

$p_1^{\prime }$ 至 $p_5^{\prime }$ 可线性拟合，将分母整理，使其与原模型的分母形式相似

$$\left(1-\frac{p_6^{\prime }}{{\sigma }_{x_1}}\overline{x_1}-\frac{p_7^{\prime }}{{\sigma }_{x_2}}\overline{x_2}-\frac{p_8^{\prime }}{{\sigma }_{x_3}}\overline{x_3}-\frac{p_9^{\prime }}{{\sigma }_{x_4}}\overline{x_4}\right)+\frac{p_6^{\prime }}{{\sigma }_{x_1}}{x}_1+\frac{p_7^{\prime }}{{\sigma }_{x_2}}{x}_2+\frac{p_8^{\prime }}{{\sigma }_{x_3}}{x}_3+\frac{p_9^{\prime }}{{\sigma }_{x_4}}{x}_4$$

为了使该分母相对于 $x_1$ 至 $x_4$ 的常数项也化为 $1$，记 $e_1=1-\frac{p_6^{\prime }}{{\sigma }_{x_1}}\overline{x_1}-\frac{p_7^{\prime }}{{\sigma }_{x_2}}\overline{x_2}-\frac{p_8^{\prime }}{{\sigma }_{x_3}}\overline{x_3}-\frac{p_9^{\prime }}{{\sigma }_{x_4}}\overline{x_4}$，标准分模型的分子分母同除 $e_1$，再比较两模型分母得

$$p_i=\frac{p_i^{\prime }}{e_1{\sigma }_{x_i-5}},i=6,7,8,9$$

可以使用向量符号来简化记号。设标准分模型参数向量 $\vec{s}=(p_6^{\prime },\cdots ,p_9^{\prime })$，平均值向量 $\vec{\overline{x}}=(\overline{x_1},\cdots ,\overline{x_4})$，标准差向量 $\vec{\sigma }=({\sigma }_{x_1},\cdots ,{\sigma }_{x_4})$，记 $\overrightarrow{e_2}=\frac{\vec{s}}{\vec{\sigma }}$，则 $e_1=1-\overrightarrow{e_2}\cdot \vec{\overline{x}}$，设原模型参数向量 $\vec{p}=(p_6,\cdots ,p_9)$，则 $\vec{p}=\frac{\overrightarrow{e_2}}{e_1}=\frac{\overrightarrow{e_2}}{1-\overrightarrow{e_2}\cdot \vec{\overline{x}}}$。