如何在lmfit最小二乘最小化中包含我的数据错误,以及lmfit中conf_interval2d函数的这个错误是什么?
作者:互联网
我是python的新手,并试图使用lmfit包检查我自己的计算,但是我不确定(1)如何包含错误的以下测试(和2)的数据(sig)的错误得到conf_interval2d如下所示):
import numpy as np
from lmfit import Parameters, Minimizer, conf_interval, conf_interval2d, minimize, printfuncs
x=np.array([ 0.18, 0.26, 1.14, 0.63, 0.3 , 0.22, 1.16, 0.62, 0.84,0.44, 1.24, 0.89, 1.2 , 0.62, 0.86, 0.45, 1.17, 0.59, 0.85, 0.44])
data=np.array([ 68.59, 71.83, 22.52,44.587,67.474 , 55.765, 20.9,41.33783784,45.79 , 47.88, 6.935, 34.15957447,44.175, 45.89230769, 57.29230769, 60.8,24.24335594, 34.09121287, 42.21504003, 26.61161674])
sig=np.array([ 11.70309409, 11.70309409, 11.70309409, 11.70309409,11.70309409, 11.70309409, 11.70309409, 11.70309409,11.70309409, 11.70309409, 11.70309409, 11.70309409,11.70309409, 11.70309409, 11.70309409, 11.70309409,11.70309409, 11.70309409, 11.70309409, 11.70309409])
def residual(pars, x, data=None):
a=pars['a'].value
b=pars['b'].value
model = a + (b*x)
if data is None:
return model
return model-data
params=Parameters()
params.add('a', value=70.0)
params.add('b', value=40.0)
mi=minimize(residual, params, args=(x, data))
#mi=minimize(residual, params, args=(x,), kws={'data': data})#is this more correct?
ci, trace = conf_interval(mi, trace=True)
这个工作到目前为止工作正常,但如上所述,我如何包含数据的错误(sig_chla),以便我可以计算加权和减少的卡方?
第2部分:FURTHERMORE,当我尝试使用以下内容以便绘制置信区间时,
xs,ys,grid = conf_interval2d(mi,’a’,’b’,20,20)
我收到以下错误:
* ValueError:无法创建意图(缓存|隐藏)|可选数组 – 必须具有已定义的维度但得到(0,)
要么
参数4到例程DGESV不正确
DGESV中的Mac OS BLAS参数错误,参数#0,(不可用),为0
解决方法:
您必须自己在residual()函数中通过错误对数据进行加权.
从lmfit docs(虽然不是很容易找到):
Note that the calculation of chi-square and reduced chi-square assume that the returned residual function is scaled properly to the uncertainties in the data. For these statistics to be meaningful, the person writing the function to be minimized must scale them properly.
但是,这并不难.来自weighted least-square fitting的维基百科条目:
If, however, the measurements are uncorrelated but have different uncertainties, a modified approach might be adopted. Aitken showed that when a weighted sum of squared residuals is minimized, is BLUE if each weight is equal to the reciprocal of the variance of the measurement.
然而,lmfit接收残差,而不是平方残差,所以不要只是去
# This is what you do with no errorbars -> equal weights.
resids = model - data
return resids
你应该做这样的事情(scipy作为sp导入):
# Do this to include errors as weights.
resids = model - data
weighted = sp.sqrt(resids ** 2 / sig ** 2)
return weighted
这应该给你正确的加权拟合.
标签:python,least-squares,confidence-interval,uncertainty 来源: https://codeday.me/bug/20190625/1285570.html