如何做固定點的多項式擬合

Question

我一直在使用numpy（使用最小二乘法）在python中進行一些擬合。

我想知道是否有得到它，以適應數據，同時迫使它通過一些固定點的一種方式？如果不是有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 TLY通過點：

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

Answer 3

一個簡單且直接的方式是利用約束，其中約束加權一個相當大的數目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, 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() 

，併產生輸出，如：

編輯：增加了 '精確' 的解決方案：

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