在使用lmfit时如何指定Y变量的误差？

Question

我几乎是Python的新手，我正在尝试使用lmfit来适应大学的数据。 Y变量的可变误差为3％。如何将该错误添加到拟合过程中？我正在改变scipy的曲线拟合和scipy，这很容易做到，只需创建一个带有错误值的数组，并在拟合时通过添加文本“sigma = [yourarray]”指定错误这是我当前的代码：

from lmfit import Minimizer, Parameters, report_fit
import matplotlib.pyplot as plt
w1, V1, phi1, scal1 = np.loadtxt("./FiltroPasaBajo_1.txt", delimiter = "\t", unpack = True)

t = w1
eV= V1*0.03 + 0.01

def funcion(parametros, x, y):
    R = parametros['R'].value
    C = parametros['C'].value

    modelo = 4/((1+(x**2)*(R**2)*(C**2))**1/2)
    return modelo - y
parametros = Parameters()
parametros.add('R', value = 1000, min= 900, max = 1100)
parametros.add('C', value = 1E-6, min = 1E-7, max = 1E-5)

fit = Minimizer(funcion, parametros, fcn_args=(t,V1))
resultado = fit.minimize()


final = V1 + resultado.residual
report_fit(resultado)

try:
    plt.plot(t, V1, 'k+')
    plt.plot(t, final, 'r')
    plt.show()
except ImportError:
    pass

V1是我测量的值，eV是错误数组。 t是x坐标。感谢您的时间

Answer 1

minimize()函数在最小二乘意义上最小化数组，调整变量参数以最小化由目标函数返回的(resid**2).sum()数组的resid。它对您的数据中的不确定性甚至您的数据都不了解。要使用拟合中的不确定性，您需要在传递eV和t时传入数组V1，然后在计算数组时将其用于最小化。

人们通常希望最小化Sum[ (data-model)^2/epsilon^2 ]，其中epsilon是数据中的不确定性（你的eV），所以残差数组应该从data-model改为(data-model)/epsilon。为了你的健康，你会想要

def funcion(parametros, x, y, eps):
    R = parametros['R'].value
    C = parametros['C'].value

    modelo = 4/((1+(x**2)*(R**2)*(C**2))**1/2)
    return (modelo - y)/eps

然后使用它

fit = Minimizer(funcion, parametros, fcn_args=(t, V1, eV))
resultado = fit.minimize()
...

如果你使用lmfit.Model接口（专为曲线拟合而设计），那么你可以传入weights数组乘以data -model，因此1.0 / eV代表不确定性的加权（如上所述minimize）。使用lmfit.Model界面并提供不确定性将如下所示：

from lmfit import Model
# model function, to model the data
def func(t, r, c):
    return 4/((1+(t**2)*(r**2)*(c**2))**1/2)

model  = Model(func)
parametros = model.make_params(r=1000, c=1.e-6)

parametros['r'].set(min=900, max=1100)
parametros['c'].set(min=1.e-7, max=1.e-5)

resultado = model.fit(V1, parametros, t=t, weights=1.0/eV)

print(resultado.fit_report())

plt.errorbar(t, V1, eV, 'k+', label='data')
plt.plot(t, resultado.best_fit, 'r', label='fit')
plt.legend()
plt.show()

希望有所帮助....

Answer 2

我想你不能直接在fit.minimize（）中提供sigma。

但是我看到fit.minimize（）使用scipy的minimalsq方法（默认情况下），这与scipy的curve_fit使用的方法相同。

如果你研究scipy's curve_fit source，它会跟随sigma（对于1-d情况）。

 transform = 1.0 / sigma
 jac = _wrap_jac(jac, xdata, transform)
 res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)

由于fit.minimize（）允许你传递kwargs（Dfun）for leastsq，你可以像在scipy curve_fit中一样传递jac。