使用 GNU Scientific C 的多项式拟合
Polynomial Fitting using GNU Scientific C
我需要从 n+1 个数据点得到一个 n-degree 函数。通过使用以下 gnuplot 脚本，我得到了正确的拟合：
f(x) = a + b*x + c*x**2 + d*x**3 + e*x**4 + f*x**5 + g*x**6 + h*x**7 + i*x**8 + j*x**9 + k*x**10 + l*x**11 + m*x**12
# Initial values for parameters
a = 0.1
b = 0.1
c = 0.1
d = 0.1
e = 0.1
f = 0.1
g = 0.1
h = 0.1
i = 0.1
j = 0.1
k = 0.1
l = 0.1
m = 0.1
# Fit f to the following data by modifying the variables a, b, c
fit f(x) '-' via a, b, c, d, e, f, g, h, i, j, k, l, m
4.877263 45.036000
4.794907 44.421000
4.703827 43.808000
4.618065 43.251000
4.530520 42.634000
4.443111 42.070000
4.357077 41.485000
4.274298 40.913000
4.188404 40.335000
4.109381 39.795000
4.027594 39.201000
3.946413 38.650000
3.874360 38.085000
e
拟合后得到如下系数：
a = -781956
b = -2.52463e+06
c = 2.75682e+06
d = -553791
e = 693880
f = -1.51285e+06
g = 1.21157e+06
h = -522243
i = 138121
j = -23268.8
k = 2450.79
l = -147.834
m = 3.91268
然后，通过将数据和 f(x) 一起绘制，似乎给定的系数是正确的：
但是，我需要使用c 代码来获得这样的拟合。在某些情况下，多项式拟合 (as in this link) 的 GNU 科学图书馆代码的结果是正确的。但是对于上述数据（以及我数据集中的其他几个案例），我得到的结果是有缺陷的。
例如，以下代码（使用与上述示例相同的数据）：
void testOfPolynomialFit(){
double x[13] = {4.877263, 4.794907, 4.703827, 4.618065, 4.530520, 4.443111, 4.357077, 4.274298, 4.188404, 4.109381, 4.027594, 3.946413, 3.874360};
double y[13] = {45.036000, 44.421000, 43.808000, 43.251000, 42.634000, 42.070000, 41.485000, 40.913000, 40.335000, 39.795000, 39.201000, 38.650000, 38.085000};
double coefficients[13];
polynomialfit(13, 13, x, y, coefficients);
int i, n = 13;
for (i = 0; i < n; i++)
{
printf("%lf\t", coefficients[i]);
}
printf("\n");
}
结果：
-6817581083.803348 12796304366.105989 -9942834843.404181 3892080279.353104
-630964566.517794 -75914607.005088 49505072.518952 -5062100.000931
-1426228.491628 514259.312320 -70903.844354 4852.824607
-136.738756
对应于以下形式的函数：
c(x)=-6837615134.799868+12834646330.586414*x**1-9973474377.668280*x**2+3904659818.834625*x**3-633282611.288889*x**4-76066283.747942*x**5+49670960.939126*x**6-5091123.449217*x**7-1426628.818192*x**8+515175.778491*x**9-71055.177018*x**10+4863.969973*x**11-137.065848*x**12
可以在此处查看 c(x) 的样子：
在这样的图像中，a(x) 和 b(x) 是使用 "polynomialfit" 拟合的函数，仅针对少数点（4 和 7）。
那么，关于我在这里做错了什么，有什么提示吗？还有一些其他的 c 代码可以提供合适的配件吗？
您的数值不稳定。简单的线性回归证实了可以直观观察到的内容：99.98% 的变化可以仅使用线性模型来解释。
您提供的 link 中的代码做了很多 非常 不安全的事情：不检查 obs 或 degree 是肯定的，不检查内存分配是否成功，不 returning 任何有用的东西，...
我假设 gsl_multifit_linear 溢出，或者不能包含数值不稳定性，不检查 return 意味着我们不知道。
编辑：
根据 GSL Website 多项式回归会因计算大数的大幂而导致额外的数值不稳定性。尝试使用 x = (x - avg_x) / sd_x 预处理您的 x 值。这应该允许您在发生这种情况之前获得更多的多项式度数。
在漫长的运行中，你很可能会再次遇到这个问题。如果您使用 35 或 100 或更多数据点执行此分析，则您不太可能找到任何技术来克服不稳定性。
这两种解决方法之间的一个主要区别是，当您使用 gnuplot 时，您正在做 fit 来设置系数并在同一程序中绘制函数，而使用 GSL 时，您正在将数字从一个程序复制到另一个程序。
如果您使用 printf("%lf", ...) 的输出作为您的第二个 gnuplot 程序的输入，您会失去很多准确性，因为 printf 比任何任一程序的内部操作。因为这是一个数值不稳定的问题，一点点四舍五入会造成很大的伤害。
当 x 为 4.877263 时，x**12 大约为 181181603.850932 因此，如果您的 m 偏离 0.000001（printf 的默认舍入级别），则会引入错误 181.181603850932这比 x.
的实际 y 值大约 300% 的相对误差
尝试%.60lf，看看它是否会好转。
如果其中一个程序在内部使用 long double 而另一个程序没有，那么无论您做什么，您都可能无法获得良好的匹配。
我对问题陈述和示例代码有点困惑。 polynomialfit 函数通常期望 >= n + 2 个数据点来拟合 n 次多项式。在 n+1 个数据点的情况下，不是进行拟合，而是通过生成具有 n+1 行和列的矩阵以及给定的 n+1 集来生成精确解（如果浮点舍入不是问题）表示线性方程组的值：
| x[0]^n + x[0]^(n-1) + ... + x[0] + 1 | | c[ n] | | y[0] |
| x[1]^n + x[1]^(n-1) + ... + x[1] + 1 | | c[n-1] | | y[1] |
| ... = | ... |
| x[n]^n + x[n]^(n-1) + ... + x[n] + 1 | | c[ 0] | | y[n] |
所以只有 c[ ] 是变量，方程是线性的。反转 x 值的矩阵，然后将 y 值乘以反转矩阵以产生结果。如果实际多项式的次数低于 n（两个或多个方程将不是线性独立的），则可能会出现问题。如果发生这种情况，您可以使用更少的一行或多行，或者改用传统的多项式拟合算法。
如果有 >= n+2 个数据点，一种选择是最小二乘型多项式拟合。这是 .rtf 文档的 link，它使用正交和递归定义的多项式来拟合一组数据点。正交多项式消除了对矩阵求逆的需要，因此它可以更准确地处理更大次数的多项式。
示例以 degree = 1 和 degree = 3 运行。第一列是原始 x 值，第二列是原始 y 值，第三列是计算出的 y 值，第四列是（原始 y 值 - 计算出的 y 值）：
variance = 9.6720e-004
6.8488e+000 X**1 + 1.1619e+001
4.877263e+000 4.503600e+001 4.502243e+001 1.356604e-002
4.794907e+000 4.442100e+001 4.445839e+001 -3.739271e-002
4.703827e+000 4.380800e+001 4.383460e+001 -2.660237e-002
4.618065e+000 4.325100e+001 4.324723e+001 3.765950e-003
4.530520e+000 4.263400e+001 4.264765e+001 -1.365429e-002
4.443111e+000 4.207000e+001 4.204901e+001 2.099404e-002
4.357077e+000 4.148500e+001 4.145977e+001 2.522524e-002
4.274298e+000 4.091300e+001 4.089284e+001 2.016354e-002
4.188404e+000 4.033500e+001 4.030456e+001 3.043591e-002
4.109381e+000 3.979500e+001 3.976335e+001 3.165004e-002
4.027594e+000 3.920100e+001 3.920321e+001 -2.205684e-003
3.946413e+000 3.865000e+001 3.864721e+001 2.788204e-003
3.874360e+000 3.808500e+001 3.815373e+001 -6.873392e-002
variance = 2.4281e-004
8.0952e-001 X**3 + -1.0822e+001 X**2 + 5.4910e+001 X**1 + -5.9287e+001
4.877263e+000 4.503600e+001 4.502276e+001 1.324045e-002
4.794907e+000 4.442100e+001 4.444280e+001 -2.180419e-002
4.703827e+000 4.380800e+001 4.381431e+001 -6.306292e-003
4.618065e+000 4.325100e+001 4.323170e+001 1.929905e-002
4.530520e+000 4.263400e+001 4.264294e+001 -8.935141e-003
4.443111e+000 4.207000e+001 4.205786e+001 1.214442e-002
4.357077e+000 4.148500e+001 4.148153e+001 3.468503e-003
4.274298e+000 4.091300e+001 4.092369e+001 -1.069376e-002
4.188404e+000 4.033500e+001 4.033844e+001 -3.436876e-003
4.109381e+000 3.979500e+001 3.979160e+001 3.397859e-003
4.027594e+000 3.920100e+001 3.921454e+001 -1.354191e-002
3.946413e+000 3.865000e+001 3.862800e+001 2.199866e-002
3.874360e+000 3.808500e+001 3.809383e+001 -8.830768e-003
我需要从 n+1 个数据点得到一个 n-degree 函数。通过使用以下 gnuplot 脚本，我得到了正确的拟合：
f(x) = a + b*x + c*x**2 + d*x**3 + e*x**4 + f*x**5 + g*x**6 + h*x**7 + i*x**8 + j*x**9 + k*x**10 + l*x**11 + m*x**12
# Initial values for parameters
a = 0.1
b = 0.1
c = 0.1
d = 0.1
e = 0.1
f = 0.1
g = 0.1
h = 0.1
i = 0.1
j = 0.1
k = 0.1
l = 0.1
m = 0.1
# Fit f to the following data by modifying the variables a, b, c
fit f(x) '-' via a, b, c, d, e, f, g, h, i, j, k, l, m
4.877263 45.036000
4.794907 44.421000
4.703827 43.808000
4.618065 43.251000
4.530520 42.634000
4.443111 42.070000
4.357077 41.485000
4.274298 40.913000
4.188404 40.335000
4.109381 39.795000
4.027594 39.201000
3.946413 38.650000
3.874360 38.085000
e
拟合后得到如下系数：
a = -781956
b = -2.52463e+06
c = 2.75682e+06
d = -553791
e = 693880
f = -1.51285e+06
g = 1.21157e+06
h = -522243
i = 138121
j = -23268.8
k = 2450.79
l = -147.834
m = 3.91268
然后，通过将数据和 f(x) 一起绘制，似乎给定的系数是正确的：
但是，我需要使用c 代码来获得这样的拟合。在某些情况下，多项式拟合 (as in this link) 的 GNU 科学图书馆代码的结果是正确的。但是对于上述数据（以及我数据集中的其他几个案例），我得到的结果是有缺陷的。
例如，以下代码（使用与上述示例相同的数据）：
void testOfPolynomialFit(){
double x[13] = {4.877263, 4.794907, 4.703827, 4.618065, 4.530520, 4.443111, 4.357077, 4.274298, 4.188404, 4.109381, 4.027594, 3.946413, 3.874360};
double y[13] = {45.036000, 44.421000, 43.808000, 43.251000, 42.634000, 42.070000, 41.485000, 40.913000, 40.335000, 39.795000, 39.201000, 38.650000, 38.085000};
double coefficients[13];
polynomialfit(13, 13, x, y, coefficients);
int i, n = 13;
for (i = 0; i < n; i++)
{
printf("%lf\t", coefficients[i]);
}
printf("\n");
}
结果：
-6817581083.803348 12796304366.105989 -9942834843.404181 3892080279.353104
-630964566.517794 -75914607.005088 49505072.518952 -5062100.000931
-1426228.491628 514259.312320 -70903.844354 4852.824607
-136.738756
对应于以下形式的函数：
c(x)=-6837615134.799868+12834646330.586414*x**1-9973474377.668280*x**2+3904659818.834625*x**3-633282611.288889*x**4-76066283.747942*x**5+49670960.939126*x**6-5091123.449217*x**7-1426628.818192*x**8+515175.778491*x**9-71055.177018*x**10+4863.969973*x**11-137.065848*x**12
可以在此处查看 c(x) 的样子：
在这样的图像中，a(x) 和 b(x) 是使用 "polynomialfit" 拟合的函数，仅针对少数点（4 和 7）。
那么，关于我在这里做错了什么，有什么提示吗？还有一些其他的 c 代码可以提供合适的配件吗？
您的数值不稳定。简单的线性回归证实了可以直观观察到的内容：99.98% 的变化可以仅使用线性模型来解释。
您提供的 link 中的代码做了很多 非常 不安全的事情：不检查 obs 或 degree 是肯定的，不检查内存分配是否成功，不 returning 任何有用的东西，...
我假设 gsl_multifit_linear 溢出，或者不能包含数值不稳定性，不检查 return 意味着我们不知道。
编辑：
根据 GSL Website 多项式回归会因计算大数的大幂而导致额外的数值不稳定性。尝试使用 x = (x - avg_x) / sd_x 预处理您的 x 值。这应该允许您在发生这种情况之前获得更多的多项式度数。
在漫长的运行中，你很可能会再次遇到这个问题。如果您使用 35 或 100 或更多数据点执行此分析，则您不太可能找到任何技术来克服不稳定性。
这两种解决方法之间的一个主要区别是，当您使用 gnuplot 时，您正在做 fit 来设置系数并在同一程序中绘制函数，而使用 GSL 时，您正在将数字从一个程序复制到另一个程序。
如果您使用 printf("%lf", ...) 的输出作为您的第二个 gnuplot 程序的输入，您会失去很多准确性，因为 printf 比任何任一程序的内部操作。因为这是一个数值不稳定的问题，一点点四舍五入会造成很大的伤害。
当 x 为 4.877263 时，x**12 大约为 181181603.850932 因此，如果您的 m 偏离 0.000001（printf 的默认舍入级别），则会引入错误 181.181603850932这比 x.
y 值大约 300% 的相对误差
尝试%.60lf，看看它是否会好转。
如果其中一个程序在内部使用 long double 而另一个程序没有，那么无论您做什么，您都可能无法获得良好的匹配。
我对问题陈述和示例代码有点困惑。 polynomialfit 函数通常期望 >= n + 2 个数据点来拟合 n 次多项式。在 n+1 个数据点的情况下，不是进行拟合，而是通过生成具有 n+1 行和列的矩阵以及给定的 n+1 集来生成精确解（如果浮点舍入不是问题）表示线性方程组的值：
| x[0]^n + x[0]^(n-1) + ... + x[0] + 1 | | c[ n] | | y[0] |
| x[1]^n + x[1]^(n-1) + ... + x[1] + 1 | | c[n-1] | | y[1] |
| ... = | ... |
| x[n]^n + x[n]^(n-1) + ... + x[n] + 1 | | c[ 0] | | y[n] |
所以只有 c[ ] 是变量，方程是线性的。反转 x 值的矩阵，然后将 y 值乘以反转矩阵以产生结果。如果实际多项式的次数低于 n（两个或多个方程将不是线性独立的），则可能会出现问题。如果发生这种情况，您可以使用更少的一行或多行，或者改用传统的多项式拟合算法。
如果有 >= n+2 个数据点，一种选择是最小二乘型多项式拟合。这是 .rtf 文档的 link，它使用正交和递归定义的多项式来拟合一组数据点。正交多项式消除了对矩阵求逆的需要，因此它可以更准确地处理更大次数的多项式。
示例以 degree = 1 和 degree = 3 运行。第一列是原始 x 值，第二列是原始 y 值，第三列是计算出的 y 值，第四列是（原始 y 值 - 计算出的 y 值）：
variance = 9.6720e-004
6.8488e+000 X**1 + 1.1619e+001
4.877263e+000 4.503600e+001 4.502243e+001 1.356604e-002
4.794907e+000 4.442100e+001 4.445839e+001 -3.739271e-002
4.703827e+000 4.380800e+001 4.383460e+001 -2.660237e-002
4.618065e+000 4.325100e+001 4.324723e+001 3.765950e-003
4.530520e+000 4.263400e+001 4.264765e+001 -1.365429e-002
4.443111e+000 4.207000e+001 4.204901e+001 2.099404e-002
4.357077e+000 4.148500e+001 4.145977e+001 2.522524e-002
4.274298e+000 4.091300e+001 4.089284e+001 2.016354e-002
4.188404e+000 4.033500e+001 4.030456e+001 3.043591e-002
4.109381e+000 3.979500e+001 3.976335e+001 3.165004e-002
4.027594e+000 3.920100e+001 3.920321e+001 -2.205684e-003
3.946413e+000 3.865000e+001 3.864721e+001 2.788204e-003
3.874360e+000 3.808500e+001 3.815373e+001 -6.873392e-002
variance = 2.4281e-004
8.0952e-001 X**3 + -1.0822e+001 X**2 + 5.4910e+001 X**1 + -5.9287e+001
4.877263e+000 4.503600e+001 4.502276e+001 1.324045e-002
4.794907e+000 4.442100e+001 4.444280e+001 -2.180419e-002
4.703827e+000 4.380800e+001 4.381431e+001 -6.306292e-003
4.618065e+000 4.325100e+001 4.323170e+001 1.929905e-002
4.530520e+000 4.263400e+001 4.264294e+001 -8.935141e-003
4.443111e+000 4.207000e+001 4.205786e+001 1.214442e-002
4.357077e+000 4.148500e+001 4.148153e+001 3.468503e-003
4.274298e+000 4.091300e+001 4.092369e+001 -1.069376e-002
4.188404e+000 4.033500e+001 4.033844e+001 -3.436876e-003
4.109381e+000 3.979500e+001 3.979160e+001 3.397859e-003
4.027594e+000 3.920100e+001 3.921454e+001 -1.354191e-002
3.946413e+000 3.865000e+001 3.862800e+001 2.199866e-002
3.874360e+000 3.808500e+001 3.809383e+001 -8.830768e-003