使用 lmfit Minimizer.minimize() 会导致错误 ValueError:“具有多个元素的数组的真值不明确”

Using lmfit Minimizer.minimize() leads to error ValueError: "The truth value of an array with more than one element is ambiguous"

提问人:Paul Joh 提问时间:10/6/2023 最后编辑:Paul Joh 更新时间:10/10/2023 访问量:240

问:

我想使用 lmfit 来适应该函数:

func(x,region,E0,C0,R1,R3,R4,R5,R6,R7,R8,alpha,beta,rho,theta,delta,d)

并获取有关参数置信区间的信息。 将浮点数或值介于 0 和 78 之间的浮点数数组作为参数。该函数调用从 cpp 文件生成的 exe。然后从这个 exe 的外游返回一个与 x 值相对应的函数值数组。funcx

我还有一个称为数据数组的数据,我想适合它。我确信它包含正确的值并具有正确的大小。以下是我的代码的缩短版本:total_data_arraytotal_data_array

#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)

results = my_model.fit(total_data_array, params=parameters,x=xData)

工作得很好。此外,合身性看起来确实不错。但我希望能够查看参数的置信区间。据我所知,有两种方法可以实现这一目标。

  1. 使用我拥有的模型或

  2. 使用最小化器

如果我尝试模型方法,将以下行添加到我的代码中:

    ci = lmfit.conf_interval(my_model,results)
    lmfit.report_ci(ci)

我收到错误:

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

绘制产量图:fit_reports

[[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)

因此,所有值都与初始值相比发生了显着变化。 我在其他帖子中读到,问题可能是某些参数的相关性非常强。为了消除这种可能性,我设置了除两个参数之外的所有参数,这并没有改变错误。我为不同的“两个”参数集做了这件事。parameters['...'].vary = False

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 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:total_data_array

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 to as parameter in the function yields instead the error:fcn_args=args=minimize()

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 . But I don't know what to pass as keywords or why. I appreciate any help :)minimize()

Python Scipy 曲线 数据拟合 lmfit

评论

1赞 hpaulj 10/6/2023
'playing around' is not good programming practice. You should understand what you are doing, and why. Getting the number of parameters right requires paying attention to how they are passed to your . That's a matter of carefully reading the docs, and keeping track of what your needs. As for the ambiguity error, you don't provide any stack trace, so we can't say what's causing it. Uusually it means you supply an array to some line of code that expects a simple True/False value, such as a or .funcfuncifor
0赞 Paul Joh 10/7/2023
Thanks for the Input. I have edited my post
2赞 M Newville 10/7/2023
I would not necessarily dismiss "playing around" ;) But I would suggest going back to the first fit that you said "worked". Why can't you get confidence intervals for that fit? You did not show a fit report. You also did not show anything even approaching a complete example. What did you try to get confidence intervals? Your second example just shows a programming error that seems pretty clear, but without showing any code, you cannot expect anyone to be able to help you with that.
0赞 M Newville 10/7/2023
uh, well, what gets printed in your fit report? Does your fit "work " to change all the variables away from their initial values? What the heck is your meant to do? That seems sort of odd to me.linearize()
1赞 ken 10/9/2023
I guess might work. To get an answer that is not a "guess", please provide a reproducible example. Often, creating a reproducible example will help you solve your problem.ci = lmfit.conf_interval(results, results)

答:

1赞 M Newville 10/10/2023 #1

You give two attempts at using lmfit. The first one is correct, and the fit runs. The second, using lmfit.Minimize does not work because your is a model function not an objective function: it does not have the same call signature, which is what the exceptions you are seeing tell you. But there is no point in trying lmfit.Minimize anyway. The lmfit.Model runs, it just tells you that it cannot estimate uncertainties. Also FWIW, there is no point in running if the initial fit cannot estimate uncertainties. That will not "fix the problem".funclmfit.conf_interval

As you will have read in the lmfit FAQ, there are a few reasons why uncertainties sometimes cannot be estimated, but they all derive from the same root cause: small changes to the value of one or more of the parameters has no noticeable effect on the fit. When that happens, the covariance matrix cannot be determined -- the fit algorithm has decided that at least one of the parameters is not actually affecting the fit. One reason that can happen is that a parameter gets stuck at a boundary or initial value. That does not seem to be the case for your fit. Another possibility is that one or more parameters have reached such extreme values (say, nearly 0) that other parameters no longer have an effect on the fit. Imagine that you are modeling a peak with a center, width, and height. If the height value moves to zero, it doesn't really matter what the width or center value is, that model will just be zero.

I guess that is what is happening with your fit. The values for several of your parameters go to ~1e-7 level -- 5 or more orders of magnitude smaller than your initial estimate of them.

I must say that you do not describe your model at all, and are running some external program and then doing some fairly clunky parsing of a text of results. We don't know what any of that is doing. So, there is no way for anyone else to know that if, say, goes from 0.99 to 1.4-e7, the rest of the model even makes sense. Maybe you know that.
Anyway, that is where I would start looking. It seems likely that you have not run your model with initial values (or final values) and plotted the results: that is always highly recommended to check that your initial values are reasonable. With such a large value of Chi-square and a value of R so far from 1, it would be easy to just assert that you have a terrible fit. If you don't have a good fit, uncertainties would be meaningless anyway.
d

评论

0赞 Paul Joh 10/11/2023
Thank you, I see your reasoning and that is probably the best answer that can be given. The model is extremely complicated so I don't think it would make much sense to post it. But I have an idea now how to go on with the project :)
0赞 M Newville 10/12/2023
@PaulJoh it is okay that you cannot post the model of "run complex program, parse results". It's just that we have to assume that this is all being done well. And, maybe that model is so opaque that it is not obvious what "R5 goes to 0" actually means for the total result. So, you may need to spend some time (I might say "playing around" ;) ) running the model with various inputs to see if it is doing what you expect, and if the fit is actually any good.