Jak zrobić wielomianowy pas z ustalonymi punktami

Question

Zrobiłem pewne dopasowanie w python za pomocą numpy (który używa najmniejszych kwadratów).

Zastanawiam się, czy istnieje sposób na uzyskanie go dopasować dane, podczas gdy zmuszając go przez niektóre stałe punkty? Jeśli nie, to istnieje inna biblioteka w Pythonie (lub innym języku, do którego mogę się podłączyć - np. C)?

UWAGA wiem, że to możliwe, aby przeforsować jednego stałego punktu przesuwając go do pochodzenia i zmuszając stały termin do zera, jak jest tutaj zauważyć, ale zastanawiałem się bardziej ogólnie na 2 lub więcej punktów stałych:

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

Answer 1

Jeśli używasz curve_fit(), można użyć argumentu sigma dać każdym punkcie wagi. Poniższy przykład daje pierwszy, w środku, ostatni punkt bardzo małe sigma, więc wynik montaż będzie bardzo zbliżone do tych trzech punktach:

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

matematycznie poprawny sposób robi szału ze stałym punktami jest użycie Lagrange multipliers. Zasadniczo modyfikujesz funkcję celu, którą chcesz zminimalizować, która jest zwykle sumą kwadratów reszt, dodając dodatkowy parametr dla każdego ustalonego punktu. Nie udało mi się wprowadzić zmodyfikowanej funkcji celu do jednego z minimizerów scipy.Ale dla wielomianu pasuje, można dowiedzieć się szczegółów z pióra i papieru i przekształcić swój problem do rozwiązania układu równań liniowych:

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] 

Aby sprawdzić, czy to działa, spróbuj wykonać następujące, gdzie n jest liczba punktów, d stopień wielomianu i f liczba stałych punktów:

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() 

i oczywiście wyposażona wielomian idzie EXAC nio przez punkty:

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

Answer 3

Jeden prosty i bezpośredni sposób jest wykorzystanie ograniczone najmniejszych kwadratów, gdzie ograniczenia są ważone z wieloma sporej M, jak:

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

Oczywiście to nie jest tak naprawdę wszystko włącznie panaceum rozwiązanie, ale pozornie wydaje się działać dobrze z racjonalną prosty przykład (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() 

i wyjście produkcji jak:

Edit: 'dokładna' dodano roztwór:

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

i testowania dopasowanie:

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