Curve fitting with conditional equation

Question

Problem
I have created a curve fitting exercise (see functional code below), but I would like to add to the functionality.
I need to be able to define the following condition: slope at min(xdata) = 0.
(in words: I want the fitted curve to start out with horizontal gradient)

What I have tried
I have spent quite a bit of time researching scipy.optimize.curve_fit and evaluated other options (lmfit package, and scipy functions scipy.optimize.fmin_slsqp, scipy.optimize.minimize, etc.). lmfit only allows me to set a static condition on the parameters, such as p1 = 2 * p2 + 3. But it does not allow me to address min(xdata) dynamically, and I cannot make use of the derivate in the constraint.
Scipy only allows me to minimize the function (find an optimal x, but parameters p are already known). Or it can be used to define a specific range for the parameters. I was not able to define a second function that can be used to constrain the parameters during the curve fitting.

I need to be able to pass the condition directly to the curve fitting algorithm (rather than addressing the problem by bringing the condition into the cubic_fit() equation - it seems possible to eliminate e.g. p3 and define it as a combination of the other parameters and min(xdata)). My actual fitting function is much more complex and I need to run this script iteratively on a batch of data (varying min(xdata)). I cannot manually alter the fitting function each time...

I am grateful for any suggestions, maybe there are other packages out there that allow for a more complex definition of the curve fitting problem?

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
import scipy.optimize

# generate dummy data - on which I will run a curve fit below
def cubic_fit_with_noise(x, p1, p2, p3, p4):
    return p1 + p2*x + p3*x**2 + p4*x**3 + np.random.rand()

xdata = [x * 0.1 for x in range(0, 100)]
ydata = np.array( [cubic_fit_with_noise (x, 2, 0.4, -.2,0.02) for x in xdata] )

# now, run the curve-fit

#  set up the fitting function: 
def cubic_fit(x, p1, p2, p3, p4):
    return p1 + p2*x + p3*x**2 + p4*x**3

# define starting point:
s1 = 2.5
s2 = 0.2
s3 = -.2
s4 = 0.02

# scipy curve fitting:
popt, pcov = scipy.optimize.curve_fit(cubic_fit, xdata, ydata, p0=(s1,s2,s3,s4))
y_modelled = np.array([cubic_fit(x, popt[0], popt[1], popt[2], popt[3]) for x in xdata])

print(popt)     # prints out the 4 parameters p1,p2,p3,p4 defined in curve-fitting

plt.plot(xdata, ydata, 'bo')
plt.plot(xdata, y_modelled, 'r-')
plt.show()

The above code runs with Python3 (fix the print statement if you have Python2). As an addition, I want to bring in the derivative:

def cubic_fit_derivative(x, p1, p2, p3, p4):
    return p2 + 2.0 * p3 * x + 3 * p4 * x**2

and the constraint that cubic_fit_derivative(min(xdata), p1,p2,p3,p4) = 0.

quick guess fit a * (x - xmin)**2 * (x - b) or something like this — mikuszefski
Thanks for your feedback, everyone. The suggestions I see below all seem to work (at the time of writing). I have accepted the one by @Newville as the answer as it has the clearest approach and is highly versatile, so can easily be adopted for more complex tasks. — rde
I agree, though I can tell from experience: Don't fit with constraints if you don't have to. — mikuszefski
@mikuszefski in other words, if you have to use constraints, use them. ;). In this case, the constraint simply removes a variable parameter from the fit. One could have altered the objective function to impose p2 = -p3*xmin - 3*p4*xmin**2, and used only three variables (p1, p3, p4), but that's not a very general approach. In a sense, the model function always does impose a constraint on the fit. The approach described with lmfit just makes it easier to adjust what values the parameter take. — M Newville
@MNewville I never doubted that the constraints-approach is very general and powerful. Yet, depending on the type of constraint and the size of your data set, computationally it can be very painful. From that point of view you gave the right answer to the wrong question. A third order polynomial with slope zero at x0 is just a polynomial where the homogeneous part has a twin root at x0. No need for any constraints, just simple maths and just a simple fit. On the other hand, we do not know what the OP's true problem is, and it might be very different from a third order polynomial. Cheers. — mikuszefski

M Newville M Newville · Accepted Answer · 2017-10-10T22:50:54

Your condition that the derivative of your polynomial = 0 at xmin can be expressed as a simple constraint and means that the variables p2, p3, and p4 are not actually independent. The derivate condition is

p2 + 2*p3*xmin + 3*p4*xmin**2 = 0

where xmin is the minimum value of xdata. Furthermore, xmin will be known prior to the fit (if not necessarily when your script is written), you can use this to constrain one of the three parameters. Since xmin may be zero (in fact, it is for your case), the constraint should be that

p2 = - 2*p3*xmin - 3*p4*xmin**2

Using lmfit, the original, unconstrained fit would look like this (I cleaned it up a bit):

import numpy as np
from lmfit import Model
import matplotlib.pylab as plt

#  the model function:
def cubic_poly(x, p1, p2, p3, p4):
    return p1 + p2*x + p3*x**2 + p4*x**3

xdata = np.arange(100) * 0.1
ydata = cubic_poly(xdata, 2, 0.4, -.2, 0.02)
ydata = ydata + np.random.normal(size=len(xdata), scale=0.05)

# make Model, create parameters, run fit, print results
model  = Model(cubic_poly)
params = model.make_params(p1=2.5, p2=0.2, p3=-0.0, p4=0.0)
result = model.fit(ydata, params, x=xdata)
print(result.fit_report())

plt.plot(xdata, ydata, 'bo')
plt.plot(xdata, result.best_fit, 'r-')
plt.show()

which prints:

[[Model]]
    Model(cubic_poly)
[[Fit Statistics]]
    # function evals   = 13
    # data points      = 100
    # variables        = 4
    chi-square         = 0.218
    reduced chi-square = 0.002
    Akaike info crit   = -604.767
    Bayesian info crit = -594.347
[[Variables]]
    p1:   2.00924432 +/- 0.018375 (0.91%) (init= 2.5)
    p2:   0.39427207 +/- 0.016155 (4.10%) (init= 0.2)
    p3:  -0.19902928 +/- 0.003802 (1.91%) (init=-0)
    p4:   0.01993319 +/- 0.000252 (1.27%) (init= 0)
[[Correlations]] (unreported correlations are <  0.100)
    C(p3, p4)                    = -0.986 
    C(p2, p3)                    = -0.967 
    C(p2, p4)                    =  0.914 
    C(p1, p2)                    = -0.857 
    C(p1, p3)                    =  0.732 
    C(p1, p4)                    = -0.646

and produces a plot of

Now, to add your constraint condition, we will add xmin as a fixed parameter, and constrain p2 as above, replace the above with:

params = model.make_params(p1=2.5, p2=0.2, p3=-0.0, p4=0.0)

# add an extra parameter for `xmin`
params.add('xmin', min(xdata), vary=False)

# constrain p2 so that the derivative is 0 at xmin
params['p2'].expr = '-2*p3*xmin - 3*p4*xmin**2'

result = model.fit(ydata, params, x=xdata)
print(result.fit_report())

plt.plot(xdata, ydata, 'bo')
plt.plot(xdata, result.best_fit, 'r-')
plt.show()

which now prints

[[Model]]
    Model(cubic_poly)
[[Fit Statistics]]
    # function evals   = 10
    # data points      = 100
    # variables        = 3
    chi-square         = 1.329
    reduced chi-square = 0.014
    Akaike info crit   = -426.056
    Bayesian info crit = -418.241
[[Variables]]
    p1:     2.39001759 +/- 0.023239 (0.97%) (init= 2.5)
    p2:     0          +/- 0        (nan%)  == '-2*p3*xmin - 3*p4*xmin**2'
    p3:    -0.10858258 +/- 0.002372 (2.19%) (init=-0)
    p4:     0.01424411 +/- 0.000251 (1.76%) (init= 0)
    xmin:   0 (fixed)
[[Correlations]] (unreported correlations are <  0.100)
    C(p3, p4)                    = -0.986 
    C(p1, p3)                    = -0.742 
    C(p1, p4)                    =  0.658

and a plot like

If xmin had not been zero (say, xdata = np.linspace(-10, 10, 101), the value and uncertainty of p2 would not be zero.

Curve fitting with conditional equation

3 Answers