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

Question

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

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

Примечание Я знаю, что можно заставить через одну неподвижную точку, перемещая его происхождения и принуждать постоянный член к нулю, как отмечено здесь, но было интересно, в более общем плане для 2-х или более неподвижных точек:

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

Answer 1

Если вы используете 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") 

Answer 2

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

и, конечно, насаженным полином идет exac ждение через точку:

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

Answer 3

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

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, s.t. Cx= d, 
    based on the idea of weighting constraints with a largish number M.""" 
    return solve(dot(A.T, A)+ M* dot(C.T, C), dot(A.T, b)+ M* dot(C.T, 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() 

и производить выход, как:

Edit: Добавлен 'точное' решение:

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(A.T, A), dot(A.T, b)) 

def clsq(A, b, C, d): 
    """An 'exact' constrained least squared solution of Ax= b, s.t. Cx= d""" 
    p= C.shape[0] 
    Q, R= qr(C.T) 
    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.])