如何用固定点进行多项式拟合答案

【问题标题】：How to do a polynomial fit with fixed points如何用固定点进行多项式拟合
【发布时间】：2013-02-17 22:09:17
【问题描述】：

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

我想知道是否有办法让它适应数据，同时强制它通过一些固定点？如果没有，python 中是否还有另一个库（或我可以链接到的另一种语言 - 例如 c）？

注意我知道可以通过将其移动到原点并将常数项强制为零来强制通过一个固定点，如此处所述，但更普遍地想知道 2 个或更多固定点点：

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

【问题讨论】：

不确定插值在这里会有所帮助 - 如果我的多项式模型没有通过正确的点，那么任何插值都无法实现。
固定点是指 x 和 y 都固定，对吗？您可以在修复这些点时使用en.wikipedia.org/wiki/Lagrange_polynomial 进行插值。
谢谢...看起来很有趣。目前，我已经完成了一项工作，我将固定点添加到数据中，并加权它们的负载比其他负载更多......似乎工作正常......
@rcompton 拉格朗日多项式非常适合拟合精确通过某些点的多项式，但您建议如何使用它们来近似拟合其他点？
@Jaime 我将问题解释为“我如何将数据与您通过给定点的约束相匹配”。他们并没有真正做与您的答案相同的事情。但我想如果你要使用拉格朗日多项式，你不妨完全拟合所有点。

标签： python numpy scipy

【解决方案1】：

用固定点进行拟合的数学上正确的方法是使用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_points(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])

【讨论】：

如果使用此处建议的 poly1d() 构造函数：docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html，则参数切片的顺序与预期相反。简单的解决方法是将return params[:n + 1] 更改为return params[:n + 1][::-1]。
感谢您的代码，它非常整洁。更进一步，您认为可以添加条件“多项式下方的总面积必须为A”吗？其中A 是一个常数。

【解决方案2】：

如果您使用curve_fit()，您可以使用sigma 参数为每个点赋予权重。下面的例子给出了第一个、中间、最后一个点非常小的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")

【讨论】：

【解决方案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.])

【讨论】：