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.])
用不动点拟合的数学正确方法是使用Lagrange multipliers。基本上,你可以修改你想要最小化的目标函数,通常是残差平方和,为每个不动点增加一个额外的参数。我还没有成功地将一个修改过的目标函数输入到scipy的一个最小化函数中。但对于多项式拟合,可以用笔和纸计算出细节,并将问题转换为线性方程组的解:
要测试它是否有效,请尝试以下操作,其中
n
是点数,d
是多项式的次数,f
是不动点的数目:当然,拟合的多项式正好穿过这些点:
如果使用
curve_fit()
,则可以使用sigma
参数为每个点赋予权重。下面的例子给出了第一点、中间点、最后一点非常小的sigma,因此拟合结果将非常接近这三点:一种简单直接的方法是利用约束最小二乘法,其中约束用较大的数字M加权,例如:
显然,这并不是一个真正的“包罗万象的银弹”解决方案,但显然,对于一个简单的例子来说,它似乎工作得很好(
for M in [0, 4, 24, 124, 624, 3124]
):并产生如下输出:
编辑:添加了“精确”解决方案:
并测试适合性:
相关问题 更多 >
编程相关推荐