将多项式拟合到数据

给定一组值(x,f(x))是否有办法find最适合数据的给定度的多项式?

我知道多项式插值法 ,它用于在给定n+1数据点的情况下find度为n的多项式,但这里有大量的值,我们想要find一个低阶多项式(find最佳线性拟合,最好二次,最好立方等)。 这可能与最小二乘有关

更一般地说,当我们有一个多元函数 – 像(x,y,f(x,y))这样的点,并且想要find最好的多项式( p(x,y) )在variables中给定的程度。 (具体来说是一个多项式,而不是样条或傅立叶级数。)

理论和代码/库(最好在Python中,但任何语言都可以)将是有用的。

感谢大家的回复。 这是另一个总结他们的尝试。 如果我说了太多“明显”的东西,请原谅:我以前对最小的方块一无所知,所以对我来说一切都是新的。

NOT多项式插值

多项式插值是在给定n+1数据点的情况下拟合n多项式,例如find一个精确通过四个给定点的立方。 正如在问题中所说,这不是我想要的 – 我有很多的要点,并希望有一个小的多项式(除非我们运气好,否则它只会大致适合) – 但是因为一些答案坚持要说话关于它,我应该提到他们:) 拉格朗日多项式 , Vandermondematrix等

什么是最小二乘?

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

S =Σi(a + bx i + cy i + dx i 2 + ey i 2 + fx i y i – Z i2

理论

重要的思想是,如果将S看作(a,b,c,d,e,f)的函数,则S在其梯度 为0的点处被最小化 。 这意味着例如∂S/∂f= 0,也就是说

Σi 2(a + … + fx i y i – Z i )x i y i = 0

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

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

对于单variables情况(当只有variablesx – 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 

对于多variables情况,或者一般情况下的线性最小二乘,有SciPy。 正如文献中所解释的 ,它需要一个值为f jx i )的matrixA. (理论上它find了A的Moore-Penrose伪逆 。)用上面的例子(x i ,y i ,Z i ),拟合多项式意味着f j是单项式x () y () 。 下面find最好的二次(或任何其他度数最好的多项式,如果你改变“度= 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 2 + 2xy + y 2 +0.01。 [最后一个词有时是-0.01,有时候是0,这是因为我们添加的随机噪声而预期的。]

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

关注

  • 这受到如何select点的强烈影响。 当我有x=y=range(20)而不是随机点时,它总是产生1.33x 2 + 1.33xy + 1.33y 2 ,这是令人费解的…直到我意识到因为我总是有x[i]=y[i] ,多项式是相同的:x 2 + 2xy + y 2 = 4x 2 =(4/3)(x 2 + xy + y 2 )。 所以道德是仔细select点来得到“正确的”多项式是重要的。 (如果你可以select,你应该selectChebyshev节点进行多项式插值;不确定最小二乘是否也是如此。)
  • 过度拟合 :更高阶的多项式可以更好地适应数据。 如果将degree更改为3或4或5,它仍然主要识别相同的二次多项式(高阶项的系数为0),但对于较大的阶数,则开始拟合高阶多项式。 但是,即使是6度,取较大的n(更多的数据点而不是20,也就是200)仍然适合二次多项式。 所以道德是避免过度拟合,这可能有助于尽可能多地获取数据点。
  • 数值稳定性问题可能存在我不完全了解的问题。
  • 如果不需要多项式,则可以更好地拟合其他types的函数,例如样条 (分段多项式)。

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

你最好的select可能是从数字食谱开始。

R是免费的,会尽你所能地做更多的事情,但是它有一个很大的学习曲线。

如果您有权访问Mathematica,则可以使用拟合函数来进行最小二乘拟合。 我想像Matlab和它的开源对应Octave有一个类似的function。

对于(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) 

请记住,更高度的多项式总是更适合数据。 尽pipe(过度拟合),更高程度的多项式通常导致非常不可能的function(参见Occam's Razor )。 你想在简单性(多项式的程度)和拟合(如最小平方误差)之间find一个平衡点。 在数量上,有testing这个, Akaike信息准则或贝叶斯信息准则 。 这些testing给出了哪个模型是优选的分数。

如果你想把(xi,f(xi))拟合到n次多项式,那么你将build立一个线性最小二乘问题,其中的数据(1,xi,xi,xi ^ 2,…,xi ^ n,f(xi))。 这将返回一组系数(c0,c1,…,cn) ,使得最佳拟合多项式为* y = c0 + c1 * x + c2 * x ^ 2 + … + cn * x ^ n。 *

您可以通过在问题中包含y的幂和xy的组合来推广这两个以上的因variables。

拉格朗日多项式(贴上@jw)给出了一个精确的拟合点,但是如果多项式的次数大于5或6,则可能会出现数值不稳定。

最小二乘为您提供“最佳拟合”多项式,误差定义为各个误差的平方和。 (沿着Y轴的距离,你有的结果和function之间的距离,平方,总结)MATLAB polyfit函数做到这一点,并与多个返回参数,你可以让它自动照顾缩放/偏移问题(例如,如果在x = 312.1和312.3之间有100个点,并且您想要一个6次多项式,那么您将要计算u =(x-312.2)/0.1,这样u值就是分布的在-1和+ =之间)。

注意 ,最小二乘拟合的结果受x轴值分布的强烈影响。 如果x值是等间隔的,那么你会在两端得到更大的错误。 如果你有一个情况,你可以select x值,并且关心已知函数和插值多项式的最大偏差,那么使用切比雪夫多项式会给你一些接近理想最小多项式的东西(这是非常很难计算)。 这在数值食谱中有一些讨论。

编辑:从我收集,这一切运作良好的一个variables的function。 对于多元函数来说,如果程度超过2,那么可能会困难得多。我确实在Google Books上find了一个参考 。

在大学我们有这本书,我仍然觉得非常有用:Conte,de Boor; 基本数值分析; Mc Grow Hill。 相关的段落是6.2:数据拟合。
示例代码来自FORTRAN,并且列表也不太可读,但是解释同时是深刻和清楚的。 你最终理解你在做什么,而不仅仅是做(就像我的数值食谱经验)。
我通常从数字食谱开始,但对于这样的事情,我很快就要抓住Conte-de Boor。

也许更好地张贴一些代码…这是有点简单,但最相关的部分在那里。 它显然依赖于numpy!

 def Tn(n, x): if n==0: return 1.0 elif n==1: return float(x) else: return (2.0 * x * Tn(n - 1, x)) - Tn(n - 2, x) class ChebyshevFit: def __init__(self): self.Tn = Memoize(Tn) def fit(self, data, degree=None): """fit the data by a 'minimal squares' linear combination of chebyshev polinomials. cfr: Conte, de Boor; elementary numerical analysis; Mc Grow Hill (6.2: Data Fitting) """ if degree is None: degree = 5 data = sorted(data) self.range = start, end = (min(data)[0], max(data)[0]) self.halfwidth = (end - start) / 2.0 vec_x = [(x - start - self.halfwidth)/self.halfwidth for (x, y) in data] vec_f = [y for (x, y) in data] mat_phi = [numpy.array([self.Tn(i, x) for x in vec_x]) for i in range(degree+1)] mat_A = numpy.inner(mat_phi, mat_phi) vec_b = numpy.inner(vec_f, mat_phi) self.coefficients = numpy.linalg.solve(mat_A, vec_b) self.degree = degree def evaluate(self, x): """use Clenshaw algorithm http://en.wikipedia.org/wiki/Clenshaw_algorithm """ x = (x-self.range[0]-self.halfwidth) / self.halfwidth b_2 = float(self.coefficients[self.degree]) b_1 = 2 * x * b_2 + float(self.coefficients[self.degree - 1]) for i in range(2, self.degree): b_1, b_2 = 2.0 * x * b_1 + self.coefficients[self.degree - i] - b_2, b_1 else: b_0 = x*b_1 + self.coefficients[0] - b_2 return b_0 

请记住, 近似多项式和find确切的一个有很大的区别。

例如,如果我给你4分,你可以

  1. 用最小二乘法近似一条线
  2. 用最小二乘法近似抛物线
  3. 通过这四点find一个确切的三次函数。

一定要select适合您的方法!

如果您知道如何将最小二乘问题表示为线性代数问题,那么使用Excel的matrix函数来快速拟合是相当容易的。 (这取决于你如何可靠地认为Excel是一个线性代数求解器。)

拉格朗日多项式在某种意义上是适合给定数据点集合的“最简单的”内插多项式。

有时候这是有问题的,因为它可能在数据点之间变化很大。