拟合多项式到数据

Question

给定一组值(x,f(x))，是否有一种方法可以找到最适合数据的给定次数的多项式？我知道多项式插值，它是为了找到一个给定n+1个数据点的n次多项式，但这里有大量的值，我们想要找到低次多项式（最佳线性拟合、最佳二次曲线拟合、最佳三次曲线拟合等）。它可能与最小二乘法有关...

更一般地，当我们有一个多元函数--像(x,y,f(x,y))这样的点--并且想要在变量中找到给定次数的最佳多项式p(x,y)。（具体而言是多项式，不是样条或傅里叶级数。）

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

- ShreevatsaR

Answer 1

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

非多项式插值

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

什么是最小二乘法？

"最小二乘法"是多项式拟合中衡量多项式拟合程度的一种定义/准则/度量方法（还有其他方法，但这是最简单的）。假设您正在尝试将一个多项式p(x,y)=a+bx+cy+dx²+ey²+fxy拟合到一些给定的数据点(xᵢ,yᵢ,Zᵢ)上（其中“Zᵢ”在问题中是“f(xᵢ,yᵢ)”），使用最小二乘法，问题就是找到“最佳”的系数(a,b,c,d,e,f)，使得被最小化的是“平方残差和”，即

S = ∑ᵢ (a + bxᵢ + cyᵢ + dxᵢ² + eyᵢ² + fxᵢyᵢ - Zᵢ)²

理论

"

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

∑_i2(a + … + fx_iy_i - Z_i)x_iy_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中的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 nodes进行多项式插值; 不确定在最小二乘法中是否也是如此。）
• 过拟合：更高次数的多项式始终能够更好地拟合数据。如果将degree改为3、4或5，它仍然大多识别相同的二次多项式（更高次数的系数为0），但对于更大的次数，它开始拟合更高次数的多项式。但即使是6次，取更多的n（比如200个数据点而不是20个）仍然适合二次多项式。所以教训是要避免过拟合，这可能有助于尽可能多地采集数据点。
• 可能存在我不完全理解的numerical stability问题。
• 如果不需要多项式，可以使用其他类型的函数获得更好的拟合效果，例如splines（分段多项式）。

- ShreevatsaR

---

@Jason：你确定Chebyshev节点被认为是最适合选择最小二乘的点吗？似乎有一个不同的问题，即选择Chebyshev多项式本身作为逼近多项式，用于与最小二乘不同类型的拟合——“极小化极差”多项式。 - ShreevatsaR

df/dx = 0并不一定意味着f被最小化，它也可能被最大化。 - quant_dev

是的，但这不是我们在这里说的。我们说的是（假设函数具有偏导数等），任何最小值都会出现在梯度为0的点上（或边界上）。 - ShreevatsaR

关于您对数值稳定性的担忧：定义多项式（即“单项式的线性组合”）是一件危险的事情，因为（用非数学的话来说），次数大于4的单项式在0附近的区域内非常相似，然后它们就会“疯狂增长”。更好的方法是决定您要拟合数据的区间，重新定义自变量，使其实际上适合于（-1，1），并寻找良好多项式的线性组合而不是单项式。我会使用Chebyshev集合。 - mariotomo

@mariotomo：谢谢，你说了之后就很有道理了 :) 很好的观点。 - ShreevatsaR

Answer 2

是的，通常使用最小二乘法来实现这个目标。还有其他指定多项式拟合程度的方法，但最小二乘法的理论最简单。一般的理论被称为线性回归。

你最好从Numerical Recipes开始学起。

R是免费的，可以完成你想要的所有任务，但它有一个很大的学习曲线。

如果你有Mathematica的访问权限，可以使用Fit函数进行最小二乘拟合。我想Matlab及其开源对应物Octave也有类似的功能。

- John D. Cook

---

这很有帮助，但您知道其中哪些可以进行多元拟合吗？ - ShreevatsaR

最小二乘法可以是多元的。吉尔·斯特朗（Gil Strang）的《应用数学导论》中有一个非常好的、易读的讨论。 - duffymo

是的，谢谢...当我询问多元拟合的评论时，我对最小二乘法还不够了解 :-) - ShreevatsaR

Answer 3

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

- J.F. Sebastian

Answer 4

请注意，高阶多项式始终更适合数据。然而，高阶多项式通常会导致高度不可能的函数（参见奥卡姆剃刀），即过度拟合。您需要在简单性（多项式程度）和适配性（例如最小二乘误差）之间找到平衡点。定量上，有一些测试可以进行，如Akaike信息准则或贝叶斯信息准则。这些测试给出了哪个模型更好的分数。

- Fredriku73

---

是的，我后来意识到适合度和简洁性之间存在某种权衡，就像你所说的。感谢您提供有关标准的信息。 - ShreevatsaR

Answer 5

如果您想将点（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的组合来将此推广到多于一个依赖变量的情况中。

- David Norman

Answer 6

Lagrange多项式（正如@jw所发表的）可以在指定点上给出精确拟合，但是使用高于5或6次的多项式可能会遇到数值不稳定性问题。

最小二乘法可以给出“最佳拟合”多项式，并将误差定义为各个误差的平方和。（将您拥有的点与结果函数之间的沿y轴的距离取出，平方并相加）。MATLAB的polyfit函数可以实现这一点，并且通过多个返回参数，您可以让它自动处理缩放/偏移问题（例如，如果您有100个点都在x = 312.1和312.3之间，并且您需要一个6次多项式，那么您需要计算u =（x-312.2）/0.1，使得u值分布在-1和+ =之间）。

请注意，最小二乘拟合的结果受x轴值分布的影响很大。如果x值等间距，则在两端会得到更大的误差。如果您有一个情况可以选择x值，并且您关心已知函数和插值多项式的最大偏差，则使用切比雪夫多项式将给您接近完美的极小值多项式（这非常难计算）。这在Numerical Recipes中进行了详细讨论。 编辑：据我所知，这对于单变量函数都适用良好。对于多元函数，如果度数超过2，则可能会更加困难。我在Google Books上找到了一份参考资料。

- Jason S

---

顺便感谢您提供的参考资料。相关内容在几页后，第231页的4.10.4部分。同样的方法也适用于高阶多元多项式，虽然存在“过度拟合”的问题需要注意。 - ShreevatsaR

Answer 7

在大学时，我们有一本书，我仍然觉得非常有用：Conte，de Boor; elementary numerical analysis; Mc Grow Hill。相关段落是6.2：数据拟合。
示例代码采用FORTRAN编写，而且清单也不太易读，但解释既深入又清晰。你最终会理解你所做的事情，而不仅仅是做它（这是我对Numerical Recipes的经验）。
我通常从Numerical Recipes开始，但对于像这样的东西，我很快就要抓住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

- mariotomo

---

再次感谢，非常清晰。顺便问一下，为什么将范围归一化到（-1,1）是好的呢？ - ShreevatsaR

因为在这个范围内切比雪夫多项式的行为良好。事实上，在该范围内，您可以这样描述它们：T_n(x) = cos(n*acos(x))。对于不在(-1, 1)范围内的x，此公式没有意义。 - mariotomo

我已经对我的模块进行了测试，与numpy.polyfit（您指向的页面中的单个示例）进行了比较，令我有些惊讶的是，即使是在15位小数的情况下，我的拟合结果也与numpy.polyfit相匹配（甚至可以进行外推）。我应该尝试更多糟糕的条件...如果它们仍然匹配，那么可能numpy在幕后使用切比雪夫多项式并返回相应的单项式系数... - mariotomo

Answer 8

记住，近似多项式和找到一个精确的多项式之间有很大的区别。

例如，如果我给你4个点，你可以：

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

一定要选择适合你的方法！

- stalepretzel

---

是的，我知道 :-) 这就是为什么我在问题中提到了“多项式插值”，它可以通过四个点找到一个精确的三次曲线，或者通过 n+1 个点找到一个精确的 n 次多项式。 - ShreevatsaR

Answer 9

如果你知道如何将最小二乘问题表示为线性代数问题，那么使用 Excel 的矩阵函数快速得出一个适配方案是相当简单的。（这取决于你认为 Excel 作为线性代数求解器的可靠程度。）

- duffymo

Answer 10

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

有时它会出现问题，因为它可能在数据点之间变化很大。

- j w

---

此外，拉格朗日插值多项式具有n个点的n-1次度数 - 这就是我在关于多项式插值的问题中所写的内容 - 它并不能为给定度数提供最佳拟合多项式。 - ShreevatsaR