I want to use lmfit to fit the function:
func(x,region,E0,C0,R1,R3,R4,R5,R6,R7,R8,alpha,beta,rho,theta,delta,d)
and get the information on confidence intervals of the parameters. func takes as parameter x an float or an array of floats with values in between 0 and 78. The function calls an exe generated from a cpp-file. And from the outstream of this exe then returns an array of function values corresponding to the x-values.
I also have an array of data called total_data_array that I want to fit to. I am certain that total_data_array contains the correct values and has the correct size. The following is a shortened version of my code:
#Define fitting function
def func(x,region,E0,C0,R1,R3,R4,R5,R6,R7,R8,alpha,beta,rho,theta,delta,d):#x should be between 0 and 75
print("Function called")
print(E0,C0,R1,R3,R4,R5,R6,R7,R8,alpha,beta,rho,theta,delta,d)
#Make x iterable if it is only a float/int
if not hasattr(x,'__iter__'):
x = np.array([x])
result = []
output = subprocess.check_output([r'...\model.exe', str(region),str(E0),str(C0),str(R1),str(R3),str(R4),str(R5),str(R6),str(R7),str(R8),str(alpha),str(beta),str(rho),str(theta),str(delta),str(d)])
output_str = codecs.decode(output)
output_str_list = output_str.split('\r\n')
output_str_list.pop()
dataarray = []
index = 0
for word in output_str_list:
if index in range(416,624) or index in range(728,832):
if word == '-nan(ind)':
dataarray.append(0.0)
else:
dataarray.append(float(word))
index+=1
for x0 in x:
y = dataarray(int(x0*4))
result.append(y)
return result
parameters = lmfit.Parameters()
parameters.add('region',value=1)
parameters.add('E0',value=500,min=0.0,max = 10000.0)
parameters.add('C0',value=200,min=0.0,max = 10000.0)
parameters.add('R1',value=0.587,min=0.0,max=1.0)
parameters.add('R3',value=0.3125,min=0.0,max=1.0)
parameters.add('R4',value=0.1666,min=0.0,max=1.0)
parameters.add('R5',value=0.1,min=0.0,max=1.0)
parameters.add('R6',value=0.2,min=0.0,max=1.0)
parameters.add('R7',value=0.4,min=0.0,max=1.0)
parameters.add('R8',value=0.125,min=0.0,max=1.0)
parameters.add('alpha',value=0.09,min=0.0,max=1.0)
parameters.add('beta',value=0.25,min=0.0,max=1.0)
parameters.add('rho',value=0.2,min=0.0,max=1.0)
parameters.add('theta',value=0.26,min=0.0,max=1.0)
parameters.add('delta',value=0.77,min=0.0,max=1.0)
parameters.add('d',value=0.99,min=0.0,max=1.0)
parameters['region'].vary = False
xData = np.arange(0, 78, 1)
#Running this part of the code works
my_model = lmfit.Model(func)
results = my_model.fit(total_data_array,params=parameters,x=xData)
Calling
results = my_model.fit(total_data_array, params=parameters,x=xData)
works perfectly fine. Also the fit does look quite good. But I want to be able to view the confidence intervals of the parameters. As far as I know, there are two ways to achieve this.
-
Using the model that I have or
-
Using Minimizers.
If I try the model approach, adding the lines to my code:
ci = lmfit.conf_interval(my_model,results)
lmfit.report_ci(ci)
I get the error:
Traceback (most recent call last):
File "c:\Users\paul1\OneDrive\Desktop\epidemiology\coding\Secihurd_Model\src\python\regional fitting\fitting2.py", line 128, in <module>
ci = lmfit.conf_interval(my_model,results)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "...\lmfit\confidence.py", line 141, in conf_interval
ci = ConfidenceInterval(minimizer, result, p_names, prob_func, sigmas,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "...\confidence.py", line 185, in __init__
raise MinimizerException(CONF_ERR_STDERR)
lmfit.minimizer.MinimizerException: Cannot determine Confidence Intervals without sensible uncertainty estimates
Plotting the fit_reports yields:
[[Model]]
Model(func)
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 511
# data points = 78
# variables = 15
chi-square = 34979.5187
reduced chi-square = 555.230456
Akaike info crit = 506.253115
Bayesian info crit = 541.603747
R-squared = 0.62955613
## Warning: uncertainties could not be estimated:
[[Variables]]
region: 1 (fixed)
E0: 4.04296161 (init = 10)
C0: 1.26289805 (init = 5)
R1: 0.63653349 (init = 0.587)
R3: 0.60694072 (init = 0.3125)
R4: 0.11924779 (init = 0.1666)
R5: 2.7998e-07 (init = 0.1)
R6: 0.58190868 (init = 0.2)
R7: 0.49973503 (init = 0.4)
R8: 9.4685e-07 (init = 0.125)
alpha: 9.7636e-05 (init = 0.09)
beta: 0.08512620 (init = 0.25)
rho: 0.55511880 (init = 0.2)
theta: 0.08153075 (init = 0.26)
delta: 0.38812246 (init = 0.77)
d: 1.4235e-07 (init = 0.99)
So all values have changed significantly from the initial values. I read in other posts that the problem might be that some parameters correlate very strongly. To eliminate this possibility, I set for all but two parameters parameters['...'].vary = False, which did not change the error. I did this for different sets of 'two' parameters.
Now, if I instead follow the second approach: I add the lines:
min = lmfit.Minimizer(func,params=parameters,fcn_args=np.arange(0,78, 1))
min.minimize(method='leastsq')
ci = lmfit.conf_interval(min)
lmfit.report_ci(ci)
So now I only minimize my function. I will later subtract the total_data_array to minimize the difference between my function and the data and obtain a curve fit. For now, I just want to get this to work. I get the error:
Traceback (most recent call last):
File "...\fitting2.py", line 134, in <module>
min.minimize(method='leastsq')
File "...\lmfit\minimizer.py", line 2345, in minimize
return function(**kwargs)
^^^^^^^^^^^^^^^^^^
File "...\lmfit\minimizer.py", line 1651, in leastsq
lsout = scipy_leastsq(self.__residual, variables, **lskws)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "...\scipy\optimize\_minpack_py.py", line 415, in leastsq
shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "...\scipy\optimize\_minpack_py.py", line 25, in _check_func
res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:...\lmfit\minimizer.py", line 548, in __residual
out = self.userfcn(params, *self.userargs, **self.userkws)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: func() takes 17 positional arguments but 79 were given
so already the minimisation seems to fail. Changing fcn_args= to args= as parameter in the minimize() function yields instead the error:
File "...\fitting2.py", line 147, in <module>
min.minimize(method='leastsq')
File "...\lmfit\minimizer.py", line 2345, in minimize
return function(**kwargs)
^^^^^^^^^^^^^^^^^^
File "...\lmfit\minimizer.py", line 1651, in leastsq
lsout = scipy_leastsq(self.__residual, variables, **lskws)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "...\scipy\optimize\_minpack_py.py", line 415, in leastsq
shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "...\scipy\optimize\_minpack_py.py", line 25, in _check_func
res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "...\lmfit\minimizer.py", line 530, in __residual
if apply_bounds_transformation:
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
So it might be that the problem lies in some keyword-arguments that need to be passed to minimize(). But I don't know what to pass as keywords or why. I appreciate any help :)
from Using lmfit Minimizer.minimize() leads to error ValueError: "The truth value of an array with more than one element is ambiguous"
No comments:
Post a Comment