Search code examples
pythonmatplotlibscipycurve-fittingdata-fitting

How to choose the right initial fitting parameters for I'm guessing a double exponential decay?

I'm stuck with failure tentative in trying to guess the initial fitting parameters to give to scipy.optimize.curve_fit, in order to fit my data (which has logarithmic y-axis). I'm guessing it's a double exponential fitting due to the double decay (I'll may be wrong though). The code and the produced graph are below. enter image description here

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

def double_exponential(x, a, b, c, d):
    return a*np.exp(-c*(x-b))+d

fig, ax = plt.subplots()
x = np.array([7.13793103, 21.4137931, 35.68965517, 49.96551724, 64.24137931, 78.51724138, 92.79310345, 107.06896552, 121.34482759, 135.62068966, 149.89655172, 164.17241379, 178.44827586, 192.72413793, 207.00000000, 221.27586207, 235.55172414, 249.82758621, 264.10344828, 278.37931034, 292.65517241, 306.93103448, 321.20689655, 335.48275862, 349.75862069, 364.03448276, 378.31034483, 392.5862069, 406.86206897, 421.13793103])
y = np.array([169954, 20599, 7309, 3989, 2263, 1569, 1134, 1017, 788, 806, 693, 528, 502, 431, 445, 367, 277, 267, 255, 189, 171, 109, 76, 36, 18, 9, 4, 3, 2, 1])
y_pruned = np.where(y < 1, 1, y)
p0 = (1.0, 1.0, 1.0, 1.0)  # only an example (I've tried other values without success)
popt, pcov = curve_fit(double_exponential, x, np.log(y_pruned), p0=p0)

ax.plot(x, double_exponential(x, *popt), 'g--')
ax.plot(x, np.log(y_pruned), 'ro')

ax.set_yscale('log')  # I need to have the y-axis logarithmic
plt.show()
Solution

  • I might be wrong but it does not seem to me to be a double exponential. In the log-plot, it seems a fourth-order polynomial, and if you try, it fits quite well:

    enter image description here

    I just fitted x and np.log(y) using the function:

    def f(x, a, b, c, d):
    return a*x**3 +b*x**2 + c*x + d

    the result of the fit is:

    a:  -4.933e-07 +/- 4.82e-08   (-9.8%)  initial:1
b:  0.0003011 +/- 3.14e-05   (10.4%)  initial:1
c:  -0.06804 +/- 0.00579   (-8.5%)  initial:1
d:  11.406  +/- 0.286      (2.5%)  initial:1

    If you really want to fit x and y and not log(y) it is a bit more tricky since the residual coming from the first few points are way bigger and the fit function ignores the small point on the right:

    def f(x, a=-4.933e-07, b=0.0003011, c=-0.06804, d=11.406):
    return np.exp(a*x**3 +b*x**2 + c*x + d)

    enter image description here

    Of course, you can adjust this behavior weighing the points in the fit function, and it works, but it seems to me to be a bit overcomplicated with respect to just fitting the log of the data:

    enter image description here