2

我正在尝试使用以下形式的抛物线形状拟合数据: Y(x) = a(1 – (x/b)^2) for |x| < b ....... Y(x) = 0 否则 我添加了两个额外的自由度,因此抛物线可以水平和垂直移动。

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

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

欢迎所有建议。

谢谢

4

1 回答 1

2

我可能会这样做:

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

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

于 2013-04-25T09:42:15.157 回答