python - 将抛物线形状拟合到数据(是否可以将逻辑语句与 SciPy CurveFit 一起使用？)

Question

我正在尝试将数据拟合为以下形式的抛物线形状: Y(x) = a(1 – (x/b)^2) 对于 |x| 我添加了两个额外的自由度，以便抛物线可以水平和垂直移动。

from scipy.optimize import curve_fit

def parabolic(self, t, *p): 
    a, b, c, d = p        
    if abs(t) > b:
        return 0
    else:
        return a*(1-(((t-c)/b)**2)) + d

我正在尝试使用 scipy.optimize 中的 curve_fit 并创建了一个抛物线函数，然后尝试使用以下方法进行拟合:

coeff, cov = curve_fit(parabolic, data[:,0], data[:,1], p0) #Return co-effs for fit and      covariance, p0 = user given seed fitting parameters

尽管与在我的抛物线拟合定义中使用逻辑相关，但该脚本返回一个错误。当我删除逻辑语句并仅拟合整个数据(这是中心为抛物线形状的背景噪声级别)时，结果很糟糕，因为拟合试图包含背景噪声级别。

欢迎所有建议。

谢谢

Answer 1

我可能会这样做:

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 之间)并找到最佳的简化卡方值。

但是，我还注意到数据并不显示为抛物线。最初，当我看到完整的数据图时，我想到了“高斯”，但事实也并非如此。它看起来几乎像一个棚车功能。 如果你有一个很好的理论模型，它是抛物线，那么要么数据不正确，要么模型可能不正确。如果您只是在寻找描述性函数，也可以尝试一些其他函数。