math - 将多项式拟合到数据

Question

有没有办法，给定一组值 (x,f(x)) , 找到最适合数据的给定次数的多项式？

我知道 polynomial interpolation ，这是为了找到一个多项式n给定 n+1数据点，但这里有大量值，我们想找到一个低次多项式(找到最佳线性拟合、最佳二次、最佳三次等)。可能与least squares有关...

更一般地说，当我们有一个多元函数时，我想知道答案——像 (x,y,f(x,y)) 这样的点。 ，比如说 - 并且想要在变量中找到给定次数的最佳多项式( p(x,y) )。 (特别是多项式，而不是样条或傅立叶级数。)

理论和代码/库(最好是 Python，但任何语言都可以)都会很有用。

Answer 1

谢谢大家的回复。这是总结它们的另一种尝试。请原谅我说了太多“显而易见”的事情:我以前对最小二乘一无所知，所以一切对我来说都是新的。

非多项式插值

Polynomial interpolation拟合多项式 n给定 n+1数据点，例如找到一个恰好通过四个给定点的三次方。正如问题中所说，这不是我想要的——我有很多点并且想要一个小次多项式(这只会近似拟合，除非我们很幸运)——但是因为有些答案坚持要讨论关于它，我应该提到它们:) Lagrange polynomial , Vandermonde matrix ， 等等。

什么是最小二乘法？

“最小二乘法”是多项式拟合“有多好”的特定定义/标准/“度量”。 (还有其他的，但这是最简单的。)假设您正在尝试拟合多项式
p(x,y) = a + bx + cy + dx2 + ey2 + fxy
到一些给定的数据点 (xi,yi,Zi) (其中“Zi”在问题中是“f(xi,yi)”)。使用最小二乘法的问题是找到“最佳”系数(a、b、c、d、e、f)，使得最小化(保持“最小”)的是“残差平方和”，即

S = ∑i (a + bxi + cyi + dxi2 + eyi2 + fxiyi - Zi)2

理论

重要的想法是，如果您将 S 视为 (a,b,c,d,e,f) 的函数，则 S 是 minimized在其 gradient is 0 .这意味着例如∂S/∂f=0，即

∑i2(a + … + fxiyi - Zi)xiyi = 0

a、b、c、d、e 的类似方程。
请注意，这些只是 a...f 中的线性方程。所以我们可以用 Gaussian elimination 解决它们或 the usual methods 中的任何一个.

这仍然称为“线性最小二乘法”，因为虽然我们想要的函数是二次多项式，但它在参数 (a,b,c,d,e,f) 中仍然是线性的。请注意，当我们希望 p(x,y) 是任意函数 fj 的任何“线性组合”，而不仅仅是多项式(=“单项式的线性组合”)时，同样的事情也会起作用。

代码

对于单变量情况(当只有变量 x - fj 是单项式 xj)，有 Numpy 的 polyfit :

>>> import numpy
>>> xs = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
>>> ys = [1.1, 3.9, 11.2, 21.5, 34.8, 51, 70.2, 92.3, 117.4, 145.5]
>>> p = numpy.poly1d(numpy.polyfit(xs, ys, deg=2))
>>> print p
       2
1.517 x + 2.483 x + 0.4927

对于多元情况，或一般的线性最小二乘法，有 SciPy。 As explained in its documentation ，它采用值 fj(xi) 的矩阵 A。 (理论是它找到了 A 的 Moore-Penrose pseudoinverse。)在我们上面涉及 (xi,yi,Zi) 的例子中，拟合多项式意味着 fj 是单项式 x()y()。以下找到最佳二次(或任何其他次数的最佳多项式，如果您更改“degree = 2”行):

from scipy import linalg
import random

n = 20
x = [100*random.random() for i in range(n)]
y = [100*random.random() for i in range(n)]
Z = [(x[i]+y[i])**2 + 0.01*random.random() for i in range(n)]

degree = 2
A = []
for i in range(n):
    A.append([])
    for xd in range(degree+1):
        for yd in range(degree+1-xd):
            A[i].append((x[i]**xd)*(y[i]**yd)) #f_j(x_i)

c,_,_,_ = linalg.lstsq(A,Z)
j = 0
for xd in range(0,degree+1):
    for yd in range(0,degree+1-xd):
        print " + (%.2f)x^%dy^%d" % (c[j], xd, yd),
        j += 1

打印

 + (0.01)x^0y^0  + (-0.00)x^0y^1  + (1.00)x^0y^2  + (-0.00)x^1y^0  + (2.00)x^1y^1  + (1.00)x^2y^0

所以它发现多项式是x2+2xy+y2+0.01。 [最后一项有时是 -0.01，有时是 0，这是可以预料的，因为我们添加了随机噪声。]

Python+Numpy/Scipy 的替代方案是 R和计算机代数系统: Sage , Mathematica, Matlab, Maple。甚至 Excel 也能做到。 Numerical Recipes讨论自己实现它的方法(在 C、Fortran 中)。

顾虑

它受到的强烈影响如何选择点 .当我有 x=y=range(20)而不是随机点，它总是产生 1.33x2+1.33xy+1.33y2，这令人费解......直到我意识到，因为我总是有 x[i]=y[i] ，多项式相同:x2+2xy+y2 = 4x2 = (4/3)(x2+xy+y2)。所以道德是仔细选择点以获得“正确”的多项式很重要。 (如果可以选择，您应该为多项式插值选择 Chebyshev nodes；不确定最小二乘法是否也是如此。)

过拟合 :更高次的多项式总是可以更好地拟合数据。如果您更改 degree到 3 或 4 或 5，它仍然主要识别相同的二次多项式(更高阶项的系数为 0)，但对于更大的阶数，它开始拟合更高阶多项式。但即使是 6 次，采用更大的 n(更多的数据点而不是 20，比如 200)仍然适合二次多项式。所以道德是避免过度拟合，为此它可能有助于获取尽可能多的数据点。

可能存在 numerical stability 的问题我不完全明白。

如果您不需要多项式，您可以获得与其他类型函数的更好拟合，例如splines (分段多项式)。