54

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

我知道多项式插值，它用于查找n给定n+1数据点的次数多项式，但是这里有大量值，我们想找到一个低次数多项式（找到最佳线性拟合、最佳二次、最佳三次等）。 ）。它可能与最小二乘有关......

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

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

4

10 回答 10

61

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

非多项式插值

多项式插值拟合n给定n+1数据点的度数多项式，例如找到精确通过四个给定点的三次方。正如问题中所说，这不是我想要的——我有很多点并且想要一个小次数多项式（除非我们很幸运，否则它只会近似拟合）——但由于一些答案坚持谈论关于它，我应该提到它们 :)拉格朗日多项式、范德蒙德矩阵等。

什么是最小二乘？

“最小二乘”是多项式拟合“有多好”的特定定义/标准/“度量”。（还有其他的，但这是最简单的。）假设您试图将多项式 p(x,y) = a + bx + cy + dx 2 + ey 2 + fxy 拟合到某些给定的数据点 (x i ,y i ,Z i )（其中“Z i ”在问题中是“f(x i ,y i )”）。使用最小二乘法的问题是找到“最佳”系数（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 的类似方程。请注意，这些只是 a...f 中的线性方程。所以我们可以用高斯消元法或任何常用方法来解决它们。

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

代码

对于单变量情况（当只有变量 x - f j是单项式 x j时），有 Numpy's 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 2 +2xy+y 2 +0.01。[最后一项有时为 -0.01，有时为 0，这是可以预料的，因为我们添加了随机噪声。]

Python+Numpy/Scipy 的替代品是R和计算机代数系统：Sage、Mathematica、Matlab、Maple。甚至 Excel 也能做到。Numerical Recipes讨论了我们自己实现它的方法（在 C、Fortran 中）。

关注点

  • 它受到如何选择点的强烈影响。当我有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 )。因此，重要的是仔细选择点以获得“正确的”多项式。（如果可以选择，您应该选择Chebyshev 节点进行多项式插值；不确定最小二乘是否也是如此。）
  • 过度拟合：更高次的多项式总是可以更好地拟合数据。如果将 更改degree为 3 或 4 或 5，它仍然主要识别相同的二次多项式（系数为 0 表示更高阶项），但对于更大的阶数，它开始拟合更高阶多项式。但即使使用 6 次，取更大的 n（更多数据点而不是 20，比如 200）仍然适合二次多项式。因此，道德是避免过度拟合，这可能有助于获取尽可能多的数据点。
  • 可能存在我不完全理解的数值稳定性问题。
  • 如果您不需要多项式，则可以更好地拟合其他类型的函数，例如样条（分段多项式）。
于 2008-12-20T08:45:51.523 回答
7

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

您最好的选择可能是从数字食谱开始。

R是免费的，可以做任何你想做的事情，但它有一个很大的学习曲线。

如果您有权访问 Mathematica，则可以使用 Fit 函数进行最小二乘拟合。我想 Matlab 和它的开源对应 Octave 有类似的功能。

于 2008-12-19T21:06:37.783 回答
5

对于 (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)
于 2008-12-19T21:33:23.440 回答
4

请记住，更高次的多项式总是更好地拟合数据。更高次的多项式通常会导致非常不可能的函数（参见奥卡姆剃刀），尽管（过度拟合）。您想在简单性（多项式次数）和拟合（例如最小二乘误差）之间找到平衡。在数量上，有这方面的测试，Akaike 信息准则或贝叶斯信息准则。这些测试给出了首选模型的分数。

于 2008-12-20T09:52:55.973 回答
2

如果您想将(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的组合来概括这两个以上的因变量。

于 2008-12-19T21:35:16.230 回答
2

拉格朗日多项式（如@jw 发布的）在您指定的点上为您提供精确拟合，但如果多项式的次数超过 5 或 6，您可能会遇到数值不稳定。

最小二乘法为您提供“最佳拟合”多项式，其中误差定义为各个误差的平方和。（取你所拥有的点和产生的函数之间沿 y 轴的距离，将它们平方并求和）MATLABpolyfit函数会执行此操作，并且使用多个返回参数，你可以让它自动处理缩放/偏移问题（例如，如果您在 x=312.1 和 312.3 之间有 100 个点，并且您想要一个 6 次多项式，您将需要计算 u = (x-312.2)/0.1，因此 u 值分布在-1 和 +=)。

请注意，最小二乘拟合的结果受到 x 轴值分布的强烈影响。如果 x 值是等距的，那么你会在末端得到更大的错误。如果您有一个可以选择x 值的情况，并且您关心与已知函数的最大偏差和插值多项式，那么使用Chebyshev 多项式将为您提供接近完美极小极大多项式的东西（这是非常难以计算）。这在数字食谱中有详细讨论。

编辑：据我所知，这一切都适用于一个变量的函数。对于多元函数，如果度数超过 2，则可能会更加困难。我确实在 Google Books 上找到了参考。

于 2008-12-19T22:38:26.273 回答
2

在大学里，我们有这本书，我仍然觉得它非常有用：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
于 2009-07-28T12:39:42.003 回答
0

请记住，近似多项式和找到一个精确的多项式之间存在很大差异。

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

  1. 用最小二乘法等方法逼近一条线
  2. 用最小二乘法等方法逼近抛物线
  3. 通过这四个点找到一个精确的三次函数。

请务必选择适合您的方法！

于 2008-12-21T16:34:25.680 回答
0

如果您知道如何将最小二乘问题表示为线性代数问题，那么使用 Excel 的矩阵函数很容易快速拟合。（这取决于您认为 Excel 作为线性代数求解器的可靠性。）

于 2008-12-28T15:57:24.423 回答
-1

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

有时会出现问题，因为数据点之间的差异可能很大。

于 2008-12-19T21:40:13.847 回答