54

一連の値(x,f(x))を指定して、データに最適な次数の多項式を見つける方法はありますか?

与えられたデータポイントの次数の多項式を見つけるための多項式補間を知っていますが、ここには多数の値があり、低次の多項式を見つけたいと考えています (最良の線形適合、最良の 2 次、最良の 3 次などを見つけます。 ）。最小二乗法に関連している可能性があります...nn+1

より一般的には、多変量関数 (例えば のような点) があり、変数内の特定の(x,y,f(x,y))次数の最適な多項式 ( p(x,y)) を見つけたい場合の答えを知りたいと思います。(具体的には、スプラインやフーリエ級数ではなく、多項式です。)

理論とコード/ライブラリ (できれば Python ですが、どの言語でもかまいません) の両方が役立ちます。

4

10 に答える 10

61

皆様からの返信ありがとうございます。それらを要約する別の試みがあります。あまりにも多くの「明白な」ことを言った場合はご容赦ください。以前は最小二乗法について何も知らなかったので、すべてが私にとって新しいものでした。

多項式補間ではありません

多項式補間nは、与えられたデータポイントの次数の多項式を近似しn+1ます。たとえば、与えられた4つのポイントを正確に通過する3次を見つけます。質問で述べたように、これは私が望んでいたことではありませんでした—私は多くのポイントを持ち、小さな次数の多項式を望んでいました（運が良ければ、ほぼ適合します）—しかし、いくつかの答えは話すことを主張したのでそれについて、私はそれらに言及する必要があります:)ラグランジュ多項式、ファンデルモンド行列など。

最小二乗とは何ですか？

「最小二乗」は、多項式が「どれだけうまく」適合するかについての特定の定義/基準/「メトリック」です。（他にもありますが、これが最も簡単です。）多項式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）の関数として見ると、勾配が0になる点でSが最小化されるということです。これは、たとえば∂S/∂f= 0、つまり

∑ i 2（a+…+fx i y i --Z i）x i y i = 0

およびa、b、c、d、eの同様の方程式。これらはa…fの単なる線形方程式であることに注意してください。したがって、ガウスの消去法または通常の方法のいずれかを使用してそれらを解決できます。

これは、「線形最小二乗」と呼ばれます。これは、必要な関数が2次多項式であったにもかかわらず、パラメーター（a、b、c、d、e、f）では線形であるためです。p（x、y）を単なる多項式（=「単項式の線形結合」）ではなく、任意の関数f jの「線形結合」にしたい場合も、同じことが機能することに注意してください。

コード

単変量の場合（変数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があります。ドキュメントで説明されているように、値f j（x i ）の行列Aを取ります。（理論は、Aのムーア-ペンローズ疑似逆行列を見つけるというものです。）（x i、y i、Z i）を含む上記の例では、多項式を当てはめると、f jが単項式x （） y （）であることを意味します。以下は、最良の2次式（または「次数= 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 + y2+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）。したがって、道徳は、「正しい」多項式を取得するためにポイントを慎重に選択することが重要であるということです。（選択できる場合は、多項式補間にチェビシェフノードを選択する必要があります。最小二乗法についても同じことが当てはまるかどうかはわかりません。）
  • 過剰適合：高次の多項式は常にデータをより適切に適合させることができます。を3、4、または5に変更しdegreeても、ほとんど同じ2次多項式（高次の項の係数は0）を認識しますが、次数が大きいほど、高次の多項式の近似を開始します。ただし、次数が6の場合でも、nを大きくすると（20ではなくより多くのデータポイント、たとえば200）、2次多項式に適合します。したがって、モラルは過剰適合を回避することです。そのため、可能な限り多くのデータポイントを取得することが役立つ場合があります。
  • 私が完全には理解していない数値安定性の問題があるかもしれません。
  • 多項式が必要ない場合は、スプライン（区分的多項式）などの他の種類の関数との適合性を高めることができます。
于 2008-12-20T08:45:51.523 に答える
7

はい、これは通常、最小二乗法を使用して行われます。多項式の適合度を指定する方法は他にもありますが、理論は最小二乗法が最も単純です。一般的な理論は線形回帰と呼ばれます。

あなたの最善の策は、おそらくNumerical Recipesから始めることです。

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

次数の高い多項式は常にデータによりよく適合することに注意してください。より高い次数の多項式は、通常、非常にありそうもない関数になります ( Occam's Razorを参照してください) が、(オーバーフィッティング)。単純さ (多項式の次数) と適合 (最小二乗誤差など) のバランスを見つける必要があります。定量的には、赤池情報量基準またはベイジアン情報量基準というテストがあります。これらのテストは、どのモデルが優先されるかのスコアを示します。

于 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の組み合わせを含めることにより、この2つの複数の従属変数を一般化できます。

于 2008-12-19T21:35:16.230 に答える
2

大学では、私が今でも非常に役立つと思うこの本を持っていました。初等数値解析; マック・グロー・ヒル。関連する段落は 6.2: データ フィッティングです。
サンプル コードは FORTRAN で提供されており、リストもあまり読みにくいですが、同時に説明は深く明確です。ただやっているだけでなく、自分が何をしているのかを理解することになります (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
于 2009-07-28T12:39:42.003 に答える
2

ラグランジュ多項式（@jwが投稿したように）は、指定した点に正確に適合しますが、5または6以上の次数の多項式を使用すると、数値が不安定になる可能性があります。

最小二乗は、個々のエラーの二乗和として定義されたエラーを持つ「最適な」多項式を提供します。(あなたが持っているポイントと結果として得られる関数の間の y 軸に沿った距離を取り、それらを二乗し、それらを合計します) MATLABpolyfit関数はこれを行い、複数の戻り引数を使用して、自動的にスケーリング/処理を行うことができます。オフセットの問題 (たとえば、x=312.1 と 312.3 の間に 100 個のポイントがあり、6 次多項式が必要な場合、u = (x-312.2)/0.1 を計算して、u 値が-1 および +=)。

最小二乗近似の結果は、x 軸の値の分布に大きく影響されることに注意してください。x 値が等間隔である場合、両端でより大きなエラーが発生します。x 値を選択でき、既知の関数と補間多項式からの最大偏差に関心がある場合、チェビシェフ多項式を使用すると、完全なミニマックス多項式に近いものが得られます (これは非常に重要です)。計算が難しい）。これについては、数値レシピで詳しく説明しています。

編集：私が集めたものから、これはすべて1つの変数の関数に対してうまく機能します。多変量関数の場合、次数がたとえば 2 を超えると、はるかに難しくなる可能性があります。Google ブックスで参考文献を見つけました。

于 2008-12-19T22:38:26.273 に答える
0

最小二乗問題を線形代数の問題として表す方法を知っていれば、Excel の行列関数を使用して簡単に当てはめることができます。(それは、Excel が線形代数ソルバーとしてどれだけ信頼できるかによって異なります。)

于 2008-12-28T15:57:24.423 に答える
0

多項式を近似することと、正確な多項式を見つけることには大きな違いがあることを覚えておいてください。

たとえば、私があなたに 4 点を与えるとしたら、あなたは

  1. 最小二乗法などで直線を近似する
  2. 最小二乗法などで放物線を近似する
  3. これらの 4 点を通る正確な3 次関数を見つけます。

必ず自分に合った方法を選んでください！

于 2008-12-21T16:34:25.680 に答える
-1

ラグランジュ多項式は、ある意味で、特定のデータポイントのセットに適合する「最も単純な」補間多項式です。

データポイント間で大きく異なる可能性があるため、問題が発生する場合があります。

于 2008-12-19T21:40:13.847 に答える