[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Levenberg-Marquardt non-linear least-squares fitting in Python [follow-on]

Below I?ve included the code I ran, reasonably (I think) commented.  Note the reference to the example.  The data actually came from a pandas data frame that was in turn filled from a 100 MB data file that included lots of other data not needed for this, which was a curve fit to a calibration run.


PS:  If you want, I can probably still find a couple of the plots of the raw data and fitted result.
import numpy as np, matplotlib.pyplot as plt
#  Inverted exponential that axymptotically approaches "a" as x gets large
def func2fit(x,a,b,c):
    return a - b * np.exp(-c * x)    

# Curve fitting below from: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html
from scipy.optimize import curve_fit
def fit(xdata, ydata, run_num):
    ll = len(xdata)
#  The next four lines shift and scale the data so that the curve fit routine can
#  do its work without needing to use 8 or 16-byte precision. After fitting, we
#  will scale everything back.
    ltemp = [ydata[i] - ydata[0] for i in range(ll)]
    ytemp = [ltemp[i] * .001 for i in range(ll)]
    ltemp = [xdata[i] - xdata[0] for i in range(ll)]
    xtemp = [ltemp[i] * .001 for i in range(ll)]
#  popt is a list of the three optimized fittine parameters [a, b, c]
#  we are interested in the value of a.
#  cov is the 3 x 3 covariance matrix, the standard deviation (error) of the fit is
#  the square root of the diagonal.
    popt,cov = curve_fit(func2fit, xtemp, ytemp)
#  Here is what the fitted line looks like for plotting
    fitted = [popt[0] - popt[1] * np.exp(-popt[2] * xtemp[i]) for i in range(ll)]
#  And now plot the results to check the fit
    fig1, ax1 = plt.subplots()
    plt.title('Normalized Data ' + str(run_num))
    color_dic = {0: "red", 1: "green", 2: "blue", 3: "red", 4: "green", 5: "blue"}
    ax1.plot(xtemp, ytemp, marker = '.', linestyle  = 'none', color = color_dic[run_num])
    ax1.plot(xtemp, fitted, linestyle = '-', color = color_dic[run_num])
    plt.savefig('Normalized ' + str(run_num))
    perr = np.sqrt(np.diag(cov))
    return popt, cov, xdata[0], ydata[0], fitted, perr[0]