Как сделать многочлен с фиксированными точками

Я делал некоторые фитинги в python, используя numpy (который использует наименьшие квадраты).

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

ПРИМЕЧАНИЕ. Я знаю, что можно принудительно выполнить одну неподвижную точку, переместив ее в начало координат и заставляя постоянный член равняться нулю, как это отмечено здесь, но задавалось более общим вопросом для 2 или более неподвижных точек:

http://www.physicsforums.com/showthread.php?t=523360

  • Scipy odeint дает предупреждение lsoda
  • Как выполнить двухдисковый односторонний t-тест с numpy / scipy
  • Лог-вычисления в Python
  • Причина, почему numpy rollaxis настолько запутанна?
  • Вычисления свертки в Numpy / Scipy
  • Ускорение scipy griddata для множественных интерполяций между двумя нерегулярными сетками
  • Как получить аргументы параметров из замороженного дистрибутива spicy.stats?
  • Создание массива в numpy / scipy путем итерации в Python?
  • 3 Solutions collect form web for “Как сделать многочлен с фиксированными точками”

    Математически правильный способ делать фитинг с неподвижными точками – использовать множители Лагранжа . В основном, вы изменяете целевую функцию, которую хотите минимизировать, что обычно представляет собой сумму квадратов остатков, добавляя дополнительный параметр для каждой фиксированной точки. Мне не удалось подавать измененную целевую функцию на один из минимизаторов Scipy. Но для полиномиального соответствия вы можете выяснить детали с помощью пера и бумаги и преобразовать свою проблему в решение линейной системы уравнений:

    def polyfit_with_fixed_points(n, x, y, xf, yf) : mat = np.empty((n + 1 + len(xf),) * 2) vec = np.empty((n + 1 + len(xf),)) x_n = x**np.arange(2 * n + 1)[:, None] yx_n = np.sum(x_n[:n + 1] * y, axis=1) x_n = np.sum(x_n, axis=1) idx = np.arange(n + 1) + np.arange(n + 1)[:, None] mat[:n + 1, :n + 1] = np.take(x_n, idx) xf_n = xf**np.arange(n + 1)[:, None] mat[:n + 1, n + 1:] = xf_n / 2 mat[n + 1:, :n + 1] = xf_n.T mat[n + 1:, n + 1:] = 0 vec[:n + 1] = yx_n vec[n + 1:] = yf params = np.linalg.solve(mat, vec) return params[:n + 1] 

    Чтобы проверить, что это работает, попробуйте следующее, где n – количество точек, d степень полинома и f – число неподвижных точек:

     n, d, f = 50, 8, 3 x = np.random.rand(n) xf = np.random.rand(f) poly = np.polynomial.Polynomial(np.random.rand(d + 1)) y = poly(x) + np.random.rand(n) - 0.5 yf = np.random.uniform(np.min(y), np.max(y), size=(f,)) params = polyfit_with_fixed(d, x , y, xf, yf) poly = np.polynomial.Polynomial(params) xx = np.linspace(0, 1, 1000) plt.plot(x, y, 'bo') plt.plot(xf, yf, 'ro') plt.plot(xx, poly(xx), '-') plt.show() 

    введите описание изображения здесь

    И, конечно, подогнанный полином проходит точно через точки:

     >>> yf array([ 1.03101335, 2.94879161, 2.87288739]) >>> poly(xf) array([ 1.03101335, 2.94879161, 2.87288739]) 

    Если вы используете curve_fit() , вы можете использовать аргумент sigma чтобы дать каждой точке вес. Следующий пример дает первую, среднюю, последнюю точку очень маленькую сигму, поэтому результат подбора будет очень близок к этим трем пунктам:

     N = 20 x = np.linspace(0, 2, N) np.random.seed(1) noise = np.random.randn(N)*0.2 sigma =np.ones(N) sigma[[0, N//2, -1]] = 0.01 pr = (-2, 3, 0, 1) y = 1+3.0*x**2-2*x**3+0.3*x**4 + noise def f(x, *p): return np.poly1d(p)(x) p1, _ = optimize.curve_fit(f, x, y, (0, 0, 0, 0, 0), sigma=sigma) p2, _ = optimize.curve_fit(f, x, y, (0, 0, 0, 0, 0)) x2 = np.linspace(0, 2, 100) y2 = np.poly1d(p)(x2) plot(x, y, "o") plot(x2, f(x2, *p1), "r", label=u"fix three points") plot(x2, f(x2, *p2), "b", label=u"no fix") legend(loc="best") 

    введите описание изображения здесь

    Одним простым и простым способом является использование ограниченных наименьших квадратов, где ограничения взвешиваются с большим числом M, например:

     from numpy import dot from numpy.linalg import solve from numpy.polynomial.polynomial import Polynomial as P, polyvander as V def clsq(A, b, C, d, M= 1e5): """A simple constrained least squared solution of Ax= b, st Cx= d, based on the idea of weighting constraints with a largish number M.""" return solve(dot(AT, A)+ M* dot(CT, C), dot(AT, b)+ M* dot(CT, d)) def cpf(x, y, x_c, y_c, n, M= 1e5): """Constrained polynomial fit based on clsq solution.""" return P(clsq(V(x, n), y, V(x_c, n), y_c, M)) 

    Очевидно, что это не совсем всеобъемлющее решение для серебряных пулей , но, по-видимому, оно хорошо работает с простым примером ( for M in [0, 4, 24, 124, 624, 3124] ):

     In []: x= linspace(-6, 6, 23) In []: y= sin(x)+ 4e-1* rand(len(x))- 2e-1 In []: x_f, y_f= linspace(-(3./ 2)* pi, (3./ 2)* pi, 4), array([1, -1, 1, -1]) In []: n, x_s= 5, linspace(-6, 6, 123) In []: plot(x, y, 'bo', x_f, y_f, 'bs', x_s, sin(x_s), 'b--') Out[]: <snip> In []: for M in 5** (arange(6))- 1: ....: plot(x_s, cpf(x, y, x_f, y_f, n, M)(x_s)) ....: Out[]: <snip> In []: ylim([-1.5, 1.5]) Out[]: <snip> In []: show() 

    и производство продукции, как: подходит с прогрессивным увеличением M

    Изменить: Добавлено «точное» решение:

     from numpy import dot from numpy.linalg import solve from numpy.polynomial.polynomial import Polynomial as P, polyvander as V from scipy.linalg import qr def solve_ns(A, b): return solve(dot(AT, A), dot(AT, b)) def clsq(A, b, C, d): """An 'exact' constrained least squared solution of Ax= b, st Cx= d""" p= C.shape[0] Q, R= qr(CT) xr, AQ= solve(R[:p].T, d), dot(A, Q) xaq= solve_ns(AQ[:, p:], b- dot(AQ[:, :p], xr)) return dot(Q[:, :p], xr)+ dot(Q[:, p:], xaq) def cpf(x, y, x_c, y_c, n): """Constrained polynomial fit based on clsq solution.""" return P(clsq(V(x, n), y, V(x_c, n), y_c)) 

    и тестирование соответствия:

     In []: x= linspace(-6, 6, 23) In []: y= sin(x)+ 4e-1* rand(len(x))- 2e-1 In []: x_f, y_f= linspace(-(3./ 2)* pi, (3./ 2)* pi, 4), array([1, -1, 1, -1]) In []: n, x_s= 5, linspace(-6, 6, 123) In []: p= cpf(x, y, x_f, y_f, n) In []: p(x_f) Out[]: array([ 1., -1., 1., -1.]) 
    Python - лучший язык программирования в мире.