Python curve fitting with constraints

Question

I have been looking for Python curve fitting with constraints. One of the option is to use lmfit module and another option is to use penalization to enforce the constraints. I have following code in which I am trying to enforce a+b=3.6 as the constraint. In other words, y=3.6 when x=1 and x is always >=1 in my case.

import numpy as np
import scipy.optimize as sio
def func(x, a, b, c):
    return a+b*x**c

x = [1, 2, 4, 8, 16]
y = [3.6, 3.96, 4.31, 5.217, 6.842]
lb = np.ones(3, dtype=float)
ub = np.ones(3, dtype=float)*10.
popt, pcov = sio.curve_fit(func, x, y)
print(popt)

Ideally, I would like to use the lmfit approach and spent good amount of trying to understand examples but could not succeed. Can someone help with an example?

Answer 1

If I understand your question correctly, you want to model some data with

def func(x, a, b, c):
    return a+b*x**c

and for a particular set of data you want to impose the constraint that a+b=3.6. You could, just "hardwire that in", changing the function to be

def func2(x, b, c):
    a = 3.6 - b
    return a+b*x**c

and now you have a model function with only two variables: b, and c.

That would not be very flexible but would work.

Using lmfit gives back some of that flexibility. To do a completely unconstrained fit, you would say

from lmfit import Model

mymodel = Model(func)
params = mymodel.make_params(a=2, b=1.6, c=0.5)
result = mymodel.fit(y, params,  x=x)

(Just as an aside: scipy.optimize.curve_fit permits you to not specify initial values for the parameters and implicitly sets them to 1 without telling you. This is a terrible mis-feature - always give initial values).

If you do want to impose the constraint a+b=3.6, you could then do

params['a'].expr = '3.6-b'

result2 = mymodel.fit(y, params, x=x)
print(result2.fit_report())

When I do that with the data you provided, this prints (note that it reports 2 variables, not 3):

[[Model]]
    Model(func)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 34
    # data points      = 5
    # variables        = 2
    chi-square         = 0.01066525
    reduced chi-square = 0.00355508
    Akaike info crit   = -26.7510142
    Bayesian info crit = -27.5321384
[[Variables]]
    a:  3.28044833 +/- 0.04900625 (1.49%) == '3.6-b'
    b:  0.31955167 +/- 0.04900626 (15.34%) (init = 1.6)
    c:  0.86901253 +/- 0.05281279 (6.08%) (init = 0.5)
[[Correlations]] (unreported correlations are < 0.100)
    C(b, c) = -0.994

Your code hinted at using (but did not actually use) upper and lower bounds for the parameter values. Those are also possible with lmfit, as with

params['b'].min = 1
params['b'].min = 10

and so forth. I'm not sure you need them here, and would caution against trying to set bounds too tightly.