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

Есть ли способ, заданный набором значений (x,f(x)) , найти многочлен данной степени, который лучше всего подходит для данных?

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

В общем, я хотел бы знать ответ, когда мы имеем многомерную функцию – точки, такие как (x,y,f(x,y)) , скажем – и хотим найти лучший многочлен ( p(x,y) ) данной степени по переменным. (В частности, полином, а не сплайнами или рядами Фурье).

И теория, и код / ​​библиотеки (желательно на Python, но любой язык в порядке) были бы полезны.

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

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

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

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

«Наименьшие квадраты» – это конкретное определение / критерий / «метрика» «насколько хорошо» подходит многочлен. (Есть и другие, но это проще всего.) Скажем, вы пытаетесь сопоставить полином 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), то S минимизируется в точке, где его gradleиент равен 0 . Это означает, что, например, ∂S / ∂f = 0, т.е.

Σ i 2 (a + … + fx i y i – Z i ) x i y i = 0

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

Это все еще называется «линейными наименьшими квадратами», потому что хотя функция, которую мы хотели, была квадратичным многочленом, она по-прежнему линейна по параметрам (a, b, c, d, e, f). Заметим, что то же самое работает, когда мы хотим, чтобы p (x, y) была любой «линейной комбинацией» произвольных функций f j вместо простого многочлена (= «линейная комбинация одночленов»).

Код

Для одномерного случая (когда существует только переменная x – fj – одночлены x j ), существует многообразие 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. Как поясняется в его документации , matrix A принимает значения f j ( x i ). (Теория состоит в том, что она находит псевдоверсию Мура-Пенроуза А.). В нашем примере, включающем (x i , y i , Z i ), подгонка полинома означает, что f j – monoмы 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 2 + 2xy + y 2 +0.01. [Последний термин иногда -0.01, а иногда и 0, что следует ожидать из-за случайного шума, который мы добавили.]

Альтернативой Python + Numpy / Scipy являются R и компьютерные алгебраические системы: Sage , Mathematica, Matlab, Maple. Даже Excel может это сделать. В Numerical Recipes обсуждаются методы его реализации (в C, Fortran).

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

  • Это сильно зависит от того, как выбираются точки . Когда вместо случайных точек x=y=range(20) меня был 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 ). Итак, мораль состоит в том, что важно тщательно выбирать точки, чтобы получить «правильный» многочлен. (Если вы можете выбрать, вы должны выбрать узлы Чебышева для полиномиальной интерполяции, не уверен, что то же самое верно и для наименьших квадратов.)
  • Overfitting : полиномы более высокой степени могут всегда лучше соответствовать данным. Если вы меняете degree на 3 или 4 или 5, она по-прежнему в основном распознает один и тот же квадратичный многочлен (коэффициенты равны 0 для более высоких степеней), но для больших степеней он начинает устанавливать полиномы более высокой степени. Но даже со степенью 6, принимая больше n (больше данных вместо 20, скажем 200), по-прежнему подходит квадратичный многочлен. Таким образом, мораль заключается в том, чтобы избегать переобучения, для чего это могло бы помочь как можно большему количеству точек данных.
  • Могут быть проблемы с численной стабильностью, которые я не совсем понимаю.
  • Если вам не нужен полином, вы можете получить лучшие подходы к другим функциям, например сплайны (кусочно-полиномиальные).

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

Лучше всего начать с Numerical Recipes .

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

Если у вас есть доступ к 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) 

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

Если вы хотите сопоставить (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 в задаче.

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

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

Учтите, что результаты присадок наименьших квадратов сильно зависят от распределения значений по оси x. Если значения x равномерно распределены, тогда вы получите больше ошибок на концах. Если у вас есть случай, когда вы можете выбрать значения x, и вы заботитесь о максимальном отклонении от вашей известной функции и интерполирующем многочлене, то использование полиномов Чебышева даст вам нечто, близкое к идеальному минимаксному многочлену (что очень трудно вычислить). Это подробно обсуждается в Numerical Recipes.

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

в колледже у нас была эта книга, которую я до сих пор считаю чрезвычайно полезной: Конте, де Бур; элементарный численный анализ; Mc Grow Hill. Соответствующий параграф 6.2: Фильтрация данных.
пример кода поставляется в FORTRAN, и списки также не очень читаемы, но объяснения являются глубокими и ясными в одно и то же время. вы в конечном итоге понимаете, что вы делаете, а не просто делаете это (как мой опыт Numerical 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 

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

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

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

Не забудьте выбрать подходящий для вас метод!

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

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

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

  • R tick: объединение даты и времени в один объект
  • Роллинг-срединный алгоритм в C
  • Как сделать паузу выполнения, спать, ждать X секунд в R?
  • 3D-плоскость с наименьшими квадратами
  • Вычисление среднего арифметического (один тип среднего) в Python
  • Давайте будем гением компьютера.