Dopasowanie wielomianów do danych

Tak, zazwyczaj odbywa się to za pomocą najmniejszych kwadratów. Istnieją inne sposoby określania, jak dobrze pasuje wielomian, ale teoria jest najprostsza dla najmniejszych kwadratów. Ogólna teoria nazywa się regresją liniową.

Najlepiej zacząć od numerycznych przepisów .

R jest darmowy i zrobi wszystko, co chcesz i więcej, ale ma dużą krzywą uczenia się.

Jeśli masz dostęp do Mathematica, możesz użyć funkcji Fit, aby najmniej kwadratów pasuje. Domyślam się, że Matlab i jego open source odpowiednik Octave mają podobną funkcję.

For (x, f(x)) case:

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)

Należy pamiętać, że wielomian wyższego stopnia zawsze lepiej pasuje do danych. Wielomiany wyższego stopnia zazwyczaj prowadzą do wysoce nieprawdopodobnych funkcji (Zobacz brzytwa Occama), choć (overfitting). Chcesz znaleźć równowagę między prostotą (stopień wielomianu) a dopasowaniem (np. błąd najmniejszego kwadratu). Ilościowe, istnieją testy na to, Akaike Information Criterion lub Bayessian Information Criterion . Te testy dają wynik, który model ma być lepiej.

Jeśli chcesz dopasować (xi, f (xi)) do wielomianu stopnia n wtedy ustawisz liniowy problem najmniejszych kwadratów z danymi (1, xi, xi, xi^2,..., xi^n, f (xi)). zwróci to zbiór współczynników (c0, c1, ..., cn) tak, że najlepiej dopasowanym wielomianem jest *y = C0 + c1 * x + c2 * x^2+... + cn * x^n. *

Można uogólnić te dwie więcej niż jedną zmienną zależną przez uwzględnienie potęg y i kombinacji x i y w problemie.

Wielomiany Lagrange ' a (jak napisał @j w) dają dokładne dopasowanie do punktów, które określisz, ale z wielomianami stopnia większego niż powiedzmy 5 lub 6, możesz napotkać niestabilność numeryczną.

Najmniejsze kwadraty daje wielomian" najlepiej dopasowany " z błędem zdefiniowanym jako suma kwadratów poszczególnych błędów. (weź odległość wzdłuż osi y między punktami, które masz, a funkcją, która daje wynik, kwadrat je i zsumuj je) funkcja MATLAB polyfit robi to i z wieloma jeśli masz 100 punktów między x = 312.1 i 312.3, i chcesz wielomian 6 stopnia, będziesz chciał obliczyć u = (x-312.2) / 0.1, więc wartości u są rozłożone między -1 i+=).

Zauważ , że wyniki pasowań najmniejszych kwadratów są silnie pod wpływem rozkładu wartości osi X. Jeśli wartości x są równomiernie rozmieszczone, wtedy pojawią się większe błędy w koniec. Jeśli masz przypadek, w którym możesz wybrać wartości x i zależy ci na maksymalnym odchyleniu od znanej funkcji i wielomianu interpolującego, to użycie wielomianów Czebyszewa da ci coś, co jest bliskie idealnemu wielomianowi minimax (który jest bardzo trudny do obliczenia). Jest to omawiane w pewnej mierze w numerycznych recepturach.

Edit: z tego co wiem, to wszystko działa dobrze dla funkcji jednej zmiennej. Dla wielowymiarowych funkcje może to być znacznie trudniejsze, jeśli stopień jest większy niż, powiedzmy, 2. Znalazłem odnośnik w Google Books.

Na studiach mieliśmy tę książkę, która nadal uważam za niezwykle przydatna: Conte, de Boor; elementarna Analiza numeryczna; Mc Grow Hill. Odpowiedni punkt to 6.2: dopasowanie danych.
przykładowy kod pojawia się w FORTRAN, a listy również nie są zbyt czytelne, ale wyjaśnienia są głębokie i jasne jednocześnie. w końcu rozumiesz, co robisz, a nie tylko to robisz (tak jak moje doświadczenie z numerycznymi przepisami).
Zazwyczaj zaczynam od numerycznych przepisów, ale dla takich rzeczy szybko trzeba złapać Conte-de Boor.

Może lepiej wrzucić jakiś kod... jest trochę rozebrany, ale najważniejsze części są tam. oczywiście opiera się na 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

Pamiętaj, jest duża różnica między przybliżeniem wielomianu a znalezieniem dokładnego .

Na przykład, jeśli dam ci 4 punkty, możesz

1. przybliżenie linii metodą jak najmniejszych kwadratów
2. przybliżenie paraboli metodą jak najmniejszych kwadratów
3. Znajdź dokładną funkcję sześcienną przez te cztery punkty.

Pamiętaj, aby wybrać metodę, która jest odpowiednia dla Ciebie!

Łatwo jest przestraszyć szybkie dopasowanie za pomocą funkcji macierzy Excela, jeśli wiesz, jak przedstawić problem najmniejszych kwadratów jako problem algebry liniowej. (To zależy od tego, jak niezawodny jest Excel jako rozwiązywanie algebry liniowej.)

Wielomian lagrange ' a jest w pewnym sensie "najprostszym" wielomianem interpolacyjnym, który pasuje do danego zbioru punktów danych.

Jest to czasami problematyczne, ponieważ może się bardzo różnić między punktami danych.