Подгонка полиномов к данным

Да, обычно это делается методом наименьших квадратов. Есть и другие способы определить, насколько хорошо подходит многочлен, но теория наиболее проста для наименьших квадратов. Общая теория называется линейной регрессией.

Лучше всего начать с Числовые рецепты.

р бесплатен и будет делать все, что вы хотите, и даже больше, но для этого потребуется большая кривая обучения.

Если у вас есть доступ к системе Mathematica, вы можете использовать функцию Fit, чтобы выполнить аппроксимацию по методу наименьших квадратов. Я полагаю, что Matlab и его аналог Octave с открытым исходным кодом имеют аналогичную функцию.

Для случая (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)

Если вы хотите подогнать (xi, f (xi)) к полиному степени п, тогда вы должны настроить линейную задачу наименьших квадратов с данными (1, xi, xi, xi ^ 2, ..., xi ^ n, f (xi)). Это вернет набор коэффициентов (c0, c1, ..., сп), так что наилучшим подходящим полиномом будет * y = c0 + c1 * x + c2 * x ^ 2 + ... + cn * x ^ n. *

Вы можете обобщить эти две более чем одной зависимой переменной, включив в задачу степени y и комбинации Икс и y.

полином Лагранжа в некотором смысле является «простейшим» интерполирующим полиномом, который соответствует заданному набору точек данных.

Иногда это проблематично, потому что оно может сильно различаться между точками данных.

Полиномы Лагранжа (как опубликовано в @j w) дают вам точное соответствие в указанных вами точках, но с полиномами степени выше, чем, скажем, 5 или 6, вы можете столкнуться с числовой нестабильностью.

Метод наименьших квадратов дает вам полином «наилучшего соответствия» с ошибкой, определяемой как сумма квадратов отдельных ошибок. (возьмите расстояние по оси Y между точками, которые у вас есть, и полученной функцией, возведите их в квадрат и просуммируйте) Функция MATLAB polyfit делает это, и с несколькими возвращаемыми аргументами вы можете автоматически позаботиться о масштабировании / смещение (например, если у вас есть 100 точек между x = 312,1 и 312,3, и вам нужен полином 6-й степени, вам нужно будет вычислить u = (x-312.2) /0.1, чтобы значения u были распределены от -1 до + =).

ПРИМЕЧАНИЕ, что результаты аппроксимации методом наименьших квадратов являются сильно под влиянием распределения значений оси x. Если значения x расположены на одинаковом расстоянии, на концах будут большие ошибки. Если у вас есть случай, когда вы можете выберите значения x, и вы заботитесь о максимальном отклонении от вашей известной функции и интерполирующего полинома, то использование Полиномы Чебышева даст вам что-то, что близко к идеальному минимаксному полиному (что очень сложно вычислять). Это подробно обсуждается в «Числовых рецептах».

Редактировать: Насколько я понимаю, все это хорошо работает для функций одной переменной. Для многомерных функций это, вероятно, будет намного сложнее, если степень будет больше, скажем, 2. Я нашел ссылка на Google Книги.

Помните, что полином более высокой степени ВСЕГДА лучше соответствует данным. Полиномы более высокой степени обычно приводят к очень маловероятным функциям (см. Бритва Оккама), хотя (переоснащение). Вы хотите найти баланс между простотой (степенью полинома) и соответствием (например, ошибкой наименьших квадратов). Количественно для этого существуют тесты Информационный критерий Акаике или Байесовский информационный критерий. Эти тесты дают оценку, какой модели следует отдать предпочтение.

Помните, есть большая разница между полиномом приблизительный и поиском полинома точный.

Например, если я дам вам 4 балла, вы можете

1. Составьте приблизительную линию с помощью метода наименьших квадратов
2. Приближаем параболу методом наименьших квадратов
3. Найдите кубическую функцию точный через эти четыре точки.

Обязательно выберите метод, который подходит именно вам!

Довольно легко напугать быстрого подбора, используя матричные функции Excel, если вы знаете, как представить задачу наименьших квадратов как задачу линейной алгебры. (Это зависит от того, насколько надежен, по вашему мнению, Excel как решатель линейной алгебры.)

в колледже у нас была эта книга, которую я до сих пор считаю чрезвычайно полезной: Conte, de Boor; элементарный численный анализ; Макроу Хилл. Соответствующий абзац - 6.2: Data Fitting.
. Пример кода представлен на ФОРТРАНЕ, и листинги тоже не очень читабельны, но в то же время объяснения глубокие и ясные. в конечном итоге вы понимаете, что делаете, а не просто делаете это (как мой опыт работы с числовыми рецептами). Обычно я начинаю с числовых рецептов, но для подобных вещей мне быстро приходится брать Конте-де-Бура.

может быть, лучше опубликовать какой-нибудь код ... он немного урезан, но там есть самые важные части. очевидно, он полагается на 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