Error in Scipy curve fit for more than two parameters

Question

I am quite new to Scipy. I have a data file (https://www.dropbox.com/s/mwz8s2kap2mnwo0/data.dat?dl=0) and want to fit the function aexp(bx^c). The problem is when I give manually the value of c (say c = 0.75), the code works perfectly, but if I want to find the 'a', 'b' and 'c' from the fit, the code does not work and producing a flat line. Sorry if the problem is too silly. The code reads as:

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

filename = sys.argv[1]

data = np.loadtxt(filename)

x = np.array(data[:,0])
y = np.array(data[:,1])

def func(x, a, b, c):
    return a*np.exp(b*x**c)

params = curve_fit(func, x, y)

[a, b, c] = params[0]

perr = np.sqrt(np.diag(params[1]))

x_new = []
y_new = []

for i in np.linspace(1.00003e-05, 0.10303175629999914, num=1000):
    j = func(i, a, b, c)
    x_new.append(i)
    y_new.append(j)

x1 = np.array(x_new)
y1 = np.array(y_new)

print ("a = ", a, "error = ", perr[0], "error % = ", (perr[0]/a)*100, '\t' "b = ", b, "error = ", perr[1], "error % = ", (perr[1]/b)*100), '\t' "c = ", c, "error = ", perr[2], "error % = ", (perr[2]/c)*100,

#np.savetxt('fit.dat', np.c_[x1, y1])

plt.plot(x, y, label='data')
plt.plot(x1, y1, label = 'a*np.exp(b*x**c)')
plt.xlabel('Time(s)')
plt.ylabel('SRO')
plt.legend()
plt.show()

Answer 1

Exponential equations can be quite sensitive to the non-linear solver's initial parameter estimates. By default, many non-linear solvers - including scipy's curve_fit - use default initial parameter values of 1.0 for these initial parameter estimates if none are supplied, and in this particular case those values were not good initial estimates for your combination of data and equation. Scipy does include a genetic algorithm which can be used to determine the initial parameter estimates, and their implementation requires bounds within which to search. Here is an example graphical solver using the scipy differential_evolution genetic algorithm module for this purpose, note the ranges that I have used for the genetic algorithm to search within. It is much easier to give ranges for the parameters in this way rather than explicit values, though this is not always true it worked here. You will need to change the file path that I used to load the data.

import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution
import warnings

filename = '/home/zunzun/Downloads/data.dat'

data = numpy.loadtxt(filename)

xData = numpy.array(data[:,0])
yData = numpy.array(data[:,1])


def func(x, a, b, c):
    return a*numpy.exp(b*x**c)


# function for genetic algorithm to minimize (sum of squared error)
def sumOfSquaredError(parameterTuple):
    warnings.filterwarnings("ignore") # do not print warnings by genetic algorithm
    val = func(xData, *parameterTuple)
    return numpy.sum((yData - val) ** 2.0)


def generate_Initial_Parameters():
    # min and max used for bounds
    maxX = max(xData)
    minX = min(xData)
    maxY = max(yData)
    minY = min(yData)

    minData = min(minX, minY)
    maxData = min(maxX, maxY)

    parameterBounds = []
    parameterBounds.append([-maxData * 10.0, maxData * 10.0]) # search bounds for a
    parameterBounds.append([-maxData * 10.0, maxData * 10.0]) # search bounds for b
    parameterBounds.append([-maxData * 10.0, maxData * 10.0]) # search bounds for c

    # "seed" the numpy random number generator for repeatable results
    result = differential_evolution(sumOfSquaredError, parameterBounds, seed=3)
    return result.x

# by default, differential_evolution completes by calling curve_fit() using parameter bounds
geneticParameters = generate_Initial_Parameters()

# now call curve_fit without passing bounds from the genetic algorithm,
# just in case the best fit parameters are aoutside those bounds
fittedParameters, pcov = curve_fit(func, xData, yData, geneticParameters)
print('Fitted parameters:', fittedParameters)
print()

modelPredictions = func(xData, *fittedParameters) 

absError = modelPredictions - yData

SE = numpy.square(absError) # squared errors
MSE = numpy.mean(SE) # mean squared errors
RMSE = numpy.sqrt(MSE) # Root Mean Squared Error, RMSE
Rsquared = 1.0 - (numpy.var(absError) / numpy.var(yData))

print()
print('RMSE:', RMSE)
print('R-squared:', Rsquared)

print()


##########################################################
# graphics output section
def ModelAndScatterPlot(graphWidth, graphHeight):
    f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
    axes = f.add_subplot(111)

    # first the raw data as a scatter plot
    axes.plot(xData, yData,  'D')

    # create data for the fitted equation plot
    xModel = numpy.linspace(min(xData), max(xData))
    yModel = func(xModel, *fittedParameters)

    # now the model as a line plot
    axes.plot(xModel, yModel)

    axes.set_xlabel('X Data') # X axis data label
    axes.set_ylabel('Y Data') # Y axis data label

    plt.show()
    plt.close('all') # clean up after using pyplot

graphWidth = 800
graphHeight = 600
ModelAndScatterPlot(graphWidth, graphHeight)