Спасибо за все ответы. Вот еще одна попытка их обобщения. Простите, если я скажу слишком много «очевидных» вещей: раньше я ничего не знал о наименьших квадратах, так что все было для меня новым.
НЕ полиномиальная интерполяция
Полиномиальная интерполяция соответствует полиному степени 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), чтобы минимизировать (сохранить «наименьшее») «сумму квадратов невязок», а именно: 1035 *
S = & sum; i (a + bx i + cy i + dx i 2 + эй я 2 + фх я у я - Z я ) 2
Теория
Важная идея заключается в том, что если вы посмотрите на 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 у я = 0
и аналогичные уравнения для a, b, c, d, e.
Обратите внимание, что это просто линейные уравнения в a & hellip; f. Таким образом, мы можем решить их с помощью Гауссова исключения или любого из обычных методов .
Это все еще называется "линейным наименьшим квадратом", потому что хотя функция, которую мы хотели, была квадратичным полиномом, она все еще линейна в параметрах (a, b, c, d, e, f) , Обратите внимание, что то же самое работает, когда мы хотим, чтобы p (x, y) была любой «линейной комбинацией» из произвольных функций f j вместо просто полиномиальной (= "линейной комбинации мономов ").
Код
Для одномерного случая (когда есть только переменная x - f j являются мономами x j ), есть Numpy's 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 i ). (Теория состоит в том, что он находит псевдообратную Мура-Пенроуза of A.) В нашем приведенном выше примере (x i , y i , Z i ), подгонка полинома означает, что f j являются мономами x () y () . Следующее находит лучший квадратик (или лучший многочлен любой другой степени, если вы измените строку «степень = 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 может сделать это. Числовые рецепты обсуждает методы, чтобы реализовать его самостоятельно (на 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 * * + одна тысяча сто семьдесят одна х + у * 2 * одна тысяча сто семьдесят два * +1173 *). Поэтому мораль заключается в том, что важно тщательно выбирать точки, чтобы получить «правильный» многочлен. (Если вы можете выбрать, вы должны выбрать чебышевские узлы для полиномиальной интерполяции; не уверены, что то же самое верно и для наименьших квадратов.)
- Переоснащение : полиномы более высокой степени всегда могут лучше соответствовать данным. Если вы измените
degree на 3, 4 или 5, он по-прежнему в основном распознает тот же квадратичный полином (коэффициенты равны 0 для членов более высокой степени), но для больших степеней он начинает подгонять полиномы более высокой степени. Но даже с 6-й степенью взятие большего n (больше точек данных вместо 20, скажем, 200) по-прежнему соответствует квадратичному полиному. Таким образом, мораль состоит в том, чтобы избежать переобучения, для чего это может помочь собрать как можно больше данных.
- Могут быть проблемы числовой стабильности Я не до конца понимаю.
- Если вам не нужен многочлен, вы можете получить лучшее согласование с другими видами функций, например, сплайны (кусочно-полиномы).