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

Математически правильный способ делать фитинг с неподвижными точками – использовать множители Лагранжа . В основном, вы изменяете целевую функцию, которую хотите минимизировать, что обычно представляет собой сумму квадратов остатков, добавляя дополнительный параметр для каждой фиксированной точки. Мне не удалось подавать измененную целевую функцию на один из минимизаторов 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() 

и производство продукции, как:

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

 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.])