Присоединение полиномов к данным

Ответ 1

Спасибо за ответы. Вот еще одна попытка их обобщения. Простите, если я говорю слишком много "очевидных" вещей: раньше я ничего не знал о наименьших квадратах, поэтому все было для меня новичком.

НЕ полиномиальная интерполяция

Полиномиальная интерполяция устанавливает полином степени n с учетом n+1 точек данных, например. находя кубику, которая проходит точно через четыре заданные точки. Как сказал в этом вопросе, это было не так, как хотелось бы, у меня было много очков, и мне нужен был многочлен с малой степенью (который будет только приблизительно соответствовать, если только нам не повезет), но поскольку некоторые из ответов настаивали на разговоре об этом, я должен упомянуть их:) полином Лагранжа, матрица Вандермонда и т.д.

Что такое наименьшие квадраты?

"Наименьшие квадраты" - это конкретное определение/критерий/ "метрика" "насколько хорошо" подходит многочлен. (Есть другие, но это проще всего.) Скажем, вы пытаетесь подогнать многочлен p (x, y) = a + bx + cy + dx ² + ey ² + fxy к некоторым заданным точкам данных (x _i, y _i, Z _i) (где "Z _i" было "f (x _i, y _i)" в вопросе). С наименьшими квадратами проблема состоит в том, чтобы найти "лучшие" коэффициенты (a, b, c, d, e, f), такие, что минимизированное ( "наименьшее" ) означает "сумму квадратов остатков", а именно

S = & sum; _i (a + bx _i + cy _i + dx _i ² + ey _i ² + fx _i y _i - Z _i) ²

Теория

Важная идея состоит в том, что если вы посмотрите на S как функцию (a, b, c, d, e, f), то S будет сведено к минимуму в точке, где ее градиент равен 0. Это означает, что, например, & part; S/& part; f = 0, то есть, что

& sum; _i 2 (a + & hellip; + fx _i y _i - Z _i) x _i y _i= 0

и аналогичные уравнения для a, b, c, d, e. Заметим, что это просто линейные уравнения в & hellip; f. Поэтому мы можем решить их с помощью Gaussian elim или любого из обычных методов.

Это все еще называется "линейными наименьшими квадратами", потому что хотя функция, которую мы хотели, была квадратичным многочленом, она по-прежнему линейна по параметрам (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. Как поясняется в его документации, он принимает матрицу A значений f _j (x <суб > явспом > ). (Теория состоит в том, что она находит псевдо-обратную сторону Мура-Пенроуза). В нашем примере, в котором используется (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 может это сделать. Numericical Recipes обсуждает методы его реализации (в C, Fortran).

Обеспокоенность

• На него сильно влияет способ выбора точек. Когда вместо случайных точек у меня была x=y=range(20), она всегда создавала 1.33x ² + 1.33xy + 1.33y ² что было озадачивающе... пока я не понял, что потому что у меня всегда было x[i]=y[i], многочлены были одинаковыми: x ² + 2xy + y ²= 4x ²= (4/3 ) (х ² + х + у ²). Итак, мораль состоит в том, что важно тщательно выбирать точки, чтобы получить "правильный" многочлен. (Если вы можете выбрать, вы должны выбрать узлы Чебышева для полиномиальной интерполяции, не уверены, что то же самое верно и для наименьших квадратов.)
  
• Overfitting: полиномы более высокой степени могут всегда лучше соответствовать данным. Если вы измените значение degree на 3 или 4 или 5, он по-прежнему в основном распознает один и тот же квадратичный многочлен (коэффициенты равны 0 для более высоких степеней), но для больших степеней он начинает устанавливать полиномы более высокой степени. Но даже со степенью 6, принимая больше n (больше данных вместо 20, скажем 200), по-прежнему подходит квадратичный многочлен. Таким образом, мораль заключается в том, чтобы избегать переобучения, для чего это могло бы помочь как можно большему количеству точек данных.
• Могут возникнуть проблемы численная стабильность Я не совсем понимаю.
• Если вам не нужен полином, вы можете получить лучшие приемы с другими видами функций, например. сплайны (кусочные полиномы).

Ответ 2

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

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

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

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

Ответ 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)

Ответ 4

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

Ответ 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.

Ответ 6

Полиномы Лагранжа (как отправлено @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 Книгах.

Ответ 7

в колледже у нас была эта книга, которую я до сих пор считаю чрезвычайно полезной: Конте, де Бур; элементарный численный анализ; Mc Grow Hill. Соответствующий параграф - 6.2: Установка данных.
пример кода поставляется в FORTRAN, и списки также не очень читаемы, но объяснения являются глубокими и ясными в одно и то же время. вы в конечном итоге понимаете, что вы делаете, а не просто делаете это (как и мой опыт Numericical Recipes).
Я обычно начинаю с Numerical Recipes, но для таких вещей мне быстро нужно схватить Конте-де-Бура.

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

Ответ 8

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

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

• Приблизительная строка с методом, подобным наименьшему квадрату
• Приблизительная парабола с методом, подобным наименьшим квадратам
• Найти точную кубическую функцию через эти четыре точки.

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

Ответ 9

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

Ответ 10

lagrange polyomial в некотором смысле является "самым простым" интерполирующим полиномом, который подходит для заданного набора точек данных.

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