Wie man eine polynomische Anpassung mit Fixpunkten macht

Question

Ich habe einige Anpassung in Python mit numpy (die kleinsten Quadrate verwendet).

Ich fragte mich, ob es eine Möglichkeit gab, die Daten anzupassen, während sie durch einige feste Punkte gezwungen wurden. Wenn nicht, gibt es eine andere Bibliothek in Python (oder eine andere Sprache, zu der ich verlinken kann - zB c)?

HINWEIS Ich weiß, es ist möglich, einen Punkt zu zwingen, durch fixiert, indem sie auf den Ursprung und zwingt den konstanten Term auf Null bewegt, wie hier festgestellt wird, wurde aber allgemein für 2 oder mehr Fixpunkten fragen:

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

Answer 1

Wenn Sie curve_fit() verwenden, können Sie sigma Argument jedem Punkt ein Gewicht zu geben. Das folgende Beispiel gibt die erste, in der Mitte letzten Punkt sehr kleines Sigma, so dass das Anpassungser wird auf diese drei Punkte sehr nahe:

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

Die mathematisch korrekte Art und Weise eine Passung mit festen tun Punkte ist Lagrange multipliers zu verwenden. Im Grunde modifizieren Sie die Zielfunktion, die Sie minimieren möchten, die normalerweise die Summe der Quadrate der Residuen ist, indem Sie für jeden Fixpunkt einen zusätzlichen Parameter hinzufügen. Es ist mir nicht gelungen, eine modifizierte Zielfunktion einem von Scipys Minimierern zuzuführen.Aber für eine Polynomfit, können Sie die Details mit Stift und Papier heraus und wandeln Ihr Problem in die Lösung eines linearen Gleichungssystems:

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] 

Um zu testen, dass es funktioniert, versuchen Sie die folgenden, wo n ist die Anzahl der Punkte, d der Grad des Polynoms und f die Anzahl der Fixpunkte:

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

und natürlich das angepasste Polynom geht Exac tly durch die Punkte:

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

Answer 3

Eine einfache und unkomplizierte Art und Weise ist der kleinsten Quadrate zu verwenden, beschränkt, wo Einschränkungen mit einer ziemlich großen Zahl M gewichtet werden, wie:

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

Offensichtlich ist dies nicht wirklich ein alle inklusive Königsweg Lösung, aber anscheinend scheint es vernünftig, gut mit einem einfachen Beispiel (for M in [0, 4, 24, 124, 624, 3124]) zu arbeiten:

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

und wie Erzeugen von Ausgangs:

Edit: Added 'genau' Lösung:

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

und Prüfung der Passform:

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