如何对固定点进行多项式拟合 [英] How to do a polynomial fit with fixed points

问题描述

我一直在使用 numpy(使用最小二乘法)在 python 中进行一些拟合.

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

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

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

解决方案

数学上正确的固定点拟合方法是使用 拉格朗日乘数.基本上，你修改你想要最小化的目标函数，它通常是残差的平方和，为每个固定点添加一个额外的参数.我没有成功地将修改后的目标函数提供给 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, 轴=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/2mat[n + 1:, :n + 1] = xf_n.Tmat[n + 1:, n + 1:] = 0vec[:n + 1] = yx_nvec[n + 1:] = yf参数 = np.linalg.solve(mat, vec)返回参数[:n + 1]

要测试它是否有效，请尝试以下操作，其中 n 是点数，d 多项式的次数和 f固定点数:

n, d, f = 50, 8, 3x = 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.5yf = np.random.uniform(np.min(y), np.max(y), size=(f,))参数 = 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数组([1.03101335，2.94879161，2.87288739])>>>聚(xf)数组([1.03101335，2.94879161，2.87288739])

I have been doing some fitting in python using numpy (which uses least squares).

I was wondering if there was a way of getting it to fit data while forcing it through some fixed points? If not is there another library in python (or another language i can link to - eg c)?

NOTE I know it's possible to force through one fixed point by moving it to the origin and forcing the constant term to zero, as is noted here, but was wondering more generally for 2 or more fixed points :

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

解决方案

The mathematically correct way of doing a fit with fixed points is to use Lagrange multipliers. Basically, you modify the objective function you want to minimize, which is normally the sum of squares of the residuals, adding an extra parameter for every fixed point. I have not succeeded in feeding a modified objective function to one of scipy's minimizers. But for a polynomial fit, you can figure out the details with pen and paper and convert your problem into the solution of a linear system of equations:

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]

To test that it works, try the following, where n is the number of points, d the degree of the polynomial and f the number of fixed points:

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

And of course the fitted polynomial goes exactly through the points:

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

这篇关于如何对固定点进行多项式拟合的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！