A robust way to fit an exponential function in loglog axis

Question

I have the following script to fit an exponential function in loglog axis:

start_exp = 21  # Start of fit
end_exp   = 40  # end of fit

# fix the parameters of  exponential Y = exponent1 *a*exp(exponent2 *b*t)
exponent1    = 3.      
exponent2    = 2.    


### The original data to perform the fit ###

bin_centers = np.array([9.31939514e-01, 1.71773809e+00, 2.58227857e+00, 3.16611120e+00,
       3.88194374e+00, 4.75962032e+00, 5.83573259e+00, 7.15514527e+00,
       8.77286667e+00, 1.07563420e+01, 1.31882653e+01, 1.61700271e+01,
       1.98259415e+01, 2.43084291e+01, 2.98043714e+01, 3.65429025e+01,
       4.48049619e+01, 5.49350072e+01, 6.73553750e+01, 8.25838892e+01,
       1.01255449e+02, 1.24148500e+02, 1.52217488e+02, 1.86632650e+02,
       2.28828807e+02, 2.80565179e+02, 3.43998732e+02, 4.21774108e+02,
       5.17133876e+02, 6.34053729e+02, 7.77408231e+02, 9.53174043e+02,
       1.16867911e+03])

norm_counts = np.array([7.47683145e-02, 2.44802278e-02, 6.40790461e-02, 4.24335863e-02,
       2.91339699e-02, 2.34427016e-02, 1.63884469e-02, 2.19590994e-02,
       1.42326906e-02, 1.64066907e-02, 1.03597068e-02, 1.03504798e-02,
       9.22668847e-03, 6.22943476e-03, 6.37955519e-03, 5.30701726e-03,
       3.61688252e-03, 3.48879003e-03, 2.46230484e-03, 1.91634352e-03,
       1.45802119e-03, 9.72116428e-04, 8.17790259e-04, 4.57538385e-04,
       3.66533958e-04, 2.17783386e-04, 1.36803586e-04, 7.28849304e-05,
       4.84365929e-05, 2.93293649e-05, 2.24564621e-05, 1.19448771e-05,
       4.87111765e-06])

def expo(x, a, b) :
    return a*np.exp(-b*x)

x_exp = bin_centers[start_exp:end_exp]
y_exp = norm_counts[start_exp:end_exp]

a_exp,b_exp = sc.optimize.curve_fit(lambda t,a,b: exponent1*a*np.exp(-b*t*(exponent2)),  x_exp,  y_exp,  p0=(0.0005, 0.0058))

expon       = (expo(bin_centers[start_exp:end_exp],a_exp[0],a_exp[1]))

plt.plot(bin_centers[start_exp:end_exp], expon,'-.',color='black',linewidth='0.8',label='${\propto} $ $e^{%.4f}$' %(-a_exp[1] ))#, 'r-',label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt))
plt.legend()

plt.plot(bin_centers, norm_counts)
plt.yscale('log')
plt.xscale('log')

The resut I a getting is prety much good.

However, I have to trie really hard to find the appropirate parameters every time. And I have to do ths fit for many curves.

Is there a way to automate ths process? A more robust algorithm that will allow me to estimate the best fit?

Answer 1

In order to clarify the notations on the graphs one define the next new variables :

One see that the quite-linearisation occurs when Y is drawn as a function of X^2 .

Thus a correct linear model might be : Y=a+b.X^2

The least mean square fitting is very simple :

As a conclusion you should use the software not for fitting the data (x,y) but for fitting the transformed data (X=log(x), Y=log(y)).

Of course instead of logarithm base 10 you could use natural logarithm on the same way.