Using scipy.optimize.curve_fit to fit a piecewise function

Question

I am attempting to write a program that reads two sets data from a .csv file into arrays, and then fits it to a piecewise function. What's most important to me is that these fits are done simultaneously because they have the same parameters. This piecewise function is my attempt to do so, though if you know a better way to fit them simultaneously I'd also greatly appreciate advice regarding that.

To avoid having to upload the csv files I've added the data directly into the arrays.

import numpy
import csv
import matplotlib
from scipy import optimize

xdata = [2.0, 10.0, 30.0, 50.0, 70.0, 90.0, 110.0, 130.0, 150.0, 250.0, 400.0, 1002.0, 1010.0, 1030.0, 1050.0, 1070.0, 1090.0, 1110.0, 1130.0, 1150.0, 1250.0, 1400.0]
ydata = [0.013833958803215633, 0.024273268442992078, 0.08792766000711709, 0.23477725658012044, 0.31997367288103884, 0.3822895295625711, 0.46037063893452784, 0.5531831477605121, 0.559757863748663, 0.6443036770720387, 0.7344601382896991, 2.6773979205076136e-09, 9.297289736857164e-10, 0.10915332214935693, 0.1345307163724643, 0.1230161681870127, 0.11286094974672768, 0.09186485171688986, 0.06609131137369342, 0.052616358869021135, 0.034629686697483314, 0.03993853791147095]

The first 11 points I want to fit to the function labeled 'SSdecay', and the second 11 points I want to fit to the function labeled 'SUdecay'. My attempt at doing this simultaneously was making the piecewise function labeled 'fitfunciton'.

#defines functions to be used in fitting

#to fit the first half of data
def SSdecay(x, lam1, lam2, norm, xoff):
    return norm*(1 + lam2/(lam1 - lam2)*numpy.exp(-lam1*(x - xoff)) -
        lam1/(lam1 - lam2)*numpy.exp(-lam2*(x - xoff)))

#to fit the second half of data
def SUdecay(x, lam1, lam2, norm, xoff):
    return norm*(lam1/(lam1 - lam2))*(-numpy.exp(-lam1*(x - xoff)) +
        numpy.exp(-lam2*(x - xoff)))

#piecewise function combining SS and SU functions to fit the whole data set
def fitfunction(x, lam1, lam2, norm, xoff):
    y = numpy.piecewise(x,[x < 1000, x >= 1000],[SSdecay(x, lam1, lam2, norm, xoff),SUdecay(x, lam1, lam2, norm, xoff)])
    return y

#fits the piecewise function with initial guesses for parameters
p0=[0.01,0.02,1,0]
popt, pcov = optimize.curve_fit(fitfunction, xdata, ydata, p0)

print(popt)
print(pcov)

After running this I get the error:

ValueError: NumPy boolean array indexing assignment cannot assign 22 input values to the 11 output values where the mask is true

It seems as though curve_fit does not like that I'm using a piecewise function but I am unsure why or if it is a fixable kind of problem.

I have a suggestion: try fitting each piece individually first. When I tried this for the first 11 data points and SSDecay() the fit was not good, which might be the initial parameters(?). If you can get good fits for each individual function to its respective portion of the data, then combining them into a single piecewise model will be much easier. — James Phillips, Jun 12 '18 at 21:50
Hey thanks for your suggestion, I realized I was using the wrong, unnormalized data set for the ydata. The data should now actually resemble the function and produce a fit. However when I try to fit the piecewise it still throws the error. — Aaron Rea, Jun 13 '18 at 18:18

score 0 · Answer 1 · answered Jun 13 '18 at 22:22

Here are my results for separately fitting the two functions using the normalized data. It looks unlikely that these will work as a single piecewise equation, please see the plot image and source code below. I also have very different fitted parameters for the two equations:

SS parameters: [ 0.0110936, 0.09560932, 0.72929264, 6.82520026]

SU parameters: [ 3.46853883e-02, 9.54208972e-03, 1.99877873e-01, 1.00465563e+03]

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

xdata = [2.0, 10.0, 30.0, 50.0, 70.0, 90.0, 110.0, 130.0, 150.0, 250.0, 400.0, 1002.0, 1010.0, 1030.0, 1050.0, 1070.0, 1090.0, 1110.0, 1130.0, 1150.0, 1250.0, 1400.0]
ydata = [0.013833958803215633, 0.024273268442992078, 0.08792766000711709, 0.23477725658012044, 0.31997367288103884, 0.3822895295625711, 0.46037063893452784, 0.5531831477605121, 0.559757863748663, 0.6443036770720387, 0.7344601382896991, 2.6773979205076136e-09, 9.297289736857164e-10, 0.10915332214935693, 0.1345307163724643, 0.1230161681870127, 0.11286094974672768, 0.09186485171688986, 0.06609131137369342, 0.052616358869021135, 0.034629686697483314, 0.03993853791147095]

#to fit the first half of data
def SSdecay(x, lam1, lam2, norm, xoff):
    return norm*(1 + lam2/(lam1 - lam2)*numpy.exp(-lam1*(x - xoff)) -
        lam1/(lam1 - lam2)*numpy.exp(-lam2*(x - xoff)))

#to fit the second half of data
def SUdecay(x, lam1, lam2, norm, xoff):
    return norm*(lam1/(lam1 - lam2))*(-numpy.exp(-lam1*(x - xoff)) +
        numpy.exp(-lam2*(x - xoff)))

# some initial parameter values
initialParameters_ss = numpy.array([0.01, 0.02, 1.0, 0.0])
initialParameters_su = initialParameters_ss # same values for this example

# curve fit the equations individually to their respective data
ssParameters, pcov = curve_fit(SSdecay, xdata[:11], ydata[:11], initialParameters_ss)
suParameters, pcov = curve_fit(SUdecay, xdata[11:], ydata[11:], initialParameters_su)

# values for display of fitted function
lam1_ss, lam2_ss, norm_ss, xoff_ss = ssParameters
lam1_su, lam2_su, norm_su, xoff_su = suParameters

# for plotting the fitting results
y_fit_ss = SSdecay(xdata[:11], lam1_ss, lam2_ss, norm_ss, xoff_ss) # first data set, first equation
y_fit_su = SUdecay(xdata[11:], lam1_su, lam2_su, norm_su, xoff_su) # second data set, second equation

plt.plot(xdata, ydata, 'D') # plot the raw data as a scatterplot
plt.plot(xdata[:11], y_fit_ss) # plot the SS equation using the fitted parameters
plt.plot(xdata[11:], y_fit_su) # plot the SU equation using the fitted parameters
plt.show()

print('SS parameters:', ssParameters)
print('SU parameters:', suParameters)

Using scipy.optimize.curve_fit to fit a piecewise function

1 Answers1