我想用lmfit来拟合函数:
func(x,region,E0,C0,R1,R3,R4,R5,R6,R7,R8,alpha,beta,rho,theta,delta,d)
并获取参数置信区间的信息。
func
将浮点数或浮点数数组作为参数,其值在 0 到 78 之间。该函数调用从 cpp 文件生成的 exe。然后从该 exe 的输出返回与 x 值相对应的函数值数组。我还有一个名为 x
的数据数组,我想要适应它。我确信
total_data_array
包含正确的值并且具有正确的大小。以下是我的代码的缩短版本:total_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)
我收到错误:
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
所以现在我只是最小化我的功能。稍后我将减去
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)
以最小化函数和数据之间的差异并获得曲线拟合。现在,我只想让它发挥作用。我收到错误:
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
更改为
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()
的关键字参数上。但我不知道要传递什么作为关键字或为什么。 我很感激任何帮助:)
minimize()
是模型函数而不是目标函数:它没有相同的调用签名,这是您看到的异常告诉您的。但尝试 lmfit.Minimize 无论如何都是没有意义的。 lmfit.Model 运行,它只是告诉您它无法估计不确定性。另外,FWIW,如果初始拟合无法估计不确定性,则运行
func
是没有意义的。这不会“解决问题”。正如您在 lmfit FAQ中所读到的,有时无法估计不确定性有几个原因,但它们都源于相同的根本原因:一个或多个参数值的微小变化不会引起明显的变化对贴合度的影响。发生这种情况时,无法确定协方差矩阵 - 拟合算法已确定至少其中一个参数实际上不会影响拟合。 可能发生的原因之一是参数卡在边界或初始值处。这似乎不适合您的情况。另一种可能性是一个或多个参数已达到极值(例如接近 0),以致其他参数不再对拟合产生影响。想象一下,您正在对一个具有中心、宽度和高度的峰进行建模。如果高度值移至零,那么宽度或中心值是多少并不重要,该模型将为零。 我想这就是你的健康状况所发生的情况。您的几个参数的值达到 ~1e-7 级别——比您对它们的初始估计小 5 个或更多数量级。
我必须说,您根本没有描述您的模型,而是运行一些外部程序,然后对结果文本进行一些相当笨拙的解析。我们不知道它们在做什么。因此,其他人无法知道,如果
lmfit.conf_interval
从 0.99 变为 1.4-e7,模型的其余部分甚至是有意义的。也许你知道这一点。
无论如何,这就是我要开始寻找的地方。您似乎没有使用初始值(或最终值)运行模型并绘制结果:始终强烈建议您检查初始值是否合理。由于卡方值如此大且 R 值远离 1,因此很容易断言您的拟合效果很差。如果你不合适,那么不确定性无论如何都是没有意义的。