将多项式拟合到数据
给定一组值(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 i ） 2
理论
重要的思想是，如果将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 j （ x 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的幂和x和y的组合来推广这两个以上的因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分，你可以
- 用最小二乘法近似一条线
- 用最小二乘法近似抛物线
- 通过这四点find一个确切的三次函数。
一定要select适合您的方法！
如果您知道如何将最小二乘问题表示为线性代数问题，那么使用Excel的matrix函数来快速拟合是相当容易的。 （这取决于你如何可靠地认为Excel是一个线性代数求解器。）
拉格朗日多项式在某种意义上是适合给定数据点集合的“最简单的”内插多项式。
有时候这是有问题的，因为它可能在数据点之间变化很大。
