Curve fitting for a complex function with 4 parameters

Question

I try to fit the coefficients in my curve equation to the experimental data using curve_fit from scipy.optimize and the code works fine with various equations (straight line, quadratic, cubic function), but I get an error when I try to fit my equation to a curve:

RuntimeError: Optimal parameters not found: Number of calls to function has reached maxfev = 10000.`

the equation looks like this coefficients a, b, c, d need to be matched

my code:

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

x = np.array([1.8652, 1.8966, 1.9125, 1.9230, 1.9313, 1.9359, 1.9435, 1.9523, 1.9585, 1.9649])
y = np.array([-0.0597, -0.0596, -0.0607, -0.0617, -0.0628, -0.0637, -0.0646, -0.0658, -0.0668, -0.0677])

fig=plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(x, y, 'ko')
plt.xlabel('x') 
plt.ylabel('y')
plt.scatter(x, y)
plt.grid() 
plt.show()


def f(x, a, b, c, d):
    return (np.log10((((2*(((a*b*(m.exp(1)**((-0.5)*38.94336875*x)))*(1-(m.exp(1)**(2*38.94336875*x))))/((a+b)+(c+d)*(m.exp(1)**(38.94336875*x)))))*96485)/(0.000035))))
args, _ = curve_fit(f, x, y, maxfev=10000)

print('args (a, b, c, d):', args)


args, _ = curve_fit(f, x, y, maxfev=10000) 
a, b, c, d = args[0], args[1], args[2], args[3]
y_fit1 = (np.log10((((2*(((a*b*(m.exp(1)**((-0.5)*38.94336875*x)))*(1-(m.exp(1)**(2*38.94336875*x))))/((a+b)+(c+d)*(m.exp(1)**(38.94336875*x)))))*96485)/(0.000035))))
plt.plot(x, y, 'bo', label="y - experiment") 
plt.plot(x, y_fit1, label="function")
plt.xlabel('x') 
plt.ylabel('y') 
plt.legend(loc = 'best', fancybox = True, shadow = True) 
plt.grid(True) 
plt.show()

I tried increasing the number of iterations to 10,000,000 but I get the same error. I also got rid of the logarithm in the equation, but still got the same error.

Answer 1

First try not to increse maxfev as reaching the iteration limit typically tells you that something is going wrong with your fit. Typical cases are that either your fit function cannot reproduce your lineshape at all or the starting values for the fit does not allow the fit to reach a nice fitting of the curve, either because they shoot in the wrong direction or get trapped in a local minimum.

Normaly you play around with good satring values and give them as a starting point to the fit:

myp0 = [a0, b0, c0, d0] # put good sarting values here
args, _ = curve_fit(f, x, y, p0 = myp0)

Looking at your formula, the first thing you should notice is, that c and d are not inependent! So first thing is to get rid of d.

Notice that you do several calculation in your fit function, you can pre compute it, use np.exp and simplify a bit:

def f(x, a, b, c):
    return np.log10(2756714285.714286*a*b*
                    np.exp(-19.471684375*x)*
                    (1-np.exp(77.8867375*x))/
                    ((a+b)+(c*np.exp(38.94336875*x))))

Now, note that due to the log, you can extract the constant in the beginning out of the log:

def f(x, a, b, c):
    return 9.44039 * np.log10(a*b* np.exp(-19.471684375*x)*
                              (1-np.exp(77.8867375*x))/
                              ((a+b)+(c*np.exp(38.94336875*x))))

This helps you to work out some good starting values. As we know the first data point is (1.8652, -0.0597), this tells you that at this point, the log should be approximatly -0,0063, so the inner part of the log10 must be approx. 0.9855.

With that we know that 0.9855 = a*b*1.6867e-16 * -1.2355e+63/((a+b)+(c*3.515e+31))

Here you already see, that the exp functions are exploding and everything is dominated by the -1.2355e+63 part. So it seems that something is very fishy with your fitting function and you should double check if it indeed makes sense for your dataset. If so continue to find good starting values for a, b, c.