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

问题描述

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

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

我想知道是否存在一种使数据适合数据同时又迫使它们经过一些固定点的方法?如果没有，是否存在python中的另一个库(或我可以链接到的另一种语言-例如c)?

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

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

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

推荐答案

在数学上正确地进行定点拟合的方法是使用拉格朗日乘数.基本上，您可以修改要最小化的目标函数，该目标函数通常是残差的平方和，并为每个固定点添加一个额外的参数.我没有成功地将修改后的目标函数提供给scipy的最小化器之一.但是对于多项式拟合，您可以用纸和笔找出细节，然后将问题转换为线性方程组的解:

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]

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

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屋！