在给定一组值(x,f(x)) 的情况下，有没有办法找到最适合数据的给定次数的多项式？
我知道polynomial interpolation，它是用于在给定n+1 数据点的情况下找到次数为n 的多项式，但是这里有大量的值，我们想找到一个低次数多项式（找到最佳线性拟合、最佳二次方、最佳三次方等）。可能与least squares有关...
更一般地说，当我们有一个多元函数时，我想知道答案——比如像(x,y,f(x,y)) 这样的点——并且想在变量中找到给定度数的最佳多项式 (p(x,y))。 （特别是多项式，不是样条或傅里叶级数。）
理论和代码/库（最好是 Python，但任何语言都可以）都会很有用。
【问题讨论】：
标签： math statistics
感谢大家的回复。这是总结它们的另一种尝试。如果我说了太多“显而易见”的事情，请原谅：我以前对最小二乘一无所知，所以一切对我来说都是新的。
Polynomial interpolation 在给定n+1 数据点的情况下拟合度为n 的多项式，例如找到一个恰好通过四个给定点的立方。正如问题中所说，这不是我想要的——我有很多点并且想要一个小次数多项式（它只会近似拟合，除非我们很幸运）——但是因为有些答案坚持要讲，我应该提一下：）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 在其 gradient is 0 的某个点上是 minimized。这意味着例如∂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 Recipes 开头。
R 是免费的，可以做任何你想做的事情，但它有一个很大的学习曲线。
如果您可以访问 Mathematica，则可以使用 Fit 函数进行最小二乘拟合。我想 Matlab 和它的开源对应 Octave 也有类似的功能。
【讨论】：
在大学时，我们有这本书，我仍然觉得它非常有用：Conte, de Boor；基本数值分析；麦格罗希尔。相关段落是 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
【讨论】：
如果您知道如何将最小二乘问题表示为线性代数问题，那么使用 Excel 的矩阵函数很容易快速拟合。 （这取决于您认为 Excel 作为线性代数求解器的可靠性。）
【讨论】：
请记住，近似多项式和找到一个精确项之间存在很大差异。
例如，如果我给你 4 分，你可以
一定要选择适合你的方法！
【讨论】：
请记住，更高次的多项式总是更适合数据。更高次的多项式通常会导致非常不可能的函数（请参阅Occam's Razor），尽管（过度拟合）。您想在简单性（多项式次数）和拟合（例如最小二乘误差）之间找到平衡。在数量上，有针对此的测试，Akaike Information Criterion 或Bayesian Information Criterion。这些测试给出了首选模型的分数。
【讨论】：
拉格朗日多项式（如@j w 发布的）在您指定的点处为您提供精确拟合，但如果多项式的次数超过 5 或 6，您可能会遇到数值不稳定。
最小二乘法为您提供“最佳拟合”多项式，其中误差定义为各个误差的平方和。 （取你所拥有的点和产生的函数之间沿 y 轴的距离，将它们平方，然后求和）MATLAB polyfit 函数会执行此操作，并且使用多个返回参数，你可以让它自动处理缩放/偏移问题（例如，如果您在 x=312.1 和 312.3 之间有 100 个点，并且您想要一个 6 次多项式，那么您将需要计算 u = (x-312.2)/0.1 所以 u 值分布在 -1 和 += 之间）。
注意最小二乘拟合的结果强烈受 x 轴值分布的影响。如果 x 值是等距的，那么你会在末端得到更大的错误。如果您有一个可以选择 x 值的情况，并且您关心与已知函数的最大偏差和插值多项式，那么使用Chebyshev polynomials 将为您提供接近于完美的极小极大多项式（很难计算）。这在数值食谱中有详细的讨论。
编辑：据我所知，这一切都适用于一个变量的函数。对于多变量函数，如果度数大于 2，则可能会更加困难。我确实找到了 reference on Google Books。
【讨论】：
lagrange polynomial 在某种意义上是适合给定数据点集的“最简单”插值多项式。
有时会出现问题，因为数据点之间的差异可能很大。
【讨论】：
如果您想将 (xi, f(xi)) 拟合到 n 次多项式，那么您将使用数据设置线性最小二乘问题(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 的组合来概括这两个以上的因变量。
【讨论】：
对于 (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)
【讨论】：