我可能会这样做:
from __future__ import division
from __future__ import print_function
from scipy.optimize import curve_fit
import numpy as np
def parabola(t, *p):
a, b, c, d = p
y = np.zeros(t.shape)
indices = np.abs(t) < b
y[indices] = (a*(1-(((t[indices]-c)/b)**2)) + d)
return y
p0 = [1, 2, 3, 4]
x = np.linspace(-10, 10, 20)
y = parabola(x, *p0)
coeff, cov = curve_fit(parabola, x, y, p0)
print(coeff)
更新
使用提供的数据(请参阅评论中的链接),事情开始变得更加清晰:
这是关于其行为的棘手数据。许多“零”数据无助于拟合函数,因为随着零数据的增多,实际函数的权重变得越来越小。通常,最好的办法是减少一些并专注于相关数据(一旦你有了这些,你可以尝试用你最好的参数来拟合所有数据)。
标准不正确:它以零为中心,而不是围绕数据的峰值 ( c
)。
b
在函数和索引标准中都有参数c
会使事情变得更难:函数的行为越来越不线性。我通过使用固定标准找到了最初的最佳参数,这是下面注释掉的行。
提供良好的启动参数。[1, 2, 3, 4]
非常通用,在非线性最小二乘法的情况下,很难拟合。
因此,考虑到以上所有因素,我想出了这个:
from __future__ import division
from __future__ import print_function
from scipy.optimize import curve_fit
import numpy as np
from matplotlib import pyplot as plt
def parabola(t, *p):
a, b, c, d = p
y = np.zeros(t.shape)
# The indices criterion was first fixed to values that appeared reasonably;
# otherwise the fit would completely fail.
# Once decent parameters were found, I replaced 28 and 0.3 with the center `c`
# and the width `b`.
#indices = np.abs(t-28) < 0.3
indices = np.abs(t-c) < b
y[indices] = (a*(1-(((t[indices]-c)/b)**2)) + d)
return y
out = np.loadtxt('data.dat')
# Limit the data to only the interesting part of the data
# Once we have the fit correct, we can always attempt a fit to all data with
# good starting parameters
xdata = out[...,0][450:550]
ydata = out[...,1][450:550]
# These starting parameters are either from trial fitting, or from theory
p0 = [2, 0.2, 28, 6.6]
coeff, cov = curve_fit(parabola, xdata, ydata, p0)
plt.plot(xdata, ydata, '.')
xfit = np.linspace(min(xdata), max(xdata))
yfit = parabola(xfit, *coeff)
plt.plot(xfit, yfit, '-')
plt.show()
请注意,生成的b
参数仍然表示不合适。我猜这些数据和这个函数的组合很棘手。一种选择是迭代各种合理的值b
(例如,在 0.2 和 0.3 之间)并找到最佳的约简卡方。
但是,我也注意到数据不会显示为抛物线。最初,当我看到完整的数据图片时,我以为是“高斯”,但也不是这样。它看起来几乎像一个棚车功能。如果你有一个很好的理论模型表明它是抛物线,那么要么数据不正确,要么模型可能不正确。如果您只是在寻找描述性功能,请尝试其他一些功能。