Scipy curve_fit fails with many-parameter fit. Is there a way to improve results?

Question

I'm trying to use scipy curve_fit method to fit to an oscillating data. Unfortunately I have 8 parameters, and the dimension can't be reduced (or I don't see a way). This is the function to fit:

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)

The first 3 parameters are irrelevant, I need the last 5 to proceed calculations.

I have a function which reads the input and does the fitting with all the options I can provide(initial parameters, bounds):

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)

Even if I generate synthetic data with the function I get wrong results.

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]

I know that curve_fit is struggling with too many parameters, that's why I needed to set maxfev higher.

And this is not even close to the real-world-like dataset, which can be noisy, etc.

Is there anything I'm doing wrong? Maybe I should search for another algorithm? The fitted function must be in that from I defined above, because this way the dispersion coefficients(which I need to find) are related to b1,b2.. I really appreciate any help/improvement on the code.

EDIT:

After disabling the referencePoint it fits perfectly, but only the generated datasets. Fitting to real dataset still leads to wrong results. Here is the updated function:

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)

I provide dataset generator there:

# 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])

You can see now with copy-pasting, curve_fit fails to produce good results.

Answer 1

I have a suggestion as a possible path out of the difficulty. It should be easier to fit a small subset of the data than to fit the entire data set - and when that works, those parameters can be used as the initial parameter estimates for a larger data subset, and so on. Here is your code modified as follows to use the first 50 data points:

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

with the following result: