【问题标题】:Scipy curve_fit fails with many-parameter fit. Is there a way to improve results?Scipy curve_fit 因多参数拟合而失败。有没有办法改善结果?
【发布时间】:2019-08-24 13:51:06
【问题描述】:

我正在尝试使用 scipy curve_fit 方法来拟合振荡数据。不幸的是,我有 8 个参数,并且无法减小维度(或者我看不到方法)。 这是适合的函数:

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

def cosFitForCFF(x,c0, c1, b0, b1, b2, b3, b4, b5):
    """
    Function for fit.
    c0, c1, b0 irrelevant parameters
    b1, b2, b3, b4, b5 are the important parameters
    """
    return c0 + c1*np.cos(b0+b1*x+b2*x**2+b3*x**3+b4*x**4+b5*x**5)

前3个参数无关紧要,我需要后5个进行计算。

我有一个函数可以读取输入并使用我可以提供的所有选项(初始参数、边界)进行拟合:

def CFFMethod(initSpectrumX, initSpectrumY, referenceArmY, sampleArmY, 
p0=[1, 1, 1, 1, 1, 1, 1, 1], referencePoint = 2.5):
    """
    Phase modulated cosine function fit method. p0 is the array containing inital parameters for fitting.
    referencePoint is some point in initSpectrumX
    """

    #best bounds I can  provide
    bounds=((-1, -1, -1, -np.inf, -np.inf, -np.inf, -np.inf, -np.inf), 
                (1, 1, 1, np.inf, np.inf, np.inf, np.inf, np.inf))

    #reading inputs
    if len(initSpectrumY) > 0 and len(referenceArmY) > 0 and len(sampleArmY) > 0:
        Ydata = (initSpectrumY-referenceArmY-sampleArmY)/(2*np.sqrt(referenceArmY*sampleArmY))
        Ydata = np.asarray(Ydata)
    elif len(initSpectrumY) == 0:
        raise ValueError('Please load the spectrum!\n')
    elif len(referenceArmY) == 0 or len(sampleArmY) == 0:
        Ydata = np.asarray(initSpectrumY)

    Xdata = initSpectrumX-referencePoint
    Xdata = np.asarray(Xdata)
    #end of reading inputs

    #fitting
    try:
        popt, pcov = curve_fit(cosFitForCFF, Xdata, Ydata, p0, maxfev = 5000, bounds = bounds)
        #plot
        fig1 = plt.figure()
        fig1.canvas.set_window_title('Cosine function fit method')
        plt.plot(Xdata, Ydata,'r-',label = 'dataset')
        plt.plot(Xdata, cosFitForCFF(Xdata, *popt),'k*', label = 'fitted')
        plt.legend()
        plt.grid()
        plt.show()
        return popt
    except Exception as e:
        print(e)

即使我使用该函数生成合成数据,我也会得到错误的结果。

xs = np.linspace(2.5, 3, 1000)
ys = cosFitForCFF(xs, 0, 1, 0, 0, 50, 0, 0, 0)
params = [0, 1, 0, 0, 50, 0, 0, 0] #exact same that was generated
reference = 2.7 # some point in the data, irrelevant

result = CFFMethod(xs, ys, [],[], p0 = params, referencePoint = reference)
print(result)
#outputs to: 
#[-5.12643151e-01  1.00000000e+00  9.99999995e-01  2.05339754e-01
 # 1.01356470e+01 -3.83963354e+01 -3.53998314e+02  1.33074662e+03]

我知道curve_fit 参数太多,所以我需要将maxfev 设置得更高。

而且这甚至与现实世界中的数据集相差甚远,可能会有噪音等。

我做错了什么吗?也许我应该寻找另一种算法?拟合函数必须是我上面定义的函数,因为这样分散系数(我需要找到)与 b1,b2.. 我非常感谢对代码的任何帮助/改进。

编辑:

禁用参考点后,它完全适合,但只适合生成的数据集。拟合真实数据集仍然会导致错误的结果。 这是更新后的功能:

def CFFMethod(initSpectrumX, initSpectrumY, referenceArmY, sampleArmY, p0=[1, 1, 1, 1, 1, 1, 1, 1]):
    """
    Phase modulated cosine function fit method. p0 is the array containing inital parameters for fitting.
    referencePoint is some point in initSpectrumX
    """

    #provided bounds
    bounds=((-1, -1, -1, -np.inf, -np.inf, -np.inf, -np.inf, -np.inf), (1, 1, 1, np.inf, np.inf, np.inf, np.inf, np.inf))

    #reading inputs
    if len(initSpectrumY) > 0 and len(referenceArmY) > 0 and len(sampleArmY) > 0:
        Ydata = (initSpectrumY-referenceArmY-sampleArmY)/(2*np.sqrt(referenceArmY*sampleArmY))
        Ydata = np.asarray(Ydata)
    elif len(initSpectrumY) == 0:
        raise ValueError('Please load the spectrum!\n')
    elif len(referenceArmY) == 0 or len(sampleArmY) == 0:
        Ydata = np.asarray(initSpectrumY)

    Xdata = np.asarray(initSpectrumX)
    #end of reading inputs

    #fitting
    try:
        popt, pcov = curve_fit(cosFitForCFF, Xdata, Ydata, p0, maxfev = 5000, bounds = bounds)
        #plot
        fig1 = plt.figure()
        fig1.canvas.set_window_title('Cosine function fit method')
        plt.plot(Xdata, Ydata,'r-',label = 'dataset')
        plt.plot(Xdata, cosFitForCFF(Xdata, *popt),'k*', label = 'fitted')
        plt.legend()
        plt.grid()
        plt.show()
        return popt
    except Exception as e:
        print(e)

我在那里提供数据集生成器:

# GENERATOR FUNCTIONS 
def _ensureInput(start, stop, center, resolution):
    if start >= stop:
        raise ValueError('start value must be less than stop')
    if center < start or center > stop:
        raise ValueError('center must be between start and  stop')  
    if resolution > (stop-start):
        raise ValueError('resolution is too big')
    else:
        pass

def _disp(x ,GD=0, GDD=0, TOD=0, FOD=0, QOD=0):
    return x*GD+(GDD/2)*x**2+(TOD/6)*x**3+(FOD/24)*x**4+(QOD/120)*x**5 


def generatorFreq(start, stop, center ,delay, GD=0, GDD=0, TOD=0, FOD=0, QOD=0, resolution = 0.1,
                     delimiter =',',pulseWidth = 0.02, includeArms = False):
    _ensureInput(start, stop, center, resolution)
    c = 299.793 
    deltaL = delay 
    omega0 = center 
    window = pulseWidth
    lamend = (2*np.pi*c)/start
    lamstart = (2*np.pi*c)/stop
    lam = np.arange(lamstart, lamend+resolution, resolution) 
    omega = (2*np.pi*c)/lam 
    relom = omega-omega0
    i1 = np.exp(-(relom)**2/(window))
    i2 = np.exp(-(relom)**2/(window))
    i = i1 + i2 + 2*np.cos(_disp(relom,GD=GD, GDD= GDD, TOD=TOD, FOD=FOD, QOD=QOD)+2*deltaL*omega/c)*np.sqrt(i1*i2) 
    if includeArms:
        return omega, i, i1, i2
    else:
        return omega, i, np.array([]), np.array([])

## using the generator to make a dataset
a,b,c,d = generatorFreq(2 ,3, 2.5, 0, GD = 0, GDD = 200, TOD = 4000, FOD = 0, QOD = 0, resolution = 0.1, delimiter = ',', pulseWidth = 0.02, includeArms = True)
# fit to data 
result = CFFMethod(a, b, c,d , p0 = [0, 1, 0, 0, 200, 4000, 0, 0])

您现在可以看到,通过复制粘贴,curve_fit 无法产生良好的效果。

【问题讨论】:

    标签: python numpy scipy curve-fitting


    【解决方案1】:

    我有一个建议作为摆脱困境的可能途径。拟合数据的一小部分应该比拟合整个数据集更容易——当这种方法有效时,这些参数可以用作更大数据子集的初始参数估计,依此类推。这是您的代码修改如下以使用前 50 个数据点:

    #fitting
    Xdata = Xdata[:50]
    Ydata = Ydata[:50]
    try:
        popt, pcov = curve_fit(cosFitForCFF, Xdata, Ydata, p0, maxfev = 5000, bounds = bounds)
    

    结果如下:

    【讨论】:

    • 我一定会尝试并发布结果。
    猜你喜欢
    • 2012-11-02
    • 1970-01-01
    • 1970-01-01
    • 2014-07-12
    • 2019-06-09
    • 1970-01-01
    • 2016-02-25
    • 2016-11-29
    • 2020-04-04
    相关资源
    最近更新 更多