将多项式拟合到数据

在给定一组值的情况下(x,f(x)),有没有办法找到最适合数据的给定度数的多项式？

我知道多项式插值,它用于找到n给定n+1数据点的度数多项式,但这里有大量的值,我们想找到一个低次多项式(找到最佳线性拟合,最佳二次,最佳立方等. ).它可能与最小二乘有关 ...

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

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

1> ShreevatsaR..：

---

感谢大家的回复.这是另一种总结它们的尝试.请原谅,如果我说太多"明显"的事情:我之前对最小二乘一无所知,所以一切对我来说都是新的.

非多项式插值

多项式插值拟合n给定n+1数据点的多项式,例如找到精确通过四个给定点的立方.正如在问题中所说,这不是我想要的 - 我有很多分数并且想要一个小程度多项式(除非我们很幸运,它们只能大致适合) - 但是因为一些答案坚持谈论关于它,我应该提到它们:拉格朗日多项式,Vandermonde矩阵等.

什么是最小二乘？

"最小二乘法"是多项式拟合"有多好"的特定定义/标准/"度量".(还有其他的,但这是最简单的.)假设你试图将多项式p(x,y)= a + bx + cy + dx ² + ey ² + fxy拟合到某些给定的数据点(x _i,y _i),Z _i)(其中"Z _i "在问题中是"f(x _i,y _i)").对于最小二乘问题,问题是找到"最佳"系数(a,b,c,d,e,f),使得最小化(保持"最小")的是"残差平方和",即

S =Σ _我(A + BX _我 + CY _我 + DX _我² + EY _我² + FX _我 ÿ _我 - z _我)²

理论

重要的想法是,如果将S视为(a,b,c,d,e,f)的函数,则S 在其梯度为0的点处被最小化.这意味着例如∂S/∂f= 0,即

Σ _我图2(a + ... + FX _我 ÿ _我 - z _我)X _我 Ŷ _我 = 0

和a,b,c,d,e的类似方程式.请注意,这些只是... f中的线性方程式.所以我们可以用高斯消元法或任何常用方法来解决它们.

这仍称为"线性最小二乘法",因为虽然我们想要的函数是二次多项式,但它在参数(a,b,c,d,e,f)中仍然是线性的.注意,当我们希望p(x,y)是任意函数f _j的任何"线性组合"时,同样的事情是有效的,而不仅仅是多项式(="单项式的线性组合").

码

对于单变量情况(当只有变量x - f _j是单项式x ^j时),有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.如其文档中所解释的,它采用值f _j(x _i)的矩阵A. (理论是它找到了A 的Moore-Penrose伪逆.)在上面的例子中涉及(x _i,y _i,Z _i),拟合多项式意味着f _j是单项式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

所以它发现多项式是x ² + 2xy + y ² +0.01.[最后一个术语有时是-0.01,有时是0,由于我们添加的随机噪声,这是预期的.

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

关注

它受如何选择点的影响很大.当我有x=y=range(20),而不是随机点,它总是产生1.33X ² + 1.33xy + 1.33y ²,这是令人费解...直到我意识到,因为我总是有x[i]=y[i],多项式都是一样的:X ² + 2XY + Y ² = 4x ² =(4/3)(x ² + xy + y ²).因此,道德是仔细选择要点以获得"正确的"多项式是很重要的.(如果你可以选择,你应该选择Chebyshev节点进行多项式插值;不确定最小二乘方是否也是如此.)

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

可能存在我不完全理解的数值稳定性问题.

如果不需要多项式,则可以更好地拟合其他类型的函数,例如样条(分段多项式).

---

关于你对数值稳定性的关注:定义一个多项式(="单项式的线性组合")是一件危险的事情,因为(用非数学的话来说)4级以上的单项式(比如说)是非常相似的.在0左右的区域然后他们只是"变得疯狂".更好的方法是决定你试图在哪个区间拟合数据,重新定义你的自变量,使你实际适合(-1,1),并寻找良好多项式的线性组合而不是单项式.我会用Chebyshev套装.

2> John D. Cook..：

---

是的,通常这样做的方法是使用最小二乘法.还有其他方法可以指定多项式的拟合程度,但对于最小二乘法,理论最简单.一般理论称为线性回归.

你最好的选择可能是从Numerical Recipes开始.

R是免费的,可以做你想要的一切,但它有一个很大的学习曲线.

如果您可以访问Mathematica,则可以使用"拟合"功能进行最小二乘拟合.我想Matlab及其开源对应物Octave具有类似的功能.

3> jfs..：

---

对于(x,f(x))案例:

import numpy

x = numpy.arange(10)
y = x**2

coeffs = numpy.polyfit(x, y, deg=2)
poly = numpy.poly1d(coeffs)
print poly
yp = numpy.polyval(poly, x)
print (yp-y)