I'm using lmfit to fit a "complex" model onto experimental data.
I define the parameters of the model like this:
def model(t, y0, h, A1, x1, s1, A2, x2, s2, A3, x3, s3, A4, x4, s4):
y = y0 + h * (
A1 * (0.5 + 0.5 * erf((t - x1) / (s1 * 2**0.5))) +
A2 * (0.5 + 0.5 * erf((t - x2) / (s2 * 2**0.5))) +
A3 * (0.5 + 0.5 * erf((t - x3) / (s3 * 2**0.5))) +
A4 * (0.5 + 0.5 * erf((t - x4) / (s4 * 2**0.5)))
)
return y
mod = Model(model)
# sys.setrecursionlimit(150000)
# Set the parameters, define their bounds and their relations
mod.set_param_hint('y0', value=self.spin_y0.value())
mod.set_param_hint('h', value=self.spin_h.value())
mod.set_param_hint('A1', value=self.spin_A1.value(), min=0, max=1),
mod.set_param_hint('x1', value=self.spin_x1.value(), min=0, max=self.x[-1])
mod.set_param_hint('s1', value=self.spin_s1.value(), min=0)
mod.set_param_hint('A2', value=self.spin_A2.value(), min=0, max=1),
mod.set_param_hint('x2', value=self.spin_x2.value(), min=0, max=self.x[-1])
mod.set_param_hint('s2', value=self.spin_s2.value(), min=0)
mod.set_param_hint('A3', value=self.spin_A3.value(), min=0, max=1),
mod.set_param_hint('x3', value=self.spin_x3.value(), min=0, max=self.x[-1])
mod.set_param_hint('s3', value=self.spin_s3.value(), min=0)
mod.set_param_hint('A4', value=self.spin_A4.value(), min=0, max=1),
mod.set_param_hint('x4', value=self.spin_x4.value(), min=0, max=self.x[-1])
mod.set_param_hint('s4', value=self.spin_s4.value(), min=0)
mod.set_param_hint('A1', expr='1-A2-A3-A4')
mod.set_param_hint('A2', expr='1-A1-A3-A4')
mod.set_param_hint('A3', expr='1-A1-A2-A4')
mod.set_param_hint('A4', expr='1-A1-A2-A3')
mod.make_params()
# Fit the data !
result = mod.fit(self.y, t=self.x)
I use spinBoxes to get the values which will be used as initial values (my program is a GUI).
I have a problem with A1, A2, A3 and A4. They are dependent, and their sum must be equal to 1.
If I comment
mod.set_param_hint('A1', expr='1-A2-A3-A4')
mod.set_param_hint('A2', expr='1-A1-A3-A4')
mod.set_param_hint('A3', expr='1-A1-A2-A4')
mod.set_param_hint('A4', expr='1-A1-A2-A3
The parameters vary between 0 and 1 (I define this range when I define the parameters), but their sum is never equal to 1.
When these lines are not commented, I have the following exception:
File "/usr/lib/python3.4/site-packages/lmfit/minimizer.py", line 255, in __update_paramval
self.__update_paramval(dep)
File "/usr/lib/python3.4/site-packages/lmfit/minimizer.py", line 250, in __update_paramval
if getattr(par, 'expr', None) is not None:
RuntimeError: maximum recursion depth exceeded
I don't really know what it means. I tried to increase the maximum recursion depth, but it changed nothing.
Do you have any idea about what is the problem ?
Ok, for future me and other persons looking for a solution, I found one. I used scimin.fmin_slsqp to do the job. Here is a sample of code:
def model(t, params):
y0, h, A1, x1, s1, A2, x2, s2, A3, x3, s3, A4, x4, s4 = [value for value in params]
y = y0 + h * (
A1 * (0.5 + 0.5 * erf((t - x1) / (s1 * 2**0.5))) +
A2 * (0.5 + 0.5 * erf((t - x2) / (s2 * 2**0.5))) +
A3 * (0.5 + 0.5 * erf((t - x3) / (s3 * 2**0.5))) +
A4 * (0.5 + 0.5 * erf((t - x4) / (s4 * 2**0.5)))
)
return y
def sum_residuals(params): # the function we want to minimize
y0, h, A1, x1, s1, A2, x2, s2, A3, x3, s3, A4, x4, s4 = [value for value in params]
return sum((list_y - model(list_x, params))**2)
def constraints(params):
y0, h, A1, x1, s1, A2, x2, s2, A3, x3, s3, A4, x4, s4 = [value for value in params]
return 1 - (A1 + A2 + A3 + A4)
pfit = scimin.fmin_slsqp(sum_residuals,
p0,
f_eqcons=constraints,
bounds=bounds,
iprint=0,
iter=300
)
Where p0 is a sequence of the initial parameters. However, fmin_slsqp doesn't return a covariance matrix, I couldn't calculate confidence intervals on the fitted parameters. For that I used the bootstrap method.